Journal of Petrology Advance Access originally published online on June 3, 2005
Journal of Petrology 2005 46(11):2167-2195; doi:10.1093/petrology/egi049
© The Author 2005. Published by Oxford University Press. All rights reserved. For Permissions, please e-mail: journals.permissions@oxfordjournals.org
Lower Crustal Magma Genesis and Preservation: a Stochastic Framework for the Evaluation of Basalt–Crust Interaction
J. DUFEK* and
G. W. BERGANTZ
DEPARTMENT OF EARTH AND SPACE SCIENCE, UNIVERSITY OF WASHINGTON, BOX 351310, SEATTLE, WA 98195, USA
RECEIVED
SEPTEMBER 4, 2004;
ACCEPTED
APRIL 13, 2005
 |
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 K
2O, typical of mature arc settings, without the associated
mass balance problem.
KEY WORDS: crustal anatexis; delamination; lower crust; magma mixing; thermal model
|
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 x 10–4 to 1·0 x 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.
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) |
where
is the reported melt volume (or length reported in one-dimensional simulations) and
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) |
where
Conservation of mass:
 | (4) |
Momentum:
 | (5) |
Species conservation:
 | (6) |
Side
boundary conditions:
 | (7) |
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.
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).
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
).
View larger version (24K):
[in this window]
[in a new window]
|
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) |
where
 | (9) |
 | (10) |
 | (11) |
and
 | (12) |
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).
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
).
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.
View larger version (35K):
[in this window]
[in a new window]
|
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/m
2 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/m
2 (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/m
3 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/m
3 and the mantle
region had a heat production of 0·02 mW/m
3 (Chapman &
Furlong, 1992
).
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) |
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.
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
x 10
3 years) and long (10
7 years) time scales examined magma
mixing and ductile creep, respectively.
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 m
3/m
2 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/m
3. 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.
View larger version (46K):
[in this window]
[in a new window]
|
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 <10
6 years. At lower fluxes of basalt
(<0·001 m
3/m
2 per year), individual intrusions result
in localized, transient thermal anomalies that are primarily
responsible for melt production.
View larger version (21K):
[in this window]
[in a new window]
|
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
, 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.
View larger version (29K):
[in this window]
[in a new window]
|
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 m
3/m
2 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 m
3/m
2 per year and
0·005 m
3/m
2 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 10
6 years following
the initiation of basalt intrusion, only scenarios (b) and (c)
with a basalt flux of 0·005 m
3/m
2 per year produce crustal
melts. All model realizations with basalt fluxes of 0·0001
m
3/m
2 per year create thermal anomalies over 10
6 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 x 107 years to reach 0·1. Even after 5·0 x 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 x 107 years and 5·0 x 107 years (Dimalanta et al., 2002). In contrast, mean crustal melt fractions exceeding 0·4 at 50 km depth occur after 1·5 x 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.
View larger version (22K):
[in this window]
[in a new window]
|
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.
View larger version (24K):
[in this window]
[in a new window]
|
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 m
3/m
2 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 10
6 years.
View larger version (23K):
[in this window]
[in a new window]
|
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 10
6 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 .
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.
View larger version (24K):
[in this window]
[in a new window]
|
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
m
3/m
2 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.
Similar results are produced after 10
7 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.
View larger version (43K):
[in this window]
[in a new window]
|
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/m
2 surface heat flow, very little to no crustal melting can be
expected for geologically constrained basalt fluxes. At 10
6 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 10
6 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 % SiO
2.
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 x 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 x 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.
View larger version (15K):
[in this window]
[in a new window]
|
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 x 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
TOP
ABSTRACT
INTRODUCTION
