## Abstract

We present a quantitative assessment of the thermal and dynamic response of an amphibolitic lower crust to the intrusion of basaltic dike swarms in an arc setting. We consider the effect of variable intrusion geometry, depth of intrusion, and basalt flux on the production, persistence, and interaction of basaltic and crustal melt in a stochastic computational framework. Distinct melting and mixing environments are predicted as a result of the crustal thickness and age of the arc system. Shallow crustal (∼30 km) environments and arc settings with low fluxes of mantle-derived basalt are likely repositories of isolated pods of mantle and crustal melts in the lower crust, both converging on dacitic to rhyodacitic composition. These may be preferentially rejuvenated in subsequent intrusive episodes. Mature arc systems with thicker crust (∼50 km) produce higher crustal and residual basaltic melt fractions, reaching ∼0·4 for geologically reasonable basalt fluxes. The basaltic to basaltic andesite composition of both crustal and mantle melts will facilitate mixing as the network of dikes collapses, and Reynolds numbers reach 10−4–1·0 in the interiors of dikes that have been breached by ascending crustal melts. This may provide one mechanism for melting, assimilation, storage and homogenization (MASH)-like processes. Residual mineral assemblages of crust thickened by repeated intrusion are predicted to be garnet pyroxenitic, which are denser than mantle peridotite and also generate convective instabilities where some of the crustal material is lost to the mantle. This reconciles the thinner than predicted crust in regions that have undergone a large flux of mantle basalt for a prolonged period of time, and helps explain the enrichment of incompatible elements such as K2O, typical of mature arc settings, without the associated mass balance problem.

## INTRODUCTION

The genesis of new continental crust in arc settings is ultimately driven by mantle melting and injection of basalt into the crust; however, the controls on basalt–crust interaction remain poorly understood. Basaltic magma transports both mass and enthalpy, and the petrological diversity that occurs when basalt reaches crustal depths has been postulated to originate from closed-system fractionation of the primitive basalt (Grove et al., 2003), crustal melting (Fornelli et al., 2002; Saleeby et al., 2003) and intermediate mixtures of the two processes (Anderson, 1976; DePaolo et al., 1992; Feeley et al., 2002). The range of processes reflects a continuum of basalt–crust interactions occurring at different pressure–temperature conditions, and lithological variations in the crust. However, quantification of the importance of crustal melting as a result of the variable style, flux, location and temporal response of basalt intrusion remains elusive.

Geological and geophysical data provide an incomplete picture of the geological expressions of the interaction of basalt with the lower crust (at depths greater than 25 km). Datasets that provide information about basaltic interaction with the lower crust include: (1) seismic velocity changes indicating a mafic lower crust (Furlong & Fountain, 1986; Rudnick, 1990; Holbrook et al., 1992; Ducea et al., 2003); (2) mafic xenoliths with mineral assemblages indicative of higher pressure (DeBari et al., 1987; Rudnick & Taylor, 1987; Ducea & Saleeby, 1998; Lee et al., 2001); (3) outcrops of lower crustal terrains in the Sierra Nevada, California; Fiordland, New Zealand; Kohistan, Pakistan; the Chipman Dikes, Saskatchewan; and in the Peninsular Terrane, Alaska (Jan & Howie, 1981; DeBari & Coleman, 1989; Pickett & Saleeby, 1993; Williams et al., 1995). Additionally, the geochemical characteristics of erupted magmas, such as trace element concentrations attributed to residual garnet (Hildreth & Moorbath, 1988), isotopic data (Griffin et al., 2002; Hart et al., 2002), and major element geochemical trends indicative of the suppression of plagioclase crystallization (Green, 1982; Grove et al., 2003), have all been interpreted to imply either crystallization or melting processes occurring at lower crustal pressures. In the light of geochemical data from the southern volcanic zone of the Andes, Hildreth & Moorbath (1988) suggested that the lower crust may be a region of enhanced melting and mixing of mantle-derived and crustally derived magmas. The paradigm of lower crustal melting, assimilation, storage and homogenization (MASH) was proposed to describe the along-arc variability they observed as a function of crustal thickness, and this concept has since been applied to other arc settings (Hopson & Mattinson, 1994; Kobayashi & Nakamura, 2001; Hart et al., 2002). Additional support for enhanced melt production in the lower crust is provided by the predictions of steady-state geotherms that approach the solidus of many lithologies in the lower crust (Chapman & Furlong, 1992), making these regions prone to melting with the addition of enthalpy supplied by mafic magmas.

Attempts to generalize basalt–crust interaction are usually presented in terms of two end-members. The first involves intrusion of large, spatially coherent bodies of basaltic magma, commonly represented in the geological record as gabbro–diorite–norite plutons. The process of crustal melting and contact metamorphism associated with the stages of intrusion and assembly of this type has been well described (Barboza & Bergantz, 2000). Notably, the most significant thermal impact of this process can be limited to a modest contact aureole (Barboza & Bergantz, 2000), despite the large volumes of basaltic material involved. Hence there is little geological support for the notion that the assembly of mafic complexes into a contiguous body leads to substantial crustal melting or regionally extensive metamorphism (Barboza et al., 1999). Alternatively, basaltic input as dike swarms, overlapping in space and time, may provide for a more efficient means of crustal melting and magma mingling (Hopson & Mattinson, 1994; Bergantz, 1995). However, a comprehensive, quantitative assessment of this form of basalt–crust interaction has not been made.

The objective of this study is to illuminate some of the thermal and dynamic consequences of basalt intrusion as dike swarms. To address a variety of possible intrusion configurations, we consider the random intrusion of basaltic dikes, and the effect of dike geometry, depth of intrusion and basalt flux on the crustal melting efficiency and the persistence and interaction of both crustal and basaltic melt through time. From this analysis, major element geochemical trends can be predicted for the developing thermal environments. We will further consider the criteria for crustal-scale density instability and stratification in the context of a crust that is actively growing through mafic addition. Lastly, we will evaluate the length and time scales associated with mixing and mingling that result from different thermal environments in the lower crust, to quantitatively assess MASH-like processes.

### Basalt flux in arc environments

Estimates of basalt flux in arc systems provide a global constraint for models aimed at understanding basalt–crust interaction. Gravity and seismic data have been used to estimate crustal thickness in some island arcs (Crisp, 1984; Dimalanta et al., 2002). These estimates can be divided by the age of the arc system since the initiation of subduction to yield an average volume flux into the crust. These data are presented as volume flux through area per unit time in Table 1. A second method of calculating basalt flux utilizes estimates of the amount of crustal assimilant in erupted magmas and the enthalpy associated with basalt intrusion required to melt this amount of crust (Grunder, 1995). The two methods of estimating the basalt flux yield results that are within an order of magnitude of each other from ∼1·0 × 10−4 to ∼1·0 × 10−3 m3/m2 per year.

Table 1:

Estimates of basalt flux into the lower crust

Location

Estimate of basalt flux (m3/m2 per year)

References

Gravity–seismic method
Marianas 4·93 × 10−4 Dimalanta et al. (2002)
Marianas 1·92 × 10−4 Crisp (1984)
Izu–Bonin 4·89 × 10−4 Dimalanta et al. (2002)
Aleutians 5·46 × 10−4 Dimalanta et al. (2002)
Aleutians 3·40 × 10−4 Crisp (1984)
Tonga 7·41 × 10−4 Dimalanta et al. (2002)
New Hebrides 1·04 × 10−3 Dimalanta et al. (2002)
Kuril 4·72 × 10−4 Crisp (1984)
Geochemical–thermal method
Eastern Nevada 4·0 × 10−4 Grunder (1995)
Location

Estimate of basalt flux (m3/m2 per year)

References

Gravity–seismic method
Marianas 4·93 × 10−4 Dimalanta et al. (2002)
Marianas 1·92 × 10−4 Crisp (1984)
Izu–Bonin 4·89 × 10−4 Dimalanta et al. (2002)
Aleutians 5·46 × 10−4 Dimalanta et al. (2002)
Aleutians 3·40 × 10−4 Crisp (1984)
Tonga 7·41 × 10−4 Dimalanta et al. (2002)
New Hebrides 1·04 × 10−3 Dimalanta et al. (2002)
Kuril 4·72 × 10−4 Crisp (1984)
Geochemical–thermal method
Eastern Nevada 4·0 × 10−4 Grunder (1995)

Basalt flux range employed in this study was 1·0 × 10−4 to 5·0 × 10−3 m3/m2 per year.

Both the seismic and gravity estimates, and the enthalpy balance calculations are probably underestimates of the amount of basalt flux. The seismic and gravity studies do not incorporate the loss of mass caused by erosion, or the lateral flow of material in the upper mantle–lower crust. Likewise, the enthalpy calculation assumes near-perfect efficiency in the transfer of enthalpy, and as such represents a minimum end-member for the amount of basalt required for melting.

If the range of estimated basalt fluxes is approximately constant in all arcs, the volume of material predicted to accumulate in mature arcs is substantial. Erosion (Montgomery et al., 2001), lower crustal flow (Meissner & Mooney, 1998), delamination (Kay et al., 1992; Lee et al., 2001; Saleeby et al., 2003), and Rayleigh–Taylor type density instability (Jull & Kelemen, 2001) have all been called upon as mechanisms to remove crustal material, and several of these mechanisms may operate simultaneously. The observation that melts leaving the mantle wedge are probably basaltic in composition, and that the bulk crust is andesitic (Kay et al., 1992), adds the further constraint that some removal of material from the crust must be weighted toward the more mafic components. Although the details of the mechanisms differ, Bird (1979), Kay & Kay (1993), and Jull & Kelemen (2001) have suggested that garnet-rich assemblages can develop greater densities than mantle material, and Jull & Kelemen (2001) demonstrated that ductile dripping instabilities can form on a time scale of 107 years, provided that the underlying mantle temperature exceeds 700°C. Thus any crustal-scale model of basalt–crust interaction must address the mass balance relationship between basaltic input and crustal thickness.

### Developing a rationalizing framework for assessing the thermal state of the lower crust

The quantitative assessment of basaltic underplating of the crust, and the resulting melting and mixing in the lower crust, have been considered in a number of numerical studies (Table 2). Analysis of the variety of results associated with these simulations provides motivation for the more general, stochastic approach applied in this study. The models can be divided into two groups: two-dimensional simulations of isolated sections of the crust, which allow for convective transport, typically with constant temperature boundary conditions (Bittner & Schmeling, 1995; Barboza & Bergantz, 1996; Raia & Spera, 1997), and one-dimensional conduction simulations (Younker & Vogel, 1976; Wells, 1980; Bergantz, 1989; Petford & Gallagher, 2001; Annen & Sparks, 2002). Other hybrid approaches, such as that of Huppert & Sparks (1988), use a one-dimensional parameterized convection model.

Table 2:

Thermal models of the lower crust

Study

Model type*

Magma intrusion style

Intrusion rate (m/year)

Total intruded magma (km)

Tinit (°C)

Depth intruded (km)

Lithology

Liquidus, solidus (°C)

Super-heat (K)

Max. crustal melt§

Efficiency (%)

Younker & Vogel (1976) 1-D, conduction, no bottom heat loss Single intrusion NA 2·0 500 NA Basalt L: 1200, S: 1100 200·0 H = 375 m 32
Biotite granite L: 1100, S: 800  @ 40 kyr
Wells (1980) 1-D, conduction, over-accretion Multiple intrusion 0·004 40·0 200 10·0 Tonalite L: 1050, S: 800 –25·0 H = 1237·5 m
@ 1·6 Myr
Huppert & Sparks (1988) 1-D, parameterized convection, no bottom heat loss Single intrusion NA 0·5 500 NA Basalt L: 1200, S: 1091 0·0 H = 220 m 44
Granodiorite L: 1000, S: 850  @ 100 years
Bergantz (1989) 1-D, conduction, no bottom heat loss Single intrusion NA 16·6 700 30–35 Basalt L: 1250, S: 980 0·0 V = 500 km3 38
Pelite L: 1200, S: 725
Bittner & Schmeling (1995) 2-D, convection Single intrusion NA 5·0 756 27·0 Basalt L: 1100, S: 950 0·0 NA NA
Granite L: 1050, S: 760
Barboza & Bergantz (1996) 2-D, convection Fixed T bottom boundary NA NA 600 30–35 Pelite L: 1200, S: 750 0·0 H = 1250 m NA
@ 79·3 kyr
Tb/Tint = 1·66
Raia & Spera (1997) 2-D, convection Fixed T bottom boundary NA NA 1195 NA (CaAl2Si2O8–CaMgSi2O6L: 1547 NA NA NA
S: 1277
Tb/Tint = 1·10
Pedersen et al. (1998) 1-D, conduction, over-accretion Multiple intrusion 0·002 10·0 650 35 Basalt L: 1250, S: 1100 0·0 H = 650 m
Granodiorite L: 1000, S: 710  @ 5·0 Myr
Petford & Gallagher (2001) 1-D, conduction, over-accretion Multiple intrusion 1·0 1·0 650 35·5 Basalt L: 1250, S: 1050 50·0 H = 50 m
Amphibolite L: 1075, S: 1010  @ 1 kyr
Annen & Sparks (2002) 1-D, conduction, over-accretion 50 m sills/10 kyr Multiple intrusion 0·005 8·0 600 30·0 Basalt L: 1300, S: 620 0·0 H = 400 m
Amphibolite L: 1075, S: 1010  @ 1·6 Myr
Study

Model type*

Magma intrusion style

Intrusion rate (m/year)

Total intruded magma (km)

Tinit (°C)

Depth intruded (km)

Lithology

Liquidus, solidus (°C)

Super-heat (K)

Max. crustal melt§

Efficiency (%)

Younker & Vogel (1976) 1-D, conduction, no bottom heat loss Single intrusion NA 2·0 500 NA Basalt L: 1200, S: 1100 200·0 H = 375 m 32
Biotite granite L: 1100, S: 800  @ 40 kyr
Wells (1980) 1-D, conduction, over-accretion Multiple intrusion 0·004 40·0 200 10·0 Tonalite L: 1050, S: 800 –25·0 H = 1237·5 m
@ 1·6 Myr
Huppert & Sparks (1988) 1-D, parameterized convection, no bottom heat loss Single intrusion NA 0·5 500 NA Basalt L: 1200, S: 1091 0·0 H = 220 m 44
Granodiorite L: 1000, S: 850  @ 100 years
Bergantz (1989) 1-D, conduction, no bottom heat loss Single intrusion NA 16·6 700 30–35 Basalt L: 1250, S: 980 0·0 V = 500 km3 38
Pelite L: 1200, S: 725
Bittner & Schmeling (1995) 2-D, convection Single intrusion NA 5·0 756 27·0 Basalt L: 1100, S: 950 0·0 NA NA
Granite L: 1050, S: 760
Barboza & Bergantz (1996) 2-D, convection Fixed T bottom boundary NA NA 600 30–35 Pelite L: 1200, S: 750 0·0 H = 1250 m NA
@ 79·3 kyr
Tb/Tint = 1·66
Raia & Spera (1997) 2-D, convection Fixed T bottom boundary NA NA 1195 NA (CaAl2Si2O8–CaMgSi2O6L: 1547 NA NA NA
S: 1277
Tb/Tint = 1·10
Pedersen et al. (1998) 1-D, conduction, over-accretion Multiple intrusion 0·002 10·0 650 35 Basalt L: 1250, S: 1100 0·0 H = 650 m
Granodiorite L: 1000, S: 710  @ 5·0 Myr
Petford & Gallagher (2001) 1-D, conduction, over-accretion Multiple intrusion 1·0 1·0 650 35·5 Basalt L: 1250, S: 1050 50·0 H = 50 m
Amphibolite L: 1075, S: 1010  @ 1 kyr
Annen & Sparks (2002) 1-D, conduction, over-accretion 50 m sills/10 kyr Multiple intrusion 0·005 8·0 600 30·0 Basalt L: 1300, S: 620 0·0 H = 400 m
Amphibolite L: 1075, S: 1010  @ 1·6 Myr
*

Dimension, heat conduction or convection.

Intruded magma listed first, then country rock.

Intruding magma, physical configuration of intrusion (sill, etc.) or specified temperature boundary condition.

§

Integrated melt volume/basal area, or for 1-D models integrated melt height.

To compare models that consider different geometric configurations and lithologies, a rationalizing framework was developed to evaluate the efficiency of the melting process. A completely efficient melting process is defined such that all of the enthalpy associated with cooling and crystallizing a volume of intruded magma is used to heat and melt only the volume of crust that becomes molten. This is efficient because no heat is ‘wasted’ heating regions of the crust that do not become molten, and all of the enthalpy is instantaneously applied. A similar approach has been described in detail by Grunder (1995). We define the efficiency of the melting process in these studies with respect to the efficient end-member:

(1)
$E\%\ =\ 100\ {\times}\ V_{\mathrm{crust}}^{\mathrm{mod.}}/V_{\mathrm{crust}}^{\mathrm{eff.}}$
where
$$V_{\mathrm{crust}}^{\mathrm{mod.}}$$
is the reported melt volume (or length reported in one-dimensional simulations) and
$$V_{\mathrm{crust}}^{\mathrm{eff.}}$$
is the volume of the efficient end-member predicted using the volume of intruded basalt, and the thermal parameters (liquidus, solidus, thermal diffusivity, latent heat and heat capacities) reported in the respective studies.

The predicted efficiency of the melting process corresponds closely to the geometric configuration assumed by the models. One-dimensional, vertically stacked, over-accretion conduction models that cool on both sides of the intrusions are 4–8% efficient (Wells, 1980; Pedersen et al., 1998; Petford & Gallagher, 2001; Annen & Sparks, 2002). One-dimensional conduction simulations that cool on one side only (insulated boundary condition on the other side) are 32–38% efficient (Younker & Vogel, 1976). The parameterized convection model (Huppert & Sparks, 1988) with no bottom heat loss is 44% efficient. In the light of the predicted melting efficiencies, uncertainty still exists as to which modeling assumptions best represent the melting process in the lower crust, especially considering the more complex intrusion geometries likely in natural settings. To avoid a priori specification of an assumed intrusion configuration, we have adopted the approach of treating the intrusion process as stochastic. Multiple simulations with different dike geometries can be combined to determine the average melting behavior of the crust for a given flux of basalt.

## STOCHASTIC DIKE INTRUSION MODEL

### Conservation equations

A two-dimensional model has been developed to illuminate the possible progression of thermal and compositional heterogeneities following basalt intrusion as random dikes in the lower crust. Two sets of simulations were performed: one with only heat transfer, but no bulk flow between the basaltic dike intrusions and the crust, and another in which heat transfer, local fluid (melt plus crystals) motion, and ductile creep are considered. The two-dimensional forms of the conservation equations per unit volume are as follows.

Conservation of enthalpy:

(2)
$\frac{{\partial}H}{{\partial}t}\ +\ \frac{{\partial}}{{\partial}\mathbf{x}_{j}}\ \left(\mathbf{v}_{j}\ H\right)\ =\ \frac{{\partial}}{{\partial}\mathbf{x}_{i}}\ k_{\mathrm{mix}}\ \left(\frac{{\partial}}{{\partial}\mathbf{x}_{i}}\ T\right)$
where
Conservation of mass:
(4)
$\frac{{\partial}\mathbf{v}_{j}}{{\partial}\mathbf{x}_{j}}\ =\ 0.$
Momentum:
(5)
$\frac{{\partial}}{{\partial}t}\left({\rho}_{\mathrm{mix}}\ \mathbf{v}_{i}\right)\ +\ \frac{{\partial}}{{\partial}\mathbf{x}_{j}}\left({\rho}_{\mathrm{mix}}\ \mathbf{v}_{j}\ \mathbf{v}_{i}\right)\ =\ {-}\frac{{\partial}P}{{\partial}\mathbf{x}_{i}}\ +\ {\mu}_{\mathrm{mix}}\ \frac{{\partial}^{2}\mathbf{v}_{i}}{{\partial}\mathbf{x}_{j}^{2}}\ +\ {\rho}_{\mathrm{mix}}\mathbf{g}.$
Species conservation:
(6)
$\frac{{\partial}{\gamma}}{{\partial}t}\ +\ \frac{{\partial}}{{\partial}\mathbf{x}_{j}}\left(\mathbf{v}_{j}{\gamma}\right)\ =\ \frac{{\partial}}{{\partial}\mathbf{x}_{i}}D\left(\frac{{\partial}}{{\partial}\mathbf{x}_{i}}\ {\gamma}\right).$
Side boundary conditions:
(7)
$\frac{{\partial}T}{{\partial}\mathbf{x}_{1}}\ =\ 0;{\ }\frac{{\partial}\mathbf{v}_{2}}{{\partial}\mathbf{x}_{1}}\ =\ 0;{\ }\mathbf{v}_{1}\ =\ 0.$

Table 3 has a list of symbols and nomenclature used. Einstein summation is implied for repeated vector indices. The advective (second from the left) term in the momentum equation was retained because it is important in describing chaotic advection (at Reynolds number ∼10−2–100) that may occur in melt-dominated regions.

Table 3:

Key to nomenclature

Symbol

Parameter

Units

H enthalpy
t time
k thermal conductivity W/m K
c specific heat J/kg K
T temperature
f melt fraction
L latent heat J/kg
g gravity m/s2
M mean value of basalt melt
N number of model realizations
vi velocity m/s
CMF critical melt function
D chemical diffusivity m2/s
$${\mu}_{\mathrm{m}}^{\mathrm{i}}$$

melt dynamic viscosity Pa s
μmix mixture dynamic viscosity Pa s
μparam. viscosity parameterization for f < CMF Pa s
$${\rho}_{\mathrm{c}}^{\mathrm{l}}$$

amphibolite melt density kg/m3
$${\rho}_{\mathrm{c}}^{\mathrm{s}}$$

solid amphibolite density kg/m3
$${\rho}_{\mathrm{b}}^{\mathrm{l}}$$

basalt melt density kg/m3
$${\rho}_{\mathrm{b}}^{\mathrm{s}}$$

solidified basalt density kg/m3
γ composition variable
ζ% percent change in mean
Symbol

Parameter

Units

H enthalpy
t time
k thermal conductivity W/m K
c specific heat J/kg K
T temperature
f melt fraction
L latent heat J/kg
g gravity m/s2
M mean value of basalt melt
N number of model realizations
vi velocity m/s
CMF critical melt function
D chemical diffusivity m2/s
$${\mu}_{\mathrm{m}}^{\mathrm{i}}$$

melt dynamic viscosity Pa s
μmix mixture dynamic viscosity Pa s
μparam. viscosity parameterization for f < CMF Pa s
$${\rho}_{\mathrm{c}}^{\mathrm{l}}$$

amphibolite melt density kg/m3
$${\rho}_{\mathrm{c}}^{\mathrm{s}}$$

solid amphibolite density kg/m3
$${\rho}_{\mathrm{b}}^{\mathrm{l}}$$

basalt melt density kg/m3
$${\rho}_{\mathrm{b}}^{\mathrm{s}}$$

solidified basalt density kg/m3
γ composition variable
ζ% percent change in mean

In multiphase regions of melt and solid, local relative motion between crystals and melt is not considered, and the thermal and physical transport properties of the system are developed as volume-weighted mixtures of material properties using the local composition and melt fraction parameters. Details of this procedure and physical and thermal properties of the magma and solid are described in the Appendix.

### Thermodynamic closure of conservation equations: melt fraction to enthalpy relationship

The geological complexity of multiphase solidification and melting can be incorporated into the enthalpy model by invoking suitable functions relating the melt fraction to the local enthalpy (Bergantz, 1990; Barboza & Bergantz, 1997). Appropriate forms of the thermodynamic closure are obtained by parameterizing experimental data that link the melt fraction of a particular lithology to temperature.

#### Amphibolite melt fraction relationship

A mafic, amphibolitic composition provides an end-member proxy for the composition of the lower crust in arc systems. In response to basaltic thermal input, the amphibolite may experience a dehydration reaction if there is an absence of a free vapor phase (Sen & Dunn, 1994)—this is the so-called ‘damp melting’. Several studies have examined this reaction (Beard & Lofgren, 1991; Rapp et al., 1991; Rushmer, 1991; Sen & Dunn, 1994; Wolf & Wyllie, 1994; Rapp & Watson, 1995) and are used to parameterize the melt fraction as a function of temperature (Fig. 1). The melt fraction and the mode of the residuum determine the major element composition of the melts, their physical properties, as well as the partitioning of trace elements (Patiño-Douce & Johnston, 1991; Barboza & Bergantz, 1997; Ducea, 2002).

Fig. 1.

Melt fraction of amphibolite as a function of temperature. Parameterized melt fraction functions based on the dehydration melting experiments of Beard & Lofgren (1991), Rushmer (1991), Patiño-Douce & Beard (1994), Sen & Dunn (1994), Wolf & Wyllie (1994), and Rapp & Watson (1995).

Fig. 1.

Melt fraction of amphibolite as a function of temperature. Parameterized melt fraction functions based on the dehydration melting experiments of Beard & Lofgren (1991), Rushmer (1991), Patiño-Douce & Beard (1994), Sen & Dunn (1994), Wolf & Wyllie (1994), and Rapp & Watson (1995).

Amphibolite dehydration reactions are pressure dependent, as reflected in the melt fractions and modal abundance of residual phases (Sen & Dunn, 1994). Departure from a linear melt-fraction to temperature relationship in amphibolites during dehydration results from the incongruent melting of amphibole ± plagioclase (Wolf & Wyllie, 1994) (Fig. 2). At greater than ∼10 kbar, amphibole and plagioclase generally react in a peritectic relationship to form an increase in the garnet and clinopyroxene mode as well as a tonalitic melt (Wolf & Wyllie, 1994).

Fig. 2.

Modal abundance (by volume) of an amphibolite undergoing dehydration melting. From the experiments of Wolf & Wyllie (1994).

Fig. 2.

Modal abundance (by volume) of an amphibolite undergoing dehydration melting. From the experiments of Wolf & Wyllie (1994).

To obtain a general melt fraction vs temperature relationship appropriate for the modeling of amphibolite melting, we combined experiments based on a range of mafic compositions with both calc-alkaline and tholeiitic trends. The functional form of the melt fraction curve was determined by parameterizing the data of Sen & Dunn (1994) at 15 kbar along with the projected liquidus for their composition using the MELTS thermodynamic package (Ghiorso & Sack, 1995):

(8)
\begin{eqnarray*}&&f\ =\ {-}2{\cdot}0968\ {\times}\ 10^{{-}12}\ \left(T^{{\ast}5}\right)\ +\ 1{\cdot}09308\ {\times}\ 10^{{-}8}\ \left(T^{{\ast}4}\right)\\&&{-}\ 2{\cdot}26718\ {\times}\ 10^{{-}5}\ \left(T^{{\ast}3}\right)\ +\ 2{\cdot}33912\ {\times}\ 10^{{-}2}\ \left(T^{{\ast}2}\right)\\&&{-}12{\cdot}0048\left(T^{{\ast}}\right)\ +\ 2451{\cdot}69\end{eqnarray*}
where
(9)
$T^{{\ast}}\ =\ T\ +\ 12{\cdot}0\left(15\ \mathrm{kbar}\ {-}\ P\right)\ \left(T\ \mathrm{in\ degrees}\ C,\ P\ \mathrm{in\ kbar}\right)$

(10)
$T_{\mathrm{solidus}}\ {\geq}\ T^{{\ast}}\ {\geq}\ T_{\mathrm{liquidus}}$

(11)
$T_{\mathrm{solidus}}\ =\ 0{\cdot}1495P^{2}\ +\ 11{\cdot}309P\ +\ 697{\cdot}86$
and
(12)
$T_{\mathrm{liquidus}}\ =\ 0{\cdot}25P^{2}\ +\ 35{\cdot}0P\ +\ 1260{\cdot}0.$
The ‘stair-step’ form of the melt fraction function corresponds to three stages of reactions (Fig. 1). The initial stage of melting consumes the small amount of quartz present (∼2 wt %) along with plagioclase and amphibole incongruent reaction to produce a significant increase in melt over a 50°C temperature range. This is followed by the continuing reaction of amphibole and plagioclase to form melt, garnet, and clinopyroxene. The final sequence of reactions corresponds to the melting of clinopyroxene and garnet. We acknowledge there is uncertainty because of the lack of experimental data at melt fractions greater than ∼0·5. However, as the thermal calculations demonstrate, these higher melt fractions are unlikely to occur as a result of anatexis in the lower crust.

For internal consistency, the functional form of the melt fraction was extrapolated to lower pressure using the solidus temperature at 10 kbar (Wolf & Wyllie, 1994) and fitting the data of Rapp & Watson (1995) and the quartz amphibolite data of Patiño-Douce & Beard (1995), both at ∼12 kbar. Below 10 kbar, garnet ceases to form as a reaction product (Beard & Lofgren, 1991; Rushmer, 1991). In general, the meta-basalts examined by Beard & Lofgren (1991) at 6·9 kbar had lower melt fractions than given by equation (8) (Fig. 1), although they fall within the range of melt fractions in the other experiments at higher pressure. This implies that at lower pressures the crust may be slightly less fertile than predicted by these calculations.

#### Wet basalt melt fraction relationship

The melt fraction–temperature relationship used for the wet basalt is modified from the experiments of Muntener et al. (2001) for a high Mg-number basaltic andesite at 12 kbar and 3·8 wt % initial water content (Fig. 3).

Fig. 3.

Melt fraction of wet basalt as a function of temperature. Parameterized melt fraction function based on the experiments of Green & Ringwood (1968), Green (1972, 1982), Holloway & Burnham (1972), Muntener et al. (2001), and Grove et al. (2003).

Fig. 3.

Melt fraction of wet basalt as a function of temperature. Parameterized melt fraction function based on the experiments of Green & Ringwood (1968), Green (1972, 1982), Holloway & Burnham (1972), Muntener et al. (2001), and Grove et al. (2003).

The solidus of the basalt is constrained by the experiments of Green & Ringwood (1968). Supplemental MELTS calculations have been performed for a similar composition and produce an equivalent melt fraction diagram (Ghiorso & Sack, 1995). Likewise the phases predicted by MELTS are very similar to those reported by Muntener et al. (2001) for the period from the start of crystallization until the formation of amphibole near the furthest extent of crystallization (0·39 melt fraction) reached in these experiments. The order of crystallization and comparison with MELTS calculations is depicted in the mode, Fig. 4, with clinopyroxene being the dominant near-liquidus phase, followed by garnet then amphibole. At ∼10–12 kbar the solidus of a wet tholeiitic basalt is at a minimum of ∼620°C (Green, 1982) and varies little in the range 8–15 kbar. The liquidus increases with increasing pressure at a rate of ∼5°C/kbar (Green, 1982).

Fig. 4.

Modal abundances of a crystallizing wet basaltic andesite (3·8 wt % H2O) as a function of melt fraction: (a) based on the experiments of Muntener et al. (2001); (b) calculated using the same starting composition as in (a) with the MELTS thermodynamic algorithm (Ghiorso & Sack, 1995).

Fig. 4.

Modal abundances of a crystallizing wet basaltic andesite (3·8 wt % H2O) as a function of melt fraction: (a) based on the experiments of Muntener et al. (2001); (b) calculated using the same starting composition as in (a) with the MELTS thermodynamic algorithm (Ghiorso & Sack, 1995).

### A stochastic dike intrusion model: the framework for evaluating a range of intrusion geometries

A stochastic, dike intrusion model was developed to study the thermal, rheological and chemical consequence of basalt–mafic crust interaction that may occur near the Mohorovicic discontinuity in arc settings. The random nature of the approach provides for a wide range of possible intrusion sequences, although still being constrained by the long-term basalt flux averages (Table 1). Conduction simulations were performed to explore melt production as a function of crustal thickness and basalt flux. Advection simulations were also performed to examine magma mixing in regions of high melt fraction, the mingling phenomena in the low melt fraction creeping regime, and the formation of density instabilities at the mantle–crust interface (Rudnick, 1990; Kay & Mahlburg-Kay, 1991).

The lower crust is idealized as depicted in Fig. 5. This two-dimensional geometry assumes that there is little variability in the perpendicular dimension to the modeled plane (along arc axis). The side boundaries are ‘reflecting’ so that lateral temperature gradients do not develop. The two-dimensional model is linked to a larger-scale one-dimensional model that allows the thermal anomaly created by the injection of magma to propagate both above and below the two-dimensional domain. The smallest scale that can be resolved in the two-dimensional grid is 1 m2 unless otherwise stated, and the one-dimensional domain has a resolution of 10 m.

Fig. 5.

Schematic representation of lower crustal conduction model. A two-dimensional model is linked to a one-dimensional conduction model that extends from the surface to the base of the mantle lithosphere. Side boundaries in the two-dimensional domain are reflecting. Dikes can intrude up to the specified maximum height of ascent.

Fig. 5.

Schematic representation of lower crustal conduction model. A two-dimensional model is linked to a one-dimensional conduction model that extends from the surface to the base of the mantle lithosphere. Side boundaries in the two-dimensional domain are reflecting. Dikes can intrude up to the specified maximum height of ascent.

The initial condition for the numerical simulations is determined by calculating a steady-state geotherm for a generic arc setting. The assumed initial heat flux at the surface is 68 mW/m2 and the temperature is fixed at 0°C. For comparison, this is approximately the current average heat flow in the Washington Cascades (Touloukian et al., 1981) and approximately equal to the global average heat flow of 65 mW/m2 (Petford & Gallagher, 2001). The initial, steady-state geotherm (Fig. 6) was calculated by the method of Chapman & Furlong (1992). Heat production from radioactive elements was assumed to be 0·94 mW/m3 at the surface, with exponential decrease with depth with a characteristic length scale of 15 km. Lower crustal heat production (below 15 km) was assumed to be 0·5 mW/m3 and the mantle region had a heat production of 0·02 mW/m3 (Chapman & Furlong, 1992).

Fig. 6.

The steady-state geotherm calculated by the method of Chapman & Furlong (1992). Surface heat flux is 68 mW/m2, approximately equal to the average heat flux in the Washington Cascades (Touloukian et al., 1981). This calculation is valid if the mechanical boundary layer is significantly thicker than the crustal thickness. Also shown is the inferred solidus of amphibolite from the experiments of Rushmer (1991), Wolf & Wyllie (1994), Patiño-Douce & Beard (1995), and Rapp & Watson (1995). Symbols are as in Fig. 1. Numbers adjacent to the symbols indicate the melt fraction.

Fig. 6.

The steady-state geotherm calculated by the method of Chapman & Furlong (1992). Surface heat flux is 68 mW/m2, approximately equal to the average heat flux in the Washington Cascades (Touloukian et al., 1981). This calculation is valid if the mechanical boundary layer is significantly thicker than the crustal thickness. Also shown is the inferred solidus of amphibolite from the experiments of Rushmer (1991), Wolf & Wyllie (1994), Patiño-Douce & Beard (1995), and Rapp & Watson (1995). Symbols are as in Fig. 1. Numbers adjacent to the symbols indicate the melt fraction.

We have modeled the injection of basalt as small aspect-ratio dikes as exemplified by the Chipman dikes and the Chelan complex (Hopson & Mattinson, 1994; Williams et al., 1995) (Fig. 5). Observations of outcrops with dikes at lower crustal pressures motivate the modeled fracture-assisted intrusion of mantle-derived basalt into a mafic, amphibolitic lower crust (Williams et al., 1995). This implicitly assumes instantaneous emplacement of the basalt, as the fracture mechanism is assumed to be much faster than subsequent viscous and thermal relaxation.

We will use the term ‘realization’ to denote a distinct episode of progressive basalt injection in the lower crust with multiple diking events. Ensemble averages of many different realizations of the model can then be used to describe the average behavior expected for a given flux of basalt into the lower crust, as well as the range of variability. Dike orientation in any particular realization is random to avoid a priori designation of dike geometry. The maximum height of ascent of a dike is fixed to simulate stress or material property changes that inhibit further dike propagation. Any intrusion up to the height of maximum ascent is possible (Fig. 5). The time interval between dike intrusion events is random, up to twice the specified average interval of intrusion, with an equal likelihood of a diking event occurring at any time up to this maximum value. The average flux of basalt is fixed in sets of simulations to facilitate comparison, but any particular realization of the model will have different dike geometries and variable intrusion histories.

The criteria for assessing the number of realizations that adequately describe this variability was defined such that adding another realization changes the mean value of the volume of basalt melt present by less than 0·01%:

(13)
${\xi}_{\%}\ =\ 100\left[{\int}_{0}^{t}\ M\left(N\right)\mathrm{d}t\ {-}\ {\int}_{0}^{t}\ M\left(N\ {-}\ 1\right)\mathrm{d}t\right]/{\int}_{0}^{t}\ M\left(N\right)\mathrm{d}t.$
Here ξ% denotes the percent change between realizations, and M is the mean for a given, N, number of realizations. This is schematically depicted in Fig. 7. For the criterion ξ% < 0·01, between 25 and 50 realizations are typically required for the conditions studied.

Fig. 7.

The percent change in the mean volume of basaltic melt (ξ) with the progressive incorporation of more realizations. Typically, 25–50 realizations are required to ensure ξ% < 0·01.

Fig. 7.

The percent change in the mean volume of basaltic melt (ξ) with the progressive incorporation of more realizations. Typically, 25–50 realizations are required to ensure ξ% < 0·01.

### Dike intrusion parameters

The effect of varying width, depth of intrusion, and maximum height of ascent of a propagating dike was assessed in separate suites of simulations. Constraints on the magnitude of these values were provided by outcrop and theoretical arguments. For example, observations of basaltic dike widths commonly less than 10 m motivated our choice of using widths of 1, 5 and 10 m in the simulations. Magma viscosity, magma over-pressure, and stress state of the host rock are all factors contributing to dike width (Fialko & Rubin, 1999). In a survey of dikes in SW Japan and Peru, Wada (1994) observed a correlation between the viscosity of the intruding magma and the width of dikes. For dikes with basaltic viscosity, widths of 1–10 m were most commonly observed. [A notable exception is the larger width dikes associated with flood basalt volcanism. The 100 m+ dike width observed in these regions has been attributed to greater magma over-pressure and extensive meltback (Fialko & Rubin, 1999).] Other observations suggest that similar controls on dike width operate at lower crustal pressures. The Chipman dikes provide evidence for dikes of tholeiitic basalt 1–10 m+ width at an inferred pressure of 10 kbar (Williams et al., 1995). Dikes from the upper mantle Balmuccia massif range in size from <1 cm to 1 m (Mukasa & Shervais, 1999). Although these observations do not restrict the possibility of larger width dikes, they do demonstrate that in many settings basalt dikes are limited to <10 m.

Crustal thicknesses between 30 and 50 km were considered, to constrain the melt production in a range of arc settings. A thicker crust will influence crustal melt production in two ways: (1) the ambient temperature will increase with depth; (2) the phase diagram, and hence the melt fraction diagram, is also a function of pressure. [For instance, garnet becomes a stable phase in the dehydration of amphibole after the pressure exceeds >10 kbar (Wolf & Wyllie, 1994).] This will have an important effect on the major and trace element composition, as well as potentially altering the volume of melts produced. The crustal thickness range of 30–50 km corresponds to the thickest island arc settings (Dimalanta et al., 2002), as well as many continental arc settings.

### Evaluating melt productivity

Two principal metrics will be applied to evaluate the degree of melt production. To facilitate comparison with previous one-dimensional simulations (Petford & Gallagher, 2001; Annen & Sparks, 2002) the sum of all melt fractions is normalized by the basal length of the simulations to give a melt length. This is equivalent to the ‘melt thickness’ reported by Petford & Gallagher (2001) or the ‘compaction thickness’ of Annen & Sparks (2002). Although the total amount of melt is important in determining the overall mass balance of magma entering, residing in and potentially leaving the lower crust, the local melt fractions determine the composition of any particular melt. At each time-step, there is a range of residual basalt and crustal melt fractions distributed throughout the modeled section of the lower crust. The basalt melt fractions vary from some minimum melt fraction that reflects basalt that has cooled to ambient conditions, up to a maximum melt fraction of 1·0 when the basalt intrudes at its liquidus. Similarly, a range of crustal melt fractions can coexist in one time-step from unmelted crust up to a maximum crustal melt fraction in regions in proximity to voluminous or recent basalt intrusion. For both the basalt and crustal material, the respective mean melt fraction is calculated by averaging the melt fractions of all areas that exceed their solidus.

## RESULTS

Depth of emplacement, maximum height of ascent of the intrusions, and basalt flux were varied in the conduction simulations to elucidate the range of residual basalt and crustal melt fractions and the amount of melt produced (Table 4). A second set of simulations, also considering advective transport, was conducted using the temperatures from conduction simulations as the initial and boundary conditions. Dynamic simulations over short (5·0 × 103 years) and long (107 years) time scales examined magma mixing and ductile creep, respectively.

Table 4:

Selected conduction results (maximum height of ascent is 100 m)

Time (years)

Flux (m3/m2 per year)

Dike width (m)

Depth (km)

Max. amp. melt fraction

Mean amp. fraction

Mean residual basalt melt fraction

Min. residual basalt melt fraction

Normalized amp. melt volume (m)

Normalized basalt melt volume (m)

Ratio (crustal melts/mantle melts)

S.D. basalt melt

S.D. crustal melt

Diff. length scale

106 0·0005 10 34 0·03 0·02 14·45 2·8 470
106 0·0005 34 0·03 0·02 15·6 3·1 332
106 0·0005 34 0·04 0·03 18·9 3·2 148
106 0·001 10 34 0·04 0·02 39·5 5·8 331
106 0·0005 10 44 0·08 0·04 0·09 0·06 118·9 43·3 2·74 6·2 10 470
106 0·001 34 0·05 0·04 44·5 6·3 234
106 0·001 10 34 0·05 0·04 52·1 6·8 331
106 0·0005 44 0·19 0·09 0·18 0·13 296·2 90·15 3·28 8·7 31 148
106 0·001 10 44 0·19 0·09 0·18 0·14 293·5 179·8 1·63 13· 28 331
106 0·001 44 0·20 0·10 0·18 0·15 307·6 183·4 1·67 16· 16 234
106 0·001 44 0·23 0·12 0·19 0·16 392·2 191·3 2·05 15· 22 104
106 0·005 10 34 0·10 0·05 492·5 44· 148
106 0·005 34 0·07 0·06 0·16 0·15 43·6 795·5 0·054 48· 11 47
106 0·005 10 44 0·32 0·22 0·21 0·16 662 1042·5 0·634 72· 42 148
106 0·005 44 0·33 0·22 0·21 0·16 683 1051·5 0·649 76· 25 105
106 0·005 44 0·47 0·28 0·23 0·26 1157·9 1139·5 1·016 66· 42 47
Time (years)

Flux (m3/m2 per year)

Dike width (m)

Depth (km)

Max. amp. melt fraction

Mean amp. fraction

Mean residual basalt melt fraction

Min. residual basalt melt fraction

Normalized amp. melt volume (m)

Normalized basalt melt volume (m)

Ratio (crustal melts/mantle melts)

S.D. basalt melt

S.D. crustal melt

Diff. length scale

106 0·0005 10 34 0·03 0·02 14·45 2·8 470
106 0·0005 34 0·03 0·02 15·6 3·1 332
106 0·0005 34 0·04 0·03 18·9 3·2 148
106 0·001 10 34 0·04 0·02 39·5 5·8 331
106 0·0005 10 44 0·08 0·04 0·09 0·06 118·9 43·3 2·74 6·2 10 470
106 0·001 34 0·05 0·04 44·5 6·3 234
106 0·001 10 34 0·05 0·04 52·1 6·8 331
106 0·0005 44 0·19 0·09 0·18 0·13 296·2 90·15 3·28 8·7 31 148
106 0·001 10 44 0·19 0·09 0·18 0·14 293·5 179·8 1·63 13· 28 331
106 0·001 44 0·20 0·10 0·18 0·15 307·6 183·4 1·67 16· 16 234
106 0·001 44 0·23 0·12 0·19 0·16 392·2 191·3 2·05 15· 22 104
106 0·005 10 34 0·10 0·05 492·5 44· 148
106 0·005 34 0·07 0·06 0·16 0·15 43·6 795·5 0·054 48· 11 47
106 0·005 10 44 0·32 0·22 0·21 0·16 662 1042·5 0·634 72· 42 148
106 0·005 44 0·33 0·22 0·21 0·16 683 1051·5 0·649 76· 25 105
106 0·005 44 0·47 0·28 0·23 0·26 1157·9 1139·5 1·016 66· 42 47

amp., amphibolite; S.D., standard deviation.

Table 5:

Selected conduction results and crustal melting efficiency (maximum height of ascent is 1000 m)

Time (years)

Flux (m3/m2 per year)

Dike width (m)

Depth (km)

Max. amp. melt fraction

Mean amp. melt fraction

Mean basalt melt fraction

Min. basalt melt fraction

Normalized amp. melt volume (m)

Normalized basalt melt volume (m)

Efficiency (%)

106 0·001 10 34 0·02 0·01 17·9
106 0·005 10 34 0·06 0·03 310
106 0·001 10 44 0·13 0·07 0·12 0·09 35·8 124 10·4
106 0·005 10 44 0·17 0·10 0·14 0·11 50·3 715 3·4
5 × 106 0·001 10 34 0·08 0·03 402·5
5 × 106 0·005 10 34 0·15 0·07 0·15 0·08 242·9 3757·5 0·5
5 × 106 0·001 10 44 0·35 0·17 0·22 0·08 693·9 1124·5 7·1
5 × 106 0·005 10 44 0·45 0·22 0·30 0·19 1038·6 7437·5 2·3
107 0·001 10 34 0·10 0·06 1044
107 0·005 10 34 0·18 0·08 0·16 0·12 532·0 8115 0·6
107 0·001 10 44 0·43 0·21 0·25 0·13 1044·5 2480 5·7
107 0·005 10 44 0·52 0·27 0·33 0·25 1479·1 16565 1·6
Time (years)

Flux (m3/m2 per year)

Dike width (m)

Depth (km)

Max. amp. melt fraction

Mean amp. melt fraction

Mean basalt melt fraction

Min. basalt melt fraction

Normalized amp. melt volume (m)

Normalized basalt melt volume (m)

Efficiency (%)

106 0·001 10 34 0·02 0·01 17·9
106 0·005 10 34 0·06 0·03 310
106 0·001 10 44 0·13 0·07 0·12 0·09 35·8 124 10·4
106 0·005 10 44 0·17 0·10 0·14 0·11 50·3 715 3·4
5 × 106 0·001 10 34 0·08 0·03 402·5
5 × 106 0·005 10 34 0·15 0·07 0·15 0·08 242·9 3757·5 0·5
5 × 106 0·001 10 44 0·35 0·17 0·22 0·08 693·9 1124·5 7·1
5 × 106 0·005 10 44 0·45 0·22 0·30 0·19 1038·6 7437·5 2·3
107 0·001 10 34 0·10 0·06 1044
107 0·005 10 34 0·18 0·08 0·16 0·12 532·0 8115 0·6
107 0·001 10 44 0·43 0·21 0·25 0·13 1044·5 2480 5·7
107 0·005 10 44 0·52 0·27 0·33 0·25 1479·1 16565 1·6

amp., amphibolite.

### Single realizations and stochastic ensemble behavior

The progress of a single realization illustrates the spatial distribution of intruded basalt and crustal melt fractions that can coexist at a single instant in time (Fig. 8). In this example calculation, the zone of intrusion is 100 m thick, the average basalt flux is 0·001 m3/m2 per year, and the crustal thickness is 44 km or at ∼12·5 kbar pressure assuming an overlying crust with an average density of 2900 kg/m3. The time after initiation of basalt intrusions is 50 000 years. Each dike is injected as a separate event, and the most recent dike intrusion and associated temperature perturbation can be identified in the lower left-hand corner of the simulated domain depicted in Fig. 8.

Fig. 8.

Example of a single realization of random dike intrusion. The zone of intrusion is 100 m thick, average basalt flux is 0·001 m3/m2 per year, crustal thickness is 44 km, and the time after the initiation of intrusions is 50 000 years. Dike position, temperature, and basalt and crustal melt fractions are shown.

Fig. 8.

Example of a single realization of random dike intrusion. The zone of intrusion is 100 m thick, average basalt flux is 0·001 m3/m2 per year, crustal thickness is 44 km, and the time after the initiation of intrusions is 50 000 years. Dike position, temperature, and basalt and crustal melt fractions are shown.

The unique intrusion history of each realization yields a range of melt volumes when an ensemble of simulations are combined (Fig. 9). The mean and standard deviations of melt volumes for the basalt and crustal melts from a selected group of simulations are presented in Table 4. The melt volume standard deviations range from 5 to 40% of the mean volumes. In general, larger percentage standard deviations occur at lower basalt intrusion rates for time scales <106 years. At lower fluxes of basalt (<0·001 m3/m2 per year), individual intrusions result in localized, transient thermal anomalies that are primarily responsible for melt production.

Fig. 9.

The mean, maximum and minimum volume of crustal amphibolite and residual basalt melt normalized by the basal area of the simulations (yielding a melt length) as a function of time. Basalt flux is 0·005 m3/m2 per year, dike thickness is 1 m, crustal thickness is 34 km, and the maximum height of ascent is 100 m. Other melt volumes and standard deviations are shown in Table 4.

Fig. 9.

The mean, maximum and minimum volume of crustal amphibolite and residual basalt melt normalized by the basal area of the simulations (yielding a melt length) as a function of time. Basalt flux is 0·005 m3/m2 per year, dike thickness is 1 m, crustal thickness is 34 km, and the maximum height of ascent is 100 m. Other melt volumes and standard deviations are shown in Table 4.

### Thermal evolution of the lower crust following basalt intrusion

The progressive emplacement of basaltic dikes increases the ambient temperature in the lower crust and a temperature anomaly develops with respect to the steady-state geotherm (Fig. 10). This phenomenon has been noted in several other numerical studies of intrusion in the lower crust (Wells, 1980; Pedersen et al., 1998; Annen & Sparks, 2002) and is simply understood as the input of thermal energy associated with successive basalt injection accumulating at a faster rate than the heat can be diffused. Hence, two thermal time scales are relevant in the heating of the crust through thin intrusions: (1) the time scale of diffusion for individual intrusions; (2) the time scale of diffusion for the broader thermal anomaly associated with the amalgamation of successive intrusions. The shorter wavelength, higher amplitude, thermal perturbations associated with the former are primarily responsible for the maximum basalt and crustal melt fractions, but also decay on time scales of

$${\sim}d_{\mathrm{w}}^{2}/{\kappa}$$
, where dw is the dike width and κ is the thermal diffusivity. For 1 m wide dikes this time scale is of the order of 10–20 days. The broad, low-amplitude thermal anomalies associated with the slow accumulation of basalt through successive intrusions decay on time scales more closely associated with the zone of active intrusion. For a 100 m thick zone of active intrusion diffusive time scales on the order of 400 years are applicable.

Fig. 10.

Temperature profiles as a function of depth. The shaded region is the domain of the two-dimensional models in which active intrusion is occurring. Extending beyond this domain is the one-dimensional conduction model. In all simulations the crust is 34 km thick, and two basalt fluxes are considered: 0·0001 and 0·005. In (a) the maximum height of ascent is 1000 m and dike width is 10 m; in (b) the maximum height of ascent is 100 m with 10 m dikes; in (c) the maximum height of ascent is 100 m with 1 m dikes.

Fig. 10.

Temperature profiles as a function of depth. The shaded region is the domain of the two-dimensional models in which active intrusion is occurring. Extending beyond this domain is the one-dimensional conduction model. In all simulations the crust is 34 km thick, and two basalt fluxes are considered: 0·0001 and 0·005. In (a) the maximum height of ascent is 1000 m and dike width is 10 m; in (b) the maximum height of ascent is 100 m with 10 m dikes; in (c) the maximum height of ascent is 100 m with 1 m dikes.

For basalt fluxes of 0·0001–0·005 m3/m2 per year the temperature near the zone of intrusions increases with time, and on average, the temperature following an intrusion does not have sufficient time to relax to the steady-state profile completely before the next intrusion. Figure 10 depicts the width and amplitude of the mean thermal anomalies for several scenarios. Two basalt fluxes (0·0001 m3/m2 per year and 0·005 m3/m2 per year) are considered to intrude crust 34 km thick. In scenario (a) the maximum height of ascent is 1000 m and dike width is 10 m; in scenario (b) the maximum ascent height is 100 m with 10 m wide dikes; in (c) the maximum ascent height of the dikes is 100 m with 1 m dikes. The mean temperature increases the fastest when the intrusions are narrower and concentrated in the smallest region of the crust. After 106 years following the initiation of basalt intrusion, only scenarios (b) and (c) with a basalt flux of 0·005 m3/m2 per year produce crustal melts. All model realizations with basalt fluxes of 0·0001 m3/m2 per year create thermal anomalies over 106 years with the ambient temperature increased only 20–40°C over the steady-state geotherm.

The amplitude of the thermal anomalies and the time required to reach the crustal solidus at different crustal depths varies greatly, and is an important constraint in determining the composition of melts from young vs mature arc systems. Whether or not crustal melting can occur depends on the thickness of the crust, the basalt intrusion rate, and the length of time the basalt intrusion rate remains approximately constant. The temperature difference between the solidus of the amphibolite and steady-state geotherm decreases with increasing depth (Fig. 6). As a consequence, for a given increase in temperature above the steady-state geotherm, amphibolite crust at depth will have a higher degree of melt. This is demonstrated in Fig. 11, in which the basalt flux is held constant at 0·001 m3/m2 per year, the zone of intrusion is 100 m thick, and dike widths are 1 m thick. For a crustal thickness of 30 km, the mean crustal melt fraction takes ∼1·3 × 107 years to reach 0·1. Even after 5·0 × 107 years of basalt intrusions the mean crustal melt fraction does not exceed 0·4 at this depth. For comparison, the period of time from the initiation of subduction to the present in several Pacific island arcs varies between ∼2·5 × 107 years and ∼5·0 × 107 years (Dimalanta et al., 2002). In contrast, mean crustal melt fractions exceeding 0·4 at 50 km depth occur after ∼1·5 × 106 years of intrusion for the same basalt flux. These calculations predict that production of significant volumes of crustal melt is much more likely in the thickened crust of mature, continental arcs.

Fig. 11.

Mean crustal and basalt melt fractions as a function of time and depth. All simulations use a flux of 0·001 m3/m2 per year, 100 m maximum height of ascent and 1 m dikes. For comparison, the times since the initiation of subduction for the New Hebrides, Tonga and Marianas island arcs are shown [all these regions have crust <30 km thick (Dimalanta et al., 2002)]. The stars mark conditions that were used for the initial conditions of the advection simulations examining magma mixing and mingling.

Fig. 11.

Mean crustal and basalt melt fractions as a function of time and depth. All simulations use a flux of 0·001 m3/m2 per year, 100 m maximum height of ascent and 1 m dikes. For comparison, the times since the initiation of subduction for the New Hebrides, Tonga and Marianas island arcs are shown [all these regions have crust <30 km thick (Dimalanta et al., 2002)]. The stars mark conditions that were used for the initial conditions of the advection simulations examining magma mixing and mingling.

At depths greater than ∼40 km, immediately following the initiation of intrusions, there will be greater volumes of crustal melt than residual basalt (Fig. 12). This can simply be explained by the greater initial temperatures from the steady-state geotherm. As basalt continues to intrude at these depths the volume of residual basalt will become greater than crustal volumes. At shallow depths (<40 km) the volume of mantle melt is always greater than the crustal melt volume.

Fig. 12.

Melt volume ratio (volume of crustal melt/volume of mantle melt) as a function of crustal depth and time. Conditions are the same as in Fig. 11. At greater depths, and just following the initiation of intrusions, greater volumes of crustal melt can be produced, whereas shallow and long-lived systems favor greater volumes of residual basaltic melt.

Fig. 12.

Melt volume ratio (volume of crustal melt/volume of mantle melt) as a function of crustal depth and time. Conditions are the same as in Fig. 11. At greater depths, and just following the initiation of intrusions, greater volumes of crustal melt can be produced, whereas shallow and long-lived systems favor greater volumes of residual basaltic melt.

### Style of emplacement of basalt dikes: maximum height of ascent

The degree to which basalts can penetrate the overlying crust will influence the efficiency of crustal melting. Previous studies have focused on an end-member intrusion geometry where basalt accumulates in sill-like bodies at the base of the crust, and each successive basaltic sill is emplaced on top of the other sills. However, it is clear that primitive magmas reach the surface (Cole, 1982; Muntener et al., 2001), motivating further examination of deviations from the over-accretion assumption as used in previous models (Table 2). To assess the sensitivity of the ascent of basaltic magma to different crustal levels, the maximum height of ascent of the dikes was varied while holding the basalt flux constant at 0·001 m3/m2 per year with 5 m wide dikes. The development of crustal melt and residual basalt melt volumes is depicted in Fig. 13 after periods of intrusion lasting 106 years.

Fig. 13.

Maximum height of ascent as a function of normalized melt volume. The normalized volume of melt for both the crustal melt and basalt is given as a function of the maximum height of ascent after 106 years of intrusion, a basalt flux of 0·001 m3/m2 per year and dike width of 5 m. Two crustal thicknesses are considered: 34 km and 44 km. Over-accretion maximizes the melt volume for the conditions examined.

Fig. 13.

Maximum height of ascent as a function of normalized melt volume. The normalized volume of melt for both the crustal melt and basalt is given as a function of the maximum height of ascent after 106 years of intrusion, a basalt flux of 0·001 m3/m2 per year and dike width of 5 m. Two crustal thicknesses are considered: 34 km and 44 km. Over-accretion maximizes the melt volume for the conditions examined.

The assumption of over-accretion maximizes the volumes of both residual basalt melt and crustal melts. For thick crust (>44 km) some crustal melting will occur if the dikes are intruding up to 5 km from the base of the crust. However, the crustal melt volume produced for such long range propagation decreases three orders of magnitude compared with over-accretion (Fig. 13). Basalt injections that propagate less than 100 m from the base of the crust produce a greater volume of crustal melt than residual basalt melt at 44 km depth over 106 years of basalt intrusion. However, as the maximum height of ascent increases and the basalt reaches shallower depths, a larger ratio of residual basalt will coexist with crustal melts.

Dikes ascending to shallower levels in the crust are less likely to produce significant crustal melting for three reasons. (1) The greater length of the dikes spreads the flux of enthalpy over greater areas, and hence if the flux of basalt is constant, the cumulative probability that dikes will overlap for a given period of time decreases. (2) The average interval of time between dike events is increased to maintain a constant flux, and with this increased time the longer dikes can diffuse their enthalpy over a greater area. A large amount of energy is expended increasing the temperature of a large volume of crust by only a few degrees C. (3) Dikes that reach shallower regions of the crust also encounter a larger temperature contrast between the basalt liquidus and the surrounding crust, which promotes solidification of the dikes without crustal melting.

Thin crust (∼34 km) is associated with a much smaller degree of crustal melting, provided that the lithosphere thickness is constant as implied by the steady-state geotherm calculations (McKenzie & Bickle, 1988). For all propagation lengths, the residual basalt volumes are greater than volumes of crustal melt. For maximum heights of ascent >100 m, only isolated patches of crustal melting occur, and no long-term reservoir of crustal melt is developed over 106 years. When the basalt dikes extend greater than ∼5 km from the base of this thinner crust, no residual basalt melt accumulates beyond the time scale required to conductively cool single intrusions. These calculations do not rule out dike events extending to the top of the crust and rapidly erupting material from depth. However, for this thermal environment, if dikes stall during their ascent no melt will be retained on time scales of the order of

$${\sim}d_{\mathrm{w}}^{2}/{\kappa}$$
.

### The effect of intrusion width on crustal melting

The width of the basalt intrusions will be influenced by the basalt magma viscosity, material properties of the surrounding crust and the magma over-pressure. Dike widths of 1, 5 and 10 m were considered at different crustal depths and basalt flux (Fig. 14) to characterize the melt productivity as a function of dike width.

Fig. 14.

Normalized volume of crustal melt and basalt melt as a function of dike width. Dikes of 1, 5 and 10 m width were modeled. The smaller width dikes maximize the melt volume in the stochastic simulations.

Fig. 14.

Normalized volume of crustal melt and basalt melt as a function of dike width. Dikes of 1, 5 and 10 m width were modeled. The smaller width dikes maximize the melt volume in the stochastic simulations.

For a given flux of basalt, smaller width dikes increase the melt volumes for both the residual basalt and the crustal melt. The average time interval between dike intrusions is much less for the 1 m dikes with constant intrusion rate, and consequently the probability of overlapping thermal anomalies of large amplitude is greater. In most scenarios 1 m dikes produced ∼5–15% greater melt than 10 m dikes. This effect on crustal melting is particularly pronounced for the highest flux rate (0·005 m3/m2 per year) and 44 km depth.

The more complex intrusion geometry, at least partially, explains the difference between the dependence on dike width in these simulations compared with the results of one-dimensional over-accretion models in which the location of intrusion is prescribed. Petford & Gallagher (2001) reported that crustal melting is maximized when the period between intrusions equals the diffusive heat loss time scale of the basalt sills. In their analysis, if heat is lost faster than the time scale for diffusion less crustal melting will result. Annen & Sparks (2002) reported little difference in crustal melt production for 10, 50 and 500 m sills at 20 km depth. The crucial difference between the over-accretion condition and our two-dimensional intrusion conditions is the probability of intersection of thermal anomalies. Over-accretion specifies that successive sill intrusions be adjacent to each other, and so the analysis of Petford & Gallagher (2001) applies. However, in the two-dimensional calculations the probability of intersection of dikes and their associated thermal anomalies is less than unity. Hence dike intrusion configurations that promote the greater intersection of dikes, and hence thermal anomalies, will develop greater crustal melting.

### Basalt flux and crustal thickness: primary controls on the crustal melting process

The two most important factors controlling the growth of thermal anomalies in response to the intrusion of magma are the depth of emplacement and the flux of the intruding magma. Flux was varied in this suite of simulations between 0·0001 and 0·005 m3/m2 per year to simulate the estimated range of basalt flux into the crust in arc settings (Table 1). The crustal thickness was varied between 30 and 50 km. A summary of crustal melt fractions and the residual basalt melt fractions at 106 years for 100 m intrusion zone and 1 m dikes is presented in Fig. 15.

Fig. 15.

Melt fraction of (a) crustal and (b) basaltic material as a function of depth and flux of basalt. The maximum height of ascent is 100 m, dike width 1 m and the information was extracted 106 years following the initiation of the modeled intrusions. For the amphibolite composition both the average and maximum melt fractions are shown. Melt fractions of the amphibolite in the simulations vary from zero to the maximum melt fraction. The mean and minimum basalt melt fractions are shown. Basalt melt fractions will vary from their intruded melt fraction of unity to the minimum melt fraction. The predicted wt % SiO2 is given for the amphibolite and basalt based on the experiments of Wolf & Wyllie (1994) and Muntener et al. (2001), respectively, and from MELTS calculations.

Fig. 15.

Melt fraction of (a) crustal and (b) basaltic material as a function of depth and flux of basalt. The maximum height of ascent is 100 m, dike width 1 m and the information was extracted 106 years following the initiation of the modeled intrusions. For the amphibolite composition both the average and maximum melt fractions are shown. Melt fractions of the amphibolite in the simulations vary from zero to the maximum melt fraction. The mean and minimum basalt melt fractions are shown. Basalt melt fractions will vary from their intruded melt fraction of unity to the minimum melt fraction. The predicted wt % SiO2 is given for the amphibolite and basalt based on the experiments of Wolf & Wyllie (1994) and Muntener et al. (2001), respectively, and from MELTS calculations.

Similar results are produced after 107 years for 10 m wide dikes and a maximum ascent height of 1000 m (Fig. 16). Crustal melt fractions vary from zero to the maximum crustal melt fraction in the upper left columns of Figs 15 and 16. Likewise, basalt melt fractions are bounded by their intruded melt fraction, 1·0, and their minimum melt fraction shown in the lower left of Figs 15 and 16.

Fig. 16.

Melt fraction of (a) crustal and (b) basaltic material as a function of depth and flux of basalt. The maximum height of ascent is 1000 m, dike width is 10 m and the information was extracted 107 years following the initiation of intrusion. (For further explanation, see Fig. 15.)

Fig. 16.

Melt fraction of (a) crustal and (b) basaltic material as a function of depth and flux of basalt. The maximum height of ascent is 1000 m, dike width is 10 m and the information was extracted 107 years following the initiation of intrusion. (For further explanation, see Fig. 15.)

The simulations predict that for a thin crust (∼30 km) of mafic composition, and the specified initial condition of 68 mW/m2 surface heat flow, very little to no crustal melting can be expected for geologically constrained basalt fluxes. At 106 years the intrusive zones of 100 m and 1000 m fail to produce mean crustal melt fractions above 0·1 even at the highest basalt flux. Except for the very lowest basalt fluxes, a small amount of residual melt from the basalt will remain after reaching ambient conditions. However, the basalt ensemble average melt fractions do not exceed 0·2 after 106 years since the initiation of intrusions for shallow depths. Residual basalt magma compositions in these melt fraction ranges are generally dacitic to rhyo-dacitic, and the crustal melts are dacitic (tonalitic) with ∼65 wt % SiO2.

Melting is facilitated with increasing depth, and with a 40 km thick crust, crustal melt fractions exceed 0·2 for both the 100 and 1000 m maximum heights of ascent after 106 years of intrusion. After 107 years the crustal melt fractions reach ∼0·3. Coeval basalt melt fractions have approximately the same range as the crustal melt fractions (0·2–0·3). At 50 km depth the mean crustal melt fractions are slightly lower than the residual basalt melt fractions. Average crustal melt fractions reach ∼0·35 at these depths after 107 years of intrusion for a maximum ascent height of 1000 m. The coeval maximum basalt melt fraction is ∼0·38. After 107 years of intrusion at depths >40 km and basalt fluxes exceeding 0·001 m3/m2 per year, both the residual basalt and crustal melts are andesite to basaltic andesite in composition and can persist as long as basalt flux continues into the lower crust system.

Within the stochastic framework, the efficiency of crustal melting is slightly enhanced relative to the reported one-dimensional over-accretion values (4–8%) for greater depths (>40 km), low basalt flux (<0·0005 m3/m2 per year), and immediately following initiation of intrusion. For instance, after 1 Myr and at 44 km depth, 10 m wide dikes are ∼10% efficient at producing crustal melts (Table 4). However, because the active zone of intrusion remains fixed at the base of the crust, the efficiency of crustal melting decreases with time as more enthalpy is consumed, elevating the temperature of recently intruded basalt. After 10 Myr, 10 m dikes are only 5% efficient at producing crustal melts at 44 km depth. For the assumed initial geotherm, the injection of basalt in dikes is less efficient at all times relative to over-accretion for thin crust.

### The geotherm and surface heat flux

Although the initial steady-state geotherm was calculated with a surface heat flux of 68 mW/m2, the surface heat flux with the progressive intrusion of basalt increased with time (Fig. 17). After 5·0 × 107 years of simulated time, and a basalt flux of 0·001 m3/m2 per year, the surface heat flux from intrusions at 30 km depth exceeds 100 mW/m2. For similar conditions, but with intrusions at 50 km depth, the surface heat flux is ∼80 mW/m2 after 5·0 × 107 years of intrusions. As expected, the simulations indicate that the thicker crust produces a greater amount of crustal melting, but has a lower surface heat flux than the thinner crust with the same flux of basalt.

Fig. 17.

Average surface heat flux as a function of time. Three simulated conditions are considered: 0·001 m3/m2 per year basalt flux with 30 km, 42 km and 50 km thick crust. Basalt dikes are 10 m thick and the maximum height of ascent is 1000 m. After 5·0 × 107 years surface heat flux exceeds 100 mW/m2 for the 30 km thick crust, whereas the 50 km crust has a surface heat flux of ∼80 mW/m2 for the same period of intrusion.

Fig. 17.

Average surface heat flux as a function of time. Three simulated conditions are considered: 0·001 m3/m2 per year basalt flux with 30 km, 42 km and 50 km thick crust. Basalt dikes are 10 m thick and the maximum height of ascent is 1000 m. After 5·0 × 107 years surface heat flux exceeds 100 mW/m2 for the 30 km thick crust, whereas the 50 km crust has a surface heat flux of ∼80 mW/m2 for the same period of intrusion.

### Advection simulations: mixing and mingling in the lower crust

The intrusion of basaltic magma into the crust creates a hybrid thermal and compositional framework that may facilitate mixing of mantle-derived basalts with crustal melts. Additionally, subsolidus creep in the ductile lower crust may mingle mantle and crustal material, and may give rise to density instabilities where portions of the dense lower crustal material are transferred into the mantle (Jull & Kelemen, 2001; Lee et al., 2001; Ducea, 2002). Both hyper-solidus mixing and sub-solidus mingling are potentially important in determining the mass balance and chemical nature of the crust, although they operate on very different time scales. Two sets of simulations were performed to illuminate (1) magma mixing under different thermal conditions, and (2) subsolidus creep and mingling of basalt and crust, as well as the conditions required for a density instability and crustal delamination.

### Magma mixing

The homogenization of mantle-derived basalts with partially molten crust was modeled in a suite of simulations that used the conduction results for initial and boundary conditions. The mixing calculations simulated 5000 years of elapsed time, and the same intrusion configuration was used in all simulations with an average basalt flux of 0·001 m3/m2 per year. The initial conditions for the advection simulations were extracted from the conduction simulations for depths of 34 km, 38 km and 42 km after 107 years of simulated time (shown as stars in Fig. 11). Although this is not an exhaustive assessment of the possible outcomes, it does illustrate the potential mixing conditions for a variety of melt fractions. In these simulations, the mixture density of the magma and solid was used.

In keeping with the rheological model [Appendix equations (A4) and (A5)], an abrupt transition in the ability to mix magmas was observed at conditions in which the mean crustal melt fraction was >0·3 and maximum melt fractions of both the intruded and crustal material exceed the critical melt fraction (CMF) in close juxtaposition. For the three conditions considered, only the 42 km depth simulations produced conditions conducive to rapid mixing (Fig. 18). However, it is reasonable to extrapolate the ease of mixing to greater depths and greater flux of basalt as these conditions would probably also produce high melt fractions. Maximum crustal melt fractions exceed 0·4 in the 38 km depth simulations; however, adjacent regions of basalt are generally quenched, impeding mixing. In addition, these regions of melt above the CMF are also short-lived, typically decaying within some 100 days with length scales of 1–5 m. For both the 34 and 38 km depth simulations, isolated pods of basalt melt are separated by crustal material with melt fractions <0·2. This impedes large-scale commingling of material on these time scales.

Fig. 18.

Mixing and mingling of crustal and basaltic melts. Initial conditions for these 5000 year simulations are given in Fig. 11 (stars) and are for a basalt flux of 0·001 m3/m2 per year, 107 years following the initiation of intrusion. Depths of 34 km, 38 km and 42 km were modeled. After 5000 years only the 42 km simulation produced extensive mixing or mingling. Grey scale indicates the extent of mixing or mingling.

Fig. 18.

Mixing and mingling of crustal and basaltic melts. Initial conditions for these 5000 year simulations are given in Fig. 11 (stars) and are for a basalt flux of 0·001 m3/m2 per year, 107 years following the initiation of intrusion. Depths of 34 km, 38 km and 42 km were modeled. After 5000 years only the 42 km simulation produced extensive mixing or mingling. Grey scale indicates the extent of mixing or mingling.

In the 42 km depth simulations, mixing is driven by the crystallization of the basalt, which yields a crystal–melt mixture density greater than the partially molten surrounding crust; consequently, the basalt mixture will begin to ‘sink’ back to a level of neutral buoyancy. This flow may trap isolated patches of crustal melt, and initiate low Reynolds number mixing. Likewise, partial melting will generate tonalitic melts from the amphibolite that will be buoyant. This is exemplified when a pool of crustal melt accumulates below a dike. Provided that portions of the dike have melt fraction >0·4, the crustal melt can ‘tunnel’ into the overlying dike, and mix in the interior. The Reynolds number in the interior of the dikes ranges from 1·0 × 10−4 to 1·0. Upon reaching the interior of the dike, mixing occurs on time scales of 1–100 days, providing a mechanism for relatively rapid mixing and mingling in the vicinity of the intruded dikes. In these simulations, the buoyancy reversals that provided the potential energy driving mixing operated mostly on length scales of <10 m, which suggests that for dike swarm intrusions in the lower crust the MASH process is fast and local in the thermally mature portions of a thickened crust.

As the melt fractions of both the crustal melts and the residual mantle melts are often similar for any given thermal condition, their compositions and viscosity will also be similar, and magma mixing may be expected provided the magmas are closely juxtaposed. The simulations predict that mixed magmas will be volumetrically weighted toward the mantle contribution in shallow systems and in arc systems with a high flux of basalt (Fig. 12). A low flux of basalt, and especially thicker crust (∼50 km) produces the greatest ratio of crustal melt to mantle melt. Although the ratio of crustal melt in erupted magmas would be difficult to discern on the basis of major element geochemistry alone, the isotopic signature of the arc melts will be affected, provided there is sufficient isotopic contrast between the crustal and mantle magmas (Hart et al., 2002).

### Gravitational instability and ductile creep

Density differences between the crystallizing basalt, melting amphibolite, and mantle peridotite may also produce subsolidus creep. Jull & Kelemen (2001) demonstrated that ultramafic cumulates with densities greater than typical mantle lithologies can form and, at temperatures greater than 700°C, produce crustal density instabilities on time scales of 107 years. We expand on their model by invoking two end-member density models: for one end-member, crustal and residual basaltic melt is assumed to be completely removed from the solid matrix leaving a dense cumulate, and, at the other extreme, no melt is removed and the mixture density is calculated for the magma and solid.

Only in the simulations in which melt was removed did the lower crustal cumulates (both residual material from the intruded magma and crustal solids) become more dense than mantle peridotite (Appendix Fig. A4). When the mean crustal and intruded melt fraction is low (<0·2 melt fraction), the removal of melt was not sufficient to produce cumulates that had greater density than the mantle. The lower density of the plagioclase component, which remains a residual phase for small melt fractions, lowers the total cumulate density below the density of peridotite.

Density contrasts created by the web of basaltic dikes initiate the instabilities for melt fractions greater than ∼0·2. These instabilities coalesce as they descend into the mantle (Fig. 19). The constitutive relation used in these simulations is similar to the weak, dense layer considered by Jull & Kelemen (2001); however, the random orientation of the dikes, and the progressive addition of mass, complicates using the scaling formulated by those workers for the Rayleigh–Taylor instability of the dense layer. Our results for the initiation of density instabilities are consistent with the results of Jull & Kelemen (2001) and the conceptual model of Ducea (2002), and add additional constraints concerning the formation of instabilities in a crust where melting of amphibolite and/or the crystallization of wet basalt is occurring. The process requires: (1) the stability of garnet, which necessitates pressures >9 kbar and average melt fractions <0·5; (2) the removal of most of the plagioclase component, which requires melt fractions typically greater than ∼0·2; (3) efficient segregation of melt leaving dense cumulates. To illustrate the basalt intrusion conditions that can produce a density instability, the average cumulate density from the simulations given in Fig. 16 (10 m dikes, 1000 m maximum ascent, 107 years after the initiation of intrusion) is shown in Fig. 20. Densities that exceed ∼3300·0 kg/m3 are sufficient to initiate an instability. If melt is efficiently segregated, a lower flux (0·0005 m3/m2 per year), and a crustal thickness of 50 km are required to produce crustal density greater than the mantle. However, an order of magnitude greater flux can melt most of the plagioclase component at shallower depths, and if the melt is removed this crust may also become unstable at depths just below the garnet stability field. The trace element concentration of melts in this range should distinctively show the presence of residual garnet with little residual plagioclase. For instance, elevated La/Yb and Sr/Y ratios would be expected (Tulloch & Kimbrough, 2003; Bachmann et al., 2005). Furthermore, seismic velocity variations in the upper mantle for such a process may be discerned, such as the inferred ‘density drip’ beneath the Sierra Nevada (Saleeby & Foster, 2004).

Fig. 19.

Temporal evolution of the cumulate density (all melt has been extracted) and ductile creep at the mantle–crust interface. The depth scale is referenced to the Moho (shown as a white dashed line). The cumulate density calculations are for a flux of 0·0002 m3/m2 per year, 10 m wide dike and 50 km crustal depth, after 1 × 107 and 2 × 107 years of intrusion. The cumulate density exceeds the mantle peridotite density, creating convective instabilities. Instabilities on the scale of individual dike intrusions coalesce as they descend. (Resolution 10 m2).

Fig. 19.

Temporal evolution of the cumulate density (all melt has been extracted) and ductile creep at the mantle–crust interface. The depth scale is referenced to the Moho (shown as a white dashed line). The cumulate density calculations are for a flux of 0·0002 m3/m2 per year, 10 m wide dike and 50 km crustal depth, after 1 × 107 and 2 × 107 years of intrusion. The cumulate density exceeds the mantle peridotite density, creating convective instabilities. Instabilities on the scale of individual dike intrusions coalesce as they descend. (Resolution 10 m2).

Fig. 20.

Average basalt cumulate density as a function of crustal depth and flux of basalt. (Simulation conditions are the same as in Fig. 16: 10 m dikes, 1000 m maximum ascent, 107 years after the initiation of intrusion). The cumulate density becomes greater than that of mantle peridotite (∼3300 kg/m3) if the melt fraction is sufficiently high that little plagioclase is in the solid assemblage (>0·2), pressures are great enough that garnet is stable (>9 kbar), the melt fraction is below ∼0·5 so that garnet is a residual phase, and the melt is efficiently segregated from the crystalline residuum.

Fig. 20.

Average basalt cumulate density as a function of crustal depth and flux of basalt. (Simulation conditions are the same as in Fig. 16: 10 m dikes, 1000 m maximum ascent, 107 years after the initiation of intrusion). The cumulate density becomes greater than that of mantle peridotite (∼3300 kg/m3) if the melt fraction is sufficiently high that little plagioclase is in the solid assemblage (>0·2), pressures are great enough that garnet is stable (>9 kbar), the melt fraction is below ∼0·5 so that garnet is a residual phase, and the melt is efficiently segregated from the crystalline residuum.

## DISCUSSION

### Predictions of melt compositions based on thermal modeling

The genetic relationship between tonalitic–andesitic melt and amphibolite–pyroxenite residuum predicted by dehydration experiments at lower crustal pressures is expressed in the geological record. Observations of tonalite, amphibolite, and pyroxenite in close juxtaposition have been noted in several lower crustal terrains. The amphibolite Chipman dikes in Saskatchewan have dominant hornblende and plagioclase and minor clinopyroxene and garnet in dikes that have intruded a tonalitic body. The emplacement of these dikes has been interpreted to be at a depth equivalent to ∼10 kbar at the base of an Archean island arc (Williams et al., 1995). Reheating conditions have produced pods of tonalitic melt, commonly associated with strain shadows associated with garnet. The Kohistan arc, Pakistan, interpreted as the exposed base of a Jurassic island arc at 12–14 kbar, has pyroxenitic, amphibolitic and tonalitic assemblages, although multiple stages of deformation have hindered the reconstruction of the relationship between units (Jan & Howie, 1981). Similarly, pyroxenites are common in the 9·5–11 kbar Tonsina assemblage, Alaska, in yet another basal zone of an interpreted mid-Jurassic island arc (DeBari & Coleman, 1989).

The geochemical and modal mineralogy trends of distinct melting regimes established by the depth, basalt flux, and duration of intrusion in the melting zone can be determined by combining an assumed starting lithology, compositional information from melting experiments, and the melt fractions predicted in the conduction simulations. We focus on the silica and potassium contents, and the aluminum saturation index (ASI), to compare the experiments and simulations with observed arc magmas.

For most basalt flux conditions (Figs 11, 15 and 16) crustal melt fractions and residual basalt melt fractions are <0·1 for a 30 km thick crust subjected to random diking, even after several million years of intrusion. At small crustal melt fractions (<0·1), dehydration of amphibolite, as considered by Wolf & Wyllie (1994), produces a tonalitic melt with ∼65 wt % SiO2 at 10 kbar. The maximum SiO2 content in a partial melt of amphibolite is dependent on the mode of the protolith. The Wolf & Wyllie (1994) experiments represent a primitive, high-Mg number end-member with mostly amphibole and high-anorthite plagioclase in the starting composition. Experiments with greater amounts of plagioclase, more albitic plagioclase, and with small amounts of quartz will all increase the SiO2 wt % of the melt. For instance, the quartz amphibolite composition examined by Patiño-Douce & Beard (1995) produced a melt with ∼76 wt % SiO2 at 10 kbar and 0·1 melt fraction. Regardless, tonalitic melts (essentially low-K2O dacitic to rhyolitic melts) are predicted to be produced at the low melt fractions predicted to form in the modeled thin crust environment. Even crust of intermediate thickness (30–40 km) will have crustal melts that are dacitic after tens of millions of years of intrusion. The SiO2 content of the amphibolite dehydration melt stays at or above 60 wt % until ∼0·3 melt fraction, or at approximately the amphibole-out phase boundary. Only after reaching ∼0·4 melt fraction are andesitic to basaltic–andesitic crustal melts produced. For the thin crust (∼30 km), more typical of island arcs, such high melt fractions were not achieved in the simulations even after 50 Myr of intrusion with a basalt flux of 0·001 m3/m2 per year. However, crustal melts up to 0·4 were produced in a 50 km thick crust in a few million years for a range of basalt fluxes, and basaltic–andesitic crustal melt was predicted to have a long residence time provided the crust continued to be fluxed with mantle-derived magmas. It should be stressed that this trend applies to the baseline compositions owing to lower crustal magma processing, and subsequent fractionation in the upper crust may be responsible for the production of some of the silicic magmas that are observed in thick crustal settings. This may be particularly true of magmas that have trace element patterns indicative of extensive plagioclase fractionation, which is largely inhibited at higher pressures.

The SiO2 trend of a crystallizing wet basalt is essentially the reverse of the melting experiments. For the 12 kbar, 3·8 wt % H2O experiments of Muntener et al. (2001), the crystallization of pyroxene dominates the phase assemblage to ∼0·4 melt fraction, and does not alter the SiO2 content significantly (Fig. 4). Extrapolation with MELTS predicts ∼65 wt % SiO2 for the crystallizing magma at 0·2 melt fraction (Appendix Fig. A2). Again, the thermal conditions in a thin crust or with a very low basalt flux will produce dacitic to rhyo-dacitic magma if these mantle basalts stall for a protracted period of time, whereas the thermal conditions produced within a thicker crust or with a greater basalt flux will drive the residual magmas toward an andesitic composition.

At low melt fractions (0·1–0·3) the Wolf & Wyllie (1994) amphibolite dehydration experiments have potassium concentrations of 0·3–0·4 wt %. Extrapolation of the experiments of Muntener et al. (2001), used to parameterize the crystallization of the intruded wet basalt, predicts potassium concentrations from 0·8 to 1·1 wt % in the same melt fraction range. At higher melt fractions, the more mafic intruded magma has an even lower potassium content (Appendix, Fig. A3). The low potassium content (<1·0 wt %) produced in several amphibolite melting experiments (Beard & Lofgren, 1991; Rapp et al., 1991; Rushmer, 1991; Wolf & Wyllie, 1994) is typical of many island-arc magmas (Bryan et al., 1979; Zellmer et al., 2003) and some tonalite–trondhjemite–granodiorite suites (Arth et al., 1978), but is too low when compared with many continental arc dacites and rhyolites (Anderson et al., 2000; Bachmann et al., 2002), suggesting that a prolonged multi-stage process involving fractionation is required to produce continental arc compositions rather than a singular melting event. The low potassium produced in these melting experiments is the result of low initial potassium composition because potassium is incompatible in the residual phases produced in the dehydration reaction (Rapp & Watson, 1995; Sisson et al., 2005). In the few experiments (Sen & Dunn, 1994; Sisson et al., 2005) with initial potassium at higher concentrations, potassium contents more typical of calc-alkaline granites are produced. Dehydration melting experiments on an amphibolite performed by Sen & Dunn (1994) at 15 kbar produced andesite–dacite composition melts with ∼2·5–5·0 wt % K2O. The experiments of Sisson et al. (2005) have demonstrated that small amounts of partial melting (0·1–0·3 melt fraction) of a mafic source at 8 kbar can produce similarly high potassium content in dacitic to rhyolitic liquids if the initial potassium content of the protolith is ∼1–2 wt % K2O. These results indicate that if enough potassium can be introduced into the basalt as it leaves the mantle wedge, or accumulates in this melt while it stalls at the base of the crust, and then crystallizes to form an amphibolite, subsequent partial melting can directly produce the elevated potassium content observed in continental arc settings.

Mantle melting processes, recycling of weathered exogenous sediments (McCarthy & Patiño-Douce, 1997), and repeated melting–solidification cycles (igneous distillation), have all been invoked to produce the high potassium content in lower crustal settings. The high potassium in melts in settings such as Costa Rica (Lidiak & Jolly, 1996; Hannah et al., 2002) and the Philippines (T. Vogel, personal communication, 2004) where a thick crust of meta-sediments is absent, appears to imply that in some settings incorporation of weathered material is not necessary for the production of high-K suites. Igneous distillation can produce elevated incompatible element concentrations; however, the accumulation of cumulates may create a mass balance problem. The model calculations demonstrate that crustal instability is a likely consequence of basalt intrusion and melt extraction, especially in thick crustal regions. A prolonged intrusion history with multiple melt–solidification cycles, coupled with the development of lower crustal instabilities that remove some of the mafic cumulates, may produce a lower crust that is a repository for the incompatible-element enriched melt, but still produces the crustal thickness observed seismically. Such an environment is predicted in these simulations for a relatively steady flux of basalt intruding the crust over millions of years in several stages of intrusion. For instance, after several million years of simulated basalt intrusion, a variety of basalt fluxes are capable of generating multiple melt–solidification cycles as well as generating density instabilities (Fig. 20, densities greater than 3300 kg/m3), although at different depth–time combinations. After removal of the garnet- and pyroxene-rich residuum, the remaining lithologies will be tonalitic–andesitic–dacitic, and subsequent melting processes will further increase the incompatible element concentration, including their potassium content.

Low melt fraction dehydration melting of amphibolite, as predicted by the thermal modeling, typically produces moderate to mildly peraluminious melts (Fig. 21). As the melt fraction increases the aluminum content decreases, and the melts become meta-aluminous near the amphibole-out boundary (Rapp et al., 1991). In the simulations these meta-aluminous conditions are achieved in thicker, and more mature crust. The ASI throughout the melting experiment of Wolf & Wyllie (1994) are at or slightly more peraluminous than the average range observed in calc-alkaline andesites, dacites and rhyolites (Ewart, 1982).

Fig. 21.

Aluminum saturation index (ASI) vs wt % SiO2 for partial melts of dehydrating amphibolite and fractionating wet basalt. Also shown is the ASI range of calc-alkaline magmas reported in the database of Ewart (1982) (shaded area). The model curves have melt fraction indicated (black for amphibolite dehydration and grey for the wet basalt).

Fig. 21.

Aluminum saturation index (ASI) vs wt % SiO2 for partial melts of dehydrating amphibolite and fractionating wet basalt. Also shown is the ASI range of calc-alkaline magmas reported in the database of Ewart (1982) (shaded area). The model curves have melt fraction indicated (black for amphibolite dehydration and grey for the wet basalt).

The shallow and young arc systems investigated in the conduction simulations produced low melt fraction crustal melts and fractionated mantle melts, both of which are predicted to be dacitic–rhyodacitic and mildly peraluminous. Residual mineral modes for amphibolite melt fractions <0·2 are amphibole > plagioclase > clinopyroxene ± garnet and orthopyroxene. Similar trends are predicted for the cumulates of crystallizing wet basalt at low melt fraction. Regions of greater crustal and intruded magma melt fraction, such as at the base of thickened crust, will produce metaluminous magmas that are generally andesitic. Residual modes will be dominated by pyroxene ± garnet depending on the pressure. In the conduction simulations, the continued input of enthalpy and mass from the basalt intrusions will create a lower crust with greater melt fractions and melt that becomes more mafic with time.

## CONCLUSIONS

The stochastic simulation of basalt diking demonstrates that partial melting in the lower crust is maximized when basaltic intrusions are confined to a relatively narrow vertical zone at the base of crust; basaltic dikes that ascend several kilometers from the base of the crust are inefficient at producing voluminous crustal melting. This suggests that if lower crustal melting is widespread, density or rheological barriers inhibiting the propagation of basaltic melt are required.

The conduction simulations indicate that shallow, young crust may be a repository for residual basalt melt fractions from ∼0·1 to 0·3 and amphibolite dehydration melts from zero to 0·1 melt fraction. The coupled cooling of intruded basalt and melting of crust can lead to a convergence in melt compositions in the dacitic to rhyolitic range. Except for very young arc systems, or very low basalt flux, the amount of residual melt from the basalt will be more voluminous than the amount of crustal melt, and linear mixtures of the two will be dominated by the crystallizing basalt. However, the dynamic simulations suggest that the low melt fractions present may inhibit the rapid mixing of these melts and that distinct pockets of isolated melt may form. At these low melt fractions, the residuum includes amphibole and plagioclase, and even with efficient segregation of melt, the density of the lower crustal material will be positively buoyant relative to mantle peridotite. Although little crustal melting occurs, surface heat flow may reach elevated values because of the thin crust.

Mature arc systems with thicker, and ultimately hotter, crust at depth yield greater melt fractions from both crustal melting and residual material from the intruded basalt at the base of the crust. For crustal depths of 50 km, melt fractions exceeding 0·4 for both the residual basalt liquids and crustal melts can coexist for thousands of years. Individually, these melts will be andesitic to basaltic andesite in composition. The long period of time that these liquids are at high melt fractions promotes homogenization, as demonstrated by the advection simulations. Coexisting with the andesitic meta-aluminous melts will be a garnet pyroxenitic residuum. If the melt can be extracted efficiently from this residuum, the density of this material will probably be great enough to initiate convective instabilities, providing a mechanism to remove the mafic material from the base of thickened crust.

## APPENDIX

### Mixture properties

Both phase changes and advection can result in mixtures of melt and solid, and magmas of different composition, residing in a control volume. A composition parameter (γ) describes the volume fraction of crustal material or intruded basalt residing at a particular location. If γ equals one, the region is all crustal material, and if it equals zero it is all mantle-derived basalt. Intermediate values indicate relative proportions, and hence locally mixed material. The variable f denotes the local melt fraction. In the following nomenclature, the subscripts b and c will be used to refer to properties of the basalt and crust, respectively. Superscripts s and l refer to solid and melt. The mixture properties are defined as follows.

Mixture density:

(A1)
${\rho}_{\mathrm{mix}}\ =\ {\gamma}f_{\mathrm{c}}\ {\rho}_{\mathrm{c}}^{\mathrm{l}}\ +\ {\gamma}\left(1\ {-}\ f_{\mathrm{c}}\right){\rho}_{\mathrm{c}}^{\mathrm{s}}\ +\ \left(1\ {-}\ {\gamma}\right)f_{\mathrm{b}}\ {\rho}_{\mathrm{b}}^{\mathrm{l}}\ +\ \left(1\ {-}\ {\gamma}\right)\left(1\ {-}\ f_{\mathrm{b}}\right){\rho}_{\mathrm{b}}^{\mathrm{s}}.$
Mixture heat capacity:
(A2)
$c_{\mathrm{mix}}\ =\ {\gamma}f_{\mathrm{c}}\ c_{\mathrm{c}}^{\mathrm{l}}\ +\ {\gamma}\left(1\ {-}\ f_{\mathrm{c}}\right)c_{\mathrm{c}}^{\mathrm{s}}\ +\ \left(1\ {-}\ {\gamma}\right)f_{\mathrm{b}}\ c_{\mathrm{b}}^{\mathrm{l}}\ +\ \left(1\ {-}\ {\gamma}\right)\left(1\ {-}\ f_{\mathrm{b}}\ \right)c_{\mathrm{b}}^{\mathrm{s}}.$
Mixture conductivity:
(A3)
$k_{\mathrm{mix}}\ =\ {\gamma}f_{\mathrm{c}}\ k_{\mathrm{c}}^{\mathrm{l}}\ +\ {\gamma}\left(1\ {-}\ f_{\mathrm{c}}\right)k_{\mathrm{c}}^{\mathrm{s}}\ +\ \left(1\ {-}\ {\gamma}\right)f_{\mathrm{b}}\ k_{\mathrm{b}}^{\mathrm{l}}\ +\ \left(1\ {-}\ {\gamma}\right)\left(1\ {-}\ f_{\mathrm{b}}\right)k_{\mathrm{b}}^{\mathrm{s}}.$
The rheology of solid–melt mixtures has been observed to deviate significantly from a simple linear mixing relationship, as depicted in Fig. A1 (Vigneresse et al., 1996; Barboza & Bergantz, 1998; Vigneresse & Tikoff, 1999; Renner et al., 2000). Solid volume fractions that form an interlocking crystalline framework may greatly increase the mixture viscosity (Marsh, 1981). The critical melt fraction (CMF) at this transition is likely to be dependent on the strain rate of the material (Renner et al., 2000) and may vary over a wide range of values (Barboza & Bergantz, 1998). A rapid drop in viscosity has not been observed in experiments with melt fractions reaching ∼0·2 (Rushmer, 1995). We assume a critical melt fraction of 0·4, to link the bulk rheology of a melt-supported crystal suspension to an interlocking crystal framework (Marsh, 1981). This is equivalent to the CMF suggested by Scaillet et al. (1998) for a monodispersed crystal size distribution.

Fig. A1.

Mixture viscosity as a function of melt fraction. The critical melt fraction (CMF) marks the transition from an interlocking crystalline framework and crystal–melt suspension. The CMF is assumed to equal 0·4 for a monodispersed crystal size distribution (Scaillet et al., 1998).

Fig. A1.

Mixture viscosity as a function of melt fraction. The critical melt fraction (CMF) marks the transition from an interlocking crystalline framework and crystal–melt suspension. The CMF is assumed to equal 0·4 for a monodispersed crystal size distribution (Scaillet et al., 1998).

When the melt fraction is below the critical melt fraction, an empirically derived viscosity was used, motivated by the amphibolite deformation experiments of Rushmer (1995). The ductile rheology used here assumes a linear relationship between the applied stress and strain rate, and is a macroscopic manifestation of a process that has been shown to be a combination of brittle and plastic flow (Rushmer, 1995). This functional form is equivalent to an assumption of diffusion creep at low melt fractions. For the subset of calculations that include the mantle interface, the diffusion creep of olivine was used as a proxy for mantle rheology (Hirth & Kohlstedt, 1995, 1996; Jull & Kelemen, 2001). Above the critical melt fraction, the mixture viscosity is modified from the melt viscosity to account for the macroscopic effect of increased drag between crystals and melt (Marsh, 1981; Dobran, 1992; Pinkerton & Stevenson, 1992; Scaillet et al., 1998). The mixture viscosity for melt fraction greater than the CMF becomes

(A4)
${\mu}_{\mathrm{mix}}\ =\ {\mu}_{\mathrm{m}}^{\mathrm{i}}\ \left[1\ +\ 0{\cdot}75\left(\frac{\left(1\ {-}\ f\right)/\left(1\ {-}\ \mathrm{CMF}\right)}{1\ {-}\ \left(1\ {-}\ f\right)/\left(1\ {-}\ \mathrm{CMF}\right)}\right)\right]^{2}.$
For melt fractions below the CMF, the mixture rheology for the crust and mantle is given by
(A5)
${\mu}_{\mathrm{mix}}\ =\ {\mu}_{\mathrm{param}.}$
and
(A6)
${\mu}_{\mathrm{mantle}}\ =\ \frac{A^{{-}1/n}\ \mathrm{exp}\left(Q/nRT\right)}{2}$
respectively. Here the parameters for the mantle rheology are
(A7)
$\mathrm{ln}\ A\ =\ 15{\cdot}4\ \mathrm{MPa}^{{-}n}/\mathrm{s},{\ }Q\ =\ 515\ \mathrm{kJ}/\mathrm{mol},\ n\ =\ 3{\cdot}5$
(Hirth & Kohlstedt, 1996; Jull & Kelemen, 2001).

### Numerical method

The conservation equations were discretized using a finite volume numerical method (Patankar, 1980). The iterative procedure SIMPLER (Semi-Implicit Method for Pressure Linked Equations, Revised) was used to ensure that the pressure and velocity simultaneously satisfy both momentum and continuity (Patankar & Spalding, 1972; Patankar, 1980; Versteeg & Malalasekera, 1995). The power-law algorithm of Patankar (1980) was used in dynamic simulations and uses a hybrid of central differencing and upwinding numerical schemes depending on the local Peclet number. Additionally, a predictor–corrector algorithm was used to compute the partitioning of energy into sensible and latent heat during phase change processes. Iteration and numerical under-relaxation was required for convergence given the non-linear melt fraction vs temperature relationships. The predictor–corrector phase change algorithm of Voller & Swaminathan (1991) has shown the best convergence and was used for the calculations presented.

### Physical and thermal properties of solid and melt

The composition of the magma and its volatile content determine the physical properties of the magma such as density and viscosity that are required to close equations (2)(5). The composition of the melts are projected from the amphibolite dehydration experiments of Wolf & Wyllie (1994) and the basaltic andesite crystallization experiments of Muntener et al. (2001) (Figs A2 and A3). The experimentally determined compositions were parameterized with second-order polynomials to provide a relationship between the melt fraction and the composition of the magma.

Fig. A2.

Major element geochemical trends for the modeled intruded basalt vs melt fraction (Muntener et al., 2001). Dashed line is an extrapolation of experimental values (•). These experiments were performed on a basaltic andesite at 12 kbar with 3·8 wt % water.

Fig. A2.

Major element geochemical trends for the modeled intruded basalt vs melt fraction (Muntener et al., 2001). Dashed line is an extrapolation of experimental values (•). These experiments were performed on a basaltic andesite at 12 kbar with 3·8 wt % water.

Fig. A3.

Major element geochemical trends for the modeled amphibolitic crust vs melt fraction (Wolf & Wyllie, 1994). These experiments were performed on a high Mg-number amphibolite at 10 kbar.

Fig. A3.

Major element geochemical trends for the modeled amphibolitic crust vs melt fraction (Wolf & Wyllie, 1994). These experiments were performed on a high Mg-number amphibolite at 10 kbar.

The lowest melt fraction reported by Muntener et al. (2001) for the crystallizing basalt was 0·39. Extrapolation of the composition to melt fractions lower than 0·39 was accomplished with MELTS, with the acknowledgement that error is introduced by the uncertainty in the mode of the crystallizing phases outside of the experimental range as a result of the decreased range of amphibole stability in the thermodynamic calculations. The composite PT diagram of a tholeiite with 5 wt % water developed by Green (1982) indicates that the first appearance of plagioclase is at about ∼800°C at 12 kbar, which corresponds to a melt fraction between ∼0·1 and 0·2. The MELTS calculations also reach plagioclase saturation at ∼0·2 melt fraction.

### Transport properties

Magmatic density, heat capacity, and viscosity can be determined using the correlations of Lange & Carmichael (1987), Lange & Navrotsky (1992), and Shaw (1972), respectively. Solid densities and heat capacities were determined using the experimental phase assemblages of Wolf & Wyllie (1994) and Muntener et al. (2001) and the thermodynamic database of Holland & Powell (1998). For the crystallizing basalt, the density was calculated at <0·39 melt fraction using MELTS (Ghiorso & Sack, 1995). These MELTS calculations reached amphibole saturation at 0·18 melt fraction, compared with the first appearance of amphibole at ∼0·39 for the same conditions in the experiments. This discrepancy is a likely consequence of the complexity of amphibole solid solution. The difference in mode appears to be accounted for by clinopyroxene.

Garnet in the experiments and MELTS calculations crystallizes at approximately the same melt fraction and has a mode of 17% in the MELTS calculation and 15% in the experiments at a melt fraction of 0·39. The densities of clinopyroxene and amphibole are very similar, and so although the mode of the thermodynamic calculations did not exactly reproduce the experiments, the density was reasoned to be a good proxy. Other specified parameters used in the simulations are included in Table 3. Mixture, melt and solid density of the residual melt from the intruding basalt and the amphibolite are compared as a function of melt fraction at 12 kbar in Fig. A4 along with the density of a mantle peridotite at 800°C.

Fig. A4.

Density of amphibolite and wet basalt as a function of melt fraction at 12 kbar (continuous lines). Calculations are based on the dataset of Holland & Powell (1998) (crystals) and the algorithm of Lange & Carmichael (1987) (melt). Also shown are the cumulate and melt density (dashed lines), and the mantle density at 800°C.

Fig. A4.

Density of amphibolite and wet basalt as a function of melt fraction at 12 kbar (continuous lines). Calculations are based on the dataset of Holland & Powell (1998) (crystals) and the algorithm of Lange & Carmichael (1987) (melt). Also shown are the cumulate and melt density (dashed lines), and the mantle density at 800°C.

Encouragement and discussions with O. Bachmann, T. Rushmer, T. Vogel, M. Williams, and T. Sisson have greatly improved this manuscript. Thoughtful reviews and suggestions by M. Ducea, M. Williams, A. Patiño-Douce, and Editor D. Geist were very useful in the preparation of this manuscript. This research has been supported by a DoD National Defense Science and Engineering Graduate Fellowship (J.D.) and NSF Grant EAR-0106441 (G.W.B.).

## REFERENCES

Anderson, A. T. (
1976
). Magma mixing—petrological process and volcanological tool.
Journal of Volcanology and Geothermal Research

1
(1),
3
–33.
Anderson, A. T., Davis, A. M. & Lu, F. Q. (
2000
). Evolution of Bishop Tuff rhyolitic magma based on melt and magnetite inclusions and zoned phenocrysts.
Journal of Petrology

41
(3),
449
–473.
Annen, C. & Sparks, R. S. J. (
2002
). Effects of repetitive emplacement of basaltic intrusions on the thermal evolution and melt generation in the crust.
Earth and Planetary Science Letters

203
(3–4),
937
–955.
Arth, J. G., Barker, F., Peterman, Z. E. & Friedman, I. (
1978
). Geochemistry of gabbro–diorite–tonalite–trondhjemite suite of southwest Finland and its implications for origin of tonalitic and trondhjemitic magmas.
Journal of Petrology

19
(2),
289
–316.
Bachmann, O., Dungan, M. A. & Lipmann, P. W. (
2002
). The Fish Canyon magma body, San Juan volcanic field, Colorado: rejuvenation and eruption of an upper-crustal batholith.
Journal of Petrology

43
(8),
1469
–1503.
Bachmann, O., Dungan, M. A. & Bussy, F. (
2005
). Insights into shallow magmatic processes in large silicic magma bodies: the trace element record in the Fish Canyon magma, Colorado. Contributions to Mineralogy and Petrology (in press).
Barboza, S. A. & Bergantz, G. W. (
1996
). Dynamic model of dehydration melting motivated by a natural analogue: applications to the Ivrea–Verbano zone, northern Italy.
Transactions of the Royal Society of Edinburgh

87
,
23
–31.
Barboza, S. A. & Bergantz, G. W. (
1997
). Melt productivity and rheology: complementary influences on the progress of melting.
Numerical Heat Transfer Part A—Applications

31A
,
375
–392.
Barboza, S. A. & Bergantz, G. W. (
1998
). Rheological transitions and the progress of melting of crustal rocks.
Earth and Planetary Science Letters

158
,
19
–29.
Barboza, S. A. & Bergantz, G. W. (
2000
). Metamorphism and anatexis in the mafic complex contact aureole, Ivrea Zone, Northern Italy.
Journal of Petrology

41
(8),
1307
–1327.
Barboza, S. A., Bergantz, G. W. & Brown, M. (
1999
). Regional granulite facies metamorphism in the Ivrea zone: is the Mafic Complex the smoking gun or a red herring?
Geology

27
,
447
–450.
Beard, J. S. & Lofgren, G. E. (
1991
). Dehydration melting and water-saturated melting of basaltic and andesitic greenstones and amphibolites at 1, 3, and 6·9 kb.
Journal of Petrology

32
,
365
–401.
Bergantz, G. W. (
1989
). Underplating and partial melting: implications for melt generation and extraction.
Science

245
,
1093
–1095.
Bergantz, G. W. (
1990
). Melt fraction diagrams: the link between chemical and transport models. In: Nicholls, J. & Russell, J. K. (eds) Modern Methods of Igneous Petrology: Understanding Magmatic Processes. Mineralogical Society of America, Reviews in Mineralogy 24, 240–257.
Bergantz, G. W. (
1995
). Changing paradigms and techniques for the evaluation of magmatic processes.
Journal of Geophysical Research

100
(B9),
17603
–17613.
Bird, P. (
1979
). Continental delamination and the Colorado Plateau.
Journal of Geophysical Research

84
(B13),
7561
–7571.
Bittner, D. & Schmeling, H. (
1995
). Numerical modelling of melting processes and induced diapirisms in the lower crust.
Geophysical Journal International

123
(1),
59
–70.
Bryan, W. B., Thompson, G. & Michael, P. J. (
1979
). Compositional variation in a steady-state zoned magma chamber: mid-Atlantic ridge at 36°50′N.
Tectonophysics

55
,
63
–85.
Chapman, D. S. & Furlong, K. P. (
1992
). Thermal state of the continental lower crust. In: Fountain, D. M., Arculus, R. & Kay, R. W. (eds) Continental Lower Crust. Amsterdam: Elsevier, pp. 179–198.
Cole, J. W. (
1982
). Tonga–Kermadec–New Zealand. In: Thorpe, R. S. (ed.) Andesites; Orogenic Andesites and Related Rocks. Chichester: John Wiley, pp. 245–258.
Crisp, J. A. (
1984
). Rates of magma emplacement and volcanic output.
Journal of Volcanology and Geothermal Research

20
,
177
–211.
DeBari, S. M. & Coleman, R. G. (
1989
). Examination of the deep levels of an island-arc: evidence from the Tonsina ultramafic–mafic assemblage, Tonsina, Alaska.
Journal of Geophysical Research

94
(4),
4373
–4391.
DeBari, S., Mahlburg Kay, S. & Kay, R. W. (
1987
). Ultramafic xenoliths from Adagdak volcano, Adak, Aleutian Islands, Alaska: deformed igneous cumulates from the Moho of an island arc.
Journal of Geology

95
,
329
–341.
DePaolo, D. J., Perry, F. V. & Baldridge, W. S. (
1992
). Crustal versus mantle sources of granitic magmas: a two-parameter model based on Nd isotopic studies.
Transactions of the Royal Society of Edinburgh

83
,
439
–446.
Dimalanta, C., Taira, A., Yumul, G. P. J., Tokuyama, H. & Mochizuki, K. (
2002
). New rates of western Pacific island arc magmatism from seismic and gravity data.
Earth and Planetary Science Letters

202
,
105
–115.
Dobran, F. (
1992
). Nonequilibrium flow in volcanic conduits and application to the eruptions of Mt. St. Helens on May 18, 1980, and Vesuvius in AD 79.
Journal of Volcanology and Geothermal Research

49
,
285
–311.
Ducea, M. N. (
2002
). Constraints on the bulk composition and root foundering rates of continental arcs: a California arc perspective.
Journal of Geophysical Research

107
(B11),
2303
.
Ducea, M. N. & Saleeby, J. B. (
1998
). The age and origin of a thick mafic–ultramafic keel from beneath the Sierra Nevada batholith.
Contributions to Mineralogy and Petrology

133
(1–2),
169
–185.
Ducea, M. N., Kidder, S. & Zandt, G. (
2003
). Arc composition at mid-crustal depths: insights from the Coast Ridge Belt, Santa Lucia Mountains, California.
Geophysical Research Letters

30
(13), article number 1703.
Ewart, A. (
1982
). The mineralogy and petrology of Tertiary–Recent orogenic volcanic rocks: with special reference to the andesite–basaltic compositional range. In: Thorpe, R. S. (ed.) Andesites. Chichester: John Wiley, pp. 25–95.
Feeley, T. C., Cosca, M. A. & Lindsay, C. R. (
2002
). Petrogenesis and implications of cryptic hybrid magmas from Washburn Volcano, Absaroka Volcanic Province, USA.
Journal of Petrology

43
,
663
–703.
Fialko, Y. A. & Rubin, A. M. (
1999
). Thermal and mechanical aspects of magma emplacement in giant dike swarms.
Journal of Geophysical Research

104
(10),
23033
–23049.
Fornelli, A., Piccarreta, G., Del Moro, A. & Acquafredda, P. (
2002
). Multi-stage melting in the lower crust of the Serre (southern Italy).
Journal of Petrology

43
(12),
2191
–2217.
Furlong, K. P. & Fountain, D. M. (
1986
). Continental crustal underplating: thermal considerations and seismic–petrologic consequences.
Journal of Geophysical Research

91
(B8),
8285
–8294.
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.
Contributions to Mineralogy and Petrology

119
,
197
–212.
Green, T. H. (
1972
). Crystallization of calc-alkaline andesite under controlled high-pressure, hydrous conditions.
Contributions to Mineralogy and Petrology

34
,
150
–166.
Green, T. H. (
1982
). Anatexis of mafic crust and high pressure crystallization of andesite. In: Thorpe, R. S. (ed.) Andesites; Orogenic Andesites and Related Rocks. Chichester: John Wiley, pp. 465–487.
Green, T. H. & Ringwood, A. E. (
1968
). Genesis of the calc-alkaline igneous rock suite.
Contributions to Mineralogy and Petrology

18
,
105
–162.
Griffin, W. L., Wang, X., Jackson, S. E., Pearson, N. J., O'Reilly, S. Y., Xu, X. & Zhou, X. (
2002
). Zircon chemistry and magma mixing, SE China: in-situ analysis of Hf isotopes, Tonglu and Pingtan igneous complexes.
Lithos

61
,
237
–269.
Grove, T. L., Elkins-Tanton, L. T., Parman, S. W., Chatterjee, N., Muntener, O. & Gaetani, G. A. (
2003
). Fractional crystallization and mantle-melting controls on calc-alkaline differentiation trends.
Contributions to Mineralogy and Petrology

145
(5),
515
–533.
Grunder, A. L. (
1995
). Material and thermal roles of basalt in crustal magmatism: case study from eastern Nevada.
Geology

23
(10),
952
–956.
Hannah, R. S., Vogel, T. A., Patino, L. C., Alvarado, G. E., Perez, W. & Smith, D. R. (
2002
). Origin of silicic volcanic rocks in central Costa Rica; a study of a chemically variable ash-flow sheet in the Tiribi Tuff.
Bulletin of Volcanology

64
(2),
117
–133.
Hart, G. L., Johnson, C. M., Shirey, S. B. & Clynne, M. A. (
2002
). Osmium isotope constraints on lower crustal recycling and pluton preservation at Lassen Volcanic Center, CA.
Earth and Planetary Science Letters

199
,
269
–285.
Hildreth, W. & Moorbath, S. (
1988
). Crustal contributions to arc magmatism in the Andes of Central Chile.
Contributions to Mineralogy and Petrology

98
,
455
–489.
Hirth, G. & Kohlstedt, D. L. (
1995
). Experimental constraints on the dynamics of the partially molten upper mantle; deformation in the diffusion creep regime.
Journal of Geophysical Research

100
(B8),
1981
–2001.
Hirth, G. & Kohlstedt, D. L. (
1996
). Water in the oceanic upper mantle; implications for rheology, melt extraction and the evolution of the lithosphere.
Earth and Planetary Science Letters

144
(1–2),
93
–108.
Holbrook, W. S., Mooney, W. D. & Christensen, N. I. (
1992
). The seismic velocity structure of the deep continental crust. In: Fountain, D. M., Arculus, R. & Kay, R. W. (eds) Continental Lower Crust. Amsterdam: Elsevier, pp. 1–34.
Holland, T. J. B. & Powell, R. (
1998
). An internally consistent thermodynamic data set for phases of petrological interest.
Journal of Metamorphic Geology

16
(3),
309
–343.
Holloway, J. R. & Burnham, C. W. (
1972
). Melting relations of basalt with equilibrium pressure less than total pressure.
Journal of Petrology

14
,
1
–29.
Hopson, C. A. & Mattinson, J. M. (
1994
). Chelan Migmatite Complex, Washington: field evidence for mafic magmatism, crustal anatexis, mixing and protodiapiric emplacement. In: Swanson, D. A. & Haugerud, R. A. (eds) Geologic Field Trips in the Pacific Northwest: 1994 Geological Society of America Annual Meeting. Seattle, WA: Department of Geological Sciences, University of Washington, pp. 1–21.
Huppert, H. E. & Sparks, R. S. J. (
1988
). The generation of granitic magmas by intrusion of basalt into continental crust.
Journal of Petrology

29
(3),
599
–624.
Jan, M. Q. & Howie, R. A. (
1981
). The mineralogy and geochemistry of the metamorphosed basic and ultrabasic rocks of the Jijal complex, Kohistan, NW Pakistan.
Journal of Petrology

22
,
85
–126.
Jull, M. & Kelemen, P. B. (
2001
). On the conditions for lower crustal convective instability.
Journal of Geophysical Research

106
(B4),
6423
–6446.
Kay, R. W. & Kay, S. M. (
1993
). Delamination and delamination magmatism.
Tectonophysics

219
,
177
–189.
Kay, R. W. & Mahlburg-Kay, S. (
1991
). Creation and destruction of lower continental crust.
Geologische Rundschau

80
,
259
–278.
Kay, R. W., Mahlburg Kay, S. & Arculus, R. J. (
1992
). Magma genesis and crustal processing. In: Fountain, D. M., Arculus, R. & Kay, R. W. (eds) Continental Lower Crust. Amsterdam: Elsevier, pp. 423–445.
Kobayashi, K. & Nakamura, E. (
2001
). Geochemical evolution of Akagi volcano, NE Japan: implications for interaction between island-arc magma and lower crust, and generation of isotopically various magmas.
Journal of Petrology

42
(12),
2303
–2331.
Lange, R. A. & Carmichael, I. S. E. (
1987
). Densities of Na2O–K2O–CaO–MgO–FeO–Fe2O3–Al2O3–TiO2–SiO2 liquids: new measurements and derived partial molar properties.
Geochimica et Cosmochimica Acta

51
,
2931
–2946.
Lange, R. A. & Navrotsky, A. (
1992
). Heat capacities of Fe2O3-bearing silicate liquids.
Contributions to Mineralogy and Petrology

110
,
311
–320.
Lee, C. T., Rudnick, R. L. & Brimhall, G. H. (
2001
). Deep lithospheric dynamics beneath the Sierra Nevada during the Mesozoic and Cenozoic as inferred from xenolith petrology.
Geochemistry, Geophysics, Geosystems

2
, article number 2001GC000152.
Lidiak, E. G. & Jolly, W. T. (
1996
). Circum-Caribbean granitoids; characteristics and origin.
International Geology Review

38
(12),
1098
–1133.
Marsh, B. D. (
1981
). On the crystallinity, probability of occurrence, and rheology of lava and magma.
Contributions to Mineralogy and Petrology

78
,
85
–98.
McCarthy, T. C. & Patiño-Douce, A. E. (
1997
). Experimental evidence for high-temperature felsic melts formed during basaltic intrusion of the deep crust.
Geology

25
(5),
463
–466.
McKenzie, D. & Bickle, M. J. (
1988
). The volume and composition of melt generated by extension of the lithosphere.
Journal of Petrology

29
(3),
625
–679.
Meissner, R. & Mooney, W. (
1998
). Weakness of the lower continental crust: a condition for delamination, uplift and escape.
Tectonophysics

296
(1–2),
47
–60.
Montgomery, D. R., Balco, G. & Willet, S. D. (
2001
). Climate, tectonics and the morphology of the Andes.
Geology

29
(7),
579
–582.
Mukasa, S. B. & Shervais, J. W. (
1999
). Growth of subcontinental lithosphere: evidence from repeated dike injections in the Balmuccia lherzolite massif, Italian Alps.
Lithos

48
,
287
–316.
Muntener, O., Kelemen, P. B. & Grove, T. L. (
2001
). The role of H2O during crystallization of primitive arc magmas under uppermost mantle conditions and genesis of igneous pyroxenites: an experimental study.
Contributions to Mineralogy and Petrology

141
(6),
643
–658.
Patankar, S. V. (
1980
). Numerical Heat Transfer and Fluid Flow. New York: Hemisphere.
Patankar, S. V. & Spalding, D. B. (
1972
). A calculation procedure for heat, mass, and momentum transfer in three-dimensional parabolic flows.
International Journal of Heat and Mass Transfer

15
,
1787
.
Patiño-Douce, A. & Beard, J. (
1994
). Dehydration-melting of biotite gneiss and quartz amphibolite from 3 to 15 kbar.
Journal of Petrology

36
(3),
707
–738.
Patiño-Douce, A. E. & Beard, J. (
1995
). Dehydration-melting of biotite gneiss and quartz amphibolite from 3 to 15 kbar.
Journal of Petrology

36
,
707
–738.
Patiño-Douce, A. E. & Johnston, A. D. (
1991
). Phase equilibria and melt productivity in the pelitic system: implications for the origin of peraluminous granitoids and aluminous granulites.
Contributions to Mineralogy and Petrology

107
,
202
–218.
Pedersen, T., Heeremens, M. & van der Beek, P. (
1998
). Models of crustal anatexis in volcanic rifts: applications to southern Finland and the Oslo Graben, southeast Norway.
Geophysical Journal International

132
(2),
239
–255.
Petford, N. & Gallagher, K. (
2001
). Partial melting of mafic (amphibolitic) lower crust by periodic influx of basaltic magma.
Earth and Planetary Science Letters

193
,
483
–499.
Pickett, D. A. & Saleeby, J. B. (
1993
). Thermobarometric constraints on the depth of exposure and conditions of plutonism and metamorphism at deep levels of the Sierra Nevada Batholith, Tehachapi Mountains, California.
Journal of Geophysical Research

98
(B1),
609
–629.
Pinkerton, H. & Stevenson, R. J. (
1992
). Methods of determining the rheological properties of magmas at sub-liquidus temperatures.
Journal of Volcanology and Geothermal Research

53
(1–4),
47
–66.
Raia, F. & Spera, F. J. (
1997
). Simulations of crustal anatexis: implications for the growth and differentiation of continental crust.
Journal of Geophysical Research

102
(B10),
22629
–22648.
Rapp, R. P. & Watson, E. B. (
1995
). Dehydration melting of metabasalt at 8–32 kbar: implications for continental growth and crust–mantle recycling.
Journal of Petrology

36
(4),
891
–931.
Rapp, R. P., Watson, E. B. & Miller, C. F. (
1991
). Partial melting of amphibolite/eclogite and the origin of Archean trondhjemites and tonalites.
Precambrian Research

51
,
1
–25.
Renner, J., Evans, B. & Hirth, G. (
2000
). On the rheologically critical melt fraction.
Earth and Planetary Science Letters

181
,
585
–594.
Rudnick, R. L. (
1990
). Growing from below.
Nature

347
,
711
–712.
Rudnick, R. L. & Taylor, S. R. (
1987
). The composition and petrogenesis of the lower crust: a xenolith study.
Journal of Geophysical Research

92
(B13),
13981
–14005.
Rushmer, T. (
1991
). Partial melting of two amphibolites: contrasting experimental results under fluid-absent conditions.
Contributions to Mineralogy and Petrology

107
,
41
–59.
Rushmer, T. (
1995
). An experimental deformation of partially molten amphibolite: application to low-melt fraction segregation.
Journal of Geophysical Research

100
,
15681
–15695.
Saleeby, J. & Foster, Z. (
2004
). Topographic response to mantle lithosphere removal in the southern Sierra Nevada region, California.
Geology

32
(3),
245
–248.
Saleeby, J., Ducea, M. & Clemens-Knott, D. (
2003
). Production and loss of high-density batholithic root, southern Sierra Nevada, California.
Tectonics

22
(6), article number 1064.
Scaillet, B., Holtz, F. & Pichavant, M. (
1998
). Phase equilibrium constraints on the viscosity of silicic magmas; 1. Volcanic–plutonic comparison.
Journal of Geophysical Research

103
(B11),
27257
–27266.
Sen, C. & Dunn, T. (
1994
). Dehydration melting of a basaltic composition amphibolite at 1·5 and 2·0 GPa: implications for the origin of adakites.
Contributions to Mineralogy and Petrology

117
,
394
–409.
Shaw, H. R. (
1972
). Viscosities of magmatic liquids: an empirical method of prediction.
American Journal of Science

272
,
870
–893.
Sisson, T. W., Ratajeski, K., Hankins, W. B. & Glazner, A. F. (
2005
). Voluminous granitic magmas from common basaltic sources.
Contributions to Mineralogy and Petrology

148
(6),
635
–661.
Touloukian, Y. S., Judd, W. R. & Roy, R. F. (eds) (
1981
). Physical Properties of Rocks and Minerals. New York: McGraw–Hill.
Tulloch, A. J. & Kimbrough, D. L. (
2003
). Paired plutonic belts in convergent margins and the development of high Sr/Y magmatism: the Penisular Ranges Batholith of California and the Median Batholith of New Zealand.
Geological Society of America, Special Papers

374
,
275
–295.
Versteeg, H. K. & Malalasekera, W. (
1995
). An Introduction to Computational Fluid Dynamics: the Finite Volume Method. Harlow: Pearson Education.
Vigneresse, J. L. & Tikoff, B. (
1999
). Strain partitioning during partial melting and crystallizing felsic magmas.
Tectonophysics

312
,
117
–132.
Vigneresse, J. L., Barbey, P. & Cuney, M. (
1996
). Rheological transitions during partial melting and crystallization with application to the felsic magma segregation and transfer.
Journal of Petrology

37
,
1579
–1600.
Voller, V. R. & Swaminathan, C. R. (
1991
). General source-based method for solidification phase change.
Numerical Heat Transfer

19B
,
175
–189.
1994
). On the relationship between dike width and magma viscosity.
Journal of Geophysical Research

99
(9),
17743
–17755.
Wells, P. R. A. (
1980
). Thermal models for the magmatic accretion and subsequent metamorphism of continental crust.
Earth and Planetary Science Letters

46
,
253
–265.
Williams, M. L., Hanmer, S., Kopf, C. & Darrach, M. (
1995
). Syntectonic generation and segregation of tonalitic melts from amphibolite dikes in the lower crust, Striding–Athabasca mylonite zone, northern Saskatchewan.
Journal of Geophysical Research

100
(8),
15717
–15734.
Wolf, M. B. & Wyllie, P. J. (
1994
). Dehydration-melting of amphibolite at 10 kbar: the effects of temperature and time.
Contributions to Mineralogy and Petrology

115
,
369
–383.
Younker, L. W. & Vogel, T. A. (
1976
). Plutonism and plate tectonics: the origin of circum-Pacific batholiths.

14
,
238
–244.
Zellmer, G. F., Hawkesworth, C. J., Sparks, R. S. J., Thomas, L. E., Harford, C. L., Brewer, T. S. & Loughlin, S. C. (
2003
). Geochemical evolution of the Soufrière Hills volcano, Montserrat, Lesser Antilles volcanic arc.
Journal of Petrology

44
(8),
1349
–1374.