Petrological geodynamic modeling of mid-ocean ridges

Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

Contents lists available at SciVerse ScienceDirect

Physics of the Earth and Planetary Interiors journal homepage: www.elsevier.com/locate/pepi

Petrological geodynamic modeling of mid-ocean ridges M. Tirone a,⇑, G. Sen b, J.P. Morgan c a

Institut für Geologie, Mineralogie und Geophysik, Bochum 44780, Germany College of Arts and Sciences, American University of Sharjah, Sharjah, United Arab Emirates c Department of Earth and Atmospheric Sciences, Cornell University, Ithaca, 14853 NY, USA b

a r t i c l e

i n f o

Article history: Received 8 April 2011 Received in revised form 24 October 2011 Accepted 28 October 2011 Available online 18 November 2011 Edited by Mark Jellinek Keywords: Mid-ocean ridge Mantle melting Thermodynamics Geodynamics Numerical modeling Abyssal peridotites

a b s t r a c t Mid-ocean ridges are the primary location where the Earth’s oceanic crust is formed. Beneath spreading ridges several processes such as dynamic melting and partial crystallization modify the petrology of the upper mantle and affect the Earth’s global geochemical evolution. A unified picture of the temporal and spatial evolution of melt and residual mantle, as well as crustal production and melt dynamics requires a comprehensive model that takes into account simultaneously the complexity of the physical processes involved and the petrological variations of the ridge system. Here we present the first results of a 2-D numerical approach applied to a spreading ridge that fully couples a two-phase flow model for melt and solid mantle and a chemical thermodynamic model which provides a spatial and temporal description of the minerals and melt abundance and composition. The most significant features found by this study are the following. (1) Accumulation of melt is observed at the base of the lithosphere in the off-axis region (<50 km from the ridge axis). (2) Crustal production (thickness) shows temporal variations which are mainly induced by periodic discharge of the melt accumulated underplate. (3) Magma waves develop between 10 and 30 km depth in proximity of the ridge axis. However to accurately resolve melt fluctuations, the grid size must be smaller than the compaction length for porous flow. Since in this study the compaction length decreases with depth, we have used a simplified 1-D melt model incorporating the two-phase flow dynamics and the thermodynamic formulation to show that the depth at which magma waves start to form increases by increasing the numerical resolution. Despite the limitation of the numerical grid resolution, we have observed that variations of the melt content do not appear to have significant influence on major elements composition of the residual solid and melt. (4) In the initial stage of the ridge evolution, a melting area detaches from the main melting region around the ridge axis. It is possible that this type of development may repeat over time beyond the duration of the simulation model of this study (15 Ma). Sluggish coupling between the dynamics of the lithosphere and the asthenospheric mantle flow suggests that accretion of the lithosphere by conductive cooling away from the ridge center involves portions of the upper mantle that not necessarily passed through the spreading ridge. (5) During the development of the spreading ridge, the asthenosphere affected by the melting process deflects downwards, creating in this way a chemical heterogeneity in the large mantle circulation. (6) Composition of major elements in the residual solid after partial melting is in agreement with the chemical pattern observed in abyssal peridotites. However, in order to explain the large variation of major elements content found in abyssal peridotite, a consistent petrological and geodynamic model of the evolution of the mid-ocean ridge, requires that partial crystallization of small amount of melt refertilizes the depleted mantle. The petrological model presented in this study accounts for the complexity of polybaric dynamic melting and the continuous reactions between the residual solid and melt, but it is limited by the assumption of local thermodynamic equilibrium within a domain defined by the numerical grid size. The interpretation of the petrological results needs to be carefully evaluated to ensure that the time and space scale of the numerical model complies with the constraints provided by solid–melt reactive experiments and the spatial scale of the petrological structures observed in mid-ocean ridges. (7) Melt distribution and thermal structure are revealed by the seismic shear wave map computed from the numerical model. Certain observations, such as the extent of the melting region, overall agree quite well with the evidences from seismic studies from various ridge settings. Ó 2011 Elsevier B.V. All rights reserved.

⇑ Corresponding author. E-mail addresses: [email protected] (M. Tirone), [email protected] (G. Sen), [email protected] (J.P. Morgan). 0031-9201/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2011.10.008

52

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

1. Introduction Current understanding of mantle processes under mid-ocean ridges relies on geophysical and geochemical observations. Interpretation of seismological data along with electrical conductivity and magnetotelluric data provide detailed information on the presence of a magma chamber, the internal structure of the crust and upper mantle and the extent of the melting region (e.g. Sinton and Detrick, 1992; Zhang et al., 1994; Cannat, 1996; Tucholke et al., 1997; Canales et al., 1998; Forsyth et al., 1998; Toomey et al., 1998; Webb and Forsyth, 1998; Evans et al., 1999; Canales et al., 2000; Dunn et al., 2000; Dunn and Forsyth, 2003; Dunn et al., 2005; Evans et al., 2005; Nedimovic´ et al., 2005; Baba et al., 2006; Yang et al., 2007). Using gravity and bathymetric data, it was also inferred that crust thickness varies with time (Pariso et al., 1995; Tucholke et al., 1997; Heinson et al., 2000; Bonatti et al., 2003; Lizarralde et al., 2004). The observation was also supported by petrological evidence (Cipriani et al., 2009). The general interpretation is that the temporal variations are related to the velocity of the spreading plate and melt production or some mechanism associated with melt extraction. Extensive information stored in the PetDB database (www.petdb.org) is available on the composition of MORB melt and abyssal peridotites collected from the bottom of the ocean floor (e.g. Dick and Fisher, 1984; Michael and Bonatti, 1985; Shibata and Thompson, 1986; Klein and Langmuir, 1987; Dick, 1989; Johnson et al., 1990; Johnson and Dick, 1992; Niu and Batiza, 1993; Baker et al., 1995; Seyler and Bonatti, 1997; Michael and Cornell, 1998; Yang et al., 1998; Niu, 2004; Bodinier and Godard, 2005; Workman and Hart, 2005; Eason and Sinton, 2006; Rubin and Sinton, 2007; Dick et al., 2008; Niu and O’Hara, 2008). Petrological studies on ophiolites, which represent displaced portions of ancient ridge structures, reveal a complex stratigraphic sequence of ultramafic and mafic rocks, ranging from peridotites that experienced various degrees of depletion to dunites and gabbros (e.g. Jousselin et al., 1998; Godard et al., 2000; Braun and Keleman, 2002; Girardeau et al., 2002; Le Mée et al., 2004). Quantitative interpretations of these observations are based either on petrological models or on physical descriptions of magma transport. Empirical parameterization of melting, petrological considerations and thermodynamic formulation (Ghiorso and Sack, 1995; Ghiorso et al., 2002) are used to interpret the chemical evolution of MORB and other mafic and ultramafic formations associated to mid-ocean ridges and peridotite melting (Klein and Langmuir, 1987; McKenzie and Bickle, 1988; Johnson et al., 1990; Kinzler and Grove, 1992a; Kinzler and Grove, 1992b; Langmuir et al., 1992; Plank and Langmuir, 1992; Kinzler and Grove, 1993; Niu et al., 1997; Hirschmann et al., 1998; Asimow, 1999; Hirschmann et al., 1999a; Hirschmann et al., 1999b; Asimow and Stolper, 1999; Asimow et al., 2001; Asimow, 2002; Asimow and Langmuir, 2003; Herzberg, 2004). Essential information are provided by phase equilibrium experimental data on peridotite melting (e.g. Jaques and Green, 1980; Takahashi, 1986; Falloon et al., 1988; Hirose and Kushiro, 1993; Takahashi et al., 1993; Walter and Presnall, 1994; Walter, 1998; Gudfinnsson and Presnall, 2000; Presnall et al., 2002). The petrological models usually consider ideal cases, for instance assuming either fractional or batch melting. However it appears that the melt evolution under mid-ocean ridges is the result of more complex physico-chemical processes. Previous studies suggested that dunites are the mantle conduits of melt moving through the oceanic mantle which are developed through continuous reactions of dissolution and precipitation with the transient melt (Kelemen, 1990; Kelemen et al., 1997; Braun and Keleman, 2002). Experimental work and field observations revealed that melting of peridotite rocks is more likely associated with kinetic reactive processes and partial equilibration between melt and the mineral phases in the residual solid (Kelemen et al., 1990; Kelemen et al., 1992; Kinzler and Grove,

1992a; Wagner and Grove, 1998; Longhi, 2002; Morgan and Liang, 2005; Lambart et al., 2009; Van Den Bleeken et al., 2010). Recent petrological models are beginning to account for this more complex interaction between melt and the solid assemblage (Lissemberg and Dick, 2008; Lambart et al., 2009; Collier and Kelemen, 2010). On the other hand, transport of melt in ridge settings has been extensively investigated by numerical methods using mainly a two phase flow formulation (e.g. Morgan, 1987; Spiegelman and McKenzie, 1987; Ribe, 1988; Scott and Stevenson, 1989; Sparks and Parmentier, 1991; Su and Buck, 1993; Braun et al., 2000; Ghods and Arkani-Hamed, 2000; Rabinowicz and Ceuleneer, 2005; Katz, 2008; Rabinowicz and Toplis, 2009). The melting process in these studies is usually based on simple models that not necessarily reproduce the petrological behavior of mantle rocks. One of the most challenging questions related to mid-ocean ridges is the temporal and spatial relation between the dynamics of the oceanic mantle, crust and melt and the petrological evolution that leads to the formation of MORBs, depleted lherzolite, dunites and the emplacement of intrusive rocks such as gabbros. A combined petrological and geodynamic approach may be able to address such problem and, more in general, to provide a comprehensive description of the evolution of mid-ocean ridges. An early attempt in this direction (Cordery and Morgan, 1993) using the parameterization of Kinzler and Grove (1992a) for the melting process, simultaneously explored the chemical and dynamical evolution of melt in the steady state. The study was based on certain simplifications. It assumed that the assemblage opx-cpx-ol-sp was always present and that the composition of the system was controlled only by the dynamic evolution of the residual mantle. By linking the melt content to the residual solid composition, the model assumed fractional melting, but at the same time the melt formed locally was added to the existing melt in a way that resembled a incremental batch melting. Despite these assumptions, the model reproduced reasonably well the general trend of major elements in abyssal peridotites and MORB melts. In alternative to a parameterized petrological model a thermodynamic formulation can be used instead. The main advantage of combining chemical thermodynamic principles with a dynamic model (Fullea et al., 2009; Tirone et al., 2009) is that petrological field observations can be related with the simulated results in a spatial and temporal setting and simultaneously the thermophysical properties that control the dynamics of the process are self-consistently incorporated in the procedure. For mid-ocean ridges the petrological observations should include MORBs, abyssal peridotites, depleted lherzolites and gabbros. A strong validation of the modelling work can be provided by a successful comparison of the petrological field evidences with the numerical computation together with a description of the thermal and dynamical evolution of the mid-ocean ridge. 2. Method The governing equations of the transport model, a brief summary of the numerical procedure and the parameters used in this study are discussed in Sections 2.1 and 2.3. Abundance and composition of mantle minerals and melt are retrieved from a thermodynamic model that minimizes the Gibbs energy for a given local bulk composition pressure and temperature. Section 2.2 presents the thermodynamic database that has been developed for this work. The relevant thermodynamic parameters for melt and solid phases described in the simplified system Na2O–CaO–MgO– FeO–Al2O3–SiO2 are summarized in Table 1. 2.1. Summary of the petrological and geodynamical numerical approach The numerical procedure used in this work consists of two parts. In the first part, a two-phase flow problem is solved by a set of

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70 Table 1 Reference state thermodynamic properties of melt components and modified thermodynamic properties of the solid components for melting of peridotite in the system Na2O–CaO–MgO–FeO–Al2O3–SiO2.

DH0f ðJ=molÞ

DS0f ðJ=mol KÞ

DV 0f ðJ=barÞa

Melt componentsb SiO2 Al4O6 Mg2SiO4 Fe2SiO4 Ca2Si2O6 NaSi0.5O1.5

921137 3434339 2154031 1490621 3311239 775898

37.3 87.6 103.5 156.3 148.33 78.93

2.44 5.81 4.96 5.72 9.27 3.55

Solid componentsc Ca–Al pyroxene Albite Anortite Spinel Hercynited Pyrope Almandine Grossular

3278767 3916618 4206249 2307313 1952681 6288548 5266502 6633859

259.359 332.927 248.150

DV0f at T = 1673 K for melt components. Mixing solution model is ideal; other thermodynamic properties of melt components from Ghiorso et al. (2002). c Other thermodynamic properties of solid components from Berman (1988); Thermodynamic properties of Hedenbergite from Ghiorso and Sack (1995). d Other properties of Hercynite from Ghiorso and Sack (1995). a

transport equations that allows us to compute the total pressure, temperature and velocity and phase abundance for each dynamic phase (melt and solid mantle). The two-phase model is similar to the one developed in previous studies (Mckenzie, 1984; Scott and Stevenson, 1984; Spiegelman, 1993a; Hewitt and Fowler, 2008) and applied with some modifications to investigate several geodynamic melting problems (e.g. Spiegelman and McKenzie, 1987; Iwamori, 1993; Ruedas et al., 2004; Schmeling, 2006; Sramek et al., 2007; Katz, 2008; Rabinowicz and Toplis, 2009). The set of dynamic equations is the following:

ð1Þ ð2Þ

 @ðqi /i v ci Þ @Bi  þ r  ðqi /i v ci v i Þ ¼ /i rP  r  ð/i si Þ þ @z c¼z @t X þ Mij

ð3Þ

j–i

X @ðq /i C p TÞ i

i

i

@t  T/i ai

þr  ðqi /i C pi T v i Þ  r  ðK i /i rTÞ

DP þ si : rv i Dt

 ¼

X

Chij

ð4Þ

jPi

An additional set of equations is required to describe the chemical transport of the basic components (oxides) in each dynamic phase.

X b @ðqi /i cbi Þ Cij þ r  ðqi /i cbi v i Þ ¼ @t j–i

53

ð6Þ ð7Þ ð8Þ

j¼l

b

X m @ðqi /i Þ þ r  ðqi /i v i Þ ¼ Cij @t j–i  X @ðq /i Þ i þ r  ðqi /i v i Þ ¼ 0 @t i

  @ v xi 2 sxx ¼ 2 m m  n þ j i i i i¼s r  v i i 3 @x  x  @ v i @ v zi sxzi ¼ mi þ @z @x  2  /j mj  ðv j  v i Þ Mij ¼   ki

ð5Þ

The subscript i; j denotes the dynamic phases (s solid, l melt), the superscript c the generic coordinate x; z and the superscript b a chemical component. The primary variables /; T; P; v ; c are volume fraction, pressure, temperature, velocity vector and chemical component, respectively. The other parameters t; q; Ch ; Cm ; Cb ; B; C p ; K are time, density, heat transfer rate, mass transfer rate, chemical transfer rate, body force, heat capacity at constant pressure and thermal conductivity. The auxiliary equations for viscous stress s and momentum exchange M are defined as follow:

where m and n are shear and bulk viscosity, ki is the (scalar) permeability defined as /3l =C k , where C k is the permeability constant. Similarly to a previous study (Ghods and Arkani-Hamed, 2000), the bulk viscosity is assumed to be equal to the shear viscosity. The numerical discretization follows the finite volume procedure and the coupling between pressure and velocity is solved using the SIMPLE method (Patankar, 1980). Numerical diffusion is reduced by a second order flux-corrected algorithm (Boris and Book, 1973; Book and Boris, 1975; Zalesak, 1979). The discretized equations are implicitly solved in dimensional form using a preconditioned biconjugate gradient method (Press et al., 1992). The system of equations is highly coupled, therefore, the solution for all the primary variables is obtained by a iterative procedure formed by two inner cycles, each cycle addressing a partial set of equations. The second part of the formulation solves a chemical equilibrium problem starting from the bulk composition, pressure and temperature determined by the dynamic model. The problem consists of finding the mineralogical assemblage, then minimize the Gibbs free energy of the local system and simultaneously satisfy the mass balance equations (Van Zeggeren and Storey, 1970; Eriksson, 1971; Eriksson and Rosén, 1973; Eriksson, 1975; Smith and Missen, 1991). The solution of the constrained minimization problem is obtained by using the method of Lagrange multipliers (Ganguly and Saxena, 1987; Smith and Missen, 1991). With respect to an earlier version of our model (Tirone et al., 2009), a preliminary assemblage is obtained using the linear programming technique (Van Zeggeren and Storey, 1970) which is also the foundation of the program Perple_X (Connolly and Kerrick, 1987; Connolly, 2005). This method provides an assemblage with only an approximate composition of the phases involved, therefore the initial set of minerals may not represent the equilibrium assemblage and a new starting assemblage must be defined using a different procedure. The output of the thermodynamic computation consists of the equilibrium minerals and melt abundance and composition at each grid point at prescribed time intervals. Density, heat capacity, thermal expansion, heat, mass and chemical transfer rate are also retrieved from the chemical equilibrium computation. These are properties that are transferred and applied to the transport model. Mass and chemical transfer rates are of particular relevance here. These quantities are computed from the difference between the amount of melt and chemical abundance obtained from the transport model and the values from the thermodynamic computation at each location in space. The difference normalized to the time interval at which the thermodynamic computation is performed (2500 a), constitutes the mass and chemical transfer rate that are transferred to the dynamic model for the duration of the next time interval. The time interval between thermodynamic computations is equal or greater than the time step applied to the dynamic model. It has been extended up to 10,000 a in few test runs but we found no significant differences in the overall results of the ridge model. Latent heat of melting is taken into account in the heat transport Eq. (4) following a procedure described in Tirone et al. (2009). 2.2. Thermodynamic model The initial bulk composition of the system used to study the petrology and geodynamics of a mid-ocean spreading ridge is that

54

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

of a homogeneous fertile peridotite. The composition has been chosen considering the uncertainty about the extent of depletion and the distribution of chemical heterogeneities in the asthenospheric upper mantle. Workman and Hart (2005) computed the composition of unmelted depleted mantle that is capable of generating MORB melt and leave a residual solid with composition similar to abyssal peridotites. However, the inversion procedure to compute the reference mantle composition is based on the somehow arbitrary assumption that the source mantle experienced a certain degree of fractional melting. Further work is needed to quantify the effect of different initial compositions or multiple lithologies on the evolution of melt in mid-ocean ridges obtained from numerical simulations. The peridotite composition used in this study is described by a system of 6 oxides (Na2O–CaO–MgO–FeO–Al2O3–SiO2). Only major components are taken into account in order to reduce the numerical and technical problems that arise by interfacing the thermodynamic and the geodynamic models. Nevertheless it has been shown that mantle melting can be effectively reproduced by describing the bulk composition of ultramafic rocks with a reduced number of components (Presnall et al., 2002). The nature and composition of the petrological assemblage in local equilibrium is controlled by the thermodynamic properties of the potential phases involved. The initial assessment of the thermodynamic database used in this study takes into account the original database developed for the pMELTS software package (Ghiorso et al., 2002). Some modifications were introduced to better reproduce the experimental petrological data on peridotite melting in simplified systems (Walter and Presnall, 1994; Walter, 1998; Gudfinnsson and Presnall, 2000; Presnall et al., 2002) and to accommodate the fact that minor components, such as Fe3+, Ti and Cr, which play a significant role in pMELTS database, have not been included in the present database. Experimental (Jaques and Green, 1980) and thermodynamic studies (Asimow et al., 2001) have shown that aluminum in spinel is largely replaced at low pressures (<10 kb) by chromium. Our current model does not incorporate chromium, therefore, to insure that spinel remains stable at low pressure, the stability field of aluminum-rich spinel has been expanded. Few selected thermodynamic properties were optimized by fitting the melt abundance and phase assemblages from the experimental studies using a genetic algorithm optimization procedure (Deb, 2001). Table 1 reports the relevant changes of the thermodynamic properties of solid and melt components. Fig. 1a shows the phase equilibrium diagram for a fertile peridotite computed with the thermodynamic model developed in this work along with the solidus obtained from experimental data (Jaques and Green, 1980; Walter and Presnall, 1994; Gudfinnsson and Presnall, 2000) and the solidus computed with MELTS and pMELTS (Ghiorso and Sack, 1995; Ghiorso et al., 2002) (bulk composition is slightly different). The most evident discrepancy with the solidus determined experimentally is the negative slope when plagioclase and spinel coexist. This is a well known problem of the thermodynamic parameterization (Asimow et al., 2001; Presnall et al., 2002), however its importance is perhaps overestimated. Since melting begins much deeper in the mantle, at lower depths melt contours rather than the solidus is the most relevant information that is needed to model the melting process during upwelling. Fig. 1b shows that the contours of constant melt fraction in equilibrium with pl-free residual solid assume a positive slope. This observation combined with the fact that plagioclase is the first phase to be assimilated in the melt, reinforce the conclusion that the negative slope of the solidus at low pressure has a minor effect on the results of the melting model. As it is shown in the following sections, the bulk composition of the system in a ridge setting changes continuously. These variations have some influence on the corresponding equilibrium phase

diagram. For instance, using the bulk composition observed at x ¼ 100 km, z ¼ 45 km after 10.8 Ma to compute the P,T phase diagram for the depleted peridotite system (see Fig. 2 and discussed in Section 3.1), it can be shown that the solidus and the contours of constant melt fraction are quite different from those of the fertile mantle (Fig. 1b). In particular, it is quite clear that the pressure range over which the slope of the solidus assumes a negative slope and the stability field of plagioclase + spinel are greatly reduced. 2.3. Model setup and parameters adopted in this study The dynamic model consists of a 2-D simulation box (100  200 km) that reproduces the half-space of a mid-ocean ridge. A prescribed spreading velocity on the surface is set to 4 cm/a. The imposed surface velocity implies that ridge spreading is mainly controlled by passive flow. However the buoyancy effect is incorporated in the dynamic model since the thermodynamic computation provides the density of the residual solid and melt at every location. Reference temperature at the bottom (depth = 100 km) is 1470 °C. Reference temperature rather than potential temperature is used hereafter. The concept of potential temperature is not quite applicable to problems where the dynamic of melt is not known a priori (i.e. thermal gradient with depth is not quantifiable in advance). Furthermore, in the presence of a cooling lithosphere, the projection of the potential temperature does not clearly define the mantle temperature everywhere. The assumption of a fixed reference basal temperature implies that, at a certain depth in the mantle, the same temperature is maintained over an horizontal distance that extend for several hundreds of kilometers from the ridge axis and is kept constant for the duration of the simulation (several millions of years). It is a constraint that cannot be removed without expanding the dimension of the ridge model to include the global mantle flow circulation. 3-D ridge models have been used to investigate the interaction of melt with transform faults (Morency et al., 2005; Gregg et al., 2009). The objective of this work is to study an ideal ridge far from the transform margins, and for this purpose, a 2-D model can provide a sufficiently accurate description. Viscosity of melt is set to 1 Pa s. Mantle viscosityis described  by the following relationship: 1 5 m ¼ 1:0214  1010 exp 310 RTJ mol , where R is the universal gas constant and T is the temperature. An upper threshold (1022 Pa s) is imposed to the viscosity model. The pre-exponential term is determined by assigning a reference viscosity value (1019 Pa s) at the reference temperature. The range of viscosity values obtained by the above expression is consistent with estimates for oceanic mantle on the basis of experimental and geophysical evidences (Hirth and Kohlstedt, 1995, 2003). Several viscosity models which included the effect of water and melt have been used to model the oceanic mantle under spreading ridges. Morgan (1987) proposed that the effect of melt is to substantially increase the viscosity of the asthenospheric mantle. Alternative models, supported by experimental studies (for a review, see Kohlstedt et al., 2000), suggest that the effect of melt is to decrease the viscosity (Su and Buck, 1993; Kelemen et al., 1997; Braun et al., 2000; Schmeling, 2000). The effect of water was also included in Braun et al. (2000) and Hall and Parmentier (2000). Nonetheless, assessing the influence of volatiles and melt on the rheology of mantle rocks in a mid-ocean ridge is not a straightforward task. Water content and distribution in a ridge setting are not clearly known. The effect of melt depends on the geometry and spatial distribution of melt (Schmeling, 2000) and grain size (Zimmerman and Kohlstedt, 2004) which are difficult to constrain in the oceanic mantle. It was also suggested that the influence of melt on the viscosity of mantle rocks becomes significant for melt content greater than 5% (Hirth and Kohlstedt, 1995; Zimmerman and Kohlstedt, 2004). As it will be shown later,

55

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

T (oC)

(a) 1150

1200

1250

1300

T (oC)

(b) 1350

1400

1450

1200

5

1300

1400

1500

1600

5 Pl+Sp

Pl+Sp Pl+Sp

10 Pl+Sp Pl+Sp

10

CMASF (experimental) Gudfinnsson & Presnall

Pl+Sp

Sp

P (kbar)

ed et

P(kbar)

pl de ot

ri d pe

rid pe

le rti fe

20

MELTS

us

us

20

id

CMASN (experimental) Walter & Presnall

l so

Sp

15

id

p-MELTS

l so

Sp

15

it e

ot

25

it e

Pyrolite(experimental) Jaques & Green

Sp+Gt Gt

25

Sp+Gt Gt

CMASNF (computed) this work

1% 5%

1% 5%

Fig. 1. (a) Solidus for a fertile peridotite computed with the thermodynamic database developed in this study. For comparison the following melting curves have been included: experimental results from Walter and Presnall (1994) (CMASN), Gudfinnsson and Presnall (2000) (CMASF) and Jaques and Green (1980) on a natural peridotite, solidus obtained with the program MELTS (Ghiorso and Sack, 1995) and pMELTS (Ghiorso et al., 2002). Composition is slightly different in all cases. (b) Comparison between solidus of a fertile peridotite and a depleted peridotites. Depleted peridotite composition is the bulk composition at x = 100 km, z = 45 km after 10.8 Ma (Fig. 2): (Na2O = 0.114, CaO = 2.614, MgO = 41.078, FeO = 8.278, Al2O3 = 3.111, SiO2 = 44.803 in wt%). Two contours are also shown for 1% and 5% melt content. Composition of the fertile peridotite is reported on caption of Fig. 5.

melt fraction greater than 5% is locally observed only near the ridge axis at low depths and at the base of the lithosphere. Further viscosity models were inferred from the axial topography assuming two descriptions of the internal structure of the ridge (Eberle et al., 1998a; Eberle and Forsyth, 1998b). However, these viscosity models pre-assume a defined ridge structure, therefore they cannot be really applied to this study. In summary, we have adopted the simplest viscosity model that retains the essential minimum requirement, i.e. to be able to describe a diffusion creep thermally activated process in the asthenospheric mantle and develops a rigid (impermeable) lithosphere. The permeability constant C k , which plays a major role in the velocity transport equations (Eq. (5)), is set to 109 m2. Lower and higher values, within the estimated range at mantle conditions (Spiegelman, 1993b; Rabinowicz and Toplis, 2009), have been used here mainly to show the effect on melt distribution along the ridge vertical axis. Initial conditions inside the model box are defined by uniform temperature (1100 °C) and uniform bulk composition (Na2O = 0.41, CaO = 3.15, MgO = 38.4, FeO = 8.3, Al2O3 = 4.54, SiO2 = 45.2 in wt%). As discussed in Section 2.2, the initial bulk compostion represents a fertile peridotite. Water is not included in the model. Although the effect on oceanic basalts has been already quantified in a petrological model (Asimow and Langmuir, 2003), a comprehensive petrological and geodynamical study of mid-ocean ridges including the presence of volatiles would require a model which incorporates a third dynamic phase, i.e. free water. The initial temperature setup in this study has no particular relevance, except that it provides a deformable lid dynamically driven by the imposed passive flow and essentially impermeable to melt flow. Grid spacing is set to 0.5 km in the vertical direction, starting from 50 km depth up to the surface, and in the horizontal direction for the first 50 km from the ridge axis. In the rest of the simulation box, spacing is equal to 1 km. A 1-D melting column is used in Appendix A to show the influence of smaller grid size on melt distribution, composition of the residual solid and extent of chemical equilibration (discussed in Section 4.1).

3. Results This study presents the initial results, in particular the spatial and temporal distribution of melt beneath ridges, crust production and variations of the bulk composition (solid + melt) of major elements. A comparison is also made between the composition of the residual solids at various depths and available data on major elements in abyssal peridotites. The model gives an indication of a physically consistent picture of the evolution of the residual mantle and, at the same time, provides the basis to understand the extent of chemical equilibration and the large chemical variations found in abyssal peridotites. The seismic structure for shear waves is computed from the simulated ridge model. This type of numerical simulations provides a substantial amount of information. A specific choice was made to discuss only certain results here and present in a forthcoming study additional findings on how melt composition and melt productivity vary as a function of the basal temperature and spreading rate. 3.1. Overview of the temporal and spatial distribution of melt and evolution of the solid mantle Fig. 2 presents the numerical results showing the melt distribution beneath a spreading ridge. The two half-space plots of panel (a) illustrate the melting region at two slightly different times, 4.45 and 4.47 Ma. Time zero is arbitrarily set with the arrival of the first melt in the near surface area. The first melt starts to form at about 95 km depth, where the upwelling mantle intersects the peridotite solidus in the garnet stability field. Moving upwards, melt productivity temporarily decreases when spinel becomes stable. Similar behavior of the melt productivity has been observed using the thermodynamic model of MELTS (Asimow, 1995). An opposite effect and an increase of melt productivity can be noticed in a narrow region at 20 km depth, which correspond to the location where clinopyroxene is exhausted. White streamlines on the left side of Fig. 2a and dark streamlines on the right side of

56

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

(a)

(b)

(d)

(c)

Fig. 2. (a) Spatial distribution of melt beneath a mid-ocean ridge. Reference temperature at the bottom is 1470 °C, spreading rate is 4 cm/a. The two sides of (a) are separated by a short time interval, 4.45 Ma, 4.47 Ma. Time is arbitrarily defined with respect to the arrival of the first melt near surface. White streamlines on the left side refer to melt flow, dark streamlines on the right side refer to solid mantle flow. The dotted line on the left side marks the limit above which no more melt is formed. Computation is performed on a half-size simulation box around the ridge axis, 200  100 km. Initial uniform temperature is 1100 °C. Grid size: Dt ¼ 500 a, x = z < 50 km, Dx ¼ Dz ¼ 0:5 km; x > 50 km, Dx ¼ 1 km. Boundary conditions are summarized on the left side of the figure. (b) Beginning of the separation of the melt region after 6.3 Ma. Streamlines for the solid mantle are also shown (c) incipient separation of a melting area away from the ridge axis and the deflection of the mantle flow after 10.8 Ma. Red contour lines in (b) and (c) indicate the melt fraction (b) and (c) include the total extension of the simulation box. (d) Melt distribution for a polybaric batch melting model ðv melt ¼ v solid Þ. Other numerical conditions are the same as for the model in (a). Dark streamlines refer to solid mantle and melt flow.

Fig. 2a indicate the instantaneous directions of melt and mantle flow, respectively. Vertical velocity of melt and solid mantle near the ridge axis progressively diverges as the melt move upwards. For instance, at 60 km depth, with 2% melt fraction, velocity of melt is about twice the velocity of the residual mantle. Different velocities for melt and residual solid imply that the conventional batch melting approach is not quite appropriate to describe the

melting process in a ridge setting. This is shown in Fig. 2d which illustrates the melt distribution for a polybaric batch melting model obtained applying the same conditions and parameters used in the main model. The only difference is that melt is forced to move with the same velocity of the solid mantle. The melting region, developed as a coherent unit, reaches the maximum width of 150–160 km after 4.88 Ma. This lateral

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

extension is reduced over time and approaches approximately a steady width of 100 km from the axis after 10 Ma (Fig. 2c), while, at greater distance from the spreading center, a certain amount of melt is still present. Lateral extension of the melting region presumably changes as a function of the influx of mantle material from deep regions, spreading rate, viscosity and possibly the reference temperature at the bottom of the simulation model. A first order comparison of the size of the melting region can be made with the East Pacific rise which has been extensively studied using geophysical methods. Even though the ridge is characterized by spreading plates with asymmetric velocities and the absolute spreading rate is slightly higher than the one imposed in this study, the overall depth and width of the melting region obtained by the numerical model appear to be consistent with those derived from seismological observations in the MELT area (Forsyth et al., 1998; Spahr and Forsyth, 1998; Toomey et al., 1998). These studies inferred that the depth of initial melting ranges between 100 and 200 km and the total lateral extension of the melting area is 200–300 km. Large width of the melting zone, like the one detected in East Pacific rise, is believed to be associated to plate spreading driven mainly by passive flow rather than active mantle upwelling, which agrees with the dynamic condition imposed in this study. As the ridge model evolves, a portion of the melting region appears to progressively detach from the main melting zone (Fig. 2b and c). The effect is created by the weak coupling between the relatively low viscous asthenospheric mantle and the highly viscous rigid lithosphere. This portion of the astenosphere that includes also small amount of melt, persists on the left side of the simulation model, away from the ridge axis. We can reasonably assume that, at some further distance from the ridge center outside the simulation model, melt eventually will freeze by heat dissipation. A similar type of evolution of the melting region has been observed in a different simulation with a maximum lateral box size of 350 km and using a uniform mesh spacing equal to 1 km. Some preliminary results suggest that the separation of a melting area from the main melting region and the lateral translation may not occur only in the early stage of the plate separation, but the process repeats itself over time. A combination of magnetotelluric and seismic data also provides evidences of the melt distribution under the southern East Pacific rise (Evans et al., 2005). A region of high conductivity located between 60 km and 130 depth and extending for over 300 km from the ridge axis was interpreted as the result of partial melting and indicative of the presence of water. Although water has not been considered in this study, persistence of melt at some distance from the spreading center is consistent with our numerical results. However, a high resistivity area above 60 km depth, extending for several hundred kilometers was also found under the East side of the Pacific rise (Evans et al., 1999; Baba et al., 2006). It was suggested that the sharp change in resistivity represents the boundary between a dry mantle above, and a wet mantle below 60 km. Experimental data shows that gabbroic rocks have higher conductivity than peridotites (Maumus et al., 2005). Therefore, the high resistivity above 60 km depth may be possibly related to the emplacement of large igneous bodies (gabbros) or underplate crystallized melt. Such alternative interpretation would be in better agreement with the numerical results of this study. Two implications that should be valid for passive ridges, follow from the off-axis dynamics of the melting area. Away from the ridge axis, the lower section of the asthenospheric mantle, after being depleted by the melting process, deflects downwards near the bottom side of the simulation model (Fig. 2). Therefore, it appears that the evolution of a spreading ridge not only modifies the chemical composition of the asthenosphere during melting, but also creates a dynamic mechanism that introduces a depleted residual solid in the large upper mantle flow. Previous numerical

57

models (Rabinowicz et al., 1984; Scott and Stevenson, 1989) share some similarities with the mantle flow pattern beneath spreading ridges observed in this study. It was shown that mantle circulation describes a close loop around the ridge axis, essentially triggering a diapiric effect and recycling of material near the ridge axis. Closed circular flow cannot be confirmed from this study since it would require to extend the model to greater depth (>100 km). However, the particular mantle flow around a spreading center observed in these earlier models could have a significant implications on the chemical evolution of melt and oceanic mantle over time and it may require further investigation. A qualitatively similar mantle evolution induced by buoyant melt instabilities was also observed in a setting describing the spreading of a persistently thick continental lithosphere (Hernlund et al., 2008a,b). The second implication that follow from our numerical model is that, the lithospheric mantle cooling away from the spreading center, not necessarily incorporates material that passed through the ridge melting area. Instead, the lithosphere would increase the thickness most likely by cooling the ambient mantle that is in proximity of the base of the growing lithosphere. Therefore, in a hypothetical vertical section, the composition of the lithosphere should not reflect higher degrees of depletion with decreasing depth, as it would be expected for a mantle column that passed entirely through a spreading ridge. According to our model, this should be the case for lithospheric material at depth greater than 50 km (Fig. 2c). An investigation of lithospheric xenoliths from the Hawaiian volcanic chain concluded that the portion of the oceanic lithosphere represented by the nodules has been involved in the melting process related to the East Pacific rise (Yang et al., 1998). This evidence apparently contradicts the conclusion above. However, geobarometry estimates indicate that the equilibration depth of the samples corresponds to a pressure varying between 7 and 21 kbar (Yang et al., 1998). Therefore, the petrographic indications provided by the xenoliths in Hawaii do not rule out the possibility that only the upper portion of the lithospheric mantle passed through a mid-ocean ridge stage. Both flow regimes, downflow of portions of the mantle after being affected by the ridge melting process and inflow of mantle unrelated to the melting process, have been described in a previous study (Scott and Stevenson, 1989). It was shown that, for a case of flow circulation driven by plates and buoyancy, the two flow patters occur in a dynamic sequence that leads to episodic interruptions of the crust production. The problem with the crust productivity may have been related to the simple melting model adopted in the previous study. By using a more realistic petrological melting model, it appears that the mantle can develop both flow patters while, as it will be shown later, the crust formation is periodic but continuous in time. 3.2. Melt focusing and distribution near the ridge axis Several physical mechanisms for focusing of melt from off-axis regions have been proposed in the past (e.g. Su and Buck, 1993; Spiegelman and McKenzie, 1987; Morgan, 1987; Sparks and Parmentier, 1991; Aharonov et al., 1995; Hall and Parmentier, 2000; Morgan and Holtzman, 2005; Katz et al., 2006). The numerical results of this work show that melt is collected towards the surface directly from the asthenospheric mantle only within 10 km from the spreading ridge. Further away, within a distance of 50 km from the ridge axis, the ascending melt tends to accumulate at the base of the lithosphere, then it is extracted to the surface by sliding along the shortest pathway under the lithosphere. At distances greater than 50 km, melt concentrating under the lithosphere, follows the mantle flow and then contributes to the formation of sub-lithospheric intrusive igneous bodies. This melt behavior is essentially similar to the one discussed in earlier studies

58

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

(Sparks and Parmentier, 1991; Spiegelman, 1993a). From the investigation of shallow depth seismic Love waves in the East Pacific rise (Dunn and Forsyth, 2003), it was suggests that melt may accumulate at the base of the lithosphere; our numerical model appears to confirm this conclusion. The presence of magma cumulates, magma lens and sills detected from seismic P and S waves data at the base of the lower crust and in the upper mantle near the ridge axis (Dunn et al., 2000; Nedimovic´ et al., 2005) could be also consistent with the melt distribution and extraction path observed in this study. The upwelling melt velocity is mainly controlled by melt buoyancy, porosity and permeability (Eqs. (3) and (5)). Unfortunately, permeability induced by melt flow in the upper mantle is not very well constrained. The main results of this study have been obtained assuming a permeability constant C k equal to 1  109 m2 , which is within the range (3  108  6  109 m2 ) estimated from experimental measurements and other considerations (Rabinowicz and Toplis, 2009). Ruedas et al. (2004) assumed a value of 6:5 108 m2 for plume melting, while Ghods and Arkani-Hamed (2000) in their study of melting under a spreading ridge, considered a permeability constant equal to 1  107 m2 . A recent experimental study concluded that permeability of melt in the asthenosphere could be two orders of magnitude higher than previously assumed (Connolly et al., 2009). If confirmed, such high permeability has significant implications on the dynamics of melt but also on the petrological evolution of mid-ocean ridges. Our model shows pulses of magma, starting at about 25–30 km depth in the region near the ridge axis (Fig. 2a). This depth corresponds approximately to the location where clinopyroxene is completely dissolved in the magma. It has been suggested that chemical transformations such as garnet-spinel or spinel-plagioclase transitions, that induce changes in the melt productivity, may have an effect on melt transport in the upwelling mantle (Asimow, 1995). It is shown here that, solid-melt dissolution reactions may also induce variations in the melt distribution. Porous waves create oscillations of the melt content as large as  1% around the average melt of 4%. These fluctuations are essentially controlled by the deformation of the solid matrix which is allowed to compact/decompact (Mckenzie, 1984; Spiegelman, 1993a). Variations of the magma content do not seem to form in numerical simulations with grid spacing greater than 0.5 km. This is approximately on the same order of the compaction length in the upper part of the melting area. In Appendix A the effect of using more refined numerical mesh is investigated by simulating the magma formation and upwelling in a pseudo 1-D melting model. The important observation from the 1-D simulations is that melt instabilities start to form much deeper in the melting column with the decrease of the grid size. This is because the compaction length decreases at increasing depth and only a mesh with smaller size is able to resolve porous melt variations at the scale defined by the compaction length. However, despite the fact that the depth range and amplitude of the melt fluctuations are not precisely reproduced in the 2-D model, no detectable variations of the composition of melt and residual solid have been observed at the locations where melt fraction varies between its maximum and minimum (see also Appendix A). The amount of melt present in the upwelling mantle depends, among other factors, on the permeability of the solid matrix. Fig. 3 illustrates the melt distribution along the vertical axis of the ridge from three different numerical models that assume a permeability constant of 7  108 ; 1  109 and 5  109 m2 (Fig. 3a–c). Fig. 3b refers to the same simulation shown in Fig. 2. Except for the case with a permeability constant of 5  109 m2 , in which no fluctuations of the melt content are observed (Fig. 3c), for the other two cases shown in Fig. 3a and b, melt develops oscillations

over time with different periodicity in response to a different compaction length. 3.3. Crustal production and underplate crystallized magma The calculated crustal thickness as a function of time is reported in Fig. 4a. The crustal production is obtained normalizing the solidified melt, within a rectangular area (10  7 km) next to the ridge axis with the spreading distance per unit of time. Similarly to a previous study (Ghods and Arkani-Hamed, 2000), the rectangular area represents the location where the extraction of melt to the surface takes place. The amount of melt that solidifies is retrieved from the thermodynamic chemical equilibrium computation which is performed at determined time intervals. Frozen melt is not extracted from the model but is simply transformed into a solid phase that is merged with the existing solid mantle already present in the location where melt crystallizes. The mechanism is quite efficient to remove melt from the model. However, it is unlikely that compositional changes of the crustal material are accurately reproduced since the composition of the new solid phase is a weighted average of the composition of the pre-existing residual solid and the frozen melt at a particular location. A period of large crustal production is observed initially when the first asthenospheric mantle reaches the near surface area. After this initial period, the amount of crust formed over time begins to develop temporal oscillations with a periodicity of about 0.2– 0.3 Ma. Similar temporal fluctuations, but with a much shorter period, were observed in an earlier model of ridge melting (Ghods and Arkani-Hamed, 2000). It appears that the mechanism responsible for the variations of the crust thickness is related to the melt production away from the ridge axis. The process of accumulation under the lithosphere off-axis is continuous in time, but only when melt reaches a sufficiently large amount (>10%), melt begins to move towards the near surface central area, following a path defined by the shape of the rigid lid (Fig. 4b). This result is in agreement with an earlier study suggesting that melt, accumulated in the mantle, is quasiperiodically discharged through the rift-axis crust (Tucholke et al., 1997). The propagation of magma waves in the upper portion of the model near the ridge axis, instead, does not appear to influence the instantaneous crust production. This is explained, in part, because the variation of the melt content induced by the magma waves is rather small, but also because melt ascending from below accumulates in a magma pool (melt abundance 10–15%) at 10 km depth, then it is extracted and solidifies at a different rate than the rate of melt input from below (Fig. 4b). The magma pool may be interpreted as the numerical proxy of a magma chamber that is usually believed to exist under intermediate and fast spreading ridges (e.g. Sinton and Detrick, 1992). The average amount of melt that contributes to the formation of the crust on the surface is 60% of the total solidified melt (Fig. 4). The rest of the melt freezes for the most part at the base of the lithosphere outside the rectangular box which defines the extraction area where crustal material is formed. Ghods and Arkani-Hamed (2000) found in their model that the amount of crust is slightly lower than in our study, about 50% of the total solidified material. The average crustal thickness of about 6 km is within the inferred range for the given spreading rate (Bown and White, 1994; Cannat, 1996). Oscillations of the crust thickness with a period of 2–7 Ma have been suggested for slow spreading ridges (<2 cm/a) based on gravity measurements on the Mid-Atlantic ridge (Maia and Gente, 1998; Bonatti et al., 2003; Cipriani et al., 2009). However from the seismic structure of a portion of the Mid-Atlantic ridge, it was inferred that variations of the crust depth are instead much more rapid (0.4–0.8 Ma) (Canales et al., 2000). No evidences have been presented yet on the variations of the crustal production for

59

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

(b)

(c)

10

10

15

15

15

20

25

30 0.02

t = 3.36 Ma t + 0.02 Ma t + 0.03 Ma t + 0.05 Ma Ck = 7x10 0.03

0.04

depth (km)

10

depth (km)

depth (km)

(a)

20

25

t = 4.41 Ma t + 0.02 Ma t + 0.04 Ma t + 0.06 Ma

8

Ck = 1x10

30 0.03

melt fraction

0.04

t = 2.75 Ma

25

9

0.05

melt fraction

20

30 0.05

Ck =5x10 0.06

9

0.07

melt fraction

Fig. 3. Vertical distribution of melt along the ridge axis at different time intervals for models that assume three permeability constants C k : 7  108 ; 1  109 and 5  109 m2 . Melt distribution at the initial time (t) is illustrated by dotted lines in (a) and (b), solid, dotted, dash-dotted and dashed lines indicate the melt fraction at different times. (c). Melt fraction at one arbitrary time, no periodic oscillations are observed in this case. Time (t) is different in (a), (b) and (c).

intermediate spreading ridges (i.e. 2–6 cm/a) more similar to the spreading ridge in our study. Nevertheless, it is interesting to note that the amplitude of the crustal thickness variations obtained in our model is remarkably similar to the range observed in the Mid-Atlantic Ridge (0.5–1.0 km) (Cipriani et al., 2009). 3.4. Lateral variations of bulk composition, melt and temperature The evolution of the mantle, away from the spreading center, is not only determined by the local mantle flow, but also by the transport of melt which influences, in particular, the variations of the chemical composition with depth and distance from the spreading center. The dotted line on the left side of Fig. 2 identifies the limit between a lower area where melt is formed locally and an area above where supply of melt from below does not compensate for the partial solidification and the net local melt content decreases. The depth defining the threshold between local solidification and melting increases away from the ridge spreading center, which means that partial crystallization begins deeper in the mantle. Although this representation is affected by the imposed local chemical equilibrium between melt and solid within a volume defined by the numerical grid size, it provides an indication of the general relationship that interconnects the solid mantle and melt abundance in a spreading ridge. The idea that partial solidification of melt may take place in the mantle beneath a spreading ridge is not new. Several studies, from various petrological considerations, pointed out that partial crystallization begins below the lower crust in the mantle (Grove et al., 1992; Michael and Cornell, 1998; Herzberg, 2004). For example the analysis of cpx-rich harzburgites at the base of the mantle section of Oman ophiolites suggested that the suite was affected by partial freezing of melt at the boundary between the lithosphere and the asthenosphere (Godard et al., 2000). Bulk composition (melt + solid mantle) along the ridge axis and away from the spreading center changes in response to different flow patterns of melt and solid mantle. This implies that a single P, T phase diagram of a peridotite with fixed bulk composition may not be sufficient to describe the petrology of mantle melting.

The bulk chemistry provides also information on the composition of the mantle material that, as discussed earlier, would deflect downwards and is introduced into the global mantle circulation. Fig. 5 reports the bulk composition, melt fraction, rate of mass transfer between melt and solid and temperature as a function of the distance from the ridge axis at two different times (4.45 and 10.08 Ma). Bulk composition is simply computed by adding the composition of the solid mantle and melt in proportion to their abundance at every location in the simulation model. Because, in general, the degree of melting increases moving upwards, two depths, 30 km and 45 km, are chosen to illustrate the effect of the melting process on the bulk composition and the other quantities. After 4.45 Ma the maximum extent of the melting area is about 150 km. The thermal boundary between the lithosphere and the spreading asthenospheric mantle marks also the sudden compositional change which occurs at 120 from the ridge axis at 30 km depth, and 150 km from the ridge axis at 45 km depth (Fig. 5a–f). After 10.08 Ma, the area covered by the spreading mantle is much wider. Melt content and bulk composition show an inversion in their pattern at 120 km from the ridge axis (Fig. 5a–g), which reflects the position dividing approximately two melting regions, as also shown in Fig. 2a and b. The decrease of the melt content with distance is almost linear from the ridge axis, extending up to 60 km after 4.45 Ma and 120 km after 10.08 km (Fig. 5g). After 4.45 Ma, between 60 km and 130 km from the ridge axis, the slope of the melt content decreases in response to a nearly flat temperature profile (Fig. 5j), which reflects the change of direction of the mantle flow from predominantly vertical to almost horizontal (Fig. 2a). Fig. 5h, shows the local rate of melting and crystallization that were defined in Section 2.1. After 4.45 Ma, local melt formation ceases at about 130 km from the ridge axis. This distance corresponds, approximately, to the maximum width over which melt is found at this particular time. After 10.08 Ma, melt keeps forming at distance as great as 90–100 km, although melt can be found at even greater distance. Beyond 100 km from the axis, melt is in a

60

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

(a)

8

7

6

thickness (km)

5

4

3

instantaneous crust average crust

2

average solidified melt underplate 1

0

0 0

10

20

30

40 km

5 10 15 20

2

4.30 Ma

3

0

depth (km)

depth (km)

(b)

1

4

5 time (Ma)

10

20

6

7

30

40

8

9

0

0

0

5

5

depth (km)

0

10 15 20

4.45 Ma

10

10

20

30

40

% >8 0.05

0.025

10

0

15 20

4.65 Ma

Fig. 4. (a) Instantaneous temporal variation of the crustal thickness computed inside a squared area (10  7.5 km) near the surface around the ridge axis. Average crustal thickness over time and average amount of melt solidified underplate outside the squared area are also illustrated in this figure. Square black dots refer to the temporal location of the frames shown in (b). (b) Three frames illustrating the periodic discharge of melt accumulated at the bottom of the lithosphere in proximity of the ridge axis. Black streamlines refer to the flow of melt in this area of the model.

quasi steady state and only small amount of crystallization is observed. This melt presumably will survive for some time further away from the ridge center and eventually will solidify in the asthenospheric mantle. It is interesting to observe that starting from 40 km from the spreading center, the amount of local melt production is essentially the same at 30 and 45 km depth despite the fact that the temperature at these two depths differs by about 30 °C, which clearly indicates the effect of the bulk composition on the melting process. The general behavior of the major elements in the bulk system is mainly controlled by the residual mantle composition and it should be in agreement with the results from experimental data on peridotite melting (e.g. Takahashi, 1986, 1993; Baker et al., 1995; Walter and Presnall, 1994; Walter, 1998) (earlier studies are summarized in Klein and Langmuir, 1987). At 45 km depth, the higher content of Na2O, CaO, Al2O3 and SiO2 and lower content of MgO and FeO with respect to the chemical content at 30 km depth (Fig. 5a), is consistent with the pressure effect on melting. Decrease of Na2O, CaO, Al2O3 and SiO2 and increase of MgO and FeO with distance from the spreading center are related to the incompatible/compatible behavior of these elements with the degree of melting which, overall, decreases away from the ridge

center (Fig. 5h). The general trend for the chemical components in the bulk system, however, shows some variations that reflect specific changes of the melt abundance and melt production. This can be observed, for instance, after 4.45 Ma, between 60 and 130 km from the ridge axis. Closer to the ridge axis, where the mantle flow is almost sub-vertical, the bulk composition at a particular location is the result of the local melt production and the input of melt from below. In the off-axis region, when the mantle flow becomes sub-horizontal, temperature does not vary significantly with the distance from the axis, and the bulk composition is mainly modified by the contribution of the melt that ascends from deeper regions. After 10.08 Ma, at 120 km from the spreading center, the change of slope of the melt content marks also a rapid variation of the composition of the system, in particular Na2O, CaO, MgO and Al2O3 (Fig. 5a–c, and e). Fig. 5j illustrates the temperature variation with distance in the asthenospheric mantle at 30 and 45 km depth. The overall temperature difference between these two depths is determined by the adiabatic thermal gradient. At 30 and 45 km depth, the variation of the temperature with the distance from the ridge axis is consistent with a previous estimate (0.05–0.1 °C/km) (Cochran, 1986). However Fig. 5j also shows that the rate of the temperature change is not the

61

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

depth = 30 km, t = 4.45 Ma depth = 45 km, t = 4.45 Ma depth = 30 km, t = 10.08 Ma depth = 45 km, t = 10.08 Ma

0.2

(b)

CaO (wt%)

Na2O (wt%)

0.25

0.15 0.1

(c) 3.0

47.0

2.5

45.0

MgO (wt%)

(a)

2.0 1.5 1.0

43.0 41.0

0.5

0.05

39.0

0.0 50

150

0

8.9

3.5

8.8

3.0

8.7 8.6 8.5 8.4

100 km

150

0

2.5 2.0

150

50

(h) -1

3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0

100 km

50

100 km

150

0

50

(j)

100 km

150

1420

solidification melting

-0.05 -0.10 -0.15 -0.20 -0.25 -0.30

0

43.0

150

0.05 0.00

44.0

42.0 0

rate of mass transfer (Ma )

(g)

100 km

150

1400

o

50

temperature ( C)

0

100 km

45.0

1.0

8.2

50

(f)

1.5

8.3

melt wt%

50

(e) Al2O3 (wt%)

FeO (wt%)

(d)

100 km

SiO2 (wt%)

0

1380 1360 1340 1320

0

50

100 km

150

0

50

100 km

150

Fig. 5. (a)–(f) Variation of the bulk chemical composition (solid + melt) with distance from the ridge axis. (g), (h) and (j) Variation of melt content, solid–melt transfer rate and temperature, respectively. Lateral variations ar shown at 30 km and 45 km depth and two different times, 4.45 and 10.08 Ma. Initial bulk composition is (wt%): Na2O = 0.41, CaO = 3.15, MgO = 38.4, FeO = 8.3, Al2O3 = 4.54, SiO2 = 45.2.

same at these two depths. In fact, the temperature gradient is smaller at greater depth. The different thermal behavior with distance observed at 45 and 30 km depth can be related to the extent of the local melting production and the effect of the latent heat of melting. 3.5. Residual solid composition along the ridge axis and relation with abyssal peridotites Abyssal peridotites are commonly found in proximity of sections of the lower crust and upper mantle tectonically exposed at the bottom of the ocean floor (e.g. Michael and Bonatti, 1985; Shibata and Thompson, 1986; Dick, 1989; Seyler and Bonatti, 1997; Niu, 2004). The direction of the mantle flow beneath spreading ridges and the location of the peridotite in relation to the lower crust suggest that this portion of the mantle must have followed an upwelling path in close proximity of the axis of the oceanic ridge, regardless of the velocity of the spreading ridge or the petrological and geochemical processes involved. This observation should be valid, unless the ridge axis translated with respect to a fixed reference frame, for instance when plates spread with different velocities. Experimental studies on subsolidus phase equilibrium show that plagioclase is stable up to 8–11 kb (Presnall et al., 2002;

Borghini et al., 2010), but the large majority of abyssal peridotites do not include plagioclase. This appears to contradict the common perception that chemical equilibration is achieved at any depth during the evolution of the mantle beneath a spreading ridge. For a peridotite initially in the spinel facies, considering a scenario in which the spreading rate is 10 cm/a, and direct upwelling along the ridge axis from 25 km depth, it would require about 0.25 Ma for the solid mantle to reach the near surface level, a time sufficient to re-equilibrate the peridotite in the plagioclase facies. To explain the petrological and geochemical evidences associated to abyssal peridotites, several petrological models have been put forward using different approaches. However, these models necessarily ignore the spatial and temporal variations of melt production as well as the interaction of melts from different depths that can only be described by a combined petrological and dynamical approach. Niu et al. (1997) proposed that abyssal peridotites are the end-product of a melting process that depleted the mantle, followed by a stage in which the residual peridotite assemblage is modified by precipitation of olivine. Although the model reproduces the major elements composition of abyssal peridotites, it does not take into account the equilibrium and kinetic experimental data on melt-rock interaction (e.g. Wagner and Grove, 1998; Morgan and Liang, 2005), or considers the consistency of the

62

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

petrological model with the thermal evolution of the mantle and melt transport inferred from geodynamic studies. Using a thermodynamic approach (Asimow, 1999), it was suggested that, in deep regions, fractional melting develops a melt distribution that resemble an inverse-tree-type structure. At lower depths, instead, magma focuses in large channels and melting is better described by a polybaric batch process. While this model, as the previous one, is consistent with the major elements pattern of abyssal peridotite, it implies a certain type of melting flow that has not been observed in dynamical models. The tree-type ramifications of the melt distribution and the sequence fractional-batch melting from the bottom to the top would require very particular assumptions on the permeability regime. More specifically near fractional melting is associated to very fast melt transport that usually can be obtained when the porosity or melt fraction is very high. It is unlikely that this requirement would be satisfied in the deep regions of the melting zone at conditions similar to the one imposed in our model. In the model presented by Lundstrom et al. (2000), melt is initially segregated from enriched veins in an heterogeneous mantle. Melt is continuously re-equilibrating by dissolution and precipitation with the surrounding peridotite which turns, at some point, into dunites. Highly porous dunites become then the main conduit for melt. At this stage, the melting process becomes fractional and the residual lerzholite acquires the geochemical signature of abyssal peridotites from diffusive interaction with melt and dunites. In the last stage, melt moves rapidly towards the surface in complete chemical disequilibrium. This elaborate model relies on several assumptions that aim to accommodate a large variety of geochemical information. Assessing whether the scenario is physically plausible would require a numerical geodynamic model that is able to develop lherzolites and dunites with different porosity (e.g. Liang and Parmentier, 2010) and incorporates a full description of the complex petrological interactions in a heterogeneous mantle. Presnall et al. (2002) proposed that MORB are formed by melting within a short pressure interval, 10–15 kb, near the plagioclase-spinel transition. At the end of the process, some portion of the residual mantle will eventually become part of the abyssal peridotites that are sampled on the ocean floor. Composition of abyssal peridotite is consistent with a deep localized melting process, but it is unclear why the melting process ceases at the spinel-plagioclase transition, and how melt is transported near surface from at least 30 km depth without interacting with the surrounding asthenospheric mantle and without modifying the residual mantle composition. The model presented in this study determines the composition of the mantle residual solid and melt self-consistently with the geodynamical evolution of the mid-ocean ridge. Fig. 6 shows a plot of major elements in the residual solid with respect to MgO content along the ridge axis after 4.45 Ma. The chosen time is not particularly relevant. After the thermal upwelling has reached the near surface area and the first melt is added to the new forming crust, the composition of the residual mantle at any depth along the ridge axis does not appear to change significantly over time. The red square symbols connecting the continuous lines in Fig. 6 represent the computed mantle composition at various depths; higher MgO content is indicative of lower depth, for instance MgO = 46.3 wt% corresponds approximately to a depth of 23 km. For comparison, selected data on the composition of abyssal peridotites, mostly in the spinel facies, are also included (Dick and Fisher, 1984; Michael and Bonatti, 1985; Dick, 1989; Johnson and Dick, 1992; Snow and Dick, 1995; Seyler and Bonatti, 1997; Baker and Beckett, 1999; Niu, 2004). A more extensive compilation of data for abyssal, ophiolitic and orogenic peridotites can be found in Bodinier and Godard, 2005. The thermodynamic model used in this study predicts that, for a fertile peridotite, plagioclase coexists with spinel on the solidus starting at 12 kb (Fig. 1). Experimental data on (batch) melting of natural

peridotite have shown that plagioclase is the first phase to dissolve in the melt (e.g. Jaques and Green, 1980; Takahashi, 1986; Falloon et al., 1988). In the numerical model the bulk composition of the system changes continuously during the melting process, therefore, there are certain compositional differences between the ridge model and the experimental studies. Nevertheless, the chemical variations observed in Fig. 6 are consistent with the experimental evidences in particular regarding the residual solid assemblage which does not includes plagioclase but only spinel and/or garnet. Along the ridge axis, spinel begins to be stable at 24 kbar (depth = 76 km, MgO = 39 wt%) in coexistence with garnet within a small pressure interval (<0.1 kbar). At lower pressure, spinel remains in the residual solid assemblage up to 6 kbar (depth = 18 km, MgO = 48.7 wt%). The composition of the residual solid obtained from this study is the result of a polybaric dynamic melting process (Klein and Langmuir, 1987). It essentially reflects the coexistence with a melt whose composition and amount at one particular depth is the result of the local melt production and the input of magma from deeper regions. The process shares some similarities with the incremental melting model that has been used to model the trace elements evolution observed in diopsides from abyssal peridotites (Johnson et al., 1990). The compositional trend at various depths from the numerical model appears to be in good agreement with the observed chemical data of major elements in abyssal peridotite. The model shows that the lowest MgO value consistent with abyssal peridotite corresponds to a depth of about 23 km. Nevertheless, few issues remain unclear. For instance, how the residual peridotite, shown in Fig. 6, is brought from at least 23 km depth up to the near surface area without being modified by any further process, and why the composition of abyssal peridotites collected on the bottom of the ocean floor is not restricted to a small range, but instead spans a large chemical interval, i.e. MgO range variation is 39–47 wt%. Moreover, from the modeling perspective, further melting is observed along the ridge axis from 23 km up to 10 km depth. Therefore, either some other mechanism must be in place and the dynamic melting does not entirely describe the chemical evolution of the oceanic mantle, or the geodynamic model should be revised in order to incorporate the idea that negligible amount of melt is produced above 23 km depth, and melt transport is so rapid that the geochemistry of the residual solid is not significantly modified. Partial crystallization between 3 and 6 kb has been invoked in previous studies to explain the composition of basalts in slow spreading mid-ocean ridge (Michael and Cornell, 1998), high-Al N-MORB (Eason and Sinton, 2006) and melting beneath oceanic transform faults (Gregg et al., 2009). It was also suggested that crystallization may begin at even greater pressure (11 kb) beneath oceanic ridges (Herzberg, 2004). A recent study using thermodynamic calculations concluded that reactive crystallization may be more appropriate to explain MORB signature (Collier and Kelemen, 2010). Alternatively, assimilation-fractional crystallization has been also used to correlate petrological evidences found in cumulates gabbros and migrating melt (Lissemberg and Dick, 2008). In general, crystallization at low pressure has been associated to the evolution of MORB melts (Langmuir et al., 1992), but the connection with the geochemical and geodynamical evolution of the sublithospheric mantle is still not quite established. Furthermore it remains to be clarified whether partial crystallization of oceanic melt (e.g. Collier and Kelemen, 2010), and the alternative hypothesis, suggesting reactive interaction of oceanic melt with dunites (e.g. Lambart et al., 2009) can coexist in a physical model of midocean ridges, or the two processes are mutually exclusive. 3.6. Refertilization of the residual mantle at low depth The simplest scenario compatible with the dynamic of a spreading ridge observed in this study requires that the depleted mantle

63

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

(a)

(b)

0.40

3.5

0.35 3.0

0.25 0.20 0.15

2.0 1.5 1.0

0.05

0.5 40.0

42.0 44.0 MgO (wt%)

46.0

0.0 38.0

48.0

40.0

42.0 44.0 MgO (wt%)

46.0

48.0

40.0

42.0 44.0 MgO (wt%)

46.0

48.0

(d)

10.0

4.0

9.5 9.0

Al2O3 (wt%)

FeO (wt%)

2.5

0.10

0.00 38.0

(c)

CaO (wt%)

Na2O (wt%)

0.30

8.5 8.0

3.0

2.0

7.5 1.0

7.0 38.0

(e)

40.0

42.0 44.0 MgO (wt%)

46.0

48.0

38.0

47.0

SiO2 (wt%)

46.0 45.0

z = 20 km

44.0 43.0 42.0 41.0 38.0

40.0

42.0 44.0 MgO (wt%)

46.0

48.0

Fig. 6. Variation of major elements composition in the residual solid as a function of MgO content (wt%) along the ridge axis after 4.45 Ma. Small square symbols connect the chemical composition of the residual solid at different depths computed from the numerical model. If not mentioned otherwise, field data refers to a spinel-bearing assemblage. Circles show data for abyssal peridotite corrected for alteration and averaged for each location from the Pacific and Indian ridges (Niu, 2004); triangles are abyssal peridotite compositions from the southwestern Indian and American–Antartic ridges (Dick and Fisher, 1984; Dick, 1989; Johnson and Dick, 1992; Snow and Dick, 1995); inverted triangles are samples from the mid-Atlantic ridge (Michael and Bonatti, 1985; Seyler and Bonatti, 1997); square symbols are peridotite samples with variable amount of plagioclase from the mid-Atlantic ridge (Seyler and Bonatti, 1997); diamond symbols are reconstructed compositions with 2–3 phases analyzed (Baker and Beckett, 1999). Black lines connecting empty symbols refer to the composition of the solid mantle after refertilization (see Section 3.6 for details).

is refertilized (Elthon, 1992; Seyler and Bonatti, 1997; Niu, 2004; Le Roux et al., 2007). Essentially the residual solid re-equilibrate, to some extent, with intermediate melt formed at some depth during mantle upwelling. As mentioned earlier, the model removes melt quite efficiently near the spreading center. However, it is unlikely that it is able to capture the detailed geochemical variations created by the solidification process because upwelling of melt at shallow depths is quite rapid, hence, local chemical equilibrium is most likely maintained only on a scale smaller than the numerical grid size (see discussion in Sections 4.1 and Fig. 8). In addition, the model does not permit multiple solids to coexist in the same location, therefore, partial solidification creates a solid with a composition that is a weighted average between the existing solid

and the newly crystallized assemblage. For these reasons, to understand any process that falls outside the limit of local chemical equilibrium on a scale smaller than the space grid size, we must introduce an additional computational procedure that makes use and elaborates the information provided by the dynamic model. To test the idea that residual peridotite is partially refertilized by certain amount of melt, two similar approaches have been taken. The two procedures aim to explore the possibility that either small pockets of melt remain trapped at the grain scale with the residual solid and follow the dynamics of the solid matrix, or that a small fraction of melt percolates upwards through portions of the mantle after melting of the peridotite is terminated. This small fraction of melt may have not necessarily been homogenized with

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

3.7. Computed seismic shear velocity map

the main body of melt moving through hypothetical dunite channels. A previous study on major and trace element compositions in pyroxenes and plagioclases proposed a similar idea in which melts with different compositions may coexist beneath mid-ocean ridges (Coogan et al., 2000). In the first approach, the solid and the melt found at 20 km depth from the dynamic model are combined together in proportion 99:1. The composition of the solid and melt are the following: initial solid (wt%) Na2O = 0.001, CaO = 0.66, MgO = 46.65, FeO = 8.52, Al2O3 = 0.82, SiO2 = 43.34, melt (wt%) Na2O = 1.45, CaO = 13.97, MgO = 5.67, FeO = 7.11, Al2O3 = 18.4, SiO2 = 53.4. The new solid is then recombined with the same melt but now in proportion 97:3. This operation is repeated few times to give at least an approximated indication that the process may occur over a certain spatial and temporal interval. In Fig. 6, the composition of the solid at each stage is illustrated by a black line connecting the circle symbols. In the second case, the same initial solid and melt compositions from the previous case have been used, but now, the assemblage is combined in proportion 90:10. A thermodynamic computation is performed at 6 kbar and at the particular temperature that leaves only 7% of melt in the equilibrium assemblage (final solid:melt proportion is 93:7). The new solid is then recombined with the same melt and the procedure is repeated three more times. The evolution of this solid is shown in Fig. 6 by the square symbols and the black line. The two models for a refertilization are clearly a simplification of a necessarily more complex process. The solid:melt proportion and the procedure used in these two cases are somehow arbitrary. Nevertheless, they appear to reproduce reasonably well the observed compositional range in abyssal peridotites which could be explained by a certain degree of refertilization. In summary, major elements in abyssal peridotites are consistent with a dynamic melting model that is followed by a process of partial refertilization. The overall scenario contains elements of previous models, recrystallization and refertilization of the depleted mantle (Niu et al., 1997; Niu and O’Hara, 2008), polybaric melting (Asimow, 1999) which takes place for the most part at pressure greater than 8–10 kb (Presnall et al., 2002). Above 20 km depth melt moves quite rapidly, therefore, reactive kinetics and solid–solid and solid–melt re-equilibration, may play a significant role in modifying the geochemical signature. The extent of such modifications as well as the proportion of the solid and liquid in equilibrium have been investigated experimentally (Wagner and Grove, 1998; Morgan and Liang, 2005; Van Den Bleeken et al., 2010). It is expected that a more detailed understanding of the interconnection between the mechanism of melt transport and chemical interaction between different solids and melt could be provided by combining the experimental data with new models that account for the evolution of multiple lithologies with different porosity (Liang and Parmentier, 2010).

180

160

140

120

100

km 80

0

60

20

4. Discussion The aim of the model presented in this study is to create a unified picture of the petrological and geodynamical evolution of

20

40

60

80

km 100

120

140

160

180 0

3.0 3.6

4.2

3.2

10

3.75 3.8

20 30

km

4 cm/a

0

3.8

4.0

10

4 cm/a

40

Fig. 7 reports the shear waves seismic velocity map obtained from the numerical model shown in Fig. 2b. Seismic velocities are computed following a method similar to the one outlined by Connolly and Kerrick (2002) with no correction for anelasticity and anharmonicity. Shear modulus of the mineral phases relevant for the model is obtained from various experimental sources (Bass, 1995; Liebermann, 2000; Murakami and Yoshioka, 2001). Density and adiabatic bulk modulus of the system at each location are retrieved from the chemical equilibrium computation which provides the minerals and melt proportion in the assemblage. The seismic model shown on the left side of Fig. 7, incorporates the effect of pressure, temperature, mineralogy and composition, but the contribution of melt is not included. In the regions of the mantle where melt crystallized underplate, for instance at the bottom of the lithosphere, the thermodynamic computation provides a mineralogical assemblage for a solid that is the weighted average between the solidified melt and a residual peridotite. Therefore, as the plate grows and moves away from the spreading center, the seismic map in these areas does not take into account the presence of different lithologies, like underplate gabbro and solidified basalt in the crustal layer. The effect induced by the presence of melt has been incorporated in the right plot of Fig. 7. Experimental results indicate that the decrease of the velocity per percent of melt can vary between 1.0% and 1.8% for P waves, and between 2.3% and 3.3% for S waves, depending on the geometry and interconnectivity of the melt pockets (Faul et al., 1994). The detailed melt distribution on a grain scale in mid-ocean ridges is poorly understood. Therefore, to roughly estimate the effect of melt, we have considered a fixed 3% reduction per percent of melt of the shear waves velocity. Several interesting features are shown by the seismic map which incorporates the effect of melt. The incipient separation of the melting area at approximately 120 km from the ridge axis is marked by an increase of the S waves velocity. The magma waves, shown in Figs. 3 and 4, around the ridge axis at 10–30 km depth, are quite distinguishable in the computed seismic map (Vs difference 2.5%). Attenuation of the shear waves velocity observed at the base of the lithosphere at distance less than 40 km from the ridge axis is related to the accumulation of melt in this section of the mantle. The presence of a small area with low velocity on the ridge axis at 4–6 km depth is consistent with the tomographic image of a low velocity (‘‘bull’s eye’’) zone beneath the center of the ridge at the transition between the lower crust and upper mantle (Dunn et al., 2005).

20 4.0

4.27

40

30

4.1 4.2

50

50

60

4.3

60

70

70 4.35

80 90

40

km

64

Vs w/o melt effect (km/s)

4.3

Vs with melt effect (km/s)

80 90

Fig. 7. Computed seismic image for shear waves beneath a mid-ocean ridge after 10.8 Ma. The effect of melt is included in the half-space model on the right side.

65

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

4.1. Extent of local chemical equilibrium One requirement of the model presented in this study is that local chemical equilibrium should be maintained at the time and space scale defined by the numerical grid size. When this is not the case,

for instance when melt moves very rapidly, then, kinetics or smaller equilibrium distance should be considered instead in order to model chemical variations during melting. A test on the local equilibrium approximation has been performed using the experimental results on basalt-lherzolite reactive dissolution (Morgan and Liang, 2005). The extent of the equilibration is expressed in the experimental study by a relation between the distance over which reactive dissolution is effective (equilibration distance) and the square root of the time. This relation has been applied on a melting column along the vertical axis of the ridge and at 20 km off the axis. Knowing the upwelling velocity of melt and the grid spacing (0.5–1 km) from the numerical simulation of this study, it is possible to retrieve the time required for melt to move through a grid point. Given the residence time, the equilibration distance can be inferred from the experimental data (AB Dun + Harz, Fig. 3 in Morgan and Liang, 2005). Vertical velocity and equilibration distance are shown in Fig. 8. The equilibration distance is described by discontinuous lines because the grid size changes from 0.5 km to 1 km at 50 km depth. Larger grid size means longer residence time, hence, chemical equilibration becomes effective over a greater distance. The penetration distance of the reactive equilibration between melt and solid ranges between 1 and 2 m over a significant portion of the melting column above 50 km depth, and is 3 m below 50 km depth on the ridge axis and 20 km off-axis. This implies that, if melt transport is controlled by porous flow, local equilibrium is easily achieved. Local equilibrium may also be a reasonable assumption in case of channeling flow, if the spacing between the channels is on the order of few meters. The model results suggest that local chemical equilibrium is effectively achieved though most of the melting column, however this may not be entirely the case in the upper 10–20 km in close proximity of the ridge axis, where upwelling velocity of melt increases quite rapidly. The grid size determines the length over which chemical equilibrium is achieved within the same domain. One may suggest to use coarser space grid so that chemical equilibrium can be extended over greater distance. However certain problems would arise, for instance, as it is well known, the numerical error increases with

0.5

equilibration distance (m) 1.5 2 2.5 3

1

10

3.5

4

velocity

20 30 depth (km)

mantle melting, and to show the interconnection between the physical and chemical processes beneath a mid-ocean ridge. The objective has been carried on by using a coupled thermodynamic and geodynamic approach. With this type of dynamic modeling it is possible to simulate, for instance, a polybaric dynamic melting process and to assess the lateral variation of the bulk composition and how the composition of the residual solid changes with depth. Magma waves are formed In the upper portion of the mantle melting region near the ridge axis. These magma pulses appear to be only a physical phenomenon with no influence on the crust production, mainly because a magma pool, few kilometers below the surface, buffers the melt coming from deeper regions. In addition, it seems that in the mantle region where these magma oscillations are observed, the geochemistry and petrology of the residual solid and melt are not significantly affected. A similar conclusion was drawn in an earlier study showing that solitary waves generated by an heterogeneous source do not carry a specific geochemical signature of the source from which they originated (Richter and Daly, 1989). In our numerical results there are no evidences of spontaneous localization of melt (porous channels) which previous studies described as induced by chemically reactive flow or artificially induced melt perturbations (Aharonov et al., 1995; Spiegelman et al., 2001; Spiegelman and Kelemen, 2003; Hewitt and Fowler, 2009; Hewitt, 2010). It has been proposed that reactive infiltration instabilities create porous channels with spacing on the order of 1–200 m (Spiegelman et al., 2001). The formation of porous channels in a numerical model of a mid-ocean ridge would help to explain the presence of dunite formations that have been found in opholitic complexes (Braun and Keleman, 2002). As shown in a previous study (Hewitt, 2010), the local equilibrium approximation does not prevent the development of porous channels since the dissolution rate during melt upwelling can be induced either by disequilibrium reactive kinetics or equilibrium melting processes. It is possible that channelized flow would appear with our model by adopting a smaller grid spacing. However, a preliminary numerical test using a numerical grid size of 250 m and a 100  100 km model box does not reveal new evidences. Another possibility is that, assuming a much greater permeability, the mantle may develop laterally distinctive geochemical regions. Melt will converge in specific areas where the melting rate would be higher (compositional effect), creating in this way a feedback process. Alternatively, high porosity of localized bands along shear deformation planes (Katz et al., 2006) may induce the chemical transformation required to form dunite channels. The effect of pyroxenite veins on melt productivity and thermal structure of the mantle has been investigated using a thermodynamic formulation (Morgan, 2001). Perhaps to explain the development of channels flow, the mantle should not be considered entirely formed by peridotite rocks. Coupling between melting of a chemical heterogeneous mantle and transport of melt through differently porous regions may enhance chemical variations to the extent that the mantle would naturally develop localized dunite formations. More work may also be needed to better understand whether the ideal mass transfer model used to study the reactive porous channels (i.e. feedback between melt flux and dissolution rate) (Aharonov et al., 1995; Spiegelman et al., 2001), which indeed appears to explain the trace element composition of melt in ridges (Spiegelman and Kelemen, 2003), realistically describes melt-solid reactions in the mantle beneath mid-ocean ridges.

40 50

distance

60

x=0 x = 20 km

70 80 90 -0.3

-0.2 vz (m/a)

-0.1

0

Fig. 8. Upwelling melt velocity along the ridge axis ðx ¼ 0Þ (solid line) and at 20 km from the ridge axis (dotted line) after 10.08 Ma (velocity negative upwards). Dashed and dash-dotted lines show the equilibration distance of a portion of the mantle that has reacted with melt during the time melt has moved through a distance equal to the space grid size. The dashed and dash-dotted lines is not continuous since above 50 km depth the grid size is 0.5 km, while below 50 km depth the grid size is 1 km.

66

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

coarse grid spacing. Another problem would be that certain information will be lost. Magma waves observed at depth less than 30 km form only when the grid size is less or equal to 0.5 km. Furthermore, certain petrological features are developed on a relatively small scale in comparison with the numerical model. Geochemical and field evidences in opholitic complexes show that at some depth in the oceanic mantle, melt reacts with dunite and harzburgite channels (Kelemen et al., 1997; Braun and Keleman, 2002). These channels, whose width varies between 3 mm and 100 m (Braun and Keleman, 2002), are not reproduced by the current numerical model. One possible reason is that heterogeneous structures smaller than the numerical grid size are only represented as an average geochemical information. The general conclusion is that, in order to provide a detailed description of the chemical processes under spreading ridges, the grid spacing should not be greater than the size of the observed petrological heterogeneities that are believed to affect the geochemical evolution of oceanic melt and the solid mantle. However, even when limited computational resources prevent the simulation of detailed petrological structures, it is still possible to extract certain geochemical results from the dynamic model and use additional petrological and thermodynamic calculations to interpret the field data. This approach was applied in this work to explain the compositional evolution of abyssal peridotites.

asthenospheric mantle unrelated to the ridge setting. The deeper portion of the asthenosphere that experienced various degrees of melting under the spreading ridge, shows a downflow pattern. This observation suggests that eventually a portion of depleted asthenospheric mantle will become part of the global mantle circulation, contributing in this way to the formation of chemical heterogeneities in the upper mantle. 8. The dynamic melt process that takes place over a continuous depth interval, determines a trend for major elements composition of the residual solid in the spinel stability field that is consistent with the range observed in abyssal peridotites. From a geodynamic point of view, however, it is difficult to assume that, during mantle upwelling, in particular at 20–30 km depth and above, only partial melting is responsible for the geochemical signature that characterizes abyssal peridotites. A picture more consistent with the spreading ridge model developed in this study requires that abyssal peridotites are the result of a dynamic melting process followed by refertilization at low depths (<20–30 km). 9. Computed shear waves map reveals several features of the melt distribution in the ridge area such as, width and depth of the melting region, melt accumulation at the lithosphere–asthenosphere boundary at some distance from the ridge axis, lateral convergence of melt at 4–6 km depth near the spreading center and magma waves in proximity of the ridge axis.

4.2. Summary The main results of this study can be summarized as follow: 1. The numerical simulation indicates a continuous melt region that reaches the maximum distance of 150 km from the spreading axis after approximately 4.5 Ma. The melting area extends over a steady distance of 100 km from the ridge axis after 10 Ma. 2. Melting of the upwelling mantle is observed up to 16 km depth under the ridge axis. This depth increases moving away from the ridge axis. In the shallower mantle, the amount of melt slightly decreases moving upward until melt reaches the bottom of the lithosphere or the area near surface where it rapidly freezes. 3. Melt focuses toward the surface within a region extending up to 40–50 km from the ridge axis following either the shortest path at the base of the lithosphere or directly from the asthenospheric mantle (distance <10 km from the spreading ridge). At distances greater than 50 km, melt moves towards the lithospheric–asthenospheric interface where it crystallizes into intrusive cumulates. 4. In the near surface region around the ridge axis, melt tends to form porous waves with periodic oscillations that are affected by the permeability of the mantle. Since compaction length decreases with depth, only models with sufficiently high spatial resolution can resolve detailed variations of the melt content in the lower part of the melting region (see Appendix A). 5. Crustal thickness show temporal fluctuations. Variations appear to be mainly controlled by the rapid discharge of the accumulated melt at the base of the cold lithosphere. The average crust thickness is 6.0 km, that is about 60% of the total solidified melt. 6. As the ridge evolves, a portion of the melting region separates from the main melting area around the ridge axis. The duration of the simulation does not allow us to determine the fate of the segregated portion of the melting area. It is expected that eventually it will slowly cool far from the ridge axis (outside the model box) and the residual melt will freeze in the asthenosphere. 7. Aside from the ridge axis, the depleted mantle, previously involved in the melting process, not necessarily becomes part of the cooling lithosphere that, instead, may grow incorporating

Acknowledgments This work is a partial contribution from a research developed over several years which would not have been possible without the financial support from the following grants: DFG Grant No. Schm872/13-1, NSF Grant No. DMR-0231291, NSF-OCE Grant No. 0452867. We are grateful to two anonymous reviewers and the editor for their critical comments. Y. Niu and A. Ghods provided helpful reviews that improved an earlier version of this manuscript. Thanks to J. Ganguly for the discussions during the development of the thermodynamic model and H. Schmeling for sharing

depth (km) vz=0.6 m/yr 0 10 20 30 40 50

Δ z=0.0625km

60 70 80 90

Δz=0.5km

100 0

0.01

0.02 0.03

melt fraction Fig. 9. Steady state melt distribution for a simple melt flow model assuming a fixed melt source at 97 km depth. Melt oscillations are well resolved using a grid size ðDzÞ of 0.0625 km, while no melt variations are detected when the grid size is set to 0.5 km. Initial melt content is 3%. The viscosity of the matrix is constant and equal to 1  1019 Pa s.

67

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

melt fraction 0

vz = 0.6 m/yr

0.05

melt fraction

0.1

0.05

melt fraction

0.15

0.05

10

10

10

10

20

20

20

20

30

30

30

30

a1 0 km

melt fraction

0 0.05 0.1 0.15

40 Δz = 0.5 km

b1 40 Δz = 0.25 km

40

50

50

Δz=0.125 km Δt = t+0.012 50 Δt(Gmin)=0.004

depth (km)

depth (km)

Δt = t+0.01 Δt(Gmin)=0.01

0.15

Δt = t+0.01 Δt(Gmin)=0.01

c1

40

d1 Δz=0.0625 km 50 Δt = t+0.01 Δt(Gmin)=0.003

depth (km) depth (km) equilib. dist.(m) equilib. dist.(m) equilib. dist.(m) equilib. dist.(m)

50 km

0.6 0.9 1.2 1.5

0.6 0.8

a2

1

1.2

0.4 0.6 0.8

b2

1

0.4

c2

0.6

0.8

1

d2

10

10

10

10

20

20

20

20

melt 100 km T =1470oC

eq. dis. 30

40

MgO

50

eq. dis. 30

30

40

40

MgO

50 40

45

50

MgO % (solid)

40

MgO

eq. dis.

50

45

50

MgO % (solid)

30

MgO

40

eq. dis.

50 40

45

50

MgO % (solid)

40

45

50

MgO % (solid)

Fig. 10. Left central panel. Schematic view of the pseudo 1-D model. The upper side is open to mass flow with a prescribed vertical velocity of 0.6 m/a. Bottom temperature is 1470 °C. All the other parameters and the thermodynamic model are the same of the 2-D ridge model. (a1)–(d1) Melt distribution using a vertical grid step size ðDzÞ of 0.5, 0.25, 0.125, 0.0625 km. Dt(Gmin) is the time interval separating two computations of the Gibbs energy minimization procedure. x-Scale is different for different panels. (a2)– (d2) MgO content in the residual solid and extent of chemical equilibration between melt and the solid matrix (see discussion in Section 4.1 related to Fig. 8 for a discussion on the definition of chemical equilibration). Scale for chemical equilibration is different for different panels.

his view on the relation between numerical resolution and formation of porous waves. Appendix A. Implications of the numerical resolution for the magma transport and petrological results Porous flow through a deformable matrix can generate variations in the magma content over depth (Spiegelman, 1993a,b). The nature of the magma waves is controlled by the compaction length which depends mainly on the bulk and shear viscosity of the matrix and the local porosity (Mckenzie, 1984). From the numerical standpoint, this melt behavior can be investigated only if the spatial grid size is to some extent less than the compaction length. Discretization of the 2-D ridge model is performed on a spatial grid that varies between 0.5 km near the ridge axis to 1 km away from the main melting region. Here two simple cases are used to analyzed the effect of a more fine mesh in particular on the development of magma waves and the implications on the assumption of local thermodynamic equilibrium as well as the composition of the solid mantle. In the following example it is assumed a 1-D melt column in which a melt source is located at approximately 97 km depth. Compaction length is 0.5 km when the melt supply is 3%. The setup describes a pure dynamic model in which an extraction vertical velocity of 0.6 m/a is imposed on the top side and the matrix viscosity is 1  1019 Pa s.

Fig. 9 shows that for the steady state solution of the problem melt oscillations are well developed with a grid step size of 0.0625 km, while no variations of the melt content with depth are observed using a 0.5 km grid size. A more complex scenario may arise when the melt content is allowed to vary according to the local melt production and temperature variations are taken into account. In this case the compaction length is not constant. Therefore, by increasing the spatial resolution in the vertical direction, melt instabilities may develop starting at greater depth where melt fraction and viscosity of the matrix decreases, hence compaction length is smaller. This problem has been addressed using a pseudo 1-D column in which the parameters and thermodynamic model are the same as for the 2-D model with one exception. Instead of imposing a spreading velocity (vx), a vertical velocity (vz) of 0.6 m/a is assumed on the surface. The velocity is slightly higher than the spreading velocity of 0.4 m/a considered in the 2-D model mainly to reduce the computational effort. The numerical grid consists of 10 grid points in the x direction and 200, 400, 800, 1600 points in the vertical direction which correspond to a spacing of 0.5, 0.25, 0.125, 0.0625 km, respectively. Compaction length varies from 100 m near the base of the melting column to 10,000 m in the upper portion of the model. Time step is set to insure numerical stability ðDt ¼ 50—500 aÞ. A schematic view of the pseudo 1-D model is shown in Fig. 10. The melt distribution versus depth computed with different grid spacing is illustrated in

68

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

Fig. 10a1–d1. The 1-D model in Fig. 10a with a grid size of 0.5 km is consistent with the 2-D ridge model in the main text. The main observation that can be made by comparing the results using a more refined mesh is that melt instabilities begin to develop at greater depth and progressively grow moving towards the surface. The lower panels in Fig. 10a2–d2 illustrate the variation of the MgO content in the residual solid with depth and the extent of equilibration computed following the procedure described in Section 4.1. It appears that the composition of the residual solid remains quite insensitive to the dynamic of the melt upwelling. The extent of chemical equilibration changes in accordance with the variations in the magma content mainly in response to the variations of the melt velocity. In summary, at conditions similar to the one applied to the ridge model, by increasing the grid resolution in a 1-D melting model, the accuracy of the numerical method cannot be easily assessed because melt compaction/decompaction develops inherently unstable melt oscillations deeper in the melt column. The influence of the numerical grid on the dynamics of magma does not seem to have been observed in previous studies on ridge melting (e.g. Ghods and Arkani-Hamed, 2000; Katz, 2008). The 2-D ridge model from this study develops magma oscillations that are similar to the one found by the higher resolution 1-D model, however the detailed magma transport at greater depth is not resolved accurately because the numerical discretization is not compatible with the compaction length at depths greater than 30 km. While this may have implication for understanding the physics of the magma dynamics from a theoretical point of view, practically the compaction/decompaction process does not affect the magma extraction on the surface which is mainly controlled by the lateral focusing of the melt under the lithosphere. The petrology of the residual solid is also unaffected. On the other hand, the extent of chemical equilibration may be influenced by the periodic variations of the upwelling velocity of melt. References Aharonov, E., Withehead, J.A., Kelemen, P.B., Spiegelman, M., 1995. Channeling instability of upwelling melt in the mantle. J. Geophys. Res. 100, 20433–20450. Asimow, P.D., 2002. Steady-state mantle–melt interactions in one dimension: II. Thermal interactions and irreversible terms. J. Petrol. 43, 1707–1724. Asimow, P.D., 1999. A model that reconciles major- and trace-element data from abyssal peridotites. Earth Planet. Sci. Lett. 169, 303–319. Asimow, P.D., 1995. The effect of pressure-induced solid–solid phase transitions on decompression melting of the mantle. Geochim. Cosmochim. Acta 59, 4489– 4506. Asimow, P.D., Langmuir, C.H., 2003. The importance of water to oceanic mantle melting regimes. Nature 421, 815–820. Asimow, P.D., Hirschmann, M.M., Stolper, E.M., 2001. Calculation of peridotite partial melting from thermodynamic models of minerals and melts, IV. Adiabatic decompression and the composition and mean properties of midocean ridge basalts. J. Petrol. 42, 963–998. Asimow, P.D., Stolper, E.M., 1999. Steady-state mantle–melt interactions in one dimension: I. Equilibrium transport and melt focusing. J. Petrol. 40, 475–494. Baba, K., Chave, A.D., Evans, R.L., Hirth, G., Mackie, R.L., 2006. Mantle dynamics beneath the East Pacific Rise at 17°S: insights from the mantle electromagnetic and tomography (MELT) experiment. J. Geophys. Res. 111. doi:10.1029/ 2004JB003598. Baker, M.B., Beckett, J.R., 1999. The origin of abyssal peridotites: a reinterpretation of constraints based on primary bulk compositions. Earth Planet. Sci. Lett. 171, 49–61. Baker, M.B., Hirschmann, M.M., Ghiorso, M.S., Stolper, E.M., 1995. Compositions of near-solidus peridotite melts from experiments and thermodynamic calculations. Nature 375, 308–311. Bass, J.D., 1995. Elasticity of minerals, glasses, and melts. In: Ahrens, T.J. (Ed.), Mineral physics and crystallography: a handbook of physical constants. American Geophysical Union, Washington DC, 1995, pp. 45–63. Berman, R.G., 1988. Internally-consistent thermodynamic data for minerals in the system Na2O–K2O–CaO–MgO–FeO–Fe2O3–Al2O3–SiO2–TiO2–H2O–CO2. J. Petrol. 29, 445–552. Bodinier, J.L., Godard, M., 2005. Orogenic ophiolitic and abyssal peridotites. In: Carlson, R.W. (Ed.), The Mantle and Core (Treatise on Geochemistry), vol. 2, Elsevier Ltd., United Kingdom, 2005, pp. 103–170. Bonatti, E., Ligi, M., Brunelli, D., Cipriani, A., Fabretti, P., Ferrante, V., Gasperini, L., Ottolini, L., 2003. Mantle thermal pulses below the mid-atlantic ridge and temporal variations in the formation of oceanic lithosphere. Nature 423, 499– 505.

Book, D.L., Boris, J.P., 1975. Flux-corrected transport II: generalizations of the method. J. Comput. Phys. 18, 248–283. Borghini, G., Fumagalli, P., Rampone, E., 2010. The stability of plagioclase in the upper mantle: subsolidus experiments on fertile and depleted lherzolite. J. Petrol. 51, 229–254. Boris, J.P., Book, D.L., 1973. Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works. J. Comput. Phys. 11, 38–69. Bown, J.W., White, R.S., 1994. Variation with spreading rate of oceanic crustal thickness and geochemistry. Earth Planet. Sci. Lett. 121, 435–449. Braun, M.G., Keleman, P.B., 2002. Dunite distribution in the Oman ophiolite: implications for melt flux through porous dunite conduits. Geochem. Geophys. Geosyst. 3, doi: 2001GC000289. Braun, M.G., Hirt, G., Parmentier, E.M., 2000. The effects of deep damp melting on mantle flow and melt generation beneath mid-ocean ridges. Earth Planet. Sci. Lett. 176, 339–356. Canales, J.P., Collins, J.A., Escartìn, J., Detrick, R.S., 2000. Seismic structure across the rift valley of the mid-atlantic ridge at 23°200 N (MARK Area): implications for crustal accretion processes at slow spreading ridges. J. Geophys. Res. 105, 28411–28425. Canales, J.P., Detrick, R.S., Bazin, S., Harding, A.J., Orcutt, J.A., 1998. Off-axis crustal thickness across and along the east pacific rise within the MELT area. Science 280, 1218–1221. Cannat, M., 1996. How thick is the magmatic crust at slow spreading ocean ridges? J. Geophys. Res. 101, 2847–2857. Cipriani, A., Bonatti, E., Brunelli, D., Ligi, M., 2009. 26 Million years of mantle upwelling below a segment of the mid atlantic ridge: the vema lithospheric section revised. Earth Planet. Sci. Lett. 285, 87–95. Cochran, J.R., 1986. Variations in subsidence rates along intermediate and fast spreading mid-ocean ridges. Geophys. J. Int. 87, 421–454. Collier, M.L., Kelemen, P.B., 2010. The case for reactive crystallization at mid-ocean ridges. J. Petrol. 9, 1913–1940. Connolly, J.A.D., 2005. Computation of phase equilibria by linear programming: a tool for geodynamic modeling and its application to subduction zone decarbonation. Earth Planet. Sci. Lett. 236, 524–541. Connolly, J.A.D., Schmidt, M.W., Solferino, G., Bogdassarov, N., 2009. Permeability of asthenospheric mantle and melt extraction rates at mid-ocean ridges. Nature 462, 209–212. Connolly, J.A.D., Kerrick, D.M., 2002. Metamorphic controls on seismic velocity of subducted oceanic crust at 100–250 km depth. Earth Planet. Sci. Lett. 204, 61– 74. Connolly, J.A.D., Kerrick, D.M., 1987. An algorithm and computer program for calculating composition phase diagrams. Calphad 11, 1–55. Coogan, L.A., Kempton, P.D., Saunders, A.D., Norry, M.J., 2000. Melt aggregation within the crust beneath the mid-atlantic ridge: evidence from plagioclase and clinopyroxene major and trace element compositions. Earth Planet. Sci. Lett. 176, 245–257. Cordery, M.J., Morgan, J.P., 1993. Convection and melting at mid-ocean ridges. J. Geophys. Res. 98, 19477–19503. Deb, K., 2001. Multi-Objective Optimization Using Evolutionary Algorithms. John Wiley & Sons Ltd., England, first ed., 497 pp. Dick, H.J.B., 1989. Abyssal peridotites, very slow spreading ridges and ocean ridge magmatism. In: Saunders, A.D., Norry, M.J. (Eds.), Magmatism in the Ocean Basins, 1989, pp. 71–105 (Geological Society, London, Special Publications 42). Dick, H.J.B., Tivey, M.A., Tucholke, E.B., 2008. Plutonic foundation of a slowspreading ridge segment: oceanic core complex at Kane Megamullion, 23°300 N, 45°200 W. Geochem. Geophys. Geosyst. 9. doi:10.1029/2007GC001645. Dick, H.J.B., Fisher, R.L., 1984. Mineralogic studies of the residues of mantle melting: abyssal and alpine-type peridotites. In: Kornprobst, J. (Ed.), Kimberlites II: The Mantle and Crust–Mantle Relationships. Elsevier, New York, 1984, pp. 295–308. Dunn, R.A., Lekic´, V., Detrick, R.S., Toomey, D.R., 2005. Three-dimensional seismic structure of the Mid-Atlantic Ridge (35°N): evidence for focused melt supply and lower crustal dike injection. J. Geophys. Res. 110. doi:10.1029/ 2004JB003473. Dunn, R.A., Forsyth, D.W., 2003. Imaging the transition between the region of mantle melt generation and the crustal magma chamber beneath the southern East Pacific Rise with short-period Love waves. J. Geophys. Res. 108. doi:10.1029/2002JB002217. Dunn, R.A., Toomey, D.R., Solomon, S.C., 2000. Three-dimensional seismic structure and physical properties of the crust and shallow mantle beneath the East Pacific Rise at 9°300 N. J. Geophys. Res. 105, 23537–23555. Eason, D., Sinton, J., 2006. Origin of high-Al N-MORB by fractional crystallization in the upper mantle beneath the Galápagos Spreading Center. Earth Planet. Sci. Lett. 252, 423–436. Eberle, M.A., Forsyth, D.W., 1998b. An alternative, dynamic model of the axial topographic high at fast spreading ridges. J. Geophys. Res. 103, 12309–12320. Eberle, M.A., Forsyth, D.W., Parmentier, E.M., 1998a. Constraints on a buoyant model for the formation of the axial topographic high on the East Pacific Rise. J. Geophys. Res. 103, 12291–12307. Elthon, D., 1992. Chemical trends in abyssal peridotites: refertilization of depleted suboceanic mantle. J. Geophys. Res. 97, 9015–9025. Eriksson, E., 1971. Thermodynamic studies of high temperature equilibria, III. SOLGAS a computer program for calculating the composition and heat condition of an equilibrium mixture. Acta Chem. Scand. 25, 2651–2658. Eriksson, E., 1975. Thermodynamic studies of high temperature equilibria, XII. SOLGASMIX a computer program for calculation of equilibrium compositions in multiphase systems. Chemica Scripta 8, 100–103.

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70 Eriksson, E., Rosén, E., 1973. Thermodynamic studies of high temperature equilibria, VIII. General equations for the calculation of equilibria in multiphase systems. Chemica Scripta 4, 193–194. Evans, R.L., Hirth, G., Baba, K., Forsyth, D., Chave, A., Mackie, R., 2005. Geophysical evidence from the MELT area for compositional controls on oceanic plates. Nature 437. doi:10.1038/nature04014. Evans, R.L., Tarits, P., Chave, A.D., White, A., Heinson, G., Filloux, J.H., Toh, H., Seama, N., Utada, H., Booker, J.R., Unsworth, M.J., 1999. Asymmetric electrical structure in the mantle beneath the East Pacific Rise at 17°S. Science 286, 752–756. Falloon, T.J., Green, D.H., Hatton, C.J., Harris, K.L., 1988. Anhydrous partial melting of a fertile and depleted peridotite from 2 to 30 kb and application to basalt petrogenesis. J. Petrol. 29, 1257–1282. Faul, U.H., Toomey, D.R., Waff, H.S., 1994. Intergranular basaltic melt is distributed in thin, elongated inclusions. Geophys. Res. Lett. 21, 29–32. Forsyth, D.W., Webb, S.C., Dorman, L.M., Shen, Y., 1998. Phase velocities of Rayleigh waves in the MELT experiment on the East Pacific Rise. Science 280, 1235–1238. Fullea, J., Afonso, J.C., Connolly, J.A.D., Fernàndez, M., Garcia-Castellanos, D., Zeyen, H., 2009. LitMod3D: An interactive 3-D software to model the thermal, compositional, density, seismological, and rheological structure of the lithosphere and sublithospheric upper mantle. Geochem. Geophys. Geosyst. 10, doi: 10.2009GC002391. Ganguly, J., Saxena, S.K., 1987. Mixtures and Mineral Reactions. Springer-Verlag, Berlin, New York, 291 pp. Ghiorso, M.S., Hirschmann, M.M., Reiners, P.W., Kress III, V.C., 2002. The pMELTS: A revision of MELTS for improved calculation of phase relations and major element partitioning related to partial melting of the mantle to 3 GPa. Geochem. Geophys. Geosyst. 3. doi:10.1029/2001GC000217. Ghiorso, M.S., Sack, R.O., 1995. Chemical mass transfer in magmatic processes IV. A revised and internally consistent thermodynamic model for the interpolation and extrapolation of liquid–solid equilibria in magmatic systems at elevated temperatures and pressures. Contrib. Mineral. Petrol. 119, 197–212. Ghods, A., Arkani-Hamed, J., 2000. Melt migration beneath mid-ocean ridges. Geophys. J. Int. 140, 687–697. Girardeau, J., Monnier, C., Mée, Le, Quatrevaux, F., 2002. The Wuqbah peridotite, central Oman ophiolite: petrological characteristics of the mantle in a fossil overlapping ridge setting. Mar. Geophys. Res. 23, 43–56. Godard, M., Jousselin, D., Bodinier, J.L., 2000. Relationships between geochemistry and structure beneath a paleo-spreading centre: a study of the mantle section in the Oman ophiolite. Earth Planet. Sci. Lett. 180, 133–148. Gregg, P.M., Behn, M.D., Lin, J., Grove, T.L., 2009. Melt generation, crystallization, and extraction beneath segmented oceanic transform faults. J. Geophys. Res. 114. doi:10.1029/2008JB006100. Grove, T.L., Kinzler, R.J., Bryan, W.B., 1992. Fractionation of mid-ocean ridge basalt (MORB). Mantle flow and melt generation at mid-ocean ridges. In: Morgan, P.J., Blackman, D.K., Sinton, J.M. (Eds.), Mantle Flow and Melt Generation at MidOcean Ridges. Geophysical Monograph, American Geophysical Union 71, pp. 281–310. Gudfinnsson, G.H., Presnall, D.C., 2000. Melting behavior of model lherzolite in the system CaO–MgO–Al2O3–SiO2–FeO at 0.7–2.8 GPa. J. Petrol. 41, 1241–1269. Hall, C.E., Parmentier, E.M., 2000. Spontaneous melt localization in a deforming solid with viscosity variations due to water weakening. Geophys. Res. Lett. 27, 9–12. Heinson, G., Constable, S., White, A., 2000. Episodic melt transport at mid-ocean ridges inferred from magnetotelluric sounding. Geophys. Res. Lett. 27, 2317– 2320. Hernlund, J.W.C., Stevenson, D.J., Tackley, P.J., 2008b. Buoyant melting instabilities beneath extending lithosphere: 2. Linear analysis. J. Geophys. Res. 113. doi:10.1029/2006JB004863. Hernlund, J.W.C., Tackley, P.J., Stevenson, D.J., 2008a. Buoyant melting instabilities beneath extending lithosphere: 1. Numerical models. J. Geophys. Res. 113. doi:10.1029/2006JB004862. Herzberg, C., 2004. Partial crystallization of mid-ocean ridge basalts in the crust and mantle. J. Petrol. 45, 2389–2405. Hewitt, I.J., 2010. Modelling melting rates in upwelling mantle. Earth Planet. Sci. Lett. 300, 264–274. Hewitt, I.J., Fowler, A.C., 2009. Melt channelization in ascending mantle. J. Geophys. Res. 114. doi:10.1029/2008JB006185. Hewitt, I.J., Fowler, A.C., 2008. Partial melting in an upwelling mantle column. Proc. R. Soc. A, doi:10.1098/rspa.2008.0045. Hirose, K., Kushiro, I., 1993. Partial melting of dry peridotites at high pressures: determination of compositions of melts segregated from peridotite using aggregates of diamonds. Earth Planet. Sci. Lett. 114, 477–489. Hirschmann, M.M., Ghiorso, M.S., Stolper, E.M., 1999a. Calculation of peridotite partial melting from thermodynamic models of minerals and melts. II. Isobaric variations in melts near the solidus and owing to variable source composition. J. Petrol. 40, 297–313. Hirschmann, M.M., Asimow, P.D., Ghiorso, M.S., Stolper, E.M., 1999b. Calculation of peridotite partial melting from thermodynamic models of minerals and melts. III. Controls on isobaric melt production and the effect of water on melt production. J. Petrol. 40, 831–851. Hirschmann, M.M., Ghiorso, M.S., Wasylenki, L.E., Asimow, P.D., Stolper, E.M., 1998. Calculation of peridotite partial melting from thermodynamic models of minerals and melts. I. Review of methods and comparison with experiments. J. Petrol. 39, 1091–1115. Hirth, G., Kohlstedt, D.L., 1995. Experimental constraints on the dynamics of the partially molten upper mantle: deformation in the diffusion creep regime. J. Geophys. Res. 100, 1981–2001.

69

Hirth, G., Kohlstedt, D.L., 2003. Rheology of the upper mantle and the mantle wedge: a view from the experimentalists. In Eiler, J. (Ed.), Inside the Subduction Factory, Geophysical Monograph, vol. 138, American Geophysical Union, Washington, D.C., pp. 83–105. Iwamori, H., 1993. Dynamic disequilibrium melting model with porous flow and diffusion- controlled chemical equilibration. Earth Planet. Sci. Lett. 114, 301– 313. Jaques, A.L., Green, D.H., 1980. Anhydrous melting of peridotite at 0–15 Kb pressure and the genesis of tholeiitic basalts. Contrib. Mineral. Petrol. 73, 287–310. Johnson, K.T.M., Dick, H.J.B., 1992. Open system melting and temporal and spatial variation of peridotite and basalt at the Atlantis II fracture zone. J. Geophys. Res. 97, 9219–9241. Johnson, K.T.M., Dick, H.J.B., Shimizu, N., 1990. Melting in the oceanic upper mantle: an ion microprobe study of diopsides in abyssal peridotite. J. Geophys. Res. 95, 2661–2678. Jousselin, D., Nicolas, A., Boudier, F., 1998. Detailed mapping of a mantle diapir below a paleo-spreading center in the Oman ophiolite. J. Geophys. Res. 103, 18153–18170. Katz, R.F., Spiegelman, M., Holtzman, B., 2006. The dynamics of melt and shear localization in partially molten aggregates. Nature 442, 676–679. Katz, R.F., 2008. Magma dynamics with the enthalpy method: benchmark solutions and magmatic focusing at mid-ocean ridges. J. Petrol. 49, 2099–2121. Kelemen, P.B., 1990. Reaction between ultramafic rock and fractionating basaltic magma I. Phase relations, the origin of calc-alkaline magma series, and the formation of discordant dunite. J. Petrol. 31, 51–98. Kelemen, P.B., Joyce, D.B., Webster, J.D., Holloway, J.R., 1990. Reaction between ultramafic rock and fractionating basaltic magma II. Experimental investigation of reaction between olivine tholeiite and harzburgite at 1150–1050 °C and 5 kb. J. Petrol. 31, 99–134. Kelemen, P.B., Dick, H.J.B., Quick, J.E., 1992. Formation of harzburgite by pervasive melt/rock reaction in the upper mantle. Nature 358, 635–641. Kelemen, P.B., Hirth, G., Shimizu, N., Spiegelman, M., Dick, H.J.B., 1997. A review of melt migration processes in the adiabatically upwelling mantle beneath oceanic spreading ridges. Phil. Trans. R. Soc. Lond. A 355, 283–318. Kinzler, R.J., Grove, T.L., 1993. Corrections and further discussion of the primary magmas of mid-ocean ridge basalts, 1 and 2. J. Geophys. Res. 98, 22339–22347. Kinzler, R.J., Grove, T.L., 1992b. Primary magmas of mid-ocean ridge basalts 2: applications. J. Geophys. Res. 97, 6907–6926. Kinzler, R.J., Grove, T.L., 1992a. Primary magmas of mid-ocean ridge basalts 1: experiments and methods. J. Geophys. Res. 97, 6885–6906. Klein, E.M., Langmuir, C.H., 1987. Global correlations of ocean ridge basalt chemistry with axial depth and crustal thickness. J. Geophys. Res. 92, 8089– 8115. Kohlstedt, D.L., Bai, Q., Wang, Z.C., Mei, S., 2000. Rheology of partially molten rocks. In: Bagdassarov, N., Laporte, D., Thompson, A.B. (Eds.), Physics and Chemistry of Partially Molten Rocks. Kluwer Academic Publ., Dordrecht, pp. 3–28. Lambart, S., Laporte, D., Schiano, P., 2009. An experimental study of focused magma transport and basalt-peridotite interactions beneath mid-ocean ridges: implications for the generation of primitive MORB compositions. Contrib. Mineral. Petrol. 157, 429–451. Langmuir, C.H., Klein, E.M., Plank, T., 1992. Petrological systematics of mid-ocean ridge basalts: constraints on melt generation beneath ocean ridges. In: Morgan, J.P., Blackman, D.K., Sinton, J.M. (Eds.), Mantle Flow and Melt Generation at MidOcean Ridges, Geophysical Monograph 71, American Geophysical Union, Washington DC, 2000, pp. 183–280. Le Mée, L., Girardeau, J., Monnier, C., 2004. Mantle segmentation along the Oman ophiolite fossil mid-ocean ridge. Nature 432, 167–172. Le Roux, V., Bodinier, J.L., Tommasi, A., Alard, O., Dautria, J.M., Vauchez, A., Riches, A.J.V., 2007. The Lherz spinel lherzolite: refertilized rather than pristine mantle. Earth Planet. Sci. Lett. 259, 599–612. Liang, Y., Parmentier, E.M., 2010. A Two-Porosity Double Lithology Model for Partial Melting, Melt Transport and Melt-rock Reaction in the Mantle: Mass Conservation Equations and Trace Element Transport. J. Petrol. 51, 125–152. Liebermann, R.C., 2000. Elasticity of mantle minerals: Experimental studies. In Karato, S., Stixrude, L., Liebermann, R., Masters, G., Forte, A. (Eds.), Earth’s deep interior: mineral physics and tomography from the atomic to the global scale., Geophysical Monograph 117, American Geophysical Union, Washington DC, 2000, pp. 181-199. Lissemberg, C.J., Dick, J.B., 2008. Melt-rock reaction in the lower oceanic crust and its implications for the genesis of mid-ocean ridge basalt. Earth Planet. Sci. Lett. 271, 311–325. Lizarralde, D., Gaherty, J.B., Collins, J.A., Hirth, G., Kim, S.D., 2004. Spreading-rate dependence of melt extraction at mid-ocean ridges from mantle seismic refraction data. Nature 432, 744–747. Longhi, J., 2002. Some phase equilibrium systematics of lherzolite melting: I. Geochem. Geophys. Geosyst. 3. doi:10.1029/2001GC000204. Lundstrom, C.C., Gill, J., Williams, Q., 2000. A geochemically consistent hypothesis for MORB generation. Chem. Geol. 162, 105–126. Maia, M., Gente, P., 1998. Three-dimensional gravity and bathymetry analysis of the Mid-Atlantic Ridge between 20°N and 24°N: flow geometry and temporal evolution of the segmentation. J. Geophys. Res. 103, 951–974. Maumus, J., Bagdassarov, N., Schmeling, H., 2005. Electrical conductivity and partial melting of mafic rocks under pressure. Geochim. Cosmochim. Acta 69, 4703– 4718. Mckenzie, D., 1984. The generation and compaction of partially molten rock. J. Petrol. 25, 713–765.

70

M. Tirone et al. / Physics of the Earth and Planetary Interiors 190–191 (2012) 51–70

McKenzie, D., Bickle, M.J., 1988. The volume and composition of melt generated by extension of the lithosphere. J. Petrol. 29, 625–679. Michael, P.J., Cornell, W.C., 1998. Influence of spreading rate and magma supply on crystallization and assimilation beneath mid-ocean ridges: evidence from chlorine and major element chemistry of mid-ocean ridge basalts. J. Geophys. Res. 103, 18325–18356. Michael, P.J., Bonatti, E., 1985. Peridotite composition from the North Atlantic: regional and tectonic variations and implications for partial melting. Earth Planet. Sci. Lett. 73, 91–104. Morency, C., Doin, M.P., Dumoulin, C., 2005. Three-dimensional numerical simulations of mantle flow beneath mid-ocean ridges. J. Geophys. Res. 110. doi:10.1029/2004JB003454. Morgan, J.P., 2001. Thermodynamics of pressure release melting of a veined plum pudding mantle. Geochem. Geophys. Geosyst. 2, 2000GC000049. Morgan, J.P., 1987. Melt migration beneath mid-ocean spreading centers. Geophys. Res. Lett. 14, 1238–1241. Morgan, J.P., Holtzman, B.K., 2005. Vug waves: a mechanism for coupled rock deformation and fluid migration. Geochem. Geophys. Geosyst. 6. doi:10.1029/ 2004GC000818. Morgan, Z., Liang, Y., 2005. An experimental study of the kinetics of lherzolite reactive dissolution with applications to melt channel formation. Contrib. Mineral. Petrol. 150, 369–385. Murakami, T., Yoshioka, S., 2001. The relationship between the physical properties of the assumed pyrolite composition and depth distributions of seismic velocities in the upper mantle. Phys. Earth Planet. Inter. 125, 1–17. Nedimovic´, M.R., Carbotte, S.M., Harding, J.A., Detrick, R.S., Canales, J.P., Diebold, J.B., Kent, G.M., Tischer, M., Babcock, J.M., 2005. Frozen magma lenses below the oceanic crust. Nature 436. doi:10.1038/nature03944. Niu, Y., 2004. Bulk-rock major and trace element compositions of abyssal peridotites: implications for mantle melting, melt extraction and post-melting processes beneath mid-ocean ridges. J. Petrol. 45, 2423–2458. Niu, Y., O’Hara, M.J., 2008. Global correlations of ocean ridge basalt chemistry with axial depth: a new perspective. J. Petrol. 49, 633–664. Niu, Y., Langmuir, C.H., Kinzler, R.J., 1997. The origin of abyssal peridotites: a new prospective. Earth Planet. Sci. Lett. 152, 251–265. Niu, Y., Batiza, R., 1993. Chemical variation trends at fast and slow spreading midocean ridges. J. Geophys. Res. 98, 7887–7902. Pariso, J.E., Sempéré, J.C., Rommevaux, C., 1995. Temporal and spatial variations in crustal accretion along the Mid-Atlantic Ridge (29°–31°300 N) over the last 10 m.y.: implications from a three-dimensional gravity study. J. Geophys. Res. 100, 17781–17794. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Pub. Corp., McGraw-Hill, New York, 197 pp. Plank, T., Langmuir, C.H., 1992. Effects of the melting regime on the composition of the oceanic crust. J. Geophys. Res. 97, 19749–19770. Presnall, D.C., Gudfinnsson, G.D., Walter, M.J., 2002. Generation of mid-ocean ridge basalts at pressures from 1 to 7 GPa. Geochim. Cosmochim. Acta 66, 2073–2090. Press, W.H., Flannery, B.P., Teukolosky, S.A., Vetterling, W.T., 1992. Numerical Recipes,, Cambridge Univ. Press, New York, 963 pp. Rabinowicz, M., Toplis, M.J., 2009. Melt Segregation in the lower part of the partially molten mantle zone beneath an oceanic spreading centre: numerical modelling of the combined effects of shear segregation and compaction. J. Petrol. 50, 1071–1106. Rabinowicz, M., Ceuleneer, G., 2005. The effect of sloped isotherms on melt migration in the shallow mantle: a physical and numerical model based on observations in the Oman ophiolite. Earth Planet. Sci. Lett. 229, 231–246. Rabinowicz, M., Nicolas, A., Vigneresse, J.L., 1984. A rolling mill effect in asthenosphere beneath oceanic spreading centers. Earth Planet. Sci. Lett. 67, 97–108. Ribe, N.M., 1988. On the dynamics of mid-ocean ridges. J. Geophys. Res. 93, 429– 436. Richter, F.M., Daly, S.F., 1989. Dynamical and chemical effects of melting a heterogeneous source. J. Geophys. Res. 94, 12499–12510. Rubin, K.H., Sinton, J.M., 2007. Inferences on mid-ocean ridge thermal and magmatic structure from MORB compositions. Earth Planet. Sci. Lett. 260, 257–276. Ruedas, T., Schmeling, H., Marquart, G., Kreutzmann, A., Junge, A., 2004. Temperature and melting of a ridge-centered plume with application to Iceland; part I, dynamics and crust production. Geophys. J. Int. 158, 729–743. Schmeling, H., 2006. A model of episodic melt extraction for plumes. J. Geophys. Res. 111, B03202, doi: 2004JB003423. Schmeling, H., 2000. Partial melting and melt segregation in a convecting mantle. In: Bagdassarov, N., Laporte, D., Thompson, A.B. (Eds.), Physics and Chemistry of Partially Molten Rocks. Kluwer Academic Publ., Dordrecht, pp. 141–178. Scott, D.R., Stevenson, D.J., 1989. A self-consistent model of melting, magma migration and buoyancy-driven circulation beneath mid-ocean ridges. J. Geophys. Res. 94, 2973–2988.

Scott, D.R., Stevenson, D.J., 1984. Magma solitons. Geophys. Res. Lett. 11, 1161– 1164. Seyler, M., Bonatti, E., 1997. Regional-scale melt-rock interaction in lherzolitic mantle in the Romanche fracture zone (Atlantic ocean). Earth Planet. Sci. Lett. 146, 273–287. Shibata, T., Thompson, G., 1986. Peridotites from the Mid-Atlantic Ridge at 43°N and their petrogenetic relation to abyssal tholeiites. Contrib. Mineral. Petrol. 93, 144–159. Sinton, J.M., Detrick, R.S., 1992. Mid-ocean ridge magma chamber. J. Geophys. Res. 97, 197–216. Smith, W.R., Missen, R.W., 1991. Chemical Reaction Equilibrium Analysis: Theory and Algorithms. Krieger Pub. Comp., Malabar, FL, reprint ed., 364 pp. Snow, J.E., Dick, H.J.B., 1995. Pervasive magnesium loss by marine weathering of peridotite. Geochim. Cosmochim. Acta 59, 4219–4235. Spahr, C.W., Forsyth, D.W., 1998. Structure of the upper mantle under the EPR from waveform inversion of regional events. Science 280, 1227–1229. Sparks, D.W., Parmentier, E.W., 1991. Melt extraction from the mantle beneath oceanic spreading centers. Earth Planet. Sci. Lett. 105, 368–377. Spiegelman, M., 1993b. Flow in deformable porous media, part 1 simple analysis. J. Fluid Mech. 247, 17–38. Spiegelman, M., 1993a. Physics of melt extraction: theory, implications and applications. Phil. Trans. R. Soc. A 342, 23–41. Spiegelman, M., Kelemen, 2003. Extreme chemical variability as a consequence of channelized melt transport. Geochem. Geophys. Geosyst. 4. doi:10.1029/ 2002GC000336. Spiegelman, M., Kelemen, P.B., Aharonov, E., 2001. Causes and consequences of flow organization during melt transport: the reaction infiltration instability in compactible media. J. Geophys. Res. 106, 2061–2077. Spiegelman, M., McKenzie, D., 1987. Simple 2-D models for melt extraction at midocean ridges and island arcs. Earth Planet. Sci. Lett. 83, 137–152. Sramek, O., Ricard, Y., Bercovici, D., 2007. Simultaneous melting and compaction in deformable two-phase media. Geophys. J. Int. 168, 964–982. Su, W., Buck, W.R., 1993. Buoyancy effects on mantle flow under mid-ocean ridges. J. Geophys. Res. 98, 12191–12205. Takahashi, E., 1986. Melting of a dry peridotite KLB-1 up to 14 GPa: implications on the origin of peridotitic upper mantle. J. Geophys. Res. 91, 9367–9382. Takahashi, E., Shimazaki, T., Tsuzaki, Y., Yoshida, H., 1993. Melting study of a peridotite KLB-1 to 6.5 GPa, and the origin of basaltic magmas. Phil. Trans. Phys. Sci. Eng. 342, 105–120. Tirone, M., Ganguly, J., Morgan, J.P., 2009. Modelling petrological geodynamics in the Earth’s mantle. Geochem. Geophys. Geosyst. 10, doi: 2008GC002168. Toomey, D.R., Wilcock, W.S.D., Solomon, S.C., Hammond, W.C., Orcutt, J.A., 1998. Mantle seismic structure beneath the MELT region of the East Pacific Rise from P and S wave tomography. Science 280, 1224–1227. Tucholke, B.E., Lin, J., Kleinrock, M.C., Tivey, M.A., Reed, T.B., Goff, J., Jaroslow, G.E., 1997. Segmentation and crustal structure of the western Mid-Atlantic Ridge flank, 25°250 –27°100 N and 0.29 m.y. J. Geophys. Res. 102, 10203–10223. Van Den Bleeken, G., Müntener, O., Ulmer, P., 2010. Reaction processes between tholeiitic melt and residual peridotite in the uppermost mantle: an experimental study at 08 GPa. J. Petrol. 51, 153–183. Van Zeggeren, F., Storey, S.H., 1970. The computation of Chemical Equilibria. Cambridge Univ. Press, London, 176 pp. Wagner, T.P., Grove, T.L., 1998. Melt/harzburgite reaction in the petrogenesis of tholeiitic magma from Kilauea volcano, Hawaii. Contrib. Mineral. Petrol. 131, 1– 12. Walter, M.J., 1998. Melting of garnet peridotite and the origin of komatiite and depleted lithosphere. J. Petrol. 39, 29–60. Walter, M.J., Presnall, D.C., 1994. Melting behavior of simplified lherzolite in the system CaO–MgO–Al2O3–SiO2–Na2O from 7 to 35 kbar. J. Petrol. 35, 329–359. Webb, S.C., Forsyth, D.W., 1998. Structure of the upper mantle under the EPR from waveform inversion of regional events. Science 280, 1227–1229. Workman, R.K., Hart, S.R., 2005. Major and trace element composition of the depleted MORB mantle (DMM). Earth Planet. Sci. Lett. 231, 53–72. Yang, Y., Forsyth, D.W., Weeraratne, D.S., 2007. Seismic attenuation near the East Pacific Rise and the origin of the low-velocity zone. Earth Planet. Sci. Lett. 258, 260–268. Yang, H.J., Sen, G., Shimizu, N., 1998. Mid-ocean ridge melting: constraints from lithospheric xenoliths at Oahu, Hawaii. J. Petrol. 39, 277–295. Zalesak, S.T., 1979. Fully multidimensional flux-corrected transport algorithms for fluids. J. Comput. Phys. 31, 335–362. Zhang, Y.S., Tanimoto, T., Stolper, E.M., 1994. S-Wave velocity, basalt chemistry and bathymetry along the mid-atlantic ridge. Phys. Earth Planet. Inter. 84, 79–93. Zimmerman, M.E., Kohlstedt, D.L., 2004. Rheological properties of partially molten lherzolite. J. Petrol. 45, 275–298.