Effects of melt migration on the dynamics and melt generation of diapirs ascending through asthenosphere

1. Introduction

[2] By crossing the solidus temperature, an upwelling mantle plume forms a partially molten buoyant front within itself and allows for a relative velocity to establish between melt and solid matrix, namely, melt migration. Surrounded by an impermeable viscous solid asthenosphere, melt migrates only within the partially molten part of the plume. The asthenosphere is impermeable because it does not melt or fracture. Farnetani and Richards [1995] showed that only the material within plumes, with isoviscous and mildly temperature-dependent rheology, melts. Melt can leave the front if it can hydraulically induce cracks in its surrounding impermeable rocks. Fowler and Scott [1996] investigated the mechanism of hydraulic crack propagation and concluded that the process is not viable even in a partially molten zone lying beneath the lithosphere. Therefore melt can only reach the base of the lithosphere by the ascent of the melting front because of melt retention buoyancy, chemical buoyancy (arising from density reduction of the solid matrix residue due to melting), and thermal buoyancy.

[3] The significant difference between the density of the solid matrix and melt makes the melt retention buoyancy the largest source of buoyancy even when there is no melt migration. For example, at zero pressure and 10% partial melting, melt retention and chemical changes reduce density by 1.8% and 0.36%, respectively. The density reduction due to the melt retention is equivalent to that produced by a temperature increase of ∼600 K. Therefore it is important to consider melt migration and the time-dependent behavior of melt distribution within a plume.

[4] The problem of melt generation in mantle plumes has been studied by many investigators [e.g., Campbell and Griffiths, 1990; Watson and McKenzie, 1991; Arndt and Christensen, 1992; Farnetani and Richards, 1994; Olson, 1994; Manglik and Christensen, 1997; Schmeling, 2000]. Many aspects such as the chemical buoyancy [Manglik and Christensen, 1997] and melt retention buoyancy [Olson, 1994] have been investigated. Melt is simply removed either instantaneously [Manglik and Christensen, 1997] or after a certain threshold of melting [Olson, 1994]. Schmeling [2000] has recently taken into account melt migration in an upwelling mantle plume that is restricted to occur at melt fraction <5%.

[5] Previous works on melt migration based on two-phase flow models can be grouped into three categories. The first category [e.g., Scott and Stevenson, 1984, 1986; Richter and McKenzie, 1984] was concerned with the existence and properties of solitary waves in the mantle while assuming an initial background permeability. The second category [e.g., Buck and Su, 1989; Scott and Stevenson, 1989; Su and Buck, 1993; Cordery and Phipps Morgan, 1993; Barnouin-Jha and Parmentier, 1997] investigated melt migration beneath mid-ocean ridges while neglecting the compaction of the solid matrix due to melt extraction and ignoring the transient behavior of melt migration within the upper mantle. Ito et al. [1996] and Barnouin-Jha and Parmentier [1997] considered a steady state melt migration in their models of the mid-ocean ridges. Steady state melt migration within an upwelling plume is possible only if melt is extracted or solidified in the top portion of the upwelling plume, at the same rate as the average rate of melt generation within the plume. None of these conditions occur in a partially molten upwelling plume where unsteady behavior of melt distribution is significant. The third category considered melt migration beneath mid-ocean ridges [Ghods and Arkani-Hamed, 2000] and within upper mantle plumes [Schmeling, 2000] while treating the upper mantle as a compacting two-phase flow media.

[6] In the present study, we consider the effect of melt migration on the dynamics, shape, and melt generation of an initially hot and spherical body of 20 km radius ascending in an impermeable asthenosphere using two-dimensional, two-phase flow convection models in an axisymmetric cylindrical coordinate system. Hereinafter the body is called diapir. Because of the intensive computational demand of two-phase flow models the assigned radius of diapir is much smaller than those of mantle plumes (i.e., 200–500 km). Thus the results obtained in this study do not relate directly to mantle plumes but shed light on the fundamental fluid dynamic phenomena that may occur when an upwelling mantle plume enters the partial melting region of the upper mantle. In our modeling, we first focus on the effects of melt migration by considering an initially uniform melt fraction within the diapir and allowing no further melting. The effects of melt retention and chemical and thermal buoyancy on the dynamics of the diapir are studied in detail. We then allow further melting and include all buoyancy sources. These models are compared with equivalent models which ignore differential velocity between melt and solid matrix. We also investigate the effects of melt migration on the shape and partial melting of a typical mantle plume of 200 km diameter. The melt migration produces a zone of high melt fraction at the top of the plume which ascends much faster than the bulk of the plume because of its high buoyancy.

2. Mathematical Formulation

[7] Our treatment of a two-phase flow model is based on McKenzie's [1984] derivation. A two-phase flow model consists of a compacting/decompacting continuous and permeable solid matrix and melt which flows through the matrix. In a given time step we solve equations of conservation of mass, momentum, and energy without considering further melting. At the end of the time step, melt fraction and temperature fields are corrected to include the effect of melting. The normalized governing equations are

The notations are defined in Table 1, and the normalization schemes are listed in Appendix A, where the normalized variables are indicated by asterisks. The asterisks in the above equations are omitted for clarity. The paramenters t, Π, and T are time, reduced pressure, and temperature, respectively. V^s and V_z^s are the solid matrix velocity vector and its vertical component; the z axis is taken downward. V^l, V_z^l, ϕ, ϕ_o, and K_ϕ are the melt velocity vector, its vertical component, the melt fraction, the reference melt fraction, and the permeability, respectively; ρ_s, ρ_s^o, ρ_l, C_P, and K are the density of the solid matrix, the reference density of the solid matrix, the density of the melt, the specific heat, and the thermal conductivity, respectively; δ_i3 is the Kronecker delta function. The stress tensor for the solid matrix σ is defined as

where η_s is the shear viscosity of the solid matrix, ζ is the bulk viscosity of the two-phase flow media associated with the volumetric deformation due to melt migration (i.e., compaction/decompaction of the solid matrix), and I is the unit matrix. In contrast to the considerable volumetric changes due to melting, each of the melt and solid phases has negligible volumetric changes. Therefore ζ is set to zero outside the molten zone. Equations (1) and (2) express the conservation of the solid matrix and melt, respectively, where both densities of the melt and solid are assumed to be constant. The density of the solid matrix is not constant in the momentum equation (equation (3)). Equation (4) is the momentum equation for melt which is equivalent to Darcy's law. Pressure is defined as P = Π + ρ_s^ogδ_cδ_i3, where g is the gravitational acceleration. The compaction length δ_c is the characteristic length over which compaction occurs. It is defined as

where K_ϕ^o and η_l are the reference permeability of the solid matrix and the viscosity of the melt, respectively. In the momentum equations the nondimensional number R_δ = ρ_s^ogδ_c²/η_sW_o is similar to the conventional Rayleigh number in that it describes the relative strength of mantle buoyancy versus the viscous forces [Cordery and Phipps Morgan, 1993]. W_o is the characteristic velocity of the system. The second nondimensional number A = η_s/[ϕ_o (ζ + 4η_s/3)] denotes the efficiency of compaction and melt migration. Equation (5) is the conservation of energy. The nondimensional Peclet number Pe = W_oδ_c/κ_s^o describes the ratio of the rate of thermal advection to that of thermal diffusion where κ_s^o is the thermal diffusivity of the melt and the solid matrix. The nondimensional number Di = gαδ_c/C_P is the dissipation number quantifying the importance of the adiabatic cooling/heating, where α is the thermal expansion coefficient of the solid matrix. The addition of adiabatic cooling/heating has minor effects on the results. We therefore did not further consider adding a shear heating effect since both of these effects scale similarly with Di. In writing the energy equation it is implicitly assumed that thermal equilibrium exists locally between the melt and solid matrix. No latent heat of melting, internal heating, or shear heating are considered in the energy equation. As we explain in section 3, at the end of each time step the temperature field is corrected to take into account the latent heat of melting. We also assume that the thermal conductivity of the melt and solid matrix are equal and ρ_lC_P^l = ρ_sC_P^s, where C_P^s and C_P^l are the specific heat at constant pressure of the solid matrix and melt, respectively.

[8] Melting is a chemical reaction which depletes iron content and reduces the density of the solid residue. Therefore the density of the solid matrix is a function of temperature and the depletion induced by melting:

where β is the melt depletion coefficient and ℳ is the fraction of melting. For example, ℳ = 0.1 denotes a solid residue of 10% melting. T_o is the reference temperature which is assumed to be the average of the temperature at the top and bottom of the computational domain. In the impermeable models, where no relative velocity is allowed between the melt and solid matrix, ℳ is equal to ϕ, whereas in the permeable models, ℳ is smaller than ϕ in regions of melt accumulation and larger than ϕ in regions of melt extraction.

[9] Following Richter and McKenzie [1984], permeability is considered to be a function of melt fraction:

where the grain size a is set to 1 mm. Equation (9) implicitly assumes that the porosity network is interconnected at all melt fractions which is in general agreement with experimental results [Kohlstedt, 1992].

3. Melt Generation

[10] There are two end-member melting models: batch melting where melt is retained with the solid matrix and fractional melting where melt is completely extracted. We consider both models in order to investigate effects of the type of melting on the dynamics of the impermeable upwelling diapirs. Because of melt migration, melting of the permeable diapirs is considered to be fractional.

[11] To include the clinopyroxene out criterion in the melting, a maximum melting of 25% is considered. At this melting stage the low melting temperature components of mantle rocks are lost, the melting temperature of the residual solid sharply increases, and melting effectively ceases [Cordery and Phipps Morgan, 1993].

[12] McKenzie and Bickle's [1988] batch melting model is used to calculate batch melting of the impermeable diapir. The model assumes a linear relationship between the temperature at a given depth and the fraction of melt:

where T_s^o and T_l^o are the solidus and liquidus temperatures, respectively. The solidus temperature of fertile peridotite, T_s^o, is obtained by fitting a third-order polynomial to the solution of McKenzie and Bickle [1988],

where P_H is the hydrostatic pressure. Changes in P_H due to melting are minute and are neglected in our calculations. The liquidus temperature, T_l^o, is calculated as [McKenzie and Bickle, 1988]

Solidus and liquidus temperatures in (11) and (12) are in kelvins.

[13] In permeable models, melting is dynamic so that some portion of the melt remains with the solid matrix, some portion may migrate out, or some extra melt may be added to the partially molten zone by melt migration from elsewhere. In this paper we model the dynamical melting of the permeable diapir using the solidus and liquidus curves of a fractional melting model and assume that the addition of melt to a partially molten zone by melt migration does not change the solidus temperature of the solid residue. Our fractional melting model is a simplified version of Iwamori et al.'s [1995] model and is constructed by modifying McKenzie and Bickle's [1988] batch melting model. The model is also used to calculate the volume of melt produced by the impermeable diapir through fractional melting. Melt extraction in the fractional melting models is simulated simply by setting melt fraction remaining within the solid matrix (ϕ) to zero.

[14] In the fractional melting model we replace the liquidus temperature of peridotite (equation (12)) with that of pure forsterite (i.e., T_l^o = 2164 + 4.77 P_H [Presnall and Walter, 1993]). Furthermore, to account for the dynamic nature of the problem, the ℳ field is advected in each time step using the following pure advective transport equation:

To simplify the calculations, we used bulk velocity, rather than the velocity of solid matrix, for advecting the depletion field. The bulk velocity may not substantially differ from the solid matrix velocity. This is because of smaller percentage of melt in the regions where melt migrates fast. In the high percentage melt region at the top of the plume the velocities of the melt and solid matrix are similar since no percolation of melt occurs in to the unmolten region at the top. The fraction of melt produced in a given time step, Δℳ, is computed by the following coupled equations:

where T₁ and T₂ are temperatures before and after correction, respectively, and the latent heat of melting H is considered to be constant. T_s is the instantaneous solidus temperature of the solid residue of the previous time step which is calculated using (10) and setting T_s = T to account for the increase of T_s through melting. The removal of the melt in a fractional melting model increases T_s because the solid residue becomes more refractory. The temperature field is corrected for the absorption of the latent heat of melting and Δℳ is added to the ϕ and ℳ fields.

[15] Melt fraction is not allowed to become larger than 0.25, and accordingly, the permeability of molten rocks is limited to that corresponding to a 25% molten rock. This is because the momentum equations (equations (3) and (4)) are only valid when melt flows within a continuous permeable solid matrix [McKenzie, 1984]. The equations are not valid [e.g., Spiegelman, 1993a; McKenzie, 1984] for melt fractions larger than the value at which the matrix disaggregates and forms a slurry-type two-phase fluid media, i.e., at ϕ = 0.2 ± 0.1 [Arzi, 1978]. We recognize the possibility that accumulation of melt at the top of a mantle plume may result in a very low viscosity slurry-type two-phase media. However, developing momentum equations capable of describing both continuous porous and slurry-type two-phase fluid media is a challenging task which we do not attempt to investigate in this paper.

4. Numerical Method

[16] An iterative finite volume method [Patankar, 1980; Prakash and Patankar, 1985] is adopted to solve -–(5). Discretization of momentum equations (equations (3) and (4)) are second order accurate. A hybrid scheme [Patankar, 1980] is used to solve the energy equation, and a second-order Runge-Kutta method is adopted for discretization of the time derivatives in the equations of conservation of melt (equation (2)) and energy (equation (5)). To calculate the pressure field, (1) and (2) are combined and converted to a pressure equation using (3) and (4), which is then solved using the semi-implicit method for pressure-linked equations, revised (SIMPLER) [Prakash and Patankar, 1985]. This eliminates the time derivative terms and yields a pressure equation with no explicit time dependence. Unless it is mentioned otherwise, the numerical results presented in this paper are calculated in a nonuniform rectangular grid of 71 × 201 nodes with a high-resolution area fully covering the ascending diapir path (Figure 1). The high-resolution area covers the entire computational domain along vertical axis and within 25 km from the diapir axis. The vertical grid spacing is constant at ∼0.6 km. The horizontal grid spacing in the high-resolution area is ∼0.56 km and gradually increases outward to a maximum of ∼6 km. The variable grid reasonably resolves compaction lengths larger than 1.8 km, which is always the case for ϕ > 0.048.

[17] Entrapment of melt within the diapir produces a melt fraction field with a sharp boundary between the diapir head and the surroundings, which necessitates a special treatment of equations (2) and (13). We use a Monotonic Second Order Upwind (MSOU) [Sweby, 1984] method to advect melt depletion and melt fraction fields. MSOU is one of the most accurate Eulerian-based methods [Tamamidis and Assanis, 1993]. To suppress the numerical leakage of melt fraction outside the plume, the plume shape is traced throughout its evolution by initially dying the plume with a constant concentration and then tracing it with time. This allows us to almost entirely prevent the numerical leakage of melt outside of the plume by bringing the numerically leaked melt inside the plume at each time step. Because of our very fine mesh, and thus very small time steps, and the adopted high-accuracy numerical method of MSOU, the numerical leakage is not significant in our models. The melt is always within the dyed area, indicating that it does not infiltrate into the unmolten region.

[18] To ensure convergence, four parameters are computed at each time step: (1) the absolute value of the total vertical advective heat transfer, (2) the sum of the flux of the two phases at the normalized depth of 0.5, (3) the Nusselt number at the top left corner, and (4) the Nusselt number at the bottom right corner. The convergence of Nusselt number is usually slower at the corners than any where else on the border of the computation domain. A solution is assumed converged if the ratios of the absolute value of the difference between the previous and current values of the above parameters to their current values become smaller than our assigned convergence accuracy limit of ϵ = 5 × 10⁻⁵. In addition, we check the sum of absolute values of errors in the solution of the equation of conservation of mass over the entire computational domain. This is a global constraint ensuring that mass is conserved throughout the computation domain. Because of the absolute values, even small local errors add up and can be detected by this global constraint. The residual errors for different fields are always less than the assigned convergence accuracy limit. The time step is assigned to be one fifth of that prescribed by the Courant criterion [e.g., Wendt, 1996].

[19] To test the fidelity of our numerical schemes, we use an equivalent of our code in two-dimensional (2-D) Cartesian coordinate system and calculate the shape of a single solitary wave in one dimension. The width of the solitary wave is resolved with 40 grid spacings and is advected over a distance 5 times larger than its width. Solitary waves should propagate at a constant velocity without changing shape [e.g., Spiegelman, 1993a]. The theoretical amplitude [e.g., Spiegelman, 1993a] is reproduced within ∼0.1%.

5. Diapir Dynamics

[20] We model the ascent of a partially molten diapir through the upper 120 km of an impermeable asthenosphere in order to investigate the effect of melt migration on the dynamics of upwelling diapirs. Although our diapir size is much smaller than those suggested for hot mantle plumes, we believe that it illustrates the basic physics involved with an affordable computational demand.

[21] We consider an initially spherical diapir with a radius of 20 km (Figure 1). To minimize the effect of the zero velocity bottom boundary condition, the base of the diapir is initially placed 20 km above the bottom boundary. A free-slip boundary condition is applied for both melt and solid matrix velocities at the left, right, and bottom boundaries. The velocities at the top boundary are set to zero. Melt fraction is set to zero at the top and bottom boundaries. A zero flux is specified for melt fraction and heat at the vertical side boundaries. The temperature at the top boundary is 1673 K, and on the basis of an adiabatic temperature gradient the temperature at the bottom boundary is set to 1723 K. The top boundary of the computational domain is assumed rigid, representing the base of an impermeable rigid, stagnant lithosphere of 85 km thickness. The diapir is initially 300 K hotter than its surroundings and has 10% melt.

[22] Table 2 lists the values of all parameters used in this paper. The melt density is taken to be constant at 2700 kg m⁻³, and the reference solid matrix density (ρ_s^o) is assumed to be 3300 kg m⁻³. A constant melt viscosity η_l of 1 Pa s is chosen, which is the lower limit of the viscosity range of basaltic melt (1–10 Pa s) [Kushiro, 1986]. Constant and equal shear (η_s) and bulk (ζ) viscosities of 1.29 × 10¹⁹ Pa s are considered for the solid matrix. The compaction length is 5.4 km for ϕ = 0.1. Our current grid spacing of 0.6 km resolves compaction length larger than 1.8 km (equation (7)), implying that melt fraction >0.048 is properly resolved.

[23] In the first set of models an initially uniform 10% melt distribution is assumed, and no further melting is al lowed. To illustrate the effects of melt migration, the dynamics of a permeable diapir in which melt is allowed to migrate is compared to that of an impermeable diapir in which melt and solid matrix move together with no differential velocity between them. In addition, the effects of the thermal and chemical buoyancy forces are investigated. In the second set of models, further melting is allowed. The permeable models are again compared with the corresponding impermeable models in both fractional and batch melting cases.

5.1. Nonmelting Diapirs

[24] In the first set of numerical experiments, dynamics of nonmelting impermeable and permeable diapirs are modeled while considering only melt buoyancy. In this experiment an initially constant 10% melt is assumed, and no further melting is allowed. This singles out the effect of melt migration on the dynamics of the upwelling diapir. Figure 2 shows the snapshots of melt fraction field for the permeable model. In the early stages the melt fraction increases in the top portion of the permeable diapir and decreases in the lower regions (Figure 2a), producing a high melt fraction cap at ∼0.008 m.y. (Figure 2b). Because of the nonlinear relationship between permeability and melt fraction (equation (9)), the compaction/decompaction of the partially molten diapir develops a low melt fraction zone beneath the high melt fraction cap, initiating 2-D solitary waves (Figure 2c). Melt migration from deeper parts increases the melt fraction beneath the low melt fraction zone and results in an increase of permeability there. This causes further draining of melt from regions beneath and creates a new local minimum melt fraction, where the new minimum melt fraction region obstructs the upward melt migration and initiates a second maximum and so on [Spiegelman, 1993b]. The upwelling velocity of the diapir decreases sharply in the vicinity of the upper boundary. This causes the solitary waves to coalesce and join the high melt fraction cap at the top.

[25] Figure 2b shows that melt is initially drained more along the axis. We tested this issue further by examining in detail melt distribution, temperature, velocity, and pressure fields along the axis. The axial disturbance is only seen on melt distribution field and at the very early stages of diapir ascending but not on other fields. We have checked the implementation of the boundary condition in our computer code and found no problem. We also ran the same model in Cartesian coordinates and did not find the along-axis disturbance. The axial disturbance in melt distribution field arises because of the cylindrical coordinate system. The curved boundary of the diapir at the top focuses melt toward the axis. In the case of Cartesian coordinate system (x, y, z) the diapir is actually an infinite cylinder along the y axis with a circular cross section in the (x, z) plane, and focusing is parallel to the x axis. In cylindrical coordinate system (r, θ, z), with a vertical z axis, the focusing is radial and moves all of the melt within the annulus along θ toward the axis.

[26] Comparing the shapes of impermeable and permeable diapirs (Figure 3a) shows that the general shape of the permeable diapir is cylindrical as compared to the mushroom shape of the impermeable diapir. For example, when both diapirs reach a depth of ∼15 km, the width of the permeable diapir is less than one half of the impermeable diapir. Once the melt percentage reaches a maximum of 25%, no further buoyancy polarization occurs, and the very top portion of the cylindrical diapir starts to become mushroom-shaped. Figure 3b compares temperature fields corresponding to Figure 3a. The impermeable diapir is more diffused, and its mean temperature is less than the permeable diapir by ∼100 K. Figure 4 shows the time evolution of the depth to the top of the diapir along its axis. The permeable diapir ascends ∼4 times faster than the impermeable diapir to reach at a depth of 10 km. The faster permeable diapir impinges the top boundary (the base of the lithosphere) within <0.4 m.y., whereas the impermeable diapir may never reach it. These differences are related to the polarized melt fraction distribution of the permeable model due to melt migration (see Figure 2), which produces a polarized melt retention buoyancy and leads to a more slender diapir that suffers less drag while ascending. The fast moving, slender diapir diffuses less thermal energy to its surroundings and is warmer than the impermeable model upon impinging the top boundary. This point is also reflected in the time evolution of the mean Nusselt number at the top (Figure 5). The maximum Nusselt number of the permeable diapir is much larger than that of the impermeable diapir and is attained much earlier.

[27] Figure 6 compares the time evolution of the depth to the top of the diapir along its axis for the impermeable and permeable models where thermal and chemical buoyancies are included in addition to the melt retention buoyancy. The addition of thermal and chemical buoyancies slightly increases the velocity of the diapirs. However, the velocity increase is much less pronounced for the permeable diapir than the impermeable one. This is related partly to the already larger velocity of the permeable diapir and partly to the initially uniform distribution of melt fraction and temperature within the diapir. The initial values result in almost uniform thermal and chemical buoyancies within the diapir, which reduce the ratio of the melt retention buoyancy force to the total buoyancy force and hampers the elongation of the diapir. A larger drag force is exerted on the less elongated diapir which counterbalances the buoyancy increase due to the addition of the chemical and thermal buoyancies.

5.2. Melting Diapirs

[28] In the second set of models we allow melting to take place. The initial diapir has zero melt fraction but 300 K anomalous temperature. The diapir melts as it rises through the asthenosphere. We do not restrict melting to the diapir, but it only occurs within the diapir. This is in agreement with Farnetani and Richard's [1995] conclusion that the only material that melts is within diapirs with isoviscous and mildly temperature-dependent rheology. All buoyancy forces, thermal, chemical, and melt retention, are considered in this set of models. In the first two experiments the dynamics of an upwelling impermeable diapir is investigated for cases in which the melt either is completely removed (fractional melting) or is retained within the matrix (batch melting). In the third experiment the melt is allowed to migrate within the partially molten permeable zone.

[29] Comparison of the depletion (ℳ) and melt fraction (ϕ) fields of the permeable model (Figures 7a and 7b show that the regions of melt generation are different from the melt fraction field. Depletion is maximum in the middle portion of the diapir and decreases monotonically with depth. However, the melt fraction is concentrated in three patches which are separated by low melt fraction zones. Figures 7a and 7b also indicate necking of the diapir, which occurs because both melt retention buoyancy and chemical buoyancy are significantly larger at the top portion of the diapir where the degree of melting is at maximum. The necking produces three separate regions of high melt fraction, one above the necked region, one inside and one below (Figure 7a). The latter regions form due to melting of the lower portion of the diapir at later times as the diapir rises. To illustrate the effect of melt migration on the shape of the diapir, the depletion and temperature fields of the permeable and impermeable diapirs with batch melting are compared when the diapir heads are at ∼10 km depth (Figure 8). Because the velocities of the upwelling diapirs are not equal, the times for the snapshots are different. The snapshots are taken at 2 and 0.56 m.y. for the impermeable and permeable models, respectively. The permeable diapir is elongated and lacks the mushroom shape associated with the impermeable diapir. Figure 9 compares the time evolution of the depth to the top of the diapir along its axis. The permeable diapir reaches at 10 km depth ∼4 times faster than the batch melting impermeable diapir. Furthermore, the permeable diapir results in the largest Nusselt number (Figure 10a) and produces the largest volume of melt (Figure 10b). This is because the thermal energy of the permeable diapir dissipates less due to its fast upward motion and the diapir remains hot and melts more extensively while ascending through the melting window of the upper mantle. The diapir with no melt retention buoyancy (fractional melting) is the slowest and produces the lowest volume of melt.

6. Discussion

[30] Our study indicates that melt migration causes the upwelling velocity of the melting diapir and its melt production to increase significantly. Furthermore, the initially spherical diapir evolves to an elongated cylindrical body. The melt production of our diapir model (i.e., <2%) is analogous to the possible melt production of upper mantle plumes responsible for continental flood basalts [Farnetani and Richards, 1994]. However, care should be exercised in directly extrapolating our results to the mantle plumes. One of the most important parameter is the size of our diapir (i.e., 20 km) which is considerably smaller than those of mantle plumes (i.e., 200–500 km). Our diapir models fully enter into the melting zone of the upper mantle, whereas the mantle plumes only partially enter due to their larger sizes, and they are already mushroom-shaped before their entrance. However, one can imagine that melt migrates up the slope of the upper boundary of the mushroom shape plume and forms a high melt fraction region at the vicinity of the plume axis. We investigate this scenario by considering an initially spherical plume model of 200 km diameter. The top of the plume is initially at 100 km beneath the base of the lithosphere (our top boundary). All initial values and boundary conditions are exactly similar to the rest of our models. The computation domain is increased to 600 × 600 km and has a nonuniform grid of 201 × 201 nodes with a high-resolution area in the top 80 km of the domain and within 125 km from the plume axis, fully covering the partially molten zone in the upper mantle. The vertical and horizontal grid spacings in the high-resolution area are ∼0.65 and ∼1 km, respectively, and gradually increases outward to a maximum of ∼5 km. Figures 11a–11c show three snapshots of the plume shape, and Figures 11d–11f show the corresponding melt fraction fields, for a permeable model where all buoyancy forces are considered. The melt fraction fields show only the upper part of the plume. A high melt fraction region forms at the top of the plume (Figure 11a) which with time becomes elongated and forms a nibble (Figures 11b and 11c). The nibble moves faster than the rest of the plume and delivers melt to the base of the lithosphere before the main body of the plume impinges there.

[31] The high melt percentages reported in our paper have not been observed by seismic investigations in Yellowstone [Saltzer and Humphreys, 1997] and Iceland [Wolfe et al., 1997]. One possible explanation for the low melt content of Yellowstone and Iceland plumes is that they have already impinged to the base of the lithosphere and their high percentage melt areas have been consumed through outpouring. This is similar to the case of mid-ocean ridge where the permeable zone extends to the base of the crust. The two-phase flow model of Ghods and Arkani-Hamed [2000] showed that the fast melt migration in the crust beneath mid-ocean ridge does not allow the melt fraction in the upper mantle to exceed 2–3%.

[32] Volcanism predates rifting in many continental rifting associated with flood basalts [Farnetani and Richards, 1994] such as Parana [Piccirillo et al., 1988] and Kenya Rift [Smith, 1993; Morley, 1993]. Our model is consistent with this observation and shows that melt delivery to the base of the lithosphere can occur before the arrival of the main body of the plume that exerts horizontal stress to the overlying lithosphere and causes rifting [Farnetani and Richards, 1994].

[33] The melt density is assumed to be constant in our models despite an appreciable increase of the density with pressure. The melt density of an equivalent basaltic magma devoid of iron and alkalis changes from its surface value of ∼2700 to ∼3000 kg m⁻³ at a depth of 100 km [Rigden et al., 1984]. This reduces melt retention buoyancy and hampers upward migration of melt from deeper parts of the diapir. Consequently, the formation of the compacting and decompacting regions at the bottom and top of the diapir is delayed. One can estimate from Darcy's law (equation (4)) that a reduction of melt buoyancy by half doubles the time required for the formation of the high melt fraction cap. Nevertheless, the formation time of the polarized buoyancy distribution is much shorter (∼0.01 m.y.) than the time taken for upwelling of the diapir through asthenosphere (∼0.3 m.y.). The upwelling velocities of the impermeable and permeable diapirs reduce similarly with the reduction of the total melt buoyancy. Therefore the increase of melt density with depth may not substantially change the ratio of these velocities. Results of models with 1/2 and 1/4 of density contrast used in our previous models confirm the above statements.

[34] In this study, we considered a simple permeability law (equation (9)). However, permeability may become significant only after a critical threshold of 2–3% of melt fraction [Faul, 1997]. Imposing such a threshold results in an initially slower formation of the buoyancy polarization within the diapir and thus slower speedup of the permeable diapir. Melt migration becomes significantly fast once the melt fraction exceeds the threshold. However, the time required for formation of the high melt fraction cap is very short in comparison to the time required for the upwelling of the diapir. Thus the shape and dynamics of diapir would not be significantly different from those presented in our paper. Melt/rock reaction also increases permeability [Aharonov et al., 1995]. Further investigations are needed to study effects of melt/rock interaction and a permeability threshold on the dynamics of the upwelling diapirs.

[35] Recent studies suggest that upwelling plumes may have ∼15% eclogite. Cordery et al. [1997] showed that such plumes can produce the volume and composition of melt observed in flood basalt provinces without having a high initial temperature. In our models the peridotite solidus and liquidus curves are adopted because we are primarily concerned with the effect of melt migration on the dynamics of the upwelling diapirs. The general conclusions drawn in this paper, the polarized structure of the melt distribution and the elongation of the upwelling diapirs, are expected to hold even if the melting models of Cordery et al. [1997] are used.

[36] To address the effect of the grid resolution on the dynamics of our diapir models, we examine three different grids: 51 × 101, 71 × 201, and 141 × 401. Our adopted 71 × 201 grid resolves compaction length longer than 1.8 km, which corresponds to melt fractions larger than 4.8%. Figure 12 shows the time evolution of the depth to the top of the impermeable and permeable diapir models calculated using these grids. The models show that by increasing the resolution the difference between the arrival times of the impermeable and permeable diapirs to a depth of ∼10 km increases. Furthermore, a higher resolution delays the arrival of the diapirs to the depth of 10 km, especially that of the impermeable one. However, the grid dependency of our models is not vital and does not significantly affect our major conclusions about the effects of melt migration on the dynamics of the upwelling diapirs.

[37] The solid matrix shear viscosity is an important parameter affecting its compaction and decompaction and thus the dynamics of the upwelling diapir. It is a function of melt fraction, permeability, temperature, and strain rate, among other factors. The pressure effect is insignificant for the depths considered in our models. It is expected that because of a positive feedback mechanism between the viscosity of the solid matrix and the increase of the strain rate, a strain rate dependent rheology would have a significant effect on the difference between the permeable and impermeable models. However, the dynamics of the upwelling diapir may not significantly change since the compacting and decompacting regions develop fast. The solid matrix viscosity substantially reduces after a critical value of 4% melt [Hirth and Kohlstedt, 1995]. Khodakovskii et al. [1998] showed that the reduction due to melt exceeding a threshold of 5% causes an increase in the growth rate of solitary waves and a decrease in their length scale. The solid matrix viscosity has been taken to be constant in the above models. To study the effects of the solid matrix viscosity reduction due to melt fraction, we calculated permeable and impermeable diapir models where the solid matrix viscosity is linearly reduced by a maximum of 1 order of magnitude as the met fraction increased from zero to a maximum of 25%. Figure 13 compares the time evolution of the depth to the top of the diapir along its axis for permeable diapir models with constant and melt-dependent variable solid matrix viscosity. Both models are computed on a higher-resolution nonuniform computational grid of 141 × 401 nodes. The diapir with melt-dependent variable viscosity ascends slightly faster than the diapir with constant viscosity because the viscosity reduction by melt fraction allows the portion of the diapir that enters the upper mantle melting zone to easily move upward and to enhance the buoyancy polarization and to render the diapir more slender. This result further emphasizes the main conclusions reached in this paper. The difference in upwelling velocity between variable and constant viscosity impermeable diapir models is even less important. We would like to mention that viscosity reduction due to melt fraction may be several orders of magnitude. Our 1 order of magnitude reduction is more likely a lower limit. However, we believe that simple reduction of the solid matrix viscosity enhances the processes observed in our previous models but may not result in new physical phenomena not observed in the present models. It is required to formulate the dynamics of a diapir, or a mantle plume, by a new set of equations that can simultaneously describe the dynamics of a partially molten porous medium for the region of lower melt fractions and a slurry medium for the zones with higher melt fractions. This is certainly a demanding task as described in some detail by Bergantz [1995], and we are presently investigating this scenario.

[38] Melt depletion increases the viscosity of the solid matrix [Hirth and Kohlstedt, 1996]. The large melt depletion and the low melt fraction in the middle part of the permeable diapir (Figures 7 and 8) make it more viscous relative to the top portion of the diapir where melt fraction is very high. This may hamper the upwelling of the more viscous, depleted material of the diapir. We have not yet taken into account the effect of melt depletion on the solid matrix viscosity in our models. Nor have we incorporated the effect of temperature on the viscosity. These two effects tend to cancel each other in the middle region where both depletion and temperature are high. Also, a constant bulk viscosity of solid matrix is assumed in our calculations, although it may vary with melt fraction [Schmeling, 2000].

7. Conclusions

[39] The effects of melt migration on the dynamics of a 20-km-radius diapir are investigated using two-phase flow model in a 2-D axisymmetric cylindrical coordinate system. In the first set of models we single out the effects of melt migration by modeling permeable and impermeable diapirs with a 10% melt fraction uniformly distributed initially and allowing no further melting. Melt migration to the top of the partially molten region and the corresponding melt retention buoyancy makes the permeable diapir more slender than the impermeable diapir. Although both models have equal mean buoyancy, the more elongated permeable diapir suffers less drag and ascends faster. In the second set of models the dynamics of the permeable and impermeable diapirs are investigated while melting is allowed. The permeable diapir moves faster as soon as it enters the partial melting window of the upper mantle because of the highly polarized melt retention buoyancy and its more slender shape. The fast moving diapir diffuses less thermal energy to its surrounding and generates more melt during its ascent through the upper mantle. The permeable models lack the mushroom shape observed for the impermeable models. The solid matrix viscosity reduction due to melt fraction further facilitates the processes but does not result in new physical phenomena. We also model the dynamics of a mantle plume of 200 km diameter and show that melt migration produces an axial high melt fraction nibble in the upper part of the plume that ascends much faster than the bulk of the plume.