A model of episodic melt extraction for plumes
H. Schmeling1
Received 6 September 2004; revised 10 November 2005; accepted 20 December 2005; published 23 March 2006.
[1] A model of melt segregation and extraction within rising plumes is proposed. It is based on two-phase porous flow within the partially molten region, combined with only three extraction parameters. The conservation equations of mass, momentum, and energy are solved for a two-phase melt matrix system. As a rising hot mantle region (plume) reaches the asthenosphere, decompression melting occurs, and the melt begins to percolate with respect to the matrix. Accumulation layers form, which might be the locus for the formation of buoyancy-driven propagating dikes. As dike propagation requires a
minimum dike length, melt extraction is parameterized by dex,j1, andj2. Here dexis the critical thickness of the partially molten layer in which a critical melt fraction j2 is exceeded. If this condition is met within a certain region of the melt source region, melt might be extracted from that region in the form of one or several propagating dikes, leaving behind a region of residual melt fractionj1. This simple extraction model is tested in one dimension for rising hot mantle flow. Depending on the chosen extraction parameters, multiple extraction events may be observed with a characteristic episodicity and a saw-tooth-like depth distribution. Exploring the parameter space shows that for values of dexandj2 of a few kilometers and a few melt percent, respectively, typical extraction cycles have the order of 103–104 years, and they extract melt volumes per surface area of 50 to several hundred meters each. Tentatively assuming that eruptions are tied to mantle ascent at depth, the model is applied to observed eruption frequencies and multiple extraction depths, and values forj2of about 2% and dexof 3– 5 km are derived.
Citation: Schmeling, H. (2006), A model of episodic melt extraction for plumes,J. Geophys. Res.,111, B03202, doi:10.1029/2004JB003423.
1. Introduction
[2] One of the major outstanding problems of melting processes in Earth’s mantle is developing an understanding of melt segregation and extraction. Of particular importance for melt extraction is the critical melt fraction (melt porosity) which must be exceeded to allow extraction and rapid melt ascent. From the modeling of major element chemistry of dry MORB generation, Ryabchikov [1998] inferred critical melt fractions of 6%, similar to the values proposed by Nicolas [1990] on the basis of obser- vations of the fraction of intergranular material in orogenic peridotites. Lower critical melt fractions of about 1% have been inferred for plume generated melts from greater depths [Ryabchikov, 1998].
[3] Disequilibria between short-lived U series isotopes have been used to estimate critical melt fractions and ascent velocities [e.g., Spiegelman and Elliot, 1993; Richardson and McKenzie, 1994; Sims et al., 1995, 1999; McKenzie, 2000]. These approaches invoke different models of melting and extraction within a rising column undergoing decom- pressional melting. They include models of instantaneous
melting (instantaneous melt generation and extraction), dynamic melting (continuous melt generation, but extrac- tion after reaching a critical melt porosity), and equilibrium percolation (percolating melt in equilibrium with solid before extraction; seeElliot[1997] for a discussion of these models). Because of such different model assumptions and because of the principal uncertainty of the values of the partition coefficients the inferred critical melt fractions may differ by orders of magnitudes. For example, Sims et al.
[1999] obtain maximum porosities of the melting zones of Hawaiian tholeiites between 0.3 and 3%,McKenzie [2000]
estimates Hawaiian melt porosities to be of the order of only 0.2%, and Landwehr et al. [2001] favor values of 0.001 to 0.1%, on the basis of new measurements of partition coefficients.
[4] Threshold melt-filled porosities have also been esti- mated on the basis of experimental results on melt distribution in partially molten olivine-basalt aggregates combined with percolation theory [Faul, 1997, 2001]. He argues that at the onset of melting volatile-rich low- viscosity melts are able to ascend at the small porosities of 0.1% to produce the observed U series isotopic dis- equilibria. At later melting stages the threshold porosities for basaltic melts are as high as 3%.
[5] Composition also influences threshold porosities. Be- cause of higher dihedral angles between melt and pyroxene grains, the presence of pyroxene may decrease the perme- ability and increase the threshold porosity [Toramaru and
1Geophysics Section, Department of Earth Sciences, Johann Wolfgang Goethe Universita¨t Frankfurt am Main, Frankfurt, Germany.
Copyright 2006 by the American Geophysical Union.
0148-0227/06/2004JB003423$09.00
Fujii, 1986]. This effect is reduced in the presence of water [Fujii et al., 1986], strengthening the argument that volatile- rich melts may escape at smaller melt fractions [Faul, 2001].
[6] The physical processes of melt extraction associated with the above mentioned threshold melt fraction are not well understood. As melting will be associated with weak rheology [e.g., Hirth and Kohlstedt, 1995a, 1995b] compaction waves may be triggered [e.g., Wiggins and Spiegelman, 1995] in which a complete melt-solid separation may occur [Rabinowicz et al., 2001]. While in such models minimum melt fractions as high as 5% are needed to trigger complete segregation, other models assume that as melt is generated continuously within the melting zones of plumes, it accumulates within planar or cylindrical zones, forming a drainage system of high- porosity channels [e.g., Hart, 1993]. As the melt is expected to be in thermal but not chemical equilibrium with the rock matrix, the spacing of such channels is estimated to be above 0.1m [Iwamori, 1993, Maaløe, 2003]. One key questions is: at which melt fraction do such channels form? There are essentially two mechanisms which may form melt channels within partially molten regions: two-phase reactive flows [Aharonov et al., 1995;
Kelemen et al., 1995] and two-phase flows under devia- toric stresses with porosity weakening [Stevenson, 1989].
Spiegelman et al. [2001] solved the full equations for two- phase reactive flow with compaction and obtained channel- ing at background porosities as small as 0.1 to 1%. If a partially molten rock is subjected to deviatoric stresses, porosity mediated weakening may lead to a channeling instability with channels oriented perpendicular to the direc- tion of minimum compressive stress [Stevenson, 1989]. To obtain melt channels at moderate finite strains, higher initial melt factions (1 – 3%) are required [Richardson, 1998;
Mu¨ller, 2005].
[7] Once melt channels have formed, they might provide a pathway of melt within and out of the partially molten zones. The question is whether further ascent by dikes occurs already within the partially molten zone or immedi- ately above it. WhileNicolas[1986] argued that fracturation occurs within the partially molten zone this has been questioned by Sleep[1988] who emphasized that fractura- tion is unlikely within a mantle undergoing plastic defor- mation. On the basis of geochemical arguments from Oman ophiolites a similar conclusion has been drawn for mid ocean ridges byKelemen et al.[1995]. Alternatively, stress concentrations at the tip of melt rich veins or channels may be expected to be the nucleation point for dikes [Stevenson, 1989]. Partially molten rock may exhibit ductile fracture due to melt propagating along grain boundaries [Fowler, 1990] or along planes of weakness such as boundaries between harzburgite and pyroxenite [Suhr, 1999]. This process requires the fluid overpressureDP = (rsolidrmelt) g dcto exceed the yield strength of peridotite, whereris the density, g is the gravity acceleration and dc is a critical length scale associated with the column of vertically interconnected melt. It should be noted, that dcmay differ from the actual vertical distance of interconnected melt, as a significant part of melt buoyancy is balanced by viscous stresses during porous flow. Estimates of the yield strength for peridotite are of the order 40 to 50 MPa [Nicolas, 1990],
but may possibly become as small as 1 MPa in the melting region or 10 – 15 MPa beneath a freezing boundary [Rubin, 1998;Maaløe, 1998]. Thus, within the melt source region, the yield strengths correspond to critical length scales of the order of 1 to 10 km. This estimate is in agreement with estimates of maximum dike heights of 9 to 11 km on the basis of rock strength analyses by Sato and Ryan [1994]
and Weertman [1971a, 1971b]. Once a dike has formed, buoyancy-driven vein and fracture propagation through the unmolten upper mantle and crust will be possible with high velocities if the dike length exceeds a critical length of 5 – 6 km [Dahm, 2000].
[8] This brief discussion shows that geochemical and geophysical observations and models give some general, but only partly quantitative idea about the mechanisms of melt extraction. They are not conclusive in quantifying the melt fractions necessary for extraction, nor do they provide a viable quantitative mechanism for focusing evenly dis- tributed melts into buoyantly propagating dikes, which are then able to leave the source regions. Therefore independent approaches are needed to better constrain these critical values. In this paper a different approach is presented, utilizing the episodicity of melt extraction. This simple melt extraction model is based on two-phase porous flow within the partially molten region and implicitly accounts for melt channeling and diking by introducing the three extraction parameters dex,j1, and j2.
2. Model
[9] The equations of conservation of mass, momentum and energy of a two-phase material consisting of fluid (subscript f, here melt) and solid matrix (subscript s) are given in the high Prandtl number approximation [McKenzie, 1984;Schmeling, 2000], but see also Ricard et al.[2001], who presented an alternative formulation relating compac- tion to an effective pressure difference between melt and solid)
rf @j
@tþr ~ j~vf
¼G ð1Þ
rs @ð1jÞ
@t þr ~ ðð1jÞ~vsÞ
¼ G ð2Þ
~vf~vs¼ kj hfj
rP~ þrfgdi3
ð3Þ
rgdi3rP~ þ@tij
@xj¼0 ð4Þ
rcp @T
@tþ~vrT~ þag cpvzTabs
¼r ~ krT~
þrHþyLG
ð5Þ whereris the density,jthe melt fraction, t the time,~vthe velocity, G the rate of melt generation (melt mass per
volume and time), kjthe permeability,hfthe melt viscosity, P the (fluid) pressure, g the gravity acceleration,dijthe unit matrix (Kronecker symbol),tij, the deviatoric stress, cpthe heat capacity at constant pressure, T the temperature,athe volumetric thermal expansivity, k the thermal conductivity, H the radiogenic heat generation,ythe dissipation rate, and L the latent heat. The constitutive equations are
tij¼hs @vsi
@xj þ@vsj
@xi
þdij hb2 3hs
r ~~ v ð6Þ
kj¼k0jn ð7Þ
r¼r0 1aTþrsrf rs jþcf
ð8Þ
wherehbandhsare the effective bulk and shear viscosities of the (porous) matrix, respectively, k0and n are constants of the permeability-porosity relation, c is a material constant and f is the melting degree or depletion. The compaction Boussinesq approximation (CBA) is used [see Schmeling, 2000], which essentially neglects density variations except for the buoyancy terms, and neglects compaction for the matrix flow, but accounts for compaction stresses and compaction pressure for the melt flow. The resulting nondimensional flow equations in two dimensions are (nondimensional quantities are primed) [see alsoSchmeling, 2000]
@2
@x02 @2
@z02
h0s @2Y0
@x02 @2Y0
@z02
þ4 @2
@x0@z0h0s@2Y0
@x0@z0
¼Ra@T0
@x0þRm@j
@x0Rc@f
@x0 ð9Þ
~v0f~vd0¼Rm
Rt jn1ð1jÞdi3 1 Rtjn1 @
@x0j h0s @vd0i
@x0j @vd0j
@x0i
!! 1
Rtjn1r~0h*r~0~v0s ð10Þ
where Yis the stream function, defined by
vdx ¼@Y
@z;vdz ¼ @Y
@x ð11Þ
with~vd as the divergence free part of the matrix velocity field ~vs and h* = h0b + 4h0s/3. The thermal, melt and depletion Rayleigh numbers Ra, Rm, Rc, respectively, and the melt retention number Rt are defined as
Ra¼r0gaDTh3 h0k ; Rm¼
rsrf
gh3 h0k ;
ð12Þ Rc¼r0ð1cÞgh3
h0k ; Rt¼hfh2 h0k0
h0andDT are the scaling viscosity and scaling temperature, h is the thickness of the model box, k is the thermal
diffusivity. The CBA has been tested bySchmeling[2000]
and Mu¨ller [2005], who showed that the wavelength- dependent analytical growth rates of the melt channeling instability [Stevenson, 1989] could be reproduced always better than within 10%. Other tests of the code including a comparison of velocities of analytical and numerical solitary porosity waves have been carried out bySchmeling[2000]
and show excellent agreement. Therefore one can be confident that the saw-tooth-like features developing as a consequence of extraction of finite amounts of melts (see below) are well resolved.
[10] Except for the two-dimensional (2-D) case, the heat conductivity is assumed constant and the radiogenic heating is assumed zero. The depletion f of the rock is advected using
@f
@tþ~vrf~ ¼ G
rf ð13Þ
Melting and solidification are modeled using a simplified binary solid solution model
G¼
@fm
@T P
@T
@tþ~vrT~
þ@fm
@P TrgvzforT TsþðfjÞðTlTsÞ 0 elsewhere
8<
:
ð14Þ where fm(T, P) is the temperature- and pressure-dependent melting function of peridotite, Tl and Ts are the depth- dependent liquidus and solidus temperatures of peridotite.
[11] The coupled equations (1), (5), (7), (9), (10), (11), (13) and (14) are solved by finite differences within a rectangular region to yieldj,Y,~vd,~vf,fandG. As equation (10) requires div~vs, this term is determined by solving equation (2) in which@j/@t is eliminated using (1).
[12] Usually, in numerical mantle convection or plume models, melt extraction is modeled by taking all generated melt out of the melt source region, or allowing for a certain amount of retained melt by defining a single extraction threshold value [e.g., Ruedas et al., 2004]. In order to account for the finiteness of a vertically connected melt column (i.e., to overcome the fluid overpressure needed for fracture propagation, see section 1), a new melt extraction model is formulated, which is controlled by the three parameters dex, j1, and j2 (Figure 1). The parameter dex
is defined as the critical thickness of a partially molten layer, in which a critical melt fractionj2 has to be exceeded to allow extraction. The parameter dexmay be identified with the critical length scale dc mentioned above, while the critical melt fraction j2 is that porosity at which porous flow stresses become sufficiently small so that a fluid overpressure can build up. If the actual melt fraction exceedsj2over a vertical distance of dexsomewhere within the melt source region, melt is simply extracted from that region. The natural equivalence of such an extraction event may be the final stage of formation of deep melt rich veins or dikes and further drainage and melt accumulation from the surrounding regions until the dike exceeds the a mini- mum length necessary for buoyancy-driven propagation [Dahm, 2000]. As the melt channels may not completely be drained after an extraction event, and the matrix sur-
rounding the previous melt channels may still contain some melt, a residual melt fraction of j1 is assumed to be left behind after an extraction event.It should be noted that this model approaches continuous melt extraction only in the limit dex!0 andj2!j1.
[13] An alternative extraction criterion would be to define j02instead ofj2in the following way: If the melt fraction exceedsj02at any position, melt is extracted from within a region of thickness dex below that position. This criterion leads to similar results as those shown below ifj02is chosen somewhat larger thanj2.
[14] The above equations are combined with the new extraction formalism and solved in 1-D upwelling flow and in a 2-D ridge-plume-like variable viscosity convection model. For the 1-D upwelling flow a constant divergence- free matrix velocity~vdof 10 cm/yr was assumed. As thermal boundary conditions constant temperature was assumed at the top and bottom. The thickness of the model is 100 km, either describing the depth range 50 – 150 km (first model without extraction) or 20 – 120 km (models with extraction).
The initial temperature condition was adiabatic using a potential temperature of 1350C, or following the solidus curve if the adiabat would exceed the solidus temperature.
This initial temperature distribution ensures that, initially, the models do not contain melt. At the bottom the temper- ature was assumed 300 K hotter than the normal adiabatic temperature to mimic a hot plume entering the model region. Further parameters are: rs rf= 126 kg/m3,h = 1019Pa s,hf= 1 Pa s, k0= 109m2, n = 3,hbis a function of j and varies roughly between h and 25h for melt fractions between 10 and 1% (see Schmeling [2000] for details of this law). A linear melting curve has been assumed with (@fm/@T)P= 0.9103K1.
[15] For the 2-D model the above equations are solved within a 660 km by 660 km region representing one half of an upper mantle spreading zone with a ridge centered (2-D) plume. The mechanical boundary conditions are: at the surface the velocity is prescribed to represent 1 cm/yr of spreading, the sides are symmetric, the bottom is assumed as free slip. The thermal boundary conditions are: T = 0 at the surface, T is symmetric at the sides, T = constant at the
bottom. The initial temperature is assumed adiabatic (or along the solidus as above), a temperature jump of 400 K is assumed at the bottom to initiate a thermal boundary layer.
The olivine-spinel phase transition is included as in work by Marquart et al. [2000]. The rheology is assumed to be temperature-, pressure- and stress-dependent according to laboratory measurements at dry olivine single crystals [Bai et al., 1991]. The thermal diffusivity is assumed depth- dependent (106m2/s at the top, dropping to a minimum of 0.7106m2/s between 20 and 100 km and increasing to 1.3106m2/s at 660 km depth). The model has been run sufficiently long to reach the stage of a well developed plume rising to shallow depth beneath the spreading zone.
[16] It should be mentioned that the 1-D models are intended to be constrained to the melting zone, thus a constant viscosity is assumed for simplicity. In contrast, the 2-D model assumes a variable viscosity, consistently producing an upward obstruction for the rising melt. Nev- ertheless, melt accumulation and porosity waves do not occur in the 2-D models because of efficient melt extraction.
3. Results
3.1. One-Dimensional Model Without Extraction [17] In the first model, a time-dependent 1-D upwelling mantle region is calculated allowing for melt migration due to porous flow. As an initial condition, a subsolidus adia- batic mantle temperature has been assumed as discussed in section 2 and a hot plume-like region with 300 K excess temperature enters the model from below with a velocity of 10 cm/yr (Figure 2). The different temperature and melt porosity curves show snapshots at different times (with approximately 0.2 Myr intervals). As the hot material reaches the depth associated with the asthenosphere, the temperature exceeds the solidus temperature and decom- pression melting occurs. Upon melting, latent heat is consumed, and that effect increases the thermal gradient.
Melt then begins to percolate upward with respect to the matrix. For comparison, the dashed line in Figure 2 (right) shows the melting degree or depletion. As no extraction is allowed in this model, the percolating melt accumulates near the top of the hot upwelling region near the solidifi- cation front (where the temperature intersects the solidus temperature). Near this barrier an accumulation layer forms whose thickness is controlled by the compaction length
lc¼
hbþ4 3h hf kj 0
B@
1 CA
1 2
ð15Þ
Such accumulation layers might be the locus within starting plumes in which melting/solidification processes, deviatoric stresses or other mechanisms will lead to melt focusing within subvertical melt channels and subsequently dike formation. At later stages, such accumulation zones will form immediately beneath regions which have been drained because of localized extraction events.
3.2. One-Dimensional Model With Extraction
[18] To investigate the effect of the proposed three- parameter controlled extraction mechanism, a 1-D model has been calculated similar to the previous model, but Figure 1. Schematic cartoon of proposed extraction
mechanism within an ascending mantle region of decom- pressional melting. A magma pulse is extracted oncej2and dexare simultaneously exceeded.
choosing the extraction parameters dex= 6 km,j1= 1% and j2 = 4%. Figure 3 shows snapshots of the resulting melt fraction at 6 different times, the 0 ka snapshot having been taken after sufficient time had elapsed since the transient stage. As the extraction threshold does not allow high melt fractions within layers thicker than 4 km, no single strong accumulation layer builds up as observed in the previous model. Multiple extraction events occur with a characteristic episodicity. The extraction events occur at different depths within the ascending flow. Although the ascent velocity is assumed to be constant in time and space, each of these melt extraction pulses leaves behind a region with a reduced melt fraction. This drained region acts effectively as a perme- ability barrier for melts immediately beneath. By this mechanism, new accumulation layers are formed, which subsequently are drained if the extraction parameters allow extraction. Thus melt extraction from each single accumu- lation region acts as a trigger of the next extraction event. A characteristic saw tooth distribution of the melt fraction evolves. The arrows in Figure 3 indicate those accumulation regions, in which the extraction threshold is already exceeded, but not yet the critical length, that is, those peaks which are immediately before the next flushing event. The box shaped minima near 1% melt indicate those regions which have been flushed shortly before the snapshots have been taken.
[19] From equation (3) it is possible to evaluate the pressure within the fluid phase and to determine the overpressure within the regions prior to an extraction pulse.
As an example, Figure 4 shows a snapshot of the model displaying the melt concentration (Figure 4a), percolation velocity (Figure 4b), and melt overpressure (Figure 4c). As can be seen, the overpressure has sharp maxima in the regions of maximum percolation velocity, marginally be- neath each maximum of the melt concentration. These are the regions of maximum dilatation. The pressure drops to a minimum value, marginally beneath the melt concentration
minima. The peak to peak amplitude of the overpressure is between 3 and 3.5 MPa, which lies within the range of values reported in the introduction. It should also be noted that the modeled percolation velocities is a result of the two- phase flow, and should not be confused with the (undeter- mined) ascent rates associated with the melt extraction events.
3.3. Variation of Extraction Parameters
[20] Exploring the parameter space of two of the three extraction parameters shows (Figure 5), that small values of both dex and j2 result in multiple extraction depths and short extraction cycles. Actually melt is extracted in small amounts from virtually all depths. As dexincreases, melt is accumulated at greater distances and only few depth regions are drained. Consequently, each pulse flushes huge amounts of melt, of the order of 500 m to 1 km melt layer thickness.
In fact, if dexbecomes larger than 10 km this value is above the compaction length of this model and accumulation layers hardly grow to that thickness. If both dex and j2
become large (such as 6km and 6%, respectively), extrac- tion conditions are not met anymore, and no extraction takes place. Instead a melt fraction similar to that of Figure 2 near steady stage evolves, in which melt migrates to shallow depth by percolation.
[21] It is interesting to evaluate profiles of overpressure within the melt for the different models. The profiles are similar to Figure 4c. For all models withj2> 2% the peak to peak values lie within the range of 3 to 4 MPa, that is, no correlation with dexis found. For the j2= 2% models the peak to peak amplitudes of the overpressure are of the order of 1 MPa, indicating that at small extraction threshold concentration the pressure increase due to percolation and local melt accumulation is significantly reduced. Altogether, the overpressures obtained in the models considered here are of the order of 0.5 to 4 MPa. Such values have been reported or estimated elsewhere (see introduction), indicat- Figure 2. Melt accumulation within an upwelling hot mantle region due to porous flow with
compaction and dilatation. (left) Temperature curves at different model times (approximately 0.2 m.y.
intervals). (right) Melt fractions at corresponding times. The parameters are vz = 10 cm/yr,rs rf= 126 kg/m3, h = 1019 Pa s, hf = 1 Pa s, k0 = 109m2, and n = 3;hb is a function of j and varies roughly between hand 25hfor melt fractions between 10% and 1% (seeSchmeling [2000] for details of this law). A linear melting curve has been assumed with (@fm/@T)P = 0.9 103 K1.
ing that the extraction parameters are within a realistic range of magnitudes.
[22] Interestingly, in most of the cases shown the extrac- tion events always start at a constant depth level within the upwelling melt source region. This is because the deepest extraction event always occurs at a depth at which the steady state part of the melt fraction curve of Figure 2 exceeds j2 over a distance of at least dex, and from that depth on, a similar section of the steady state melt fraction curve (shortened by the part between 0 andj1) is needed to trigger the next extraction event above, and so on. This rather regular behavior is a consequence of continuous supply of unmolten mantle from below in 1-D. In reality,
the actual extraction process may occur over a vertical continuum in growing veins and dikes arranged in both offset and overlapping arrangements. It is interesting to note that the model of production, accumulation and transport of melt byBons and van Milligen[2001] does not include this continuous influx, and quickly evolves into a state of self- organized criticality.
[23] In section 2 an alternative extraction criterion had been mentioned, using j02 instead of j2. Inspecting the models enables us to determine thej02, values which would give the same extraction behavior as the models shown: For the three depths dexof 3, 6, and 10 km, the maximum melt concentrationsj02reached are 2.26, 3.0, and 2.55% for the Figure 3. Snapshots of the melt fraction versus depth in a 1-D rising plume flow (10 cm/yr, 300 K
excess temperature) with extraction parameters dex= 6 km,j1= 1%, and j2= 4%. Extraction pulses originate from the peaks marked by arrows, as can be seen by the resulting minima (with 1% melt) in the immediately following stages.
j2= 2% models, respectively, 4.53, 5.2, and 5.92% for the j2= 4% models, respectively, and 7.1% for the j2= 6%, dex= 3 km model. This evaluation indicates, thatj2andj02
differ only by at most 1.9%.
[24] Figures 5b and 5c show the quantitative results of the extraction periods and associated amounts of melts.
4. Discussion and Application to Iceland, Hawaii, and Canary Islands
4.1. Scaling
[25] In order to apply the models to natural examples, it has to be discussed how the results depend on the particular choice of other model parameters. Quantitative scaling of the results to other scenarios is difficult, because the occurrence of the extraction events depends on (1) the shape of the steady state solution of the percolating melt (dash-dotted curve in Figure 2) and (2) the temporal evolution and thickness of melt accumulation peaks.
[26] The initial slope @j/@z of the steady state solution is given by (cp(@Ts/@z) agTabs)/(L + (TlTs)cp), thus it depends only on thermodynamic parameters of melting and the adiabatic temperature gradient, and not on dynamic
parameters such as the upwelling velocity vz. On the other hand, the flattening of the steady state j(z) curve (in other words, the vertical distance between the steady state (dash-dotted) and linear (dashed) curves of Figure 2 (right)) is given by rf (rs rf) g k0 jn/(hf G) with G = rf vz (cp (@Ts/@z) agTabs)/(L + (Tl Ts)cp). Thus, if percolation is easy (high buoyancy of the melt, low melt viscosity or high permeability) or the melting rate or upwell- ing velocity is slow, this term may get large and the steady state solution may not reach sufficiently high melt fractions for extraction. Then, for small dex, extraction is controlled by the amplitude and width of the accumulation peaks. As mentioned before, the width of these peaks scales with the compaction length (equation (15)). A small compaction length (e.g., due to smaller rock viscosities) will make it more difficult for the melt accumulation zone to exceed dex, but the amplitude of the peak will grow faster, allowing the melting process to more easily exceed the thresholdj2. If dex
is larger than the compaction length, not only the accumu- lation peak but also a sufficiently large part of the ‘‘steady state part’’ of thej(z) curve has to exceedj2.
[27] The above arguments imply, that for smaller mantle upwelling velocity the extraction behavior is roughly Figure 4. One snapshot of the 1-D model shown in Figure 3: (a) melt concentration, (b) percolation
velocity due to porous flow, (c) fluid overpressure within the melt phase, and (d) depletion.
equivalent to a decrease of j2combined with an increase of dex. For example, the model with j2 = 4% and dex = 6 km (cf. Figure 5) has been run with half the mantle ascent velocity: As a result (not shown), extraction events occur only at two depths (50 and 80 km) with an episodicty of 1 per 30,000 years.
4.2. Melt Extraction in Two Dimensions
[28] It is interesting to consider the behavior of the extraction in two dimensions. Therefore the three-parameter extraction formalism has been applied to a 2-D convection model with melting and melt segregation. In this model with nonlinear, temperature- and depth-dependent viscos- ity based on laboratory experiments, a hot plume starting from a thermal boundary layer at 660 km depth rises beneath a spreading ridge with a spreading velocity of 1 cm/yr (see Schmeling [2000] for further discussion of this kind of 2-D models with segregation). The extraction rules are now applied to vertical finite difference grid columns. As the flow and temperature field now varies with depth and horizontal coordinate, extraction events occur in a stochastic way, and the melting region of the hot upwelling becomes very heterogeneous. In Figure 6 a snapshot of this model is shown. The heterogeneous melt
fractions (Figure 6a), the anomalous temperature (TTad, where Tad is the depth-dependent adiabatic temperature with a surface potential temperature of 1350C) and the depletion are translated into seismic shear velocity anoma- lies (Figure 6b). For the temperature dependence, anhar- monic and anelastic effects are accounted for. The dependence on melt fraction is based on a finite element model applied to melt distributions observed in laboratory experiments. The method used for obtaining the seismic velocities is described in detail by Kreutzmann et al.
[2004], where also further references are given. In the lower part of the melting zone the green/blue spot indicates a region which has been flushed shortly before the snapshot was taken. If such a situation would be encountered by a seismic investigation, for example, by surface wave inversion, it might be possible that it would be interpreted as a double low-velocity zone. In fact, Kreutzmann et al.[2004] discussed seismic velocity anoma- lies of the Iceland plume, and explained the observation of a double low-velocity zone beneath central Iceland by dehydration of the melting plume. As an alternative expla- nation emerging from the present model one might consider the possibility that the melting region of the Iceland plume is strongly heterogeneous because of localized episodic Figure 5. (a) Depth-time distribution of extraction pulses for different extraction parameters dexandj2.
Each plot represents one 1-D extraction model run; the time in each plot is in m.y. and depth is in km.
Each dot represents the time and position (lowest point of a region of height dex) of one extraction pulse.
Contour lines of (b) periods in years between extraction pulses and (c) corresponding melt layer thickness of each extraction pulse as a function of dexandj2;j1has been chosen constant (1%). In the top right corner (j2= 6%, dex> 5 km) no extraction occurs.
melt extraction, and the seismologists just see a snapshot of such heterogeneities.
4.3. Comparison With Observations
[29] It is interesting to compare the results with observed extraction timescales and volumes of intraplate and rift zone magmatism, such as Hawaii and Iceland, respectively.
Eruption frequencies of the Hawaiian hot spot vary between one per year to one per 25 years with erupted volumes of 0.01 to 0.46 km3 for the Kilauea and Mauna Loa volca- noes, respectively [Ryan et al., 1981; Takada, 1999;
Maaløe, 1999]. These eruption frequencies are significantly higher than those predicted by the proposed melt extrac- tion model. In order to increase the extraction frequency to the observed Hawaiian values very low extraction param- eters of j1, j2 and dex would be needed (j1, j2 2%, dex 3 km, cf. Figure 5). On the basis of structural observation, it seems more likely that the observed Hawaiian episodicity is the result of a magma plumbing system that works under self-control by stress partitioning [Ryan, 1988;Takada, 1999]. The structure of the magmatic system beneath Hawaii has been inferred by Ryan et al.
[1981] and Ryan [1988] on the basis of seismicity inter- preted as resulting from ascending dikes and magma.
Beneath 60 km, microseismic events are typically not generated, but between 60 and 40 km depth seismic harmonic tremor has been observed over a monitoring period of 17 years indicating a generally steady process of magma ascent in a connected system of dikes [Aki and Koyanagi, 1981]. If this apparently long-term steady state nature of deep magma ascent is part of an extraction pulse as modeled in this paper, it may provide a lower bound for the extraction frequency of several decades. Such a fre- quency still corresponds to low extraction parameters of j1,j2and dex(cf. Figure 5) On the other hand, the tremor episodes and the seismicity in that region are spatially
scattered suggesting that at least some of the larger extraction events are both intermittent and spatially scat- tered as well [Ryan, 1988]. If these events are triggered by the extraction processes modeled here, very low extraction parameters of j1, j2 and dex would be predicted. At shallower depth (above 40 km) magma ascends within a well defined conduit consisting of a core of interconnected veins and dikes and occasionally activated offset dikes near the conduit margins [Ryan, 1988]. Magma partitioning into dike intrusions near the level of neutral buoyancy and into surface eruptions changes the crustal stress field, which, as a feed back mechanism, controls the magma supply of the shallow dike system [Ryan1987;Takada, 1999]. From this discussion it is clear, that any potential episodicity of melt extraction within the melting source region (60 – 120 km depth) is difficult to infer for the Hawaiian case.
[30] In Iceland, frequencies and volumes of major erup- tion and intrusion episodes vary from one per 250 years and 0.1 km3 (Krafla) to one per 2800 years and 15 km3 (Lakagigar) [Takada, 1999]. Most of the Icelandic volcanic systems consist of one central volcano, underlain by a shallow magma chamber and a system of adjacent fissures.
These fissures seem to be fed laterally by the magma chamber. Outside of these central volcano systems eruptions are rare but of large volume, fed mainly by thick, subvert- ical dikes [Gudmundsson, 1998]. Gudmundsson distin- guishes between eruptions of the central volcano with frequencies of once every few hundred years and volumes less than 0.1 km3and fissure eruptions outside the central volcano (not fed by the magma chamber of the central volcano) every several thousand years with volumes of several km3. While the temporal behavior of the central volcano eruptions are probably controlled by the crustal feedback mechanism mentioned above, the fissure swarm eruptions on Iceland not associated with a central volcano may well represent the episodicity generated in the source Figure 5. (continued)
regions [Takada, 1999]. If the Iceland plume has a vertical velocity of about 10 cm/yr within the melt generation zone [Ruedas et al., 2004], Figure 5 could explain the Icelandic periods with extraction parameters ofj2= 2 – 3% and dex= 3 – 5 km. The corresponding melt layer thickness of about 75 m corresponds well to an observed magma volume of 15 km3 for the Lakagigar eruption if distributed over an area of 200 km2.
[31] Episodic volcanism of the Iceland plume on the much larger timescale of 5 – 10 Myr has been found for the late Cretaceous and early Tertiary from dating of several
seamounts in the Rockall trough off the Atlantic coast of Scotland [O’Connor et al., 2000], and for the late Tertiary from the spacing of V-shaped bathymetry of the Reykjanes ridge [Vogt, 1971]. Obviously these timescales cannot be explained by the present model. Instead, the vigor of the plume itself has to be assumed to pulsate resulting in temporal variations in crust production. Schmeling and Marquart [1993] found that the frequency of temporal variations of melt generation in a convecting mantle with nonlinear rheology correlates with the effective astheno- spheric viscosity. From their models, magmatic recurrence Figure 6. (a) Melt fraction (%) and (b) seismic shear velocity anomaly (%) of a partially molten 2-D
plume head whose melt is extracted according to the proposed mechanism (j2= 6%,j1= 1%, and dex= 6 km). The plume rises in the central region and flows laterally out of the region shown. The actual region of calculation is larger than the region shown. The dark blue region in Figure 6b near the top represents cold lithosphere. The melt fraction, the anomalous plume temperature (locally about 180 K hotter than the reference mantle), and the depletion (maximum 18%) reduce the shear velocity locally by up to 15%,5%, and0.5%, respectively. The effect of heterogeneous distribution of melt fraction clearly dominates the anomaly and is responsible for the heterogeneous appearance of the shear velocity.
times of 5 – 10 m.y. correlate with effective asthenospheric viscosities of about 21019Pa s.
[32] Figures 3 or 5 suggest that melt may be extracted from different depths within one melt column. This is in agree- ment with major element and isotope studies of ocean island basalts from Lanzarote (Canary Islands) by Thomas et al.
[1999]. They concluded that differences in the average composition depend on the thickness of the overlying lithosphere, and the differences within single lava suites indicate different depths of extraction from the melt column.
In particular, they identified lavas which were extracted from the bottom quarter of a 28 km melting column, that is, the extraction event should have drained a 7 km section of the column. From Figure 5, extraction parameters of aboutj2= 2% and dex= 5 km could tentatively explain such an event.
Surprisingly, these values are in good agreement with the parameters derived from Icelandic eruption frequencies (see above).
[33] The extraction model proposed here implies that the residual heterogeneities within the melt source region should not vary very strongly: As a result of the 1-D model set up (constant upwelling velocity) and melting relation (equations (13) and (14)), continuous decompressional melting leads to a depletion distribution linearly increasing with z coordinate (Figure 4d). If a section of the model source region would be exposed at the surface, its distribu- tion of depletion would at most show the linear trend of Figure 4d superimposed with a frozen-in distribution of previous melt concentration (Figure 4a). Thus the model cannot explain strong repetitive compositional variations found, for example, in the Horoman peridotite (Hokkaido) [Takazawa et al., 2000].
5. Conclusions
[34] A simple model of melt segregation and extraction within rising plumes is proposed. It is based on two-phase porous flow within the partially molten region, combined with three extraction parameters, namely two porosity thresholdsj1andj2, and a finite fraction of the melt column in which thej2threshold is exceeded. This model is able to predict characteristic episodicities of melt extraction in mantle upwelling regions and multiple depths for extraction.
Thus it provides an independent method for estimating extraction parameters. Tentatively applied to natural exam- ples (Iceland, Canary Islands) an extraction threshold of at most 2 – 3% is derived, while dexshould not exceed 3 – 5 km.
This result in general agreement with isotope activity rations which demand a threshold porosity of 1% or less. It should, however, be kept in mind, that this three-parameter approach is only a very crude approximation of the complex process of melt channeling, melt accumulation and dike formation near the top of the melting zones.
[35] Acknowledgments. This work benefited from valuable discus- sions with T. Ruedas. Constructive reviews of M. P. Ryan, two anonymous reviewers, and the associate editor H. Iwamori are gratefully acknowledged.
This work was done partly under the DFG grants Schm872-6 and 11.
References
Aharonov, E., J. A. Whitehead, P. B. Kelemen, and M. Spiegelman (1995), Channeling instability of upwelling melt in the mantle,J. Geophys. Res., 100, 20,433 – 20,450.
Aki, K., and R. Koyanagi (1981), Deep volcanic tremor: An magma ascent mechanism under Kilauea, Hawaii,J. Geophys. Res.,86, 7095 – 7109.
Bai, Q., S. J. Mackwell, and D. L. Kohlstedt (1991), High temperature creep of olivine single crystals: Mechanical results for buffered samples, J. Geophys. Res.,96, 2441 – 2463.
Bons, P. D., and B. P. van Milligen (2001), New experiment to model self-organized critical transport and accumulation of melt and hydro- carbons from their source rocks,Geology,29(10), 919 – 922.
Dahm, T. (2000), On the shape and velocity of fluid-filled fractures in the Earth,Geophys. J. Int.,142, 181 – 192.
Elliot, T. (1997), Fractionation of U and Th during mantle melting: A reprise,Chem. Geol.,139, 165 – 183.
Faul, U. H. (1997), Permeability of partially molten upper mantle rocks from experiments and percolation theory, J. Geophys. Res., 102, 10,299 – 10,311.
Faul, U. H. (2001), Melt retention and segregation beneath mid-ocean ridges,Nature,410, 920 – 923.
Fowler, A. C. (1990), A compaction model for melt transport in the Earth’s asthenosphere. II. Applications, inMagma Transport and Storage, edited by M. P. Ryan, pp. 15 – 32, John Wiley, Hoboken, N. J.
Fujii, N., K. Osamura, and E. Takahashi (1986), Effect of water saturation on the distribution of partial melt in the olivine-pyroxene-plagioclase system,J. Geophys. Res.,91, 9253 – 9259.
Gudmundsson, A. (1998), Magma chambers modeled as cavities explain the formation of rift zone central volcanoes and their eruption and intru- sion statistics,J. Geophys. Res.,103, 7401 – 7412.
Hart, S. R. (1993), Equilibrium during melting: A fractal tree model,Proc.
Natl. Acad. Sci. U. S. A.,90, 11,914 – 11,918.
Hirth, G., and D. L. Kohlstedt (1995a), Experimental constraints on the dynamics of the partially molten upper mantle: 1. Deformation in the diffusion creep regime,J. Geophys. Res.,100, 1981 – 2001.
Hirth, G., and D. L. Kohlstedt (1995b), Experimental constraints on the dynamics of the partially molten upper mantle: 2. Deformation in the dislocation creep regime,J. Geophys. Res.,100, 15,441 – 15,449.
Iwamori, H. (1993), A model for disequilibrium mantle melting incorpor- ating melt transport by porous and channel flows,Nature,366, 734 – 737.
Kelemen, P. B., N. Shimizu, and V. J. M. Salters (1995), Extraction of mid-ocean-ridge basalts from the upwelling mantle by focused flow of melt in dunite channels,Nature,375, 747 – 753.
Kreutzmann, A., H. Schmeling, A. Junge, T. Ruedas, and G. Marquart (2004), Dynamics and melting of a ridge-centered plume with application to Iceland, part II: Predictions for electromagnetic and seismic observa- bles,Geophys. J. Int.,159, 1097 – 1111.
Landwehr, D., J. Blundy, E. M. Chamorro-Perez, E. Hill, and B. Wood (2001), U-series disequilibria generated by partial melting of spinel lherzolite,Earth Planet. Sci. Lett.,188, 329 – 348.
Maaløe, S. (1998), Melt dynamics of a layered mantle plume source, Contrib. Mineral. Petrol.,133(1 – 2), 83 – 95.
Maaløe, S. (1999), Magma accumulation in Hawaiian plume sources,Am.
J. Sci.,139, 139 – 156.
Maaløe, S. (2003), Melt dynamics of a partially molten mantle with ran- domly oriented veins,J. Petrol.,44, 1193 – 1210.
Marquart, G., H. Schmeling, G. Ito, and B. Schott (2000), Conditions for plumes to penetrate the mantle phase boundaries,J. Geophys. Res.,105, 5679 – 5693.
McKenzie, D. (1984), The generation and compaction of partially molten rock,J. Petrol.,25, 713 – 765.
McKenzie, D. (2000), Constraints on melt generation and transport from U-series activity ratios,Chem. Geol.,162, 81 – 94.
Mu¨ller, K. (2005), Numerische Untersuchung der spannungsangetriebenen Schmelzsegregation mit Anwendung auf einen Plume unter einem Mittel-Ozeanischen Ru¨cken, Ph.D. thesis, Johann Wolfgang Goethe Univ. Frankfurt am Main, Frankfurt, Germany.
Nicolas, A. (1986), A melt extraction model based on structural studies in mantle peridotites,J. Petrol.,27, 999 – 1022.
Nicolas, A. (1990), Melt extraction from mantle peridotites: Hydrofractur- ing or porous flow consequences on oceanic ridge activity, inMagma Transport and Storage, edited by M. P. Ryan, pp. 160 – 174, John Wiley, Hoboken, N. J.
O’Connor, J. M., P. Stoffers, J. R. Wijbrans, P. M. Shannon, and T. Morrissey (2000), Evidence from episodic seamount volcanism for pulsing of the Iceland plume in the past 70 Myr,Nature,408, 954 – 958.
Rabinowicz, M., P. Genthon, G. Ceuleneer, and M. Hillairet (2001), Com- paction in a mantle mush with high melt concentrations and the genera- tion of magna chambers,Earth Planet. Sci. Lett.,188, 313 – 328.
Ricard, Y., D. Bercovici, and G. Schubert (2001), A two phase model for compaction and damage: 2: Compaction, J. Geophys. Res., 106, 8907 – 8924.
Richardson, C. N. (1998), Melt flow in a variable viscosity matrix,Geo- phys. Res. Lett.,25, 1099 – 1102.
Richardson, C. N., and D. P. McKenzie (1994), Radioactive disequilibria from 2-D models of melt generation by plumes and ridges,Earth Planet.
Sci. Lett.,128, 425 – 437.
Rubin, A. M. (1998), Dyke ascent in a partially molten rock,J. Geophys.
Res.,103, 20,901 – 20,919.
Ruedas, T., H. Schmeling, G. Marquart, A. Kreutmann, and A. Junge (2004), Dynamics and melting of a ridge-centered plume with application to Iceland, part I: Evolution and crust production,Geophys. J. Int.,158, 729 – 743.
Ryabchikov, I. D. (1998), Geochemical modelling of magma generation in passive and active mantle plumes, inChallenges to Chemical Geology, edited by M. Novak and J. Rosenbaum, pp. 103 – 120, Czech. Geol.
Surv., Prague.
Ryan, M. P. (1987), Neutral buoyancy and the mechanical evolution of magmatic systems, inMagmatic Processes: Physiochemical Principles, edited by B. O. Mysen,Geochem. Soc. Spec. Publ.,1, 259 – 288.
Ryan, M. P. (1988), The mechanics and three-dimensional internal structure of active amagmatic systems: Kilauea Volcano, Hawaii,J. Geophys. Res., 93, 4213 – 4228.
Ryan, M. P., R. Y. Koyanagi, and R. S. Fiske (1981), Modeling the three- dimensional structure of macroscopic magma transport systems: Applica- tion to Kilauea Volcano, Hawaii,J. Geophys. Res.,86, 7111 – 7129.
Sato, H., and M. P. Ryan (1994), Generalized upper mantle thermal structure of the western United States and its relationship to seismic attenuation, heat flow, partial melt and magma ascent and emplacement, in Magmatic Systems, edited by M. P. Ryan, pp. 259 – 286, Elsevier, New York.
Schmeling, H. (2000), Partial melting and melt segregation in a convecting mantle, inPhysics and Chemistry of Partially Molten Rocks, edited by N. Bagdassarov, D. Laporte, and A. B. Thompson, pp. 141 – 178, Springer, New York.
Schmeling, H., and G. Marquart (1993), Mantle flow and the evolution of the lithosphere,Phys. Earth Planet. Inter.,79, 241 – 267.
Sims, K. W. W., D. J. DePaolo, M. T. Murrell, W. S. Baldridge, S. Goldstein, and D. Clague (1995), Mechanism of magma generation beneath Hawaii and mid-ocean ridges: Uranium/thorium and samarium/neodynium isotope evidence, Science, 267, 508 – 512.
Sims, K. W. W., D. J. DePaolo, M. T. Murrell, W. S. Baldridge, S. Goldstein, D. Clague, and M. Jull (1999), Porosity of the melting zone and variations in the solid mantle upwelling rate beneath Hawaii:
Inferences from 238U-230Th-226Ra and 235U-231Pa disequilibria, Geochim. Cosmochim. Acta,63, 4119 – 4138.
Sleep, N. H. (1988), Tapping of melt by veins and dikes,J. Geophys. Res., 93, 10,255 – 10,272.
Spiegelman, M., and T. Elliot (1993), Consequences of melt transport for uranium series disequilibrium in young lavas,Earth Planet. Sci. Lett., 118, 1 – 20. (CorrectionEarth Planet. Sci. Lett.,121, 666.)
Spiegelman, M., P. B. Kelemen, and E. Aharonov (2001), Causes and consequences of flow organization during melt transport: The reaction infiltration instability in compactible media, J. Geophys. Res., 106, 2061 – 2077.
Stevenson, D. J. (1989), Spontaneous small-scale segregation in partial melts undergoing deformation,Geophys. Res. Lett.,9, 1064 – 1070.
Suhr, G. (1999), Melt migration under oceanic ridges: Inferences from reactive transport modeling of upper mantle hosted dunites,J. Petrol., 40, 575 – 599.
Takada, A. (1999), Variations in magma supply and magma partitioning;
role of the tectonic settings,J. Volcanol. Geotherm. Res.,93, 93 – 110.
Takazawa, E., F. A. Frey, N. Shimizu, and M. Obata (2000), Whole rock compositional variations in an upper mantle peridotite (Horoman, Hokkaido, Japan): Are they consistent with a partial melting process?, Geochim. Cosmochim. Acta,64, 695 – 716.
Thomas, L. E., C. J. Hawkesworth, P. Van Calsteren, S. P. Turner, and N. W.
Rogers (1999), Melt generation beneath ocean islands: A U-Th-Ra isotope study from Lanzarote in the Canary Islands,Geochim. Cosmochim. Acta, 63, 4081 – 4099.
Toramaru, A., and N. Fujii (1986), Connectivity of melt phase in a partially molten peridotite,J. Geophys. Res.,91, 9239 – 9252.
Vogt, P. B. (1971), Asthenospheric motion recorded by the ocean floor south of Iceland,Earth Planet. Sci. Lett.,13, 153 – 160.
Weertman, J. (1971a), Theory of water-filled crevasse in glaciers applied to vertical magma transport beneath oceanic ridges,J. Geophys. Res.,76, 1171 – 1183.
Weertman, J. (1971b), Velocity at which liquid-filled cracks move in the Earth’s crust or in glaciers,J. Geophys. Res.,76, 8544 – 8553.
Wiggins, C., and M. Spiegelman (1995), Magma migration and magmatic solitary waves in 3-D,Geophys. Res. Lett.,22, 1289 – 1292.
H. Schmeling, Geophysics Section, Department of Earth Sciences, Johann Wolfgang Goethe Universita¨t Frankfurt am Main, Feldbergstr. 47, D-60323 Frankfurt, Germany. (schmeling@geophysik.uni-frankfurt.de)
