1 24 3 6 @BA - CiteSeerX

10 downloads 0 Views 1MB Size Report
process was identified by Griffiths in d! e laboratory. Щ models of rising thermals 29 , and was also. Q observed in models of the development of large f! g plume.
EPSL ELSEVIER

Earth and Planetary

Science Letters 148 (1997) 13-26

The effect of temperature dependent viscosity on the structure of new plumes in the mantle: Results of a finite element model in a spherical, axisymmetric shell Louise H. Kellogg a3*, Scott D. King b~l a Department of Geology Unir~ersityof California. DaLIis,CA 95688, USA ’ Department of Earth and Atmospheric Sciences, Purdue University West Lafayette. IN 47907, USA Received 5 July 1996; revised 18 January

1997; accepted 3 February

1997

Abstract developed a finite element model of convection in a spherical, axisymmetric shell that we use to simulate thermal plumes in the mantle. The finite element method provides the flexibility to include realistic properties

We have

upwelling

such as temperature-dependent viscosity, the focus of this paper. We used this model to investigate the effect of temperature-dependent viscosity on the structure of new plumes originating at the core-mantle boundary. Because the way in which mantle viscosity varies with temperature is not well constrained, we determined the plume structure using a variety of viscosity laws. We focus on 3 different viscosity laws: (1) constant viscosity; (2) weakly temperature-dependent viscosity, in which the viscosity increases by a factor of 10 between the hottest and the coldest material; and (3) strongly temperature-dependent viscosity, in which the viscosity varies by a factor of IOOO. In a constant viscosity fluid, the plume exhibited a spout structure without a distinctive head. The plume head and tail consisted largely of material from the hot thermal boundary layer at the base of the spherical shell that represented the mantle. When the viscosity was strongly temperature-dependent, starting plumes developed a mushroom structure with a large, slow-moving head, followed by a narrow, faster moving tail. Material from the overlying shell was assimilated into the plume head during formation of the upwelling in models with strongly temperature-dependent viscosity, while the plume tail showed little entrainment Constant viscosity models and models with weakly temperature-dependent viscosity showed almost no entrainment in the head or tail. The large plume head that formed in models with strongly temperature-dependent viscosity created and then shed “blobs” of material from the deep mantle that did not arrive at the surface near the plume but instead were deposited elsewhere in the upper mantle. Keywords: mantle plumes; geodynamics;

convection;

hot spots

1. Introduction

* Corresponding author. Tel.: + 1 916 752 3690. Fax: + 1 916 752 0951. E-mail: [email protected] ’ Tel.: + 1 317 494 3696. E-mail: [email protected]. 0012-821X/97/$17.00 Copyright PII SOOl2-821X(97)00025-3

In 1971, W.J. Morgan suggested that intraplate volcanism such as the hotspot at Hawaii derives from hot upwelling plumes in the mantle [ 11. Our current understanding of the mantle suggests that

0 1997 Elsevier Science B.V. All rights reserved.

14

L.H. Kellogg, S.D. King/Earrh

and Planetary Science titters 148 (1997) 13-26

convection takes the form of subducting plates and upwelling, approximately axisymmetric plumes. According to this model, hotspots are the surface manifestation of mantle plumes. The topography and gravity anomalies associated with hotspots have been used to determine that actively upwelling mantle plumes account for about 10% of the total heat loss of the Earth [2,3]. The source of mantle plumes remains a matter of considerable debate. A widely accepted picture is that in the Earth’s mantle, a hot thermal plume, originates as an instability at a thermal boundary layer, forming a large mushroomshaped head as it grows. This head eventually rises from the boundary layer followed by a narrow tail. These tails probably take the form of pipe-like features [4] that are less than a few hundred kilometers in diameter, probably just barely detectable by seismic methods in a few specific favorable locations

El. Geochemical studies provide another source of information about mantle plumes and show that oceanic island basalts fall into several groups according to their radiogenic isotopic signatures [61. These indicate a variety of influences on oceanic island basal&, including the presence of recycled crustal materials in some cases [7], with a possible contribution of primordial material in other locations [8], but do not provide constraints on the dynamics of the mantle source of oceanic islands. Correlations of helium and lead isotope data from mid-oceanic ridge basalts and oceanic island basalts suggest a source for both located near the mantle transition zone [9]. The noble isotope measurements have been interpreted as suggesting that plumes transport some material from the deep mantle to the surface [8] and have been used to estimate rates of mass transport across the transition zone [ 10,111; however, it is also possible that the transition zone is the main source of these anomalies. Recent observations of high 1870s/ ‘**OS in some oceanic island basalts have been cited as evidence that at least some mantle plumes originate at the core-mantle boundary [12]. No matter where they originate, as plumes rise they may, under the right circumstances, entrain significant quantities of material. Laboratory studies [ 131 and numerical simulations of plumes [14] both exhibit entrainment. Evidence for mixing between multiple sources is also seen at some ocean islands;

detailed study of the Galapagos Islands, for instance, indicates a complex pattern of mixing and entrainment between the primary source and the overlying asthenosphere [ 15,161. Material from plumes may be a significant source of material in the asthenosphere that eventually feeds into ridges [ 171. Information about mantle plumes, including their source, is indirect and subject to multiple interpretations. Careful modeling constrained by observations is essential to understanding the dynamics and structure of plumes. Numerical models of steady-state convection in an axisymmetric geometry show that upwelling plumes occur only in the presence of basal heating [ 181. The plume flow narrows and is more concentrated when temperature-dependent viscosity is introduced [ 181.Three-dimensional computer models of convection in spherical shells show that flow at low Rayleigh numbers naturally takes the form of narrow, upwelling plumes and broad, downwelling sheets [4,19]. Three-dimensional calculations in a flat layer show that boundary layer instabilities develop as long, ridge-like features [20]. Where these ridges intersect, they coalesce to form more pipe-like upwellings and downwellings. From 3D studies, phase changes seem to alter that overall pattern significantly, reducing the effect of the upwelling plumes [21]. However, in recent studies focusing on plumes with temperature-dependent viscosity, plumes can penetrate all but the strongest endothermic phase changes, with some disagreement as to the exact magnitude of the Clapeyron slope needed to stop the plume from penetrating into the upper mantle [22,23]. It appears that the form of the viscosity law and the temperature contrast across the bottom boundary layer may be the cause of the discrepancy. In addition, the choice of vertical viscosity stratification across the mantle, the ratio of internal to bottom heating, and the effects of viscous dissipation and compressibility all appear to influence the formation of plume-like structures in 3D thermal convection models (e.g., [24]). However, in all published 3D calculations the upwelling plumes are larger than estimates of mantle plumes based on surface observables and should be features that could be resolved with seismic techniques. Laboratory experiments with sucrose solutions provide striking models of the structure of plumes [25-27,13,28]. The laboratory experiments show that

15

L.H. Kellogg, S.D. King/ Earth and Planetary Science Letters 148 (1997) 13-26

SPO”I

Mushrwm

BA.Xl

Fig. I. Schematic profiles of different types of starting plumes. These plumes may take the form of a spout, a balloon or a mushroom. Only the mushroom shape involves significant entrainment of material into the plume head. After [29].

there are three idealized shapes for plumes: spoutshaped, where the tail of the plume is nearly as large as the head; mushroom-shaped; and balloon-shaped (Fig. 1). In the latter two cases, the head is much larger than the tail and takes a form implied by the name. The mushroom-shaped plumes entrain surrounding material, provided the head expands as it rises [29,30]. Several recent numerical studies have successfully produced plumes in a thermally convecting viscous fluid with similar features to those seen in the laboratory [23,14]. These calculations provide a different interpretation of entrainment of material in plumes. In the laboratory experiments, the interpretation is that a significant amount of material from the thermal boundary layer is entrained in the plume head, while the numerical experiments show that entrainment may take place at shallower levels than suggested by laboratory experiments. In this paper we present the results of finite element calculations of convection in a spherical shell with axisymmetry and temperature-dependent viscosity (Fig. 2). By imposing axisymmetry, we restrict the problem to two degrees of freedom, reducing the computational effort significantly over 3D calculations. Hence we can use high resolution calculations to determine the structure of individual plumes with variable viscosity laws. Using a temperature-dependent rheology, we calculate the structure of upwelling plumes for viscosity laws ranging from constant viscosity to strongly temperature-dependent viscosity. The structure of a starting plume depends strongly on the viscosity law. Plumes in a material with constant viscosity have a nearly constant diameter from top to bottom. In contrast, plumes in a strongly temperature-dependent viscosity material exhibit a distinctive “mushroom” shaped head ac-

companied by a narrow tail in which the flow velocities are large. When a weakly temperature-dependent viscosity law is used, starting plumes resemble the constant viscosity plumes. It has been suggested that the D” region of the Earth is the source of mantle plumes [31]. This region may also be chemically distinct from the overlying mantle, either due to a flux of material across the core-mantle boundary [32], or due to the presence of ancient subducted material, undifferentiated mantle material, or melt (see [33] for a review of the characteristics of D”). We have previously investigated the interaction of mantle plumes with negatively buoyant chemical heterogeneities at the core-mantle boundary [34]. It is not the goal of this present work to examine the effect of chemical heterogeneities on the dynamics of plume formation.

Axis of symmetry I

thehth

Mantle

Fig. 2. Geometry of the finite element model. The solution is restricted to axisymmetry about a central pole. Thus, the solution is free to vary as a function of r and 8, but is independent of 4. A patch on the inner shell is heated by maintaining it at constant temperature; the shell is cooled from above. In the figures that follow, the entire solution can be pictured by imagining that the calculated temperatures, velocities, and compositions shown are rotated around the axis of symmetry to sweep out a cone.

16

LH. Kellogg, SD. King/Earth and PlanetaryScience Letters148 (19971 13-26

However, we have included here a passive tracer (dye) in the lowermost 1/15th of the shell to determine whether the plume head that forms is composed primarily of material from the lowermost 200 km or of entrained overlying material. Here we focus on the effect of temperature-dependent viscosity on the structure of the starting plume. We will show that, in constant viscosity models, most of the material in the plume came from the base of the shell and there was little entrainment in the rising plume. In contrast, in a strongly temperature-dependent material, circulation in the plume head caused entrainment and mixing of material from the base of the shell with overlying material. A smaller plume head developed for weakly temperature-dependent viscosity materials. The size of the plume head, the velocities in the tail, and the amount of entrainment, all depended strongly on the specific rheology.

2. Governing equations and solution method Convection in the mantle can be approximated using the Boussinesq equations for flow in an incompressible fluid in a spherical shell. In order to create high resolution models within the limits of available computer power, we have restricted the flow to two degrees of freedom. While 3D calculations with temperature- and stress-dependent viscosity are now possible, 2D calculations are less computer intensive, allowing greater freedom to explore parameters and analyze results. The goal of the current work is to model the structure and evolution of mantle plumes; we therefore developed this model in spherical geometry with axisymmetry about the pole at 0 = 0 (Fig. 2). The temperature and velocity are functions of the radius, r, and the colatitude, 19, and are independent of the longitude, 4. With an axisymmetric analysis we can ensure that the plume will always rise at the pole, making analysis of the calculations easier. The equations governing the conservation of mass, momentum and energy are solved using the finite element code SCAM (Spherical Convection in an Axisymmetric Mantle). SCAM is a spherical, axisymmetric version of the finite element code ConMan [35].

2.1. Momentum and mass conservation For non-uniform viscosity, the equations for momentum and mass conservation, in dimensionless form, can be written: V.( p&) =RaTf

-VP+

(1)

and: v* u = 0

(2)

where P is the dimensionless dynamic pressure; p is the dimensionless viscosity, which can be a function of temperature, pressure, composition and/or strain rate; T is the dimensionless temperature; u is the dimensionless velocity and: &=

;(ui.j +uj,i)

is the deviatoric strain rate. All of the dimensional parameters in the single, dimensionless number, number:

Ra =

po gaAT(

(3) are contained the Rayleigh

r, - r;j3 (4)

POK

where p. is the density at T = 0; g is gravitational acceleration; a is the coefficient of thermal expansion; AT is the temperature drop across the shell (AT = Tpatch- Tsurface); r0 and r, are the outer and inner radii of the shell; K is thermal diffusivity; and p. is viscosity. In models using temperature-dependent rheology, the Rayleigh number becomes more complicated because viscosity p varies with temperature, a function of position and time. We use Eq. (4) to define an effective Rayleigh number using po, the viscosity at a dimensionless temperature of 1.0, the temperature at the base of the shell. To facilitate comparison between models, we use Ra = lo7 in all the calculations that follow. As is standard with most finite element formulations, these equations are converted to integral equations, divided into elements and integrated numerically. To enforce incompressibility, we use the penalty formulation (i.e., Eq. (2)) [36,37]. The large value of the penalty parameter leads to a stiff set of equations that is hard to solve numerically. In a

L.H. Kellogg, S.D.

King/Earth and PlanetaryScience Letters 148 (1997113-26

Cartesian geometry, we reduce the number of quadrature points used to integrate the terms including the penalty parameter (reduced integration). In spherical geometry, we use the B approach, which reduces to the reduced integration technique in Cartesian geometry [38,39]. In non-Cartesian geometries, the B approach is necessary to achieve accurate solutions to viscous flow problems using the penalty method. 2.2. Energy equation The dimensional form of the energy equation spherical axisymmetric geometry is given by:

=--

K

in a

a

Y? ar where T is the temperature; t is time; K is the thermal diffusivity; u, is the radial component of velocity; and I(@ is the azimuthal component of velocity. Eq. (5) is solved using Streamline-Upwind Petrov-Galerkin (SUPG) shape functions 1401. With bi-linear elements, this method is second-order accurate in space. Because the energy equation is transformed onto a Cartesian parent element to perform all of the integrations, the application of SUPG in spherical axisymmetric geometry is straightforward and we refer the reader to the above references for further details. 2.3. Time stepping The algorithm for time stepping the energy equation is a version of the predictor-corrector algorithms proposed by Hughes [41,42]. We have implemented both the implicit and explicit time stepping. The explicit method is second-order accurate in time, and the time step is limited by the Courant condition. The implicit method allows us to circumvent the Courant condition and take larger time steps, but it is only first-order accurate in time. All of the results presented in this paper use the explicit, second-order time-stepping scheme.

17

2.4. Rheology

In modeling mantle convection it is appropriate to incorporate realistic variations in the rheology (e.g., [43]). There are no direct measurements of the creep of (Mg,Fe)SiO, perovskite or magnesiowiistite at lower mantle pressures and temperatures. Single crystal creep experiments on analog perovskite structures agree that the creep deformation is a strong function of temperature but disagree on the mechanism of deformation. Karat0 and Li [441 find that diffusion creep is the dominant mechanism in calcium perovskite (CaTiO,), while Wright et al. [45] find a power-law mechanism dominates. Phase changes and variations in composition may also affect the rheology, but we do not consider those effects here. In this model we have assumed a diffusion creep mechanism, in which the viscosity law depends exponentially on temperature:

---

1+to (6) 1

where E * is the activation energy; R is the universal gas constant; AT is the temperature drop across the shell; T is the dimensionless temperature; t, is the surface temperature divided by the temperature drop across the shell; and p0 is the reference viscosity at T = 1.0. Since the rheological properties of the lower mantle are poorly understood [46], we choose several “ideal” rheologies: a uniform viscosity law ( E */RAT = 01, a weakly temperature-dependent viscosity law (E */RAT = 0.25328), and a stronger temperature-dependent viscosity law (E * /RAT = 3.0). The weakly temperature-dependent viscosity law corresponds to a viscosity contrast of a factor of 10 between the hottest and coldest material. The choice of this viscosity law was motivated by the observation that the interior temperature of a convective fluid adjusts so that the viscosity drop at the base of the layer is about one order of magnitude law (e.g., [471X Th e stronger temperature-dependent corresponds to a viscosity contrast of more than 1000 between the hottest and coldest material. Viscosities exceeding 1000 times pcLoare “clipped” to 1000 I_L~.This viscosity law may be appropriate if a

L.H. Kellogg, SD. King/Earth

18

and Planetary Science Letters 148 (1997) 13-26

chemical boundary layer exists at the core-mantle boundary. The three idealized rheologies allow us to examine the role of temperature dependence on the formation, structure and mixing of plumes. Calculations were also made for a wide range of viscosity laws between the weakly and strongly temperaturedependent rheologies.

with a ratio of inner to outer radius of ri/rO = 0.55 and colatitude from 13= 0 to 6 = 9r/4 (Fig. 2). This section of the shell was chosen to simulate the structure of upwelling mantle plumes originating from the core-mantle boundary. Free-slip boundary conditions were imposed on all boundaries. A patch on the inner boundary of the shell was heated by maintaining it at a constant dimensionless temperature of T = 1; elsewhere the base of the shell was insulating (U/6? = 0). At the top of the shell, a constant dimensionless temperature of T = 0 was maintained. The sides of the cone were insulating (dT/% = 0).We began each calculation with an

3. Onset of plumes in a spherical shell We modeled convection, driven by thermal buoyancy only, in a conical wedge of a spherical shell

Constant t =

0.000135

t =

Weakly

Temperature-Dependent t =

Strongly t =

0.003298

Viscosity

0.000186

0.000561

t =

0.016353

Viscosity

0.000222

Temperature-Dependent t =

t =

0.006945

Viscosity

Fig. 3. Contours of temperature (right half of each inset) and concentration of passive dye (left) for starting plumes in materials with constant viscosity (top row), weakly temperature-dependent viscosity (E * /RAT = 0.25328) (middle), and strongly temperature-dependent viscosity (E * /RAT = 3.0) (bottom row). Non-dimensional times in the right-hand column correspond to 1.5 X 10s yr (top), 2.0 X 10s yr (middle), and 4.4 X lo9 yr bottom).

19

L.H. Kellogg, S.D. King/ Earth and Planetary Science Letters 148 (1997) 13-26

initial dimensionless temperature of T = 0.25 throughout the interior of the shell, close to the average conductive temperature profile for the spherical sheil. We ran calculations on two grids: a grid consisting of 64 elements in the radial direction and 32 in the azimuthal direction, both evenly spaced in r and 8; and a higher resolution grid of 96 by 96 elements, in which elements were concentrated along the axis of symmetry and along the boundary layers. In the radial direction, half of the elements were concentrated between 1.2 < r < 1.6; in the azimuthal direction, two-thirds of the elements were concentrated between the axis of symmetry 6 = 0 and 13= r/16. Our results were not significantly different for the two grids. In the following, we show results obtained using the higher resolution grid for the three viscosity laws discussed. Additional calculations were run at lower resolution using a wide variety of viscosity laws; these calculations were used to obtain the relationship between viscosity law and plume head size, discussed at the end of this section.

..

Temperature Composition

-

I (b) Strongly T-dep. viscosity

‘. “.-

. . . . . . . . . . . . . . . ..S.............

I

0

0

n/16

46

e

Fig. 5. Profiles of temperature and composition as a function of angle for r = 1.700, approximately the narrowest part of the plume tail, for: (a) constant viscosity and (b) strongly temperature-dependent viscosity. There is little entrainment and mixing in the plume tail for all viscosity laws. Profiles are taken at the same times as the right-hand column in Fig. 3. The dots in the dotted temperature curve show the spacing of the finite element grid; 2/3rd of the elements are concentrated within 0 < 0 < n / 16.

....

5

-

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i

Temperature Composition

-

0

n/a

d16

3.1. Onset of plumes wirh constant uiscosiiy

1.. *A-x.

e 0

416

n/6

8

Fig. 4. Profiles of temperature and composition as a function of angle for r = 1.700, approximately the widest part of the plume head in Fig. 3, center, for (a) constant viscosity and (b) strongly temperature-dependent viscosity. The smooth, monotonically varying gradients of temperature and composition in the constant viscosity model show that there is little entrainment and mixing in this plume. Significantly more entrainment and mixing are shown when the vicosity is strongly temperature-dependent. Profiles are taken at the same times as the middle column in Fig. 3. The dots in the dotted temperature curve show the spacing of the finite element grid; 2/3rd of the elements are concentrated within 0