Entrainment of a dense layer by thermal plumes

Summary

The large mismatch between Earth's heat flux and that which would be produced in a uniform mantle with the composition of the mid-ocean ridge basalt source region, and several recent seismic models, suggest that the bottom part of the mantle may form a compositionally distinct layer. The long-term stability of such a dense but hot abyssal layer may hinge on the rate of entrainment of the dense material by upwelling thermal plumes that would originate from the thermal boundary layer at the compositional boundary. We have formulated high-resolution numerical models to examine the efficiency of such an entrainment. To isolate the physics of entrainment, thermal buoyancy forces that drive the flow and entrain the dense material are prescribed in both space and time in axisymmetric finite-element models with an isoviscous structure. Our models employ a marker chain method to track the evolution of the material interface. Thermal plumes entrain the dense material, forming a concentric annular structure, with a hot thermal plume surrounding an inner cylinder of dense material. The entrainment rate is controlled by two parameters: the radius of the thermal plume, r_T, and the ratio of compositional to thermal buoyancy, Ra_bt. The smaller r_T is or the larger Ra_bt, the smaller the entrainment rate will be, as expected. As Ra_bt increases, the radius of the entrained compositional plume decreases, but the vertical velocity of the compositional plume also increases. We found that the entrainment rate scales as Ra^−2.48_btr^3.80_T, while the radius of the entrained compositional plume scales as Ra^−1.43_btr^1.24_T For a mantle viscosity of 10²¹ Pas, thermal plumes with r_T∼ 100 km and temperature 600 K hotter than the ambient mantle, and with horizontal spacings that are approximately the same as the mantle thickness, and for a denser layer with a thickness of 1000 km and net negative buoyancy of approximately 1 per cent, more than 90 per cent of the dense material can survive for 4.5 Gyr.

entrainment, mantle flow, thermochemical convection

1 Introduction

Earth's mantle loses approximately 38 TW through plate tectonic processes. However, radioactive heating in a uniform mantle with the composition of the Mid-Ocean Ridge Basalt (MORB) source region would only produce 2–6 TW, because the MORB source region is depleted in heat-producing elements (O'Nions & Oxburgh 1983; Davies & Richards 1992; Kellogg et al. 1999). This large mismatch in the heat budget suggests that the mantle is compositionally stratified, with a bottom layer that has a higher heat production rate than that in the top layer (i.e. MORB sources) (O'Nions & Oxburgh 1983). This suggests that mantle convection may take place in two separate layers (Kellogg et al. 1999).

Early layered convection models suggested that the layering occurs at the 670 km seismic discontinuity (Wasserburg & DePaolo 1979; O'Nions & Oxburgh 1983). However, this depth is inconsistent with recent high-resolution seismic images of the mantle that indicate penetration of subducted lithosphere into the lower mantle (Grand et al. 1997; van der Hilst et al. 1997). However, the seismic images may be consistent with mantle layering at a depth of ∼1700 km (van der Hilst & Karason 1999). It is important to examine the dynamic stability of these models. For the layered mantle convection model, the stability of a layered mantle requires that enriched material in the bottom layer be intrinsically denser than the depleted mantle in the top layer (e.g. Kellogg et al. 1999).

In this study, we investigate the dynamic stability of a layered mantle. In general, a layered mantle is stable against major overturns as long as the bottom layer is sufficiently dense (Davies & Gurnis 1986; Tackley 1998). In this case, the long-term stability of the bottom layer should be controlled by gradual entrainment of the dense material by upwelling thermal plumes (Sleep 1988). Ascending plumes in the upper layer entrain the dense material of the bottom layer into the plumes, gradually eroding the bottom layer. This regime of stability is particularly relevant for the newly proposed layered convection model (Kellogg et al. 1999; van der Hilst & Karason 1999), because the proposed in situ density contrast must be small in order not to violate constraints from seismology.

Previously, the entrainment of dense material by upwelling thermal plumes has been studied with both laboratory studies (Davaille 1999; Gonnermann et al. 2002) and 1-D analytic models (Sleep 1988). Sleep (1988) employed a 1-D model to study the long-term stability of the D′′ layer. However, for a given compositional buoyancy for the dense material and thermal buoyancy for a plume, a 1-D analytic model cannot uniquely determine the entrainment rate. For this reason, Sleep (1988) only estimated the maximum entrainment rate, which he assumed to be the carrying capacity of thermal plumes. Sleep's studies (1988) suggest that the D′′ layer needs to be 6 per cent denser than the overriding mantle in order for the D′′ layer to maintain itself over geological history. A more accurate determination of the entrainment rate requires considering the flow at the two ends of a plume, for which numerical modelling is necessary. Indeed, numerical models suggest that a smaller density contrast is sufficient to maintain the long-term stability of the D′′ layer (Tackley 1998). However, intrinsically small-scale compositional features demand high resolution in numerical models (van Keken et al. 1997; Tackley 1998; Tackley & King 2002).

We formulate high-resolution finite-element models in a 2-D axisymmetric geometry to study how the entrainment rate depends on the thermal buoyancy of plumes and the density contrast between the two layers. Compared with previous numerical studies (van der Hilst et al. 1997; Tackley 1998), our models use much higher resolution and explore a larger parameter space in a more systematic manner. Our models also use a larger thickness for the dense layer in order to be consistent with the newly proposed layer thickness (van der Hilst & Karason 1999).

2 Physical Models and Numerical Formulation

The mantle is assumed to consist of two layers with different intrinsic densities (Fig. 1). We treat the mantle as an incompressible viscous fluid, in which case the equations governing mantle flow can be expressed as  

(1)

 

(2)

where u is the velocity, P is the pressure, η is the viscosity, g is the gravitational acceleration, e_z is the unit vector in the vertical direction and δρ is the density perturbation that is due to compositional and/or thermal effects:  

(3)

Here ρ₀ and T₀ are the reference density and the reference temperature, α is the thermal expansion coefficient, T is the temperature, C specifies the compositional field and Δρ is the compositional density difference between the two layers. In our two-layer models, C is zero and one for materials in the top and bottom layers, respectively.

Figure 1

A schematic representation of the model setup.

The transport equation for the compositional field is given by (Kellogg & King 1993; Lenardic & Kaula 1993)  

(4)

where t is time. The energy balance equation has a similar mathematical form to eq. (4) except that an additional diffusion term governs the time evolution of temperature (van Keken et al. 1997).

Previous numerical studies that have focused on the dynamics of thermochemical convection often consider the transport equations for both composition and temperature (e.g. Kellogg & King 1993; Tackley 1998). To resolve small-scale compositional features and thermal boundary layers, high numerical resolution is required. Since the entrainment is mainly controlled by thermal plumes, following Sleep (1988), we consider only the transport equation for composition and assume that the thermal structure for plumes, which we fix a priori, does not change with time (Fig. 1). This simplification enables us to have the resolution necessary to resolve both the thermal and compositional plumes. In our calculations, the horizontal resolution is ∼3 km at the centres of the plumes. Without the need to resolve the transient effects of thermal boundary layers, this simplification also makes it easier to quantify the entrainment rate.

In our model, the temperature for the thermal plume, T_p, is assumed to be the same as that for the bottom layer, whereas the top layer, except for the thermal plume, has a background temperature T₀ (Fig. 1). This approximation may be justified because the bottom layer is approximately isothermal for a layered convective system (Tackley 1998; Kellogg et al. 1999). The bottom layer temperature may be somewhat larger than T_p, but this does not significantly affect our analysis. We treat the plume temperature, T_p, and the plume radius, r_T, as independent variables that control the plume buoyancy. A thermal plume entrains dense material from the bottom layer, forming a concentric annular structure, with a hot thermal plume surrounding an inner cylinder of dense material (Fig. 1) (Sleep 1988). For different compositional (Δρ) and thermal (T_p−T₀ and r_T) buoyancy, we determine the entrainment rate (i.e. the mass flux for the compositional plume) and the radius of the compositional plume, r_c.

We substitute eq. (3) into eqs (1) and (2) and normalizing variables as  

(5)

where r denotes a spatial variable, η₀ is the background viscosity, D is the thickness of the mantle and all the variables with a prime are non-dimensional. While the resulting non-dimensional equations for the continuity equation and compositional transport equation have the same forms as eqs (1) and (4), the non-dimensional momentum equation can be written as (dropping all of the primes, for simplicity)  

(6)

where Ra_bt=Δρ/[ρ₀α(T_p−T₀)] is the ratio of compositional to the thermal buoyancy, C is 1 for the denser material in the lower layer and 0 for material in the upper layer, and T is 1 for the thermal plume and denser material and 0 for the rest of the mantle.

The momentum and continuity equations are solved using the finite-element software Citcom, which uses an Uzawa algorithm for the pressure field and a multigrid solver for the velocity (Moresi & Solomatov 1995; Moresi & Gurnis 1996). We have developed a marker chain method to solve the compositional transport equation (i.e. eq. 4). We have also extended the finite-element formulation to axisymmetric geometry. These new features are presented in. The validity of the marker chain method and our finite-element analysis are verified in the results section.

In our calculations, the mantle is assumed to be isoviscous and to have a constant coefficient of thermal expansion. All the boundaries are assumed to be free-slip. Initially the compositional boundary is set at z=z₀ and is flat. With these initial and boundary conditions, and for a given radius of a thermal plume, r_T, and the ratio of compositional to thermal buoyancy, Ra_bt, we solve the momentum, continuity and compositional transport equations to determine the mass flux and geometry for the entrained compositional plume.

A stream function φ is used to illustrate the flow field, and it can be related to the flow velocity as  

(7)

The stream function automatically satisfies the continuity equation.

3 Results

In this section, we first present a representative case with r_T= 0.05 and Ra_bt= 1.5 to show the basic physics of entrainment and the validity of our numerical techniques. We then demonstrate how the entrainment rate depends on r_T, Ra_bt and other parameters.

3.1 Entrainment For Rt= 0.05 And Rabt= 1.5

For this case the compositional boundary is initially set at z= 0.25. At the beginning of the model run, thermal buoyancy in the plume entrains the compositionally dense material into the thermal plume. A compositional plume gradually forms at the centre of the thermal plume (Fig. 2a). Upon reaching the surface, the compositional plume spreads out below the surface (Fig. 2a). The upward vertical flow near the plume indicates that the thermal plume entrains the compositionally dense material (Fig. 3a). The entrainment causes the volume of the bottom layer to decrease with time (Fig. 2b). The volume of the bottom layer, V(t), is determined with the following integral:  

(8)

where Z(r, t) is the z-coordinate for the compositional boundary for the bottom layer and is calculated with our marker chain technique. Except for the initial phase, V(t) decreases with time at a nearly constant rate (Fig. 2b; the entrainment rate can be determined from the slope of the V(t) versus t plot and is approximately 1.07 × 10⁻⁶ for this case). We also compute the average radius for the compositional plume as a function of time, r_c(t); the average is calculated between z= 0.45 and 0.65. Similar to the entrainment rate, r_c(t) also reaches a constant value shortly after the initial phase (Fig. 2b). We terminate the model run after the compositional plume fully spreads under the surface (Fig. 2a). We think that the quasi-steady-state entrainment rate is representative of the entrainment by a thermal plume in a two-layer system.

Figure 2

(a) Compositional boundary at four different non-dimensional times: 209, 572, 2590 and 6640. The origins of the coordinate system for the last three times have been shifted to better show the geometry of the compositional boundary. (b) Time evolution of the volume of the dense layer V and the radius of the compositional plume r_c (averaged over the depth range from z= 0.45 and 0.65). For this calculation, Ra_bt= 1.5 and r_T= 0.05.

Figure 3

(a) Flow field at t= 6640 for the calculation with Ra_bt= 1.5 and r_T= 0.05. A positive stream function (thin solid lines) represents clockwise circulation, while a negative stream function (thin dashed lines) is for counter-clockwise circulation. The contour levels are −0.125, −0.025, 0.025, 0.125, 0.225, 0.325 and 0.425. The thick solid and thick dashed lines are for the compositional boundary and the thermal plume, respectively. Note that only the flow field for r < 0.25 is displayed. (b) Spatial distribution of the three terms of the z-component momentum equation for the flow in Fig. 3(a). The solid and dashed lines are for the compositional boundary and the thermal plume, respectively. (c) Horizontal dependence of the dτ_rz/dr+τ_rz/r (thin line) buoyancy force (dashed line) and τ_zz,z−P_,z (thick line) that are averaged over 0.55 < z < 0.65 for the flow in Fig. 3(a).

The flow field indicates that, except for the regions where the compositionally dense material enters and exits the plume, the flow in the vicinity of the plume is nearly 1-D and only depends on the radial coordinate (Fig. 3a). This can be clearly seen in plots for different terms in the z-component of the momentum equation: ∂P/∂z, ∂τ_zz/∂z and ∂τ_rz/∂r+τ_rz/r (Fig. 3b). Among these terms, ∂P/∂z and ∂τ_zz/∂z are much smaller than ∂τ_rz/∂r+τ_rz/r near the plume, except near its ends (Fig. 3b), indicative of 1-D flow. The term ∂τ_rz/∂r+τ_rz/r is dominant near the plume, except at its ends, and is balanced by the buoyancy term T−Ra_btC (Fig. 3c for these quantities averaged over 0.55 < z < 0.65). We believe that the near cancellation of ∂τ_rz/∂r+τ_rz/r and T−Ra_btC validates our numerical method. However, the 2-D flow at the ends of the plume (Figs 3a and b) is important in determining the entrainment rate, so a 1-D approximation to the entire problem is insufficient to accurately determine the entrainment rate.

3.2 Dependence of Entrainment Rate on r_T and Ra_bt

We now present results for calculations in which the initially flat compositional boundary is at z₀= 0.25, but the radius for the thermal plume, r_T, and the ratio of compositional to thermal buoyancy, Ra_bt, are varied systematically. When r_T is fixed at 0.05, increasing Ra_bt results in a narrower compositional plume, with a smaller entrainment rate (Fig. 4a). Likewise, when Ra_bt is fixed at 1.5, decreasing r_T has a similar effect on the entrainment rate as increasing Ra_bt for a fixed r_T (Fig. 4b). This is because decreasing r_T reduces the thermal buoyancy driving the flow, while increasing Ra_bt enhances the compositional buoyancy resisting the flow, both of which hinder entrainment. For all of these calculations, except for the initial stages, the entrainment rate, Q=dV/dt, and the average compositional plume radius, r_c, are nearly constant (Fig. 4). By averaging Q and r_c over the period of the quasi-steady state, we can determine the average entrainment rate, Q, and the compositional plume radius, r_c, that we take as representative measures of the dynamics of entrainment.

Figure 4

(a) Time evolution of the volume of dense layer, V, and the averaged radius of the compositional plume, r_c, for r_T= 0.05 but different Ra_bt. (b) Time evolution of the volume of dense layer, V, and the averaged radius of the compositional plume, r_c, for Ra_bt= 1.5 but different r_T.

Our calculations show that for a given r_T, both the entrainment rate, Q, and the compositional plume radius, r_c, decrease with increasing Ra_bt and that for a given Ra_bt, Q and r_c increase with increasing r_T (Figs 5a and b). In addition to Q and r_c, we also compute the averaged vertical velocity within compositional plumes, u_zc, over the same depth ranges for which r_c is determined. For a given r_T, u_zc increases with increasing Ra_bt (Fig. 5c). This is very different from Q and r_c, which decrease with increasing Ra_bt (Figs 5a and b). The entrainment rate, Q, is the mass flux of the dense material through the compositional plume, so Q∼u_zcr_c². Therefore, Q is more sensitive to the variation of r_c with Ra_bt than to the variation of u_zc with Ra_bt. The decrease of Q with Ra_bt (Fig. 5a) arises because r_c decreases much more rapidly than u_zc increases as Ra_bt increases (Figs 5b and c).

Figure 5

Dependence of (a) the averaged entrainment rate, Q, (b) the averaged radius of the compositional plume, r_c, (c) the averaged vertical velocity for the compositional plume, u_zc, on Ra_bt and r_T, and fitting of polynomial functions to Q (d), r_c (e) and u_zc (f). While the compositional boundary is initially at z=z₀= 0.25 for most calculations, the dashed lines with squares and triangles in (a)–(c) are for calculations with z₀= 0.5 and 0.1, respectively. In (d)–(f), the horizontal and vertical axes represent results from finite-element calculations and from the fitting functions, respectively.

In all of the calculations that we have presented so far, the initially flat compositional boundary is set at z₀= 0.25. It is important to examine whether our results are sensitive to z₀, because in a dynamic mantle with two layers, the volume of the bottom layer must change with time as a result of continuous entrainment. To address this issue, we compute two sets of models with z₀= 0.1 and 0.5, respectively. For each set of models, r_T= 0.05, but Ra_bt varies systematically. Our calculations show that both Q and r_c are insensitive to z₀ (Figs 5a and b). However, compared with cases with z₀= 0.25, cases with z₀= 0.5 yield smaller u_zc, while cases with z₀= 0.1 lead to larger u_zc (Fig. 5c). This is expected because the net thermal buoyancy for the plume decreases as z₀ increases. The variations in u_zc with z₀ do not have a significant effect on Q, because of the corresponding changes in r_c that are relatively small but sufficient to offset the effects of u_zc on Q (Fig. 5b). This suggests that the entrainment rate for such a two-layer system is mainly controlled by r_T and Ra_bt, with a very weak dependence on z₀ or the volume of the bottom layer.

For models with z₀= 0.25, we fit each of Q, r_c and u_zc with a polynomial function aRa^b_btr^c_T to seek their dependence on r_T and Ra_bt. We found that Q∼Ra^−2.48_btr^3.80_T, r_c∼Ra^−1.43_btr^1.24_T and u_zc∼Ra^0.45_btr^1.28_T (see Table 1 for the fit for constants a, b and c, and their errors). Figs 5(d)–(f) and Table 1 show that these scaling laws fit our numerical results quite well with errors of less than 8 per cent, except for u_zc for which the error for constant b is ∼24 per cent. We also noticed that these scaling laws satisfy Q∼u_zcr²_c, as expected. It should be pointed out that the scaling laws for Q and r_c do not explain well the two cases with smallest Q and r_c (Figs 5d and e). We believe that these two cases with their relatively small r_T (i.e. 0.025D) and relatively large Ra_bt (i.e. Ra_bt≥ 2) are the most demanding numerically and that our numerical resolution for these two cases may not be high enough.

Table 1.

Scaling law parameters† for Q, rc, and uzc.

4 Discussion

Our numerical models have shown that a thermal plume entrains compositionally dense material, forming a compositional plume that is concentric to the thermal plume. The entrainment rate of compositionally dense material by a thermal plume is mainly dependent on the ratio of compositional to thermal buoyancy, Ra_bt, and the radius of the thermal plume, r_T. We have determined how the entrainment rate and radius of the entrained plume vary with Ra_bt and r_T. In this section, we discuss the implications of our results for previous studies on entrainment and for the long-term stability of a layered mantle.

4.1 Comparison With 1-D Entrainment Models

Before discussing their implications for the long-term stability of a compositionally stratified mantle, we compare the entrainment rate determined from our models with that from a 1-D model that is similar to that of Sleep (1988). We have extended this model to satisfy the constraints that the density contrast driving flow and the vertical velocity must have a horizontal average of zero at every depth. Our calculations demonstrate that the flow is nearly 1-D in the vicinity of a plume, except at its two ends (Fig. 4), thus justifying the approximation made by Sleep (1988). However, 1-D models with a prescribed ratio of compositional to thermal buoyancy, Ra_bt, and radius of the thermal plume, r_T, cannot uniquely determine the flow velocity, the radius of entrained compositional plume and the entrainment rate. To determine the entrainment rate, Sleep (1988) assumed that the radius of the compositional plume is such that the entrainment is maximized. It is important to examine to what extent this assumption is valid.

For a 1-D model in which a thermal plume with a radius, r_T, and density anomaly Δρ_T=ρ₀(T_p−T₀)α entrains a compositional plume with a radius, r_c, and density anomaly, Δρ_c (note Ra_bt=Δρ_c/Δρ_T), the non-dimensional averaged vertical velocity and entrainment rate are (Appendix B)  

(9)

and Q=πr_c²u_zc, where δρ_b= (r²_T−Ra_btr²_c)/(1 −r²_T) represents the background buoyancy for r_T≤r≤ 1. It is clear from eq. (9) and that unless Ra_bt, r_T and r_c are all specified, neither the vertical velocity nor the entrainment rate can be determined. Following the approach by Sleep (1988), the radius of the compositional plume r_c is obtained by maximizing the entrainment rate, Q_max, for prescribed Ra_bt and r_T.

The 1-D model predicts that for a fixed r_T, Q_max and r_c decrease with Ra_bt at a slower rate than that given by the 2-D calculations (Figs 5a, b, 6a and b). u_zc from the 1-D model is nearly constant (Fig. 6c), in contrast to the 2-D calculations which clearly show an increase of vertical velocity with Ra_bt (Fig. 5c). Depending on Ra_bt and r_T, the 1-D model can overestimate the entrainment rate by more than a factor of 4, compared with our 2-D calculations (Figs 5a and 6a). This indicates that we should exercise caution when interpreting the maximum entrainment rate predicted from 1-D entrainment models. Although the plume flow is nearly 1-D, the 2-D flow at the two ends of a plume is important in determining the entrainment rate. With numerical models, Tackley (1998) also reported smaller entrainment rates than predicted by Sleep (1988).

Figure 6

Dependence of (a) the maximum entrainment rate, Q, (b) the radius of the compositional plume, r_c, and (c) the averaged vertical velocity for the compositional plume, u_zc, on Ra_bt and r_T from the 1-D analytic models.

4.2 A Modified 1-D Model With Τzz,Z−P,Z Term

One interesting question is what causes the deviation of our 1-D entrainment model from our 2-D numerical models. As discussed in the previous section, the 1-D entrainment model includes two important assumptions (Sleep 1988): (1) r_c is determined by maximizing the entrainment rate and (2) vertical variations of deviatoric stresses and pressure, (τ_zz,z−P_,z), are ignored. Is the deviation of the 1-D entrainment models caused by the neglect of (τ_zz,z−P_,z), or the way in which r_c is determined?

To answer this question, we use the averaged compositional plume radius, r_c, from our 2-D numerical models (Fig. 5b) as the compositional plume radius in 1-D analytic models to determine the averaged vertical velocity and entrainment rate for r_T= 0.05 and different Ra_bt (i.e. we use eq. 9). The predicted entrainment rate and averaged vertical velocity (the dashed lines with squares in Figs 7a and c) are larger than those from our 2-D numerical models, although the rates at which they vary with Ra_bt become similar to those from our numerical models. This indicates that in addition to the determination of r_c, the neglect of (τ_zz,z−P_,z) may also contribute to the deviation of the 1-D models.

Figure 7

Dependence of (a) the entrainment rate, Q, (b) the radius of the compositional plume, r_c, and (c) the averaged vertical velocity for the compositional plume, u_zc, on Ra_bt for r_T= 0.05D and z₀= 0.25 from different models. In Figs 7(a)–(c), the solid lines with open circles and dashed lines with diamonds are from the 1-D analytic and 2-D numerical models, respectively (i.e. the same as those in Figs 5 and 6), the solid lines with inverted triangles are for the modified 1-D model that includes the τ_zz,z−P_,z term. In (a) and (c), the dashed lines with squares and solid lines with triangles are from the 1-D and modified 1-D models with r_c that are determined from the 2-D numerical models.

To examine the effects of (τ_zz,z−P_,z), we have modified the 1-D entrainment model to include this term (Appendix C). (τ_zz,z−P_,z) is significantly smaller than the buoyancy and shear stress terms (Figs 3b and c). However, (τ_zz,z−P_,z) reflects the 2-D nature of the entrainment problem. The 2-D numerical models show that (τ_zz,z−P_,z) may be approximated as an exponential function of the distance from the plume centre, r (Fig. 3c), that is, τ_zz,z−P_,z=−l exp(−mr). For the calculation in Fig. 3(c), l and m are found to be 0.059 and 6.47, respectively. The parameters l and m depend only weakly on Ra_bt. For 1.3 ≤Ra_bt≤ 2.25, l varies from 0.058 to 0.051, and m varies from 6.89 to 5.55.

The modified 1-D model with the (τ_zz,z−P_,z) term cannot uniquely determine the flow velocity, u_zc, the radius of the entrained compositional plume, r_c, and the entrainment rate, Q, similar to our analytic 1-D model. Again we determine r_c by maximizing the entrainment rate (Appendix C). While the predicted r_c and Q from the modified 1-D model are improved, the predicted u_zc deviates more from the 2-D numerical models, compared with the original 1-D model (in Figs 7a–c, the modified 1-D model is shown as solid lines with inverted triangles, the original 1-D model as solid lines with circles and the 2-D numerical model as dashed lines with diamonds). Although u_zc in the modified 1-D model deviates more from the 2-D numerical model, the improved r_c leads to better estimates of entrainment rate in the modified 1-D model (Figs 7a–c). This is because the entrainment rate depends more strongly on r_c than on u_zc. However, because the difference in r_c between the modified 1-D and 2-D numerical models remains significant (Fig. 7b), the entrainment rate from the modified 1-D model deviates substantially from the 2-D numerical models (Fig. 7a). This indicates that the inclusion of either r_c or (τ_zz,z−P_,z) only moderately improves the 1-D entrainment model. However, using r_c determined from the 2-D numerical model in the modified 1-D model results in an entrainment rate and a plume velocity that are nearly identical to those from the 2-D numerical model (solid lines with triangles in Figs 7a and c). This suggests that both the determination of r_c and the inclusion of (τ_zz,z−P_,z) are important for improving the 1-D entrainment model. This also confirms the accuracy of our 2-D numerical models in predicting the entrainment rate.

4.3 Implications for the Long-Term Stability of a Layered Mantle

Whether a layered mantle can survive through Earth's history depends critically on the rate at which the dense material in the lower layer is entrained into the upper layer by thermal plumes (Sleep 1988). Our models determine the entrainment rate under different mantle conditions (Fig. 5a). We find that the entrainment rate is mainly controlled by the ratio of compositional to thermal buoyancy, Ra_bt, and the radius of the thermal plume, r_T, with a functional relationship Q∼Ra^−2.48_btr^3.80_T. Either decreasing r_T or increasing Ra_bt reduces the entrainment rate (Fig. 5a). The entrainment rate is relatively insensitive to the thickness of the dense layer (Fig. 5a).

These results may be used to estimate the survival time for a layered mantle. If we assume that the mantle consists of two compositionally distinct layers, with the dense layer at the bottom entrained by thermal upwelling plumes, then the survival time of this layer is t_s=V_c/Q_t, where V_c is the volume of the dense layer and Q_t is the total entrainment rate. If the characteristic spacing between thermal plumes at mid-mantle depth is the same as the thickness of the mantle, D(2900 km), there would be ∼40 plumes, which is similar to the number of identified hotspots on the Earth's surface (e.g. Sleep 1990). If the bottom layer is ∼1000 km thick, as suggested by recent seismic models (van der Hilst & Karason 1999), then the volume of the bottom layer that is affected by entrainment for each plume is ∼0.27D³. Therefore, the survival time can be expressed as t_s= 0.27Q⁻¹t_scale, where Q is the non-dimensional rate of entrainment by a thermal plume, given by our calculations (Fig. 5a and Table 1), and t_scale=η₀/[(ρ₀gαD(T_p−T₀)]. As an example, let us first examine under what conditions the dense layer can be completely entrained in the age of the Earth. For r_T= 145 km (i.e. 0.05D) and Ra_bt= 1.5, our scaling law shows Q∼ 10⁻⁶ and t_s∼ 2.7 × 10⁵t_scale. For mantle parameters: η₀= 10²¹ Pa s, α= 2.5 × 10⁻⁵ K⁻¹, ρ₀= 3300 kg m⁻³ and T_p−T₀= 800 K, the survival time t_s is ∼4.5 × 10⁹ yr, which is equal to the age of the Earth. For these mantle conditions, with Ra_bt= 1.5, the compositional density anomaly Δρ is 3 per cent of ρ₀ and the net buoyancy (i.e. Δρ−ρ₀α(T_p−T₀)) for the dense material is only 1 per cent of ρ₀.

However, the actual survival time of the dense layer can be significantly longer than the above t_s estimate. The thermal buoyancy varies between plumes, as indicated by the estimated mass flux for different hotspots (e.g. Davies 1988; Sleep 1990). For the Hawaiian plume, possibly the strongest plume on Earth, the estimates for plume excess temperature are between 250 and 400 K (Campbell & Griffiths 1990; Zhong & Watts 2002) and the plume radius is less than 80 km (Zhong & Watts 2002). If we take r_T= 100 km (i.e. 0.0345D) and T_p−T₀= 600 K (i.e. Ra_bt= 2) in our estimate of survival time t_s, while keeping other parameters the same, we found that t_s∼ 5.56 × 10¹⁰ yr, which is more than 12 times of the age of the Earth (i.e. only ∼8 per cent of the denser layer is entrained through the age of the Earth). This is because the entrainment rate is highly sensitive to r_T and Ra_bt (Table 1). We also realize that our entrainment model includes assumptions that may affect the entrainment rates and hence the survival times. Our models assume an isoviscous mantle. A temperature-dependent rheology that would reduce the viscosity for the hot and compositionally dense material may reduce the entrainment, because the reduced viscosity may reduce mechanical coupling at the material interface (Sleep 1988). While future studies are needed to address the effects of temperature-dependent viscosity with dynamically generated thermal plumes, our estimates of the survival time suggest that a dense layer with ∼1 per cent net compositional buoyancy may survive over geological history.

5 Conclusions

The long-term stability of a compositionally layered mantle depends critically on the entrainment of dense material by upwelling thermal plumes that may originate from thermal boundary layers at the material interface. We have formulated high-resolution finite-element models in axisymmetric geometry to examine the efficiency of such an entrainment. We employ a marker chain method to track the evolution of the material interface. Our models show that a thermal plume entrains the dense material, forming a compositionally distinct plume that is concentric to the thermal plume. While the entrainment rate is relatively insensitive to the volume (or thickness) of the dense layer, it is sensitive to two parameters: the radius of the thermal plume, r_T, and the ratio of compositional to thermal buoyancy, Ra_bt. The smaller r_T or the larger Ra_bt, the smaller the entrainment rate, as expected. The entrainment rate scales with r_T and Ra_bt as Ra^−2.48_btr^3.80_T. As Ra_bt increases, the radius of the entrained compositional plume decreases, but the vertical velocity for the compositional plume increases. The functional dependence of the radius of the compositional plume is Ra^−1.43_btr^1.24_T. We have applied these results to study the long-term stability of a layered mantle. For a mantle viscosity of 10²¹ Pa s, thermal plumes with r_T∼ 100 km and temperature 600 K hotter than the ambient mantle and with horizontal spacings that are approximately the same as the mantle thickness, and for a denser layer with a thickness of 1000 km and net negative buoyancy of approximately 1 per cent, our model predicts that more than 90 per cent of the dense material can survive after 4.5 Gyr.

Acknowledgments

We thank P. Tackley for his constructive review. SZ would like to thank M. T. Zuber for her support and J. Huang for assistance with fitting the scaling laws. Some of the calculations were performed while SZ was at the University of Colorado with partial support from the David and Lucile Packard Foundation. Supported by NSF CSEDI grant EAR 9905779.

Appendix

Appendix A: A Finite-Element Analysis And A Marker Chain Method

We briefly describe the numerical methods used in our 2-D axisymmetric models. In axisymmetric cylindrical coordinates, non-zero stress components are  

(A1)

 

(A2)

 

(A3)

 

(A4)

where P is the pressure, u_r and u_z are the radial and vertical velocity components. The elemental stiffness matrix can be expressed as (Hughes 1987)  

(A5)

where e_i and e_j are the unit direction vectors, the superscript T represents the transpose of a matrix, S^e is the area of element e, D_η is a diagonal matrix with diagonal entries (2η, 2η, η, 2η), subscripts a and b are element node indices, and the B matrix can be expressed in terms of a nodal shape function, N_a, and its derivatives with respect to r and z.  

(A6)

The continuity equation in an axisymmetric cylindrical coordinate can be written as  

(A7)

The finite-element representation of this equation is (Hughes 1987),  

(A8)

Now we discuss our marker chain method for solving the advection eq. (4). For tracking a compositional boundary with a relatively simple geometry, the marker chain method is very efficient (e.g. van Keken et al. 1997). Suppose that we have a compositional field, C₀, that is defined by a set of markers with coordinates, xⁱ₀, for marker, i, and velocity, u₀, at t=t₀. Our algorithm for solving C₁ at the next time step, t=t₀+dt=t_1, can be summarized as follows.

• Using a forward Euler scheme, predict the new position for each marker i with xⁱ_1p=xⁱ₀+v₀dt and compositional field C_1p at t=t₁.

• Using the predicted C_1p, solve the Stokes equation for new velocity u_1p.

• Using a modified Euler scheme with second-order accuracy, compute the position for each marker i with xⁱ₁=xⁱ₀+ 0.5(u₀+u_1p) dt and compositional field C₁ at t=t₁. Add or delete markers to ensure that the distance between neighbouring markers satisfies certain criteria.

• Using the new compositional field, C₁, solve the Stokes equations to obtain a new velocity u₁ at t=t₁.

We have two remarks concerning our marker chain method. (1) For an element that is cut through by segments of the marker chain, the element force vector includes the contribution only from the part of element with C= 1. (2) We do not use a fourth-order accurate Runge–Kutta method in our advection scheme, which requires the velocity field at four different times. It is unclear what we would gain by using such a computationally more expensive scheme given the only second-order accurate velocity solver used in our method (Moresi & Gurnis 1996).

Appendix B: 1-D Entrainment Model

We consider a 1-D entrainment model in which a thermal plume with a radius, r_T, and density anomaly, Δρ_T, entrains a compositional plume with a radius, r_c, and density anomaly, Δρ_c (e.g. Sleep 1988). The normalized momentum equation is  

(B1)

where δρ is Ra_bt− 1, − 1 and δρ_b= (r²_T−Ra_btr²_c)/(1 −r²_T), for r≤r_c, r_c≤r≤r_T, and r_T≤r≤ 1, respectively. The buoyancy for r_T≤r≤ 1, δρ_b, is obtained from the mass balance for 0 ≤r≤ 1, and it automatically enforces zero shear stress boundary condition at r= 1. (This source of buoyancy was not included by Sleep (1988), and his model does not satisfy the condition that the horizontal average of the density perturbation must vanish.)

The vertical velocity u is  

(B2)

 

(B3)

 

(B4)

where u₀, u₁ and u₂ are velocities at r= 0, r_c and r_T, respectively. The zero net mass flux from r= 0 to 1 and continuity conditions at r_c and r_T are used to constrain u₀, u₁ and u₂.  

(B5)

 

(B6)

 

(B7)

where in eq. (B7), we have considered r_T≪ 1.

Eqs (B2) and (B5) allow the determination of the averaged vertical velocity in the compositional plume and the entrainment rate:  

(B8)

 

(B9)

Eqs (B8) and (B9) show that the plume velocity and entrainment rate cannot be determined unless Ra_bt, r_T and r_c are all given. However, for typical entrainment problems, r_c is unknown a priori. Sleep (1988) solved for the value of r_c, which leads to a maximum entrainment rate by maximizing eq. (B9) with respect to r_c.

Appendix C: A Modified 1-D Entrainment Model

The modified 1-D entrainment model considers the vertical variations of stress and pressure τ_zz,z−P_,z in the momentum equation. τ_zz,z−P_,z may be approximated as an exponential function of the distance from the plume centre r (Fig. 3c), that is τ_zz,z−P_,z=−l exp(−mr). For 1.3 ≤Ra_bt≤ 2.25, the 2-D numerical calculations suggest that l varies from 0.058 to 0.051 and m varies from 6.89 to 5.55. Similar to the 1-D entrainment model, the momentum equation is  

(C1)

Considering for r≤r_T, mr≪ 1 and, we can express the vertical velocity u as  

(C2)

 

(C3)

 

(C4)

Similar to those in the original 1-D entrainment in, zero net mass flux from r= 0 to 1 and continuity conditions at r_c and r_T are used to constrain u₀, u₁ and u₂.  

(C5)

 

(C6)

 

(C7)

where in eq. (C7), we have considered r_T≪ 1 and m∼ 6.

Similarly, we can determine the averaged vertical velocity in the compositional plume and the entrainment rate. It can be shown that these equations become identical to those in for l= 0.