## Abstract

A new fracture-mediated intrusion model resolves the sequence of magma and rock displacements generating a felsic magma system with a lower crustal source, central conduit and shallow sill pluton. Idealized intraplate conditions are assumed, to neglect regional tectonism and to focus on juvenile cracking by magma-intrinsic hydraulic and buoyant loads.

The magma source is conductively heated and develops by endothermic fluid-absent melting in approximately 10^{6} years. The idealized domical thermal anomaly and endothermic heat focusing yield a low aspect ratio source, with outer-porous and inner-permeable partial melt zones. An anatectic core region is unrealized owing to magma segregation.

Thermal stresses are readily relaxed and unimportant to source loading while crustal uplift generates tensile stress and, upon relaxation, lateral space for tensile fractures. Dilative melting generates buoyancy overpressure (Δ*P*_{B}) and a hydraulic contribution (Δ*P*_{V}) to the magma pressure (*P*_{M}). Δ*P*_{V} develops by elastic wall-rock compression as the ‘excess magma volume’, *EMV*, arises too abruptly for full relaxation by inelastic deformation, inducing a brittle response. Tensile rupture criteria are met in an effective tensile stress field with low differential stress induced by magma pore pressure and wedging by pressurized cracks, which initiate by source inflation and uplift. Preferred vein geometry reflects the starting stress field. For symmetric doming, radial vertical cracks with a central nexus form a natural conduit. A vertically extensive crack system, however, requires special explanation because wedging by Δ*P*_{V} reorients dykes to sills just above source. The solution is that volumetric crack growth accommodating non-relaxed *EMV* (*EMV**) causes Δ*P*_{V}→0. Magma transport becomes buoyancy-driven and the Δ*P*_{V} problem does not arise. The critical sill intrusion depth, *I*, is where Δ*P*_{B} exceeds the regional vertical stress curve, where columns must intrude owing to Δ*P*_{B} alone. Sill growth is mainly by floor depression, involving ductile shear of lower crust, creating sill volume, suppressing roof uplift, expelling source contents, processing protolith through the melting zone, reducing stress, σ_{H}, and widening conduits for sustained flow.

Two intrusive regimes are identified; Δ*P*_{V}>0 (hydraulic inflation) and Δ*P*_{V}=0 (buoyancy pumping). Partitioning between three sinks for *EMV* – inelastic uplift (ϕ), crack growth (η) and a non-relaxed portion (*EMV**) generating Δ*P*_{V} – defines four hydraulic subregimes. *Disequilibrium dilation* occurs during crustal relaxation prior to rupture, when η=0 and *EMV* partitions between ϕ and *EMV**. Uplift occurs readily owing to the crust's weakness in flexure, so ϕ abruptly increases while *EMV** decreases, causing abrupt variations in source failure mode, geometry and rate that smooth initial Δ*P*_{V} variations. During *equilibrium dilation* source swelling continues with ϕ dominant over *EMV**. Dilatant loading rates mean that positive Δ*P*_{V} is always maintained however, keeping the source near-isotropically inflated and prepared for rupture. *Disequilibrium cracking* begins when uplift-driven horizontal stretching initiates rupture and crack growth (η). Crack volume is initially small, but readily enlarges as dykes propagate by conversion of stored Δ*P*_{V}. *EMV* is minimized better and faster by η than by ϕ, owing to crack-tip stress concentration, giving abrupt augmentation of η and decreases in ϕ and *EMV** in a crack growth-surge until the uplift-modified stress field is balanced. In *equilibrium cracking*, once *EMV** (and Δ*P*_{V}) decrease to incipient levels, each new increment of *EMV** partitions directly into crack growth, while continued uplift maintains vertical rupture, generating a vertically extensive fracture system.

The absolute volumetric equivalence of *EMV**, at most a few tens of km^{3}, will be exhausted during dyking, sill intrusion or surface eruption. The system then becomes buoyancy-driven and, if depth *I* is reached, must intrude a sill. Relaxation of sill underburden initiates crustal decoupling and buoyancy pumping, where the downward underburden-flux drives and is balanced by upward magma flow.

## Nomenclature

*a*- sill radius
*c*- specific heat
*D*- intrusion depth
*E*- Young's modulus
*EMV*- excess magma volume
*EMV**- non-relaxed excess magma volume
*g*- acceleration due to gravity
*h*- source or magma column height
*I*- critical overburden thickness
*k*- thermal conductivity
*M*- melt proportion
*P*- ambient pressure
*P*_{M}- magma pressure
*Q*_{C}- conduit flux
*Q*_{U}- underburden flux
*r*- conduit radius
*S*- source depth
*t*- time
*T*- temperature
*V*- volume
*W*- sill half-height
*Y*- vertical shear displacement
*z*- depth (positive down)
- α
_{L} - linear coefficient of thermal expansion
- β
- thermal gradient
- γ
- heat-focusing factor
- Δ
*H*_{r} - specific heat of reaction
- Δ
*P*_{B} - buoyancy overpressure
- Δ
*P*_{E} - magma excess pressure
- Δ
*P*_{V} - volume-related hydraulic overpressure
- Δ
*T* - temperature change
- Δ
*V*_{M} - coefficient of melting-related volume increase
- ϵ
- volumetric strain
- ϵ
_{H} - horizontal strain
- ϵ
_{V} - vertical strain
- ϕ
- coefficient of EMV relaxation by inelastic deformation
- η
- coefficient of EMV relaxation by crack growth
- κ
- thermal diffusivity
- µ
_{M} - magma viscosity
- ν
- Poisson ratio
- ρ
_{R} - country rock density
- ρ
_{M} - magma density
- σ
_{1} - maximum principal stress
- σ
_{2} - intermediate principal stress
- σ
_{3} - minimum principal stress
- σ
_{d} - differential stress
- σ
_{H} - horizontal (confining) stress
- σ
_{HC} - horizontal circumferential stress
- σ
_{HR} - horizontal radial stress
- σ
_{L} - lithostatic stress
- σ
_{N} - normal stress
- σ
_{S} - shear stress
- σ
_{V} - vertical stress
- σ
_{H}^{(th)} - horizontal thermal stress
- −σ
_{H}^{(subs)} - subsidence related horizontal tensile stress
- −σ
_{H}^{(up)} - uplift related horizontal tensile stress
- τ
- tensile strength (subscripts V and H refer to vertical and horizontal strengths)
- ζ
- flexural rigidity

Most granitoid plutons are shallow (2–10 km) and tabular, suggesting their intrusion as sills or laccoliths (Corry 1988; McCaffrey & Petford 1997; Vigneresse *et al.* 1999), while the felsic magmas involved are mostly deep crustal anatectites (e.g. Kistler & Peterman 1973; Chappell & White 1974; Atherton *et al.* 1979; DePaolo 1981; McCulloch & Chappell 1982). These facts require that plutons were linked to deep-crustal magma source regions during trans-crustal magma transport. However, owing to the infrequent, dynamic and deep-seated nature of this process, its mechanism is contentious. Nevertheless, there is now increasing consensus that the intrusive displacements of magma and rock involve the propagation and exploitation of hydraulic fractures (Clemens & Mawer 1992; Vigneresse *et al.* 2000). Moreover, considerable progress has been made in understanding the putative fracture-mediated process, including field, experimental and theoretical studies of melting, deformation and magma segregation in source regions (e.g. Vielzeuf *et al*. 1990; Sawyer 1991; Petford 1995; Rushmer 1995; Rutter & Neumann 1995; Vigneresse *et al.* 1996; Petford & Koenders 1998; Daczko *et al.* 2001; Simakin & Talbot 2001; Klepeis *et al.* 2003), transport by trans-crustal dyke systems (Takada 1989; Petford *et al.* 1994; Vigneresse 1995*a*, *b*; Petford 1996) and tabular pluton emplacement (Cruden, 1998, 2006; Vigneresse *et al.* 1999; Cruden & McCaffrey 2001). Notwithstanding this, several issues are unresolved, notably the relative roles of magma buoyancy, hydraulic overpressure and regional tectonic forces in driving the process, the role of pre-existing weaknesses as magma pathways, how a fracture conduit of adequate aperture might develop and persist, and what controls the emplacement depth and growth of sills. Owing to this, a coherent mechanical account from source to post-emplacement has yet to be given. This is the present objective.

Several factors account for these uncertainties. While shallow plutons commonly are exposed, complete intrusive systems are never seen owing to their scale. Relict source regions and fracture conduits are rarer in the rock record than is attributable only to their deeper origins, probably reflecting magma drainage. The geological and tectonic diversity of crustal settings involved complicates judgement of the roles of pre-existing structure and regional stresses (Brown 1993, 1994; Vigneresse 1995*a*). This ambiguity over boundary conditions hinders efforts to model the process. Analogue models have inherent scaling problems that compromise, for example, the validity of the diapir hypothesis (Ramberg 1981) for felsic magmas (Petford 1996; Vigneresse & Clemens 2000). Numerical methods cannot, at present, holistically simulate the continuum of subprocesses, despite robust theory for each, owing to complex coupling over diverse length- and timescales, and the magma–host rock strain differential. If a mechanical description of fracture-mediated intrusion is to cohere, theoretical methods must presently suffice to address outstanding issues.

## Fracture-mediated intrusion

There are broadly two views of the fracture-mediated felsic intrusion mechanism. Fracture formation and exploitation by magma obey established principles (e.g. Jaeger & Cook 1976; Pollard 1987) so these are distinguished essentially by the origin of the cracks and nature of the driving forces. One derives from field studies in settings where tectonic forces were strong, and holds that magma exploits pre-existing faults driven by regional tectonic loads and buoyancy (e.g. Hutton 1982; Hutton *et al.* 1990; Atherton 1990; McCaffrey 1992; Tikoff & Teyssier 1992; Brown 1994). The other posits the magma to form and exploit juvenile fractures under hydraulic and buoyant loads intrinsic to the magma (Takada 1989; Clemens & Mawer 1992). Two recent sets of geological findings point to the latter as a better basis for a general model. First, there is the identification in shallow intrusive architecture of a generic structural style significantly independent of tectonic setting – plutons statistically are sill-shaped (McCaffrey & Petford 1997) while their conduit systems, conjectured to develop as systems of dykes, are characteristically centralized and pipe-like (Vigneresse *et al.* 1999). Second, there is the deduction from observations of frozen, deep-crustal, anatectic source regions of rapid magma drainage into tensile cracks formed during dilative melting (Daczko *et al.* 2001; Klepeis *et al.* 2003). Taken together, these findings suggest: (i) that juvenile cracking may occur at all structural levels, obviating any need for inherited fractures; and (ii) a crack system may propagate owing to loading by dilative melting and magma buoyancy, obviating any need for tectonic forces. Indeed, the common causative mechanical influence implied by generic structural style may, if universals like crustal weight are excluded, only be the magma. Moreover, the diversity of tectonic settings involved precludes specific crustal structures or regional stress patterns being prerequisites, irrespective of ad hoc roles. A general model must account for intrusion across all tectonic settings. An emphasis on juvenile fracture by magmatic loading thus offers more fundamental insight than studies of diverse natural systems and simplifies analysis of the key parameters – those controlling the intrusive stress field.

## Approach

Previous workers recognized the significance of the stress field in fracture-mediated intrusion. Johnson & Pollard (1973) reviewed factors controlling the buoyancy overpressure of a static magma column and analysed the deformation in sill emplacement. Subsequent work on felsic intrusive architecture has followed their static approach, a valuable example being Vigneresse *et al.* (1999), who used the insights of Jaeger & Cook (1976) and Parsons & Thompson (1991) to identify stress exchange between magma-filled fractures and host rock – magma wedging – as the key feedback between propagating cracks and the causative stress field. Approaching from a petrological perspective, Clemens *et al.* (1997) and Clemens & Droop (1998) established dilative fluid-absent melting as the likely cause of hydraulic source loading. To build on these ideas we consider theoretically how a felsic system might develop the generic style identified above from source upward. We adopt the premise that if magmatic loading is responsible, we should examine circumstances where this would be manifested most clearly, i.e. powerful, voluminous systems in structurally simple lithosphere in the absence of remote loading. We therefore consider the hypothetical case of a giant crustal felsic system intruded under idealized anorogenic conditions, isolated from the influence of mantle-derived magma. This offers the simplest basis for future efforts to encompass systems rooted in the mantle and those in complex crustal environments. As several significant subprocesses are described by robust theory, this work is in part a review. Our contribution is to address the outstanding issues and link the process in the simplest conceivable case.

A hypothetical felsic intrusive system having the generic structure above is shown in Figure 1. It comprises a partially molten lower-crustal source region, a central conduit and a shallow, sill-like pluton. We analyse how the stress field required to produce this might arise. Idealized intraplate stress conditions are assumed and reference made to a standard section of hypothetical continental lithosphere (Table 1) based on studies such as Meissner (1986) and Salisbury & Fountain (1990). We consider the geometric and structural development of the source under simple genetic circumstances, the origin and nature of the conduit, magma ascent, depth of pluton emplacement and post-emplacement processes. Thermal aspects are outside our scope (but see Petford *et al.* 1994), although findings do permit us to reassess the argument that felsic dyke conduits are thermally non-viable (Rubin 1995). Static analysis of depth–stress relations between magma and crust is employed, considering dynamical aspects only qualitatively. Nevertheless, by analysing intrusive space relations this approach identifies the complimentary roles of hydraulic and buoyancy forces, and yields a coherent sequence of events. Findings are immediately relevant to intraplate ‘supervolcano’ systems.

## Magma source evolution

The geometry, developmental timescale and internal structure of the magma source will greatly influence its deformation history. We address these factors first.

### Geometry and growth

We focus on anatectic source regions, whose occurrence depends on crustal composition, pressure, *P*, and temperature, *T* (Tuttle & Bowen 1958). Most felsic magmas are initially H_{2}O-undersaturated partial melts of hydrous, chiefly biotite- and amphibole-bearing assemblages (Eggler 1973; Thompson 1982; Clemens 1990; Clemens *et al.* 1997; Clemens & Watkins 2001), so model melting is assumed fluid-absent. Four thermodynamic features of such reactions significantly influence source evolution. First, positive entropy changes cause isograds to follow appropriate isotherms with little kinetic lag (Thompson & Tracy 1979; Vielzeuf *et al.* 1990). Source regions are thus closely spatially and temporally constrained by causative isotherm distributions. Second, they are markedly endothermic with specific heats of reaction Δ*H*_{r} of 1.0×10^{5} to 2.5×10^{5} Jkg^{−1}. By focusing heat and buffering *T* over the melting interval they exert control on the isotherm distribution (Rice & Ferry 1982). Third, for protoliths with a dominant hydrous mineral, most melting is in a narrow thermal interval 25–50 °C wide for biotite rocks, starting at 825–850 °C (Vielzeuf & Holloway 1988; Stevens *et al.* 1997), and approximately 100 °C wide for amphibolites starting at 825 to 900 °C (Rutter & Wyllie 1988; Beard & Lofgren 1991; Rushmer 1991; Wolf & Wyllie 1994; Patiño Douce & Beard 1995). Ignoring multiple hydrous phases with stepped melt production (Skjerlie and Johnston 1992), protoliths with ≥*c*. 1 wt% H_{2}O yield melt fractions (*M*)=0.4–1.0 if *T* reaches 900–1000 °C, so comprehensive melting needs only a 50–100 °C perturbation above the fluid-absent solidus (Clemens & Vielzeuf 1987). Fourth, the volume change on melting [Δ*V*_{M}=(*V*_{products}−*V*_{reactants})/*V*_{reactants}] is positive, between 0.02 and 0.20 depending on reaction stoichiometry, completion, *P* and magma compressibility (Clemens & Mawer 1992). Thus: (a) anatectites are less dense than solid counterparts and buoyant; and (b) as fluid-absent magma has low compressibility, Δ*V*_{M} may, depending on the mechanical response of enclosing rock, yield a hydraulic contribution to the magma pressure *P*_{M} (Clemens & Mawer 1992; Rushmer 1995).

For an idealized source we ignore unusual crustal structure (e.g. isolated non-refractory domains), neglect ongoing tectonism by considering only isobaric fusion (cf. Brown 1993) and invoke an idealized thermal anomaly. Copious crustal fusion requires mantle heat (but see Sandiford *et al.* 1998) transported by advection (A) by mantle-derived magma (e.g. Huppert & Sparks 1988; Petford 1995) or by conduction (C). However, as the mechanical influence of intruding mafic magma is outside our scope, so A-type or mixed A/C-type sources are also neglected despite their probable significance. Our idealized source develops only by C-type heating from below. Besides protolith structure, source geometry thus depends on the thermal anomaly structure, heat focusing by endothermic melting and how these control the isotherm distribution. Amphibolite is adopted as an idealized protolith.

An idealized lithospheric thermal anomaly is circular, forming a domical melting region as isotherms *T*≥*T*_{solidus} intersect the crust. Anomalies generating large systems like Yellowstone (cf. Smith & Braille 1994) are 600–1000 km across (Anderson 1998), so the lateral extent of an idealized melting region far exceeds the lithosphere depth. This constrains the source to have a low geometric aspect ratio. Notwithstanding the transient nature of crustal thermal perturbations, constraints on the vertical melting distribution within such an anomaly are offered by temperature (*T*) variations at the Moho for natural upper-mantle anomalies, estimated at ±200 °C (Anderson 1998), and the use of steady-state geotherms. Values of the heat flux (*Q*_{H}) reasonable for non-cratonic continents are consistent with C-type melting of appropriate lithologies in the lowest part of thick crust for a modest *T* perturbation. For example, assuming the lower crust of the amphibolite having *T*_{solidus}= 800 °C, a main melting pulse at 850–950 °C and *T*_{liquidus}=1050 °C, then, using Chapman's (1986) isotherm distributions, a 40 km-thick crust would be everywhere subsolidus at *Q*_{H}≤70 mW m^{−2} as the 800 °C isotherm is predicted to lie 5 km below the Moho (mantle solidus is>1100 °C). However, at *Q*_{H}=80 mW m^{−2} the 800 °C isotherm occurs nominally at 35 km depth, so partial melting might be expected up to 5 km above the Moho. For *Q*_{H}=80 m Wm^{−2} the 850 °C isotherm is located at 38 km, suggesting also that the lowest 2 km of the crust would be within the main melting pulse. At *Q*_{H}=90 mW m^{−2} Chapman's models predict an 800 °C isotherm at 28 km depth, with 850 °C and 950 °C isotherms at approximately 31 and 37 km, respectively, so the lowest 12 km of crust might be expected to partially melt. However, for endothermic melting by transient heating, heat focusing would cause isotherms to lie deeper (Rice & Ferry 1982), reducing the spatial melting interval as follows.

The effect of Δ*H*_{r} on the spatial melting interval corresponding to a *T* interval, Δ*T*, may be estimated from the effective specific heat, *c*′, of the melting region. This relates to *c* of solid protolith as (Carslaw & Jaeger 1959) *c*′=*c*+(Δ*H*_{r}/Δ*T*). Assuming a steady increase in *Q*_{H} with just sufficient time for isotherms to reach equilibrium positions, the retarding effect of Δ*H*_{r} on melting isotherms (*T*_{solidus}<*T*<*T*_{liquidus}) may be assessed from how *c*′ influences the effective thermal diffusivity in the melting region according to κ′=*k*′/ρ′_{P}*c*′, where *k*′ is the effective thermal conductivity and ρ′_{P} the density of the partially molten rock. Melting thus requires extra heat, which retards ascending isotherms by a heat-focusing factor γ=κ/κ′. For hypothetical amphibolite lower crust, γ may be evaluated using Δ*T*=100 K, Δ*H*_{r}=2×10^{5} J kg^{−1}, ρ=2900 kg m^{−3}, *k*=2.5 W m^{−1} K^{−1}, and *c*=900 J kg^{−1} K^{−1}. For solid protolith at *T* close to *T*_{solidus}, κ is approximately 8×10^{−7} m^{2} s^{−1}, while in the melting zone, with average properties ρ′_{P}=2700 kg m^{−3} and *k*′=2.0 W m^{−1} K^{−1}, then *c*′=2900 J kg^{−1} and κ′=2.6×10^{−7} m^{2} s^{−1}, so γ is approximately 3. For our simple model, a heat-focusing factor, γ=3, reduces the separation between the base of the melting region (Moho) and isotherm *T*_{solidus} (800 °C) from approximately 5 to about 1.7 km at *Q*_{H}=80 mW m^{−2}. This would triple the thermal gradient over the melting interval and correspondingly reduce the gradient in non-melting crust above. At *Q*_{H}=90 mW m^{−2} the 800 °C isotherm would lie at 36 km, with the main melting pulse between 37 and 39 km.

These simple considerations illustrate three points. First, the Moho is a major step in primary fertility relative to likely thermal gradients, so C-type melting will always initiate there. Few lithological and superimposed thermal profiles would permit multiple melting horizons, especially with heat focusing. Second, owing also to heat focusing, putative C-type source regions have a low aspect ratio. For the *Q*_{H} and γ values given above, maximum source heights of approximately 4 km are reasonable for thick crust, which given the surface expression of large natural systems suggests source aspect ratios of 0.001–0.01. Third, if transient heating perturbs isotherms on a timescale *t*=*z*^{2}/κ, then heat focusing yields a modified timescale *t*′ for isotherm ascent through a melting region *t*′≈*z*^{2}/κ′. For amphibolite this predicts a reduced isotherm advance rate of approximately 0.003 m year^{−1}. An idealized partial melting region would thus require approximately 8×10^{5} years to reach 2.5 km thickness. Given an appropriate thermal anomaly, the timescale for C-type source formation is thus about 10^{6} years.

### Internal structure

The chief consideration regarding how the isotherm distribution and degree of partial fusion (*M*) relate to source internal structure concerns whether the protolith is open or closed to melt escape. Closed-system fusion progressively changes solid rock to liquid magma (Arzi 1978), whereas under open-system conditions melt may segregate from solid residue (Sawyer 1991). There are two threshold *M* values during progressive fusion, a percolation threshold (PT) defining interconnection of melt pockets at grain boundaries and a subsequent rheological transition (RT) defining mechanical protolith breakdown (Vigneresse *et al.* 1996; Vigneresse & Tikoff 1999). Absolute *M* values for these transitions *M*^{PT} and *M*^{RT} depend on phases, textures, stresses, etc. (Rutter & Neumann 1995), but are less significant for source structure than how *M* varies with *T*, as illustrated in Figure 2.

Figure 2a shows a hypothetical sigmoidal melt-productivity curve for a simple hydrous protolith. Melt forms initially at the solidus *T*_{1}, attains greatest productivity over the melting pulse *T*_{2} to *T*_{3}, and ends production at the liquidus *T*_{4}. T_{3} is a ‘pseudo-liquidus’ above which residual solid persists to much higher *T*. For *M*^{PT}≈0.08 and *M*^{RT}≈0.25 (Vigneresse & Tikoff 1999) such a melting curve implies for model amphibolite a PT effect near the start of the melting pulse and an abrupt RT somewhere before the middle. The model solidus, PT and RT map onto a hypothetical, mature C-type source in Figure 2b. If deformation is neglected and equilibrium *T* and *M* values assumed, then for closed-system behaviour three concentric zones are defined. The upper (and first to form as isotherms ascend) is a porous zone (*M*<0.08) in which melt is in isolated cavities (see Rutter & Neumann 1995). Next is a permeable zone (0.08<*M*<0.25) where stress may be transmitted hydraulically by melt, which may percolate in response to pressure gradients. Downward (up-*T*) melt cavities increase in size and connectivity towards the RT, the top of a hypothetical, extensively melted inner zone, not supporting shear stress (*M*>0.25), which could be designated an *anatectic core*. However, as a source might become an unstable open system at *M*<*M*^{RT}, such a putative core region may never be realized. It is a special type of *hydraulic nucleus*, which describes any magma-filled cavity system that initiates fracture-mediated magma ascent. The nature of the hydraulic nucleus is analysed after preliminary mechanical considerations.

## Deformation of host lithosphere

In model lithosphere one of the principal stresses (σ_{1}, σ_{2}, σ_{3}) is assumed vertical (σ_{V}) while the horizontal stresses are taken initially as a confining stress (σ_{H1}=σ_{H2}=σ_{H}). The average stress at depth *z* is the pressure *P*_{(z)}. The model top is a free surface and hydrostatic stress is assumed at the Moho, while arbitrary lateral boundaries are far from any intrusion. Model lithosphere is weak in flexure and tension so vertical strain ϵ_{V}, is assumed to occur freely at low rates (Turcotte & Schubert 1982). For horizontal strain, ϵ_{H}, a non-confined state (ϵ_{H}≠0) is assumed over times permitting lateral boundary migration, whereas over shorter times effective confinement is assumed (ϵ_{H}=0). Far-field tectonic loading is neglected. The influence of crustal strength and relevant non-magmatic loads is summarized as follows.

### Crustal strength

Maximum differential stress (σ_{d}=[σ_{1}−σ_{3}]) at failure depends mainly on failure mode, rock type, *T* and loading rate (Evans & Kohlstedt 1995). Tensile shear strength curves for crustal analogues, wet quartz for weak upper crust and wet diopside for stronger lower crust, calibrated for a loading rate of 10^{−10} s^{−1} (Kohlstedt *et al.* 1995), are shown in Figure 3. This rate approximates the lower limit of the crust's short-term elastic strength, which is most relevant here as abrupt loading and elastic behaviour will accompany dilative melting as shown below. The curves illustrate that under such conditions crustal materials can indefinitely sustain σ_{d} many times their tensile strength, τ, while the strain rate-dependence of strength, which explains how nominally plastic materials may undergo brittle failure (e.g. Hibbard & Watters 1985), allows for elastic tensile failure throughout the crust if the unbalanced loads responsible require strain at rates of more than about 10^{−10} s^{−1} to relax them (cf. Pfiffner & Ramsay 1982). As strength is also composition-dependent, crustal structure also influences the stress field, because strong materials attract high σ_{d} (Hogan & Gilbert 1995) and strength contrasts concentrate stresses. To neglect these factors, vertically gradational crustal properties and an idealized initial stress distribution are assumed.

### Non-magmatic loads

Model stresses from body forces, thermal dilation and lateral stretching accompanying vertical strain are considered with reference to Figure 3. For free ϵ_{V} then σ_{V(z)} is the overburden weight ρ_{R}*gz*. Only very abrupt loading would cause σ_{V} to depart from hydrostatic. For σ_{H}, in deep, intraplate crust, a lithostatic starting state is expected (σ_{H}=σ_{V}). However, to provide an example stress field with a lower σ_{H}/σ_{V} ratio, Figure 3 also shows the horizontal elastic body force σ_{H}^{(bo)}=*v*_{(z)}σ_{V(z)}/(1−*v*_{(z)}), where *v* is Poisson's ratio (e.g. Price 1959). The possible genetic significance of such a stress field is considered subsequently.

#### Thermal stresses

For ϵ_{V}≠0 and ϵ_{H}=0 thermal dilation generates vertical displacement and a thermal stress, σ_{H}^{(th)} (tensile is negative). If a thermal gradient, β, increases to β′ yielding Δ*T*=(β′−β)*z* (Turcotte & Schubert 1982), then over the heated column ϵ_{V}=α_{L}Δ*T*, where α_{L} is the linear coefficient of thermal expansion. If ϵ_{H}=0, then the thermal stress contribution σ_{H}^{(th)}=*E*α_{L}Δ*T*, where *E* is Young's modulus. For example, the curve σ_{H1}^{(L+th)} (Fig. 3) shows the effect (summed with the lithostat) of an average Δ*T* of 3.75 °C km^{−1} for which the total displacement is 88 m (ϵ_{V}=2.2×10^{−3}). In 40 km-thick model crust this would accompany a 10 mW m^{−2} increase in *Q*_{H}. At a strain rate of 10^{−10} s^{−1} this ϵ_{V} would require only a few thousands of years, confirming that σ_{V}^{(th)} normally is essentially fully relaxed, while σ_{H}^{(th)} dissipates quickly after heating. Thermal stresses are, therefore, unlikely to be critical in crustal failure during fracture-mediated intrusion.

#### Uplift-related extension

Unconsidered previously for intrusion is that ϵ_{V} affecting layers of a spherical planet requires proportional coeval horizontal strain, ϵ_{H} (Price 1959). Uplift thus induces a tensile stress −σ_{H}^{(up)}=*E*d*L*/*L*, where d*L*/*L* is the change in length of a chord perpendicular to the uplift axis. This −σ_{H}^{(up)} becomes mechanically significant for uplifts over wide areas above low-aspect source regions. For example, for plane strain, 250 m of uplift over a 200 km-wide region of lower crust with *E*=10^{11} Pa generates −σ_{H}^{(up)}=8 MPa. Although modest (not shown on Fig. 3 for clarity) −σ_{H}^{(up)} differs importantly from thermal dilation, in that it would permanently transform initially lithostatic conditions to marginally σ_{V}-dominant. Moreover, relaxation of −σ_{H}^{(up)} would create lateral space for cracks. In the above example, 6 m of elongation would eliminate the residual stress in one horizontal direction.

## Source instability

Source instability, i.e. conduit formation and incipient magma ascent, is postulated to require (exclusively) loads and deformation effects intrinsic to the generation or presence of magma. These are considered for our idealized source.

### Magmatic driving forces

Driving forces intrinsic to volatile-undersaturated anatectic magma contribute to the excess magma pressure Δ*P*_{E} (=*P*_{M(z)}−*P*_{(z)}) and derive from the positive Δ*V*_{M} of melting. The buoyancy and hydraulic contributions to Δ*P*_{E} are written as Δ*P*_{B} and Δ*P*_{V}, respectively, where Δ*P*_{E}=Δ*P*_{B}+Δ*P*_{V}. For magma in deformable rock of greater density, ρ, then Δ*P*_{B}=*gh*Δρ, where *g* is acceleration due to gravity, *h* is magma body height (*h*=[*S*−*z*] if the base is at depth *S*) and Δρ is the height-integrated density contrast ρ_{R}−ρ_{M} (e.g. Johnson & Pollard 1973; Gudmundsson 1988). For positive Δρ, Δ*P*_{B} increases with *h*, so the maximum Δ*P*_{B}^{max} is at the body top. For sources a few kilometres thick, Δ*P*_{B}^{max} is thus less than about 15 MPa. The Δ*P*_{B(z)} curve in Figure 3 is for a magma column rooted at *z*=40 km with ρ_{M}=2300 kg m^{−3}. Implications are insensitive to small differences in ρ_{M}.

Volume created by dilative melting is termed ‘excess magma volume’, *EMV*. If *EMV* arises too abruptly for full relaxation by inelastic deformation, the non-relaxed portion, written as *EMV**, will generate hydraulic overpressure Δ*P*_{V} by elastic compression of wall rock. For a low-aspect ratio source in lithostatic crust, relaxation of *EMV* would involve overwhelmingly ϵ_{V}, the uplift depending on Δ*V*_{M}, mean melt proportion *M* and *h*, according to ϵ_{V} = Δ*V*_{M}M̄ *h*. If the crust cannot deform on an appropriately short timescale, such that ϵ_{V} approaches zero, then Δ*P*_{V} would equal the nominal σ_{V} increase in roof rocks, the elastic maximum being *E*ϵ_{V} (Price 1959). For model crust, assuming widespread ϵ_{V} (average crustal Ē = 5 × 10^{10} Pa; Table 1), then Δ*P*_{V}^{max} is 100 MPa. For very abrupt loading, however, ϵ_{V} would be absorbed within a narrower elastic zone. Assuming *E*=10^{11} Pa (lower crust) then for a 1 km-zone the increase in σ_{V} is 7500 MPa. Such an unbalanced load would, of course, never arise before deformation began. Thus: (i) melting will always provoke inelastic uplift, even if strong, abrupt deformation were initially required, except if *EMV** was accommodated another way; and (ii) if part of the potential σ increase associated with dilation evades relaxation then Δ*P*_{V}>0.

### Contrasts in loading style and rate between ΔP_{B} and ΔP_{V}

Δ*P*_{B}^{max} is at the top of any magma column, e.g. a hydraulic nucleus, so for a source developing by isotherm ascent dΔ*P*_{B}^{max}/d*t*∝d*h*/d*t*, where under stable (non-fractured) conditions d*h*/d*t* equates to the vertical progress of the isotherm for *M*^{PT}. Strains induced by Δ*P*_{B} do not relax Δ*P*_{B} except if these reduce the density gradient responsible. In contrast, owing to the hydraulic principle that a surface force is isotropic in a static fluid (e.g. Shames 1982) then for such an idealized condition Δ*P*_{V} is felt uniformly through all magma contiguous with melting rock. Even accounting for viscous dissipation by magma flow, this characteristic gives Δ*P*_{V} great mechanical efficiency over Δ*P*_{B} as an agent of deformation; pressurizing all magma effectively instantaneously up to the elastic strength of the hydraulic nucleus envelope.

The value of Δ*P*_{V} attained during dilative melting might be assessed from Δ*P*_{V}^{max} for confined dilation and the proportion relaxed, written ϕ, according to Δ*P*_{V}=(1−ϕ)Δ*P*_{V}^{max}. However, constraining ϕ is difficult, depending both on Δ*P*_{V}^{max} and the strain-rate dependence of envelope strength. Nevertheless, dΔ*P*_{V}^{max}/d*t* for a developing source will always exceed dΔ*P*_{B}^{max}/d*t* because it depends on the volumetric rather than linear melting rate. For a source with *h*=2500 m and Δρ=500 kg m^{−3} formed in 8×10^{5} years, then dΔ*P*_{B}^{max}/d*t* is 16 Pa year^{−1}. If confined, with Δ*V*_{M}=0.1 and *M*=0.3, then if the whole overburden acted elastically dΔ*P*_{V}^{max}/d*t*=125 Pa year^{−1}. If only a 1 km overburden layer absorbed this load then dΔ*P*_{V}^{max}/d*t* would be 9.4 kPa year^{−1}. Instantaneous relaxation of such Δ*P*_{V} would require strain at 10^{−6} or 10^{−4} s^{−1}, respectively. Thus, first, while perfect vertical confinement of a dilating source is unlikely, the instantaneous onset, fast augmentation and potential magnitude of Δ*P*_{V} make full relaxation by source deformation and uplift unfeasible (ϕ<1), at least temporarily, so Δ*P*_{V} will always augment *P*_{M}, at least initially. Second, as crustal materials are elastic at strain rates of 10^{−6} s^{−1} (Pfiffner & Ramsay 1982), Δ*P*_{V} will, at least temporarily, induce an elastic response from the rock envelope (Rushmer 1995).

### Magma deployment – vein system

We may now return to the nature of the hydraulic nucleus. For closed-system melting at high degrees of melting (*M*>*M*^{RT}) this could be an anatectic core (Fig. 2b). However, at lower melt fractions (0<*M*<*M*^{RT}) the linking of smaller cavities, i.e. porosity and/or structural cavities, would suffice (Turcotte 1987; Sleep 1988; Petford 1995). Deformation requires deviatoric stress, so the latter would be elongate or planar (Simakin & Talbot 2001). The failure mode producing such structural cavities depends on σ_{d} in the host rock: if σ_{d} were initially large, deviatoric stress-exchange effects and modes other than pure tensile would initially occur to reduce this, resulting in various types of shear- and transitional-tensile shear-failure (Petford & Koenders 1998). Once σ_{d} declined sufficiently, tensile rupture would be favoured to produce mode 1 fractures. Such tensile cracks formed *in situ* would be low-pressure sites into which melt would infiltrate by porous flow, i.e. veins (Sleep 1988; Stevenson 1989; Petford 1995).

The pressure gradient driving melt segregation into veins partly results from buoyancy, compacting residual matrix (Sleep 1988; Scott & Stevenson 1986), and partly from hydraulic loading by dilative melting. The influence of grain size, porosity, melt density, and melt and matrix viscosity on segregation rate are considered by Petford (1995), who has shown that for amphibolite melting under conditions similar to those considered here, rates of porous melt flow into cracks permit segregation on an appropriate timescale. The melt fraction required for segregation is also considered by Petford (1995) and Vigneresse & Tikoff (1999), who have shown that the critical melt fraction concept for melt mobilization, suggested by Wickham (1987) as *M*>0.35, is illusory. Connectivity of porosity to permit melt migration occurs at some much lower *M*=*M*^{PT} (Vigneresse *et al.* 1996).

Within the porous zone (0<*M*<*M*^{PT}) of the source (Fig. 2a), melt could infiltrate into veins only from intersecting pores. However, in the permeable zone (*M*^{PT}<*M*<*M*^{RT}) veins would drain magma from the contiguous porosity reservoir given an adequate pressure gradient. Notably, therefore, as the development of a permeable zone and vein formation within it would occur at lower *M* than required for an anatectic core (Fig. 2a), and since melt segregation into veins would represent open-system behaviour with respect to the original porosity, the formation of a hydraulic nucleus comprising one or more sets of veins will probably always predate, and thus preclude, anatectic core formation. This agrees with observations of magma deployment in formerly partially molten lower crustal rocks, which vary from weakly porous protolith cut by vein networks (Daczko *et al.* 2001; Klepeis *et al.* 2003) to gneissic rock conspicuously segregated into veins and protolithic schlieren (Sawyer 1991, 1994). The magmatic loading effects required to generate such vein networks in the absence of regional stresses are as follows.

### Magma-intrinsic deformation effects meeting source rupture criteria

A vein system provoking source instability, i.e. magma ascent, must comprise at least one set of subvertical cracks. Three criteria for vertical tensile rupture are: (i) σ_{3}=τ_{H}; (ii) σ_{d} at failure=0.3–1τ_{H}; and (iii) σ_{3}=σ_{H} (Anderson 1951; Roberts 1970). These are met by a combination of magma-intrinsic deformation effects: (a) development of an effective tensile stress field owing to magma pore pressure; and (b) reduction in σ_{d} owing to stress exchange between principal stress directions through the wedging effect of overpressured magma in planar cavities. Tensile loading may obtain under extreme extension but σ_{H} at depth is generally compressive. Mode 1 hydraulic fracture thus requires an effective tensile stress field induced by the pressure of magma in cavities (Fig. 4) analogous to how pore fluid reduces the normal stress (σ_{N}) on potential failure surfaces to permit tensile jointing (Hubbert & Rubey 1959; Secor 1965). Whereas for small, equant pores in brittle rock the effect of pore pressure is isotropic, however, the primary porosity in melting rock might initially be elongate or planar. The source matrix may also deform inelastically, permitting deviatoric stress changes prior to tensile rupture.

If σ_{d} were initially too high to permit tensile failure then deviatoric stress changes must first occur to reduce this. Without remote loading, this is explained by the orientation-dependency of stress exchange between pressurized cracks and their surroundings – the magma wedging effect (Jaeger & Cook 1976). Tensile cracks dilate selectively in the direction σ_{3} (Anderson 1951), the effect of their internal pressure (*P*_{M}) being to augment σ_{3} towards *P*_{M}, while the ambient contribution to *P*_{M} is σ_{3} (Parsons & Thompson 1991). Thus, the effect of an internal excess magma pressure Δ*P*_{E} (=*P*_{M}−σ_{3}) in a crack is first to raise σ_{3}→σ_{2} until σ_{2}≈σ_{3}, when dilation would be favoured normal to both minor principal stress directions (Fig. 2). If two orthogonal sets of cracks were present, subsequent dilation of both would cause σ_{2} and σ_{3}, now coupled, to augment together incrementally (Vigneresse *et al.* 1999). Parenthetically, this helps meet criterion (ii) given earlier, facilitating tensile rupture in the second orthogonal plane. Further dilation of both sets of cracks increases (σ_{2}≈σ_{3})→σ_{1} until σ_{1}≈σ_{2}≈σ_{3}. Under such a condition, rupture or dilation is essentially equally favoured in all three planes and, if not already present, formation of a third orthogonal set of cracks becomes favoured. In practice, if three orthogonal sets of cracks are present, each would tend to dilate alternately if overpressured, augmenting all principal stresses, until their degree of interconnection allowed the whole network to inflate isotropically. The order and significance of stress-axis switching would depend on the original stress field and boundary conditions. If criterion (ii) given earlier were not initially satisfied, then dilation of primary porosity would first induce chaotic, plastic, shear or transitional tensile deformation of the matrix more immediately favoured to reduce σ_{d}. Only after σ_{d} were reduced could organized tensile rupture and vein inflation begin.

In a porous melting zone, dilation of overpressured cracks would, if *EMV** were available to diminish σ_{d} sufficiently, result in the development of orthogonal vein sets and a pressure spike around the source. This near-isotropic stress state elevated with respect to the starting condition is designated ‘elevated – pseudo-lithostatic’. Note that the pore pressure value *P*_{M}, portrayed in Figure 4 as equal to the initial lithostat, would thus in practice evolve from an initial state (*P*_{M}=σ_{3}<σ_{1}) towards a magmatically elevated pseudo-lithostatic state (*P*_{M}=σ_{3}≈σ_{1}). However, as the reduction in σ_{N} by pressurized porosity would always reduce the effective σ_{3} back towards 0 because *P*_{M}=σ_{3}+Δ*P*_{E}, the complex evolution of the pore pressure and pressure spike is omitted from Figure 4 for clarity (but see Fig. 6 later).

Whereas σ_{d} must normally diminish to yield tensile rupture conditions, it must nevertheless be non-zero and approximate the local tensile strength, τ. Some process generating σ_{d} must therefore operate if the source is to undergo further mode 1 rupture once homogeneous inflation has begun. This is satisfied by another feature of dilative melting, which is that a low-aspect source, even homogeneously inflated, will behave essentially as a large horizontal crack. Inflation would induce ϵ_{V} (uplift) with ongoing diminution in σ_{H} owing to uplift-related horizontal stretching, and/or a temporary dynamic augmentation in σ_{V}. This provides another way to achieve criterion (iii), given earlier, by non-relaxed source dilation causing a dynamic rise in σ_{V} in rocks above the source such that σ_{V3}→σ_{V1}. These effects would maintain σ_{d} as non-zero and σ_{V}=σ_{1}, as illustrated graphically in the following subsection.

### Source rupture

Rupture scenarios are distinguished by the pre-existing stress field: (a) lithostatic or σ_{H}-dominant; and (b) σ_{V}-dominant. Triaxial rupture is discussed separately. Two type (a) cases are shown in Figures 5 and 6, where thermal stresses σ_{H}^{th} and uplift-related stretching−σ_{H}^{up} are superimposed in the sequence heating–uplift, together with the application of Δ*P*_{V} from dilative melting. In Figure 5, the initial lithostat σ_{L} [1] is the model σ_{V(z)} curve from Figure 3. The addition of a modest horizontal compressive thermal stress, σ_{H}^{(th)}, gives [2] but cooling returns this to the lithostat, showing that criterion (iii) cannot be met by thermal stress relaxation. In contrast, a modest uplift-related tensile stress contribution, −σ_{H}^{up}, reduces σ_{H} to σ_{H3}^{(−up)}, which for clarity is shown only in the accompanying Mohr plot (Fig. 5b) as [3]. Thus, if −σ_{H}^{(up)} exceeded the local tensile strength of the source rock, as is likely for the widespread uplifts considered here, then after relaxation of any σ_{H}^{(th)}, σ_{H3} becomes σ_{H3}^{(−up)} [4]. This stress distribution satisfies all criteria for vertical rupture once the effect of magma pore pressure is accounted for [5].

Figure 6 illustrates partly non-relaxed source inflation with the same initial lithostat [1] and σ_{H}^{(th)} neglected for clarity. Melting generates excess magma volume *EMV* too rapidly for complete relaxation (ϕ<1) so *EMV** generates Δ*P*_{V} in the non-fractured source porosity. This raises σ_{V} and σ_{H}, coupled owing to magma wedging, to a magmatically elevated pseudo-lithostatic state shown by the arbitrary thick grey line σ_{L}^{M} [2]. Note that the degree of magmatic stress augmentation in the pressure spike is exaggerated for clarity. If uplift were to accompany this internal overpressure development, then −σ_{H}^{(up)} would subtract from σ_{L}^{M} to give σ_{H3}^{(M−up)} [3], with σ_{V} becoming σ_{V1}^{(M)}. Conditions for vertical rupture are met at [4].

Figure 7 shows the case where the starting stress state is σ_{V}-dominant. The curve σ_{H}^{(bo)} from Figure 3, again, is only illustrative. It might represent residual stresses following tectonic extension, or derive from crustal subsidence into the source (see below). Thermal stresses are neglected. The initial σ_{V} and σ_{H} are shown as [1] and [2] (Fig. 7). Non-relaxed magma dilation in pores takes up the horizontal slack within the source, increasing σ_{H} to a magmatically elevated state σ_{H}^{(bo+M)} [3], which approaches σ_{V}. At the source top, criterion (i) is met by default owing to melting, so if σ_{H}^{(bo+M)} exceeds τ_{H}, criterion (ii) is satisfied for rupture. Vertical cracks form owing to σ_{V1}. Prograde failure as σ_{H3} approaches σ_{V1} occurs at [4], but note that cracking would first initiate in the source interior where melting begins. Note that retrograde rupture is also conceivable, where σ_{H} declines from an initially lithostatic or elevated pseudo-lithostatic state owing to lateral crustal tension.

Excluding ad hoc changes in remote loading, therefore, the simplest and most likely trigger for source rupture is uplift-related tension, with or without a dynamic elastic σ_{V} increase. Preferred crack geometry depends on the stress field as failure is approached. Rupture under σ_{V}-dominance requires horizontal inflation for prograde vertical cracking or uplift for retrograde rupture. For an initially compressive confining stress, horizontal ruptures form first. Dilation of these allows vertical inflation, equalization of σ_{H} and σ_{V}, and a dynamically elevated and low σ_{d} stress state. If the vertical load is initially σ_{V1}, the first cracks are vertical and so the source is inherently unstable; inflation could continue only if magma cannot discharge faster than *EMV** develops. Although having implications for extensional and compressive settings, we continue to analyse the idealized intraplate case.

### Three-dimensional source rupture geometry

The internal rupture pattern will control magma dynamics during segregation, dyke propagation and source drainage. For a region of symmetric doming, triaxial notation (σ_{1}≠σ_{2}≠σ_{3}) is needed to describe the stress field, with principal stresses configured as in Figure 8. Here, σ_{HR} and σ_{HC} are, respectively, the radial and circumferential horizontal axes of principal stress. σ_{HR} is the horizontal confining stress, while σ_{HC} varies only by Poisson's effect or by magma wedging, so their magnitudes are semi-independent. Under initially weakly σ_{H}-dominant conditions, as for Figure 5, triaxial stress fields with σ_{V}=σ_{3} favour horizontal cracking and are shown in Figure 8 as (a) and (b). Of these, (a) is relevant to a dilating source where σ_{HR}>σ_{HC} owing to lateral confinement. From these states, a reduction in one of the horizontal stresses below σ_{V}, or an increase in σ_{V} above the lesser σ_{H} owing to vertical source inflation, would cause (a) and (b) to reconfigure, respectively, to (c) and (d). The new σ_{3}, either σ_{HC3} (c) or σ_{HR3} (d), would induce a new preferred crack geometry (Roberts 1970; Gudmundsson *et al.* 1997); vertical, radial and pure tensile for (c), and transitional-tensile ring fractures for (d). After such a switch, however, failure would still require τ, the tensile rock strength, to be exceeded. Thus, if σ_{HR} and σ_{HC} were initially comparable, σ_{V2} might change to σ_{V1} while σ_{H1}→σ_{V2}, giving a second switch from (c) to (e) or (d) to (f). Under starting conditions where σ_{V}=σ_{1} the initial preferred crack geometry is vertical. Depending on the confining stress and relative magnitudes of σ_{HR} and σ_{HC}, this corresponds either to radial fractures (c) or (e) or ring fractures (d) or (f); radial fractures again being favoured by lateral confinement. Once an initial fracture set had developed, either horizontal or radial, the wedging effect of magma within it would be to augment σ_{3}→σ_{2} and subsequently σ_{2}→σ_{1}, as described previously. Thus, ultimately, the fracture system within a domical inflating source would come to comprise a network of interconnecting flat-lying and radial veins. The development of circumferential fractures would be strictly limited by the elastic compressibility of the rock mass and radial confining stress.

Failure within a developing source would initiate where melting begins at the base and at the centre of the uplift responsible for source rupture where −σ_{H}^{up} is greatest. Numerical simulations of internally and vertically loaded cavities in elastic media support initial central failure as the horizontal tensile stresses concentrate over the centre of a swelling low-aspect body if deeply buried (Tsuchida *et al.* 1982; Martí *et al.* 1994; Gudmundsson *et al.* 1997, 1999). Δ*P*_{B}^{max} also occurs beneath the source apex, possibly exceeding τ for a tall source, also favouring central failure. Under σ_{H1}, peripheral tensile failure would not induce source instability because the magma could not ascend in resulting sills. Rather, these would locally increase σ_{V} by wedging, stabilizing the source for further inflation. Vertical radial fractures developed in a stress field caused by doming could initiate at the inflation axis, consistent with predictions concerning the initial failure locus. Radial cracks have been simulated in experimental doming (Komuro *et al.* 1984; Martí *et al.* 1994) and observed in natural cases of magmatic tumescence (Müller & Pollard 1977). Indeed, lithospheric-scale fractures at Yellowstone form part of a huge radial array (Glen & Ponce 2002).

Tensile cracks in a flat extending layer are spaced proportional to its thickness, the total strain determining average crack width (Bai *et al.* 2000). Although such relations are more complex for radial geometries, a useful point emerges for source regions. Progressive melting will form cracks whose spacing and width increase as *h*, many thin veins localizing into few cracks at the top. Constant spacing and width would indicate instantaneous rupture.

Petford *et al.* (1994) estimated that dykes approximately 6 m wide would allow felsic magma to traverse thick crust. In contrast, Rubin (1995), applying the thermal entry length concept (Delaney & Pollard 1982), asserted that elastic felsic dykes could never achieve criticality. In an infinite elastic medium, the crack width depends on Δ*P*_{E} and *E* (Pollard 1987). However, multiple cracks interact, while tensile stresses will relax on rupture, providing additional space. Uplift, for example, generates tens of metres of horizontal extension across a wide source region. Predicted radial cracks meeting centrally will form a nexus, focusing space gained across such a region at a natural conduit, as illustrated in Figure 9. Radial dykes would also remain vertical across a domical uplift along the slightly dipping stress trajectories associated with lateral gradients in ϵ_{V}.

## Intrusive architecture

The key mechanical problem in fracture-mediated intrusion is the generation of a vertically extensive tensile fracture system and its transition to horizontal propagation at shallow depth (Vigneresse *et al.* 1999). To understand how this might occur, account must be taken of stress criteria (i)–(iii) given earlier, which continue to apply above the source. Of these, criterion (iii) is most significant; the wedging effect of overpressured magma in vertical cracks will augment σ_{H} in wall rock until the stress field favours horizontal fracture, and dykes must reorient to sills (Parsons & Thomson 1991; Vigneresse *et al.* 1999). This magma-wedging effect provides the simplest explanation for dyke–sill and sill–dyke transitions (Bradley 1965; Roberts 1970; Johnson & Pollard 1973), which must be understood if we are to explain why a transition from dyke to sill propagation does not prevent a magma fracture system from extending to shallow depth, but does occur once it reaches such a depth to generate the tabular plutons observed there. Other ideas for how dykes transform to sills have invoked neutral buoyancy, i.e. magma–rock density contrast (Gilbert 1877; Pollard & Muller 1976; Ryan 1987; Corry 1988; Lister 1991), or dyke-stopping by stress guides (Weertman 1980; Gudmundsson *et al.* 1999; Hogan & Gilbert 1995; Hogan *et al.* 1998). However, neutral buoyancy is inconsistent with density relations between felsic magmas and crustal lithologies, gravity anomalies over felsic plutons (Vigneresse & Clemens 2000) and experimental evidence (Takada 1989), while our assumption of idealized crust negates dyke-stopping by ad hoc features of regional stress and structure. Here, therefore, we adopt the simplest criteria for sill intrusion, that: (1) σ_{V}=σ_{3} (Anderson 1951; Roberts 1970); and (2) that the crack contents have sufficient overpressure to rupture and displace the rock envelope, stated by Bradley (1965) as *P*_{M}=σ_{V}+τ_{V}, where τ_{V} is the vertical tensile strength. To examine the consequences of such constraints we analyse how the intrusive stress field varies with depth above an inflating source and under what conditions such criteria would allow a vertically extensive crack system to develop and the transition to shallow sill emplacement.

### Depth–stress relations

Figure 10 shows depth–stress relations for idealized magma columns in model lithosphere under previous conditions (Figs 5, 6, 7). These are presented with three caveats. Crack propagation is effectively assumed for any Δ*P*_{E}, as justified by stress-intensification effects at crack tips (Pollard 1987; Emmerman & Marret 1990). Viscous pressure losses (Δ*P*_{losses}) are considered only qualitatively; the detailed crack geometry is unknown so the plots suggest preferred fracture geometry only after magma has come to rest and equilibrated with local stresses. Variations in magma properties owing to cooling are also neglected (but see Clemens & Mawer 1992, Petford *et al.* 1994).

#### Conditions lithostatic to σ_{H}-dominant

In Figure 10a, curve AB is the Δ*P*_{B(z)} curve for a model magma column (Fig. 3). Curve CD (Fig. 10a) is the initial lithostat (Fig. 5), while CEFD is the uplift-modified horizontal stress σ_{H3}^{(−up)}. Heavy dashed arrows (Fig. 10a) are hypothetical driving pressure gradients. Curve BD adds Δ*P*_{B(z)} to σ_{V(z)} and thus represents the pressure of a static magma column (*P*_{M(z)}) neglecting any Δ*P*_{V} and Δ*P*_{losses}. Such a column locally exceeds σ_{V(z)} and is thus unfeasible, as by magma wedging it would cause σ_{H} to become σ_{1} inducing horizontal cracking and sill injection just above the source. Incorporation of Δ*P*_{V} with (1) further augments *P*_{M(z)}, exaggerating this tendency. If Δ*P*_{V} is neglected, then Δ*P*_{losses} of 20–40 kPa m^{−1} give trajectory (2), taking *P*_{M(z)} below σ_{H3(z)}. A column with such a *P*_{M(z)} distribution could not maintain a vertical fracture open against σ_{H3(z)}, so to the left of CEFD remote extension is required to create lateral dyke space. Trajectory (3), following the reduced σ_{H3(z)}, is one of a set of viable gradients in the region CEFD, where *P*_{M(z)} balances σ_{H3(z)} without overcoming σ_{V1(z)}. However, rather precise Δ*P*_{losses} (*c*. 10 kPa m^{−1}) are needed. If Δ*P*_{V} and Δ*P*_{B} were also considered, then these additional forces must also be balanced by viscous pressure losses.

#### Conditions σ_{V}-dominant

Figure 10b is for σ_{V1} conditions (Fig. 7) caused by uplift or dynamic source inflation, or a horizontal stress drop (see below). Curves AB and CD again represent Δ*P*_{B(z)} and σ_{V(z)}, while EC is σ_{H3}^{(bo)} (Fig. 3). CFGH (Fig. 10b) is σ_{H3(z)} elevated by source dilation, such that rupture criteria are met. HB is the predicted *P*_{M(z)} for Δ*P*_{losses}=0 in lithostatic crust, derived by adding Δ*P*_{B(z)} to *P*_{M(z)}. Curve BJH adds Δ*P*_{B(z)} to the magmatically elevated σ_{H3}^{(bo+M)} distribution (CFGH) to give *P*_{M(z)} for a buoyant column (without Δ*P*_{V} or Δ*P*_{losses}). Ignoring thermal factors, a column on trajectory (2) could extend as one or more dykes to approximately 10 km depth, where Δ*P*_{B(z)} plus σ_{H3(z)} would exceed σ_{V(z)}. If Δ*P*_{losses} balanced Δ*P*_{B} exactly, the *P*_{M(z)} curve would be GFC (3), for which Δ*P*_{B(z)} crosses σ_{H1} at shallow depths, here approximately 7.5 km. The cross-over between Δ*P*_{B(z)} (minus τ) and σ_{V(z)} represents a critical depth *I*, at which static columns, i.e. those that stagnated, would be forced to intrude as sills by their buoyancy alone.

### Constraints on intrusive architecture

Together with the differential between σ_{V} and σ_{H}, magma wedging is the major constraint on crack-system development because the greater is *P*_{M(z)}, the more abruptly will dykes reconfigure the stress field, initiate failure normal to the new σ_{3}, and inject as a sill. Notably, for a deep-sourced, static magma column, even if Δ*P*_{V} is negligible, Δ*P*_{B} is sufficient on its own to induce this switch once it has reached shallow depth, which partly explains why sills are shallow. Paradoxically, viscous dissipation of Δ*P*_{E} favours dyke stability and magma transport to shallow levels. The main difficulty, however, is the wedging action of Δ*P*_{V}. Under most conditions this promotes σ_{H}-dominance just above the source, after a small vertical overshoot proportional to the rock strength.

One solution to this problem is that an overpressured magma column might develop as alternating vertical and horizontal cracks (Corry 1988; Vigneresse *et al.* 1999). Two other possibilities are identified. First, dynamic elevation of σ_{V} by source inflation and reduction of σ_{H} above a developing low aspect ratio source by uplift-related stretching would increase the σ_{d} gradient, driving magma upwards and stabilizing vertical fractures. Dynamic overburden support could only occur during source inflation, however, while its effects would be restricted to the pressure spike around the source. Second, if total crack volume came to equal *EMV**, then Δ*P*_{V} would be exhausted by crack-system enlargement, halting hydraulic crack growth. Without further *EMV**, and thus Δ*P*_{V}, the driving force for magma fracture would be limited to buoyancy. However, if melting continued after initial Δ*P*_{V} (stored as elastic energy in the source-region envelope) was exhausted, stable crack growth would be expected, with *EMV** accumulation matched volumetrically by crack growth. Given simple propagation criteria, the crack system would grow in its initial configuration until the build-up in Δ*P*_{B} with column height forced conversion to a sill at a shallow depth during a temporary period of stasis. This possibility is examined after considering the implications of sill emplacement for space relations.

## Magma emplacement

The scale of magma emplacement generating large intraplate intrusions is illustrated by two well-studied examples in Figure 11: the Lebowa Granite (LG), South Africa (Kleeman & Twist 1989; Ferré *et al.* 1999); and San Juan Batholith (SJB), USA (Plouff & Pakiser 1972; Lipman 1984). These batholithic bodies intruding thick, intraplate crust have low-aspect, sill-like cross-sections, lobate plan-forms and known or inferred central conduits, consistent with Figure 1. Their short emplacement histories agree with the proposed intrusion timescale. The LG is a Proterozoic alkali granite sill complex capping the Bushveld intrusion in 45 km-thick crust, which crops out over 66 000 km^{2} (Fig. 11a). It is a saucer shape, 1.5–3.5 km thick, 190 km in radius, dipping inward at less than 15 ° (Kleeman & Twist 1989; Ferré *et al.* 1999). The aspect ratio is 0.01 and it has a volume of approximately 120 000 km^{3}. Its 10 or so sheets were injected at approximately 5 km depth via vertical dykes over about 1.8 Ma (Walraven *et al.* 1990; Walraven & Hattingh 1993). The mildly alkaline SJB (Fig. 11b) intrudes 35 km-thick crust of the Colorado Plateau (Lipman 1984). Only apophyses are exposed, but a negative Bouguer anomaly of more than 15 000 km^{2} (Plouff & Pakiser 1972) defines three lobes approximately 100 km in radius. Gravity models suggest a domed roof, 7 km deep peripherally, less than 1 km centrally, giving a thickness of approximately 10 km, aspect ratio 0.03 and volume of approximately 30 000 km^{3}. After 1–2 km of regional doming (Atwood & Mather 1932) and central volcanism, the SJB took around 3 Ma to emplace and cool (Lipman 1984). Implications of these observations are considered below.

### Sill emplacement

Sill inflation involves vertical displacements – roof up, floor down. For an idealized circular sill with radius *a* initially small compared to its depth *D*, the simplest model is that of a mode 1 crack in an infinite medium under uniform internal Δ*P*_{E} (Gudmundsson 1990). For this case, expressions for the symmetric deflections (±*W*), radius and volume given by Sneddon (1951) are repeated in the Appendix. These are plotted in Figure 12, which shows, however, that large, tabular bodies like the LG and SJB are much thicker and less laterally extensive than elastic cracks of equivalent volume predicted using these expressions. This is not surprising, as natural intrusions are lobate and intrude beneath a free surface. The additional thickness and volume of sills is thus attributable to overburden flexure and/or floor depression.

### Roof uplift

To model laccolithic uplift, sill roofs are treated as elastic plates (Pollard & Johnson 1973), in which case the key assumptions relate to elastic plate thickness and the role of shear stresses (Timoshenko & Goodier 1987). For sills with *a*/*D*≫1, a thin-plate model can be used, while for *a*/*D*≪1, a thick plate model is required (Gudmundsson 1990), in which significant vertical and shear stresses yield smaller normalized deflections. Expressions for these cases are given by Ugural (1981) and repeated in the Appendix. Notwithstanding the appropriateness of plate flexure models for thin sills, laccolithic doming by roof flexure cannot account for the additional volumes of giant sills, however, since those with significant *a*, exceeding approximately 30 km, would generate implausibly high (>5 km) domes, not observed on Earth (Lipman 1984). Indeed, for the LG (Fig. 11a) any uplift must have been reversed as it displays the saucer shape of some smaller elastic sills, while strata above the much thicker SJB were uplifted by only about 1.5 km (Atwood & Mather 1932). The additional thickness and volume of giant sill intrusions therefore results from floor depression.

### Floor depression

All sills experience minor elastic downward deflection of the floor by Δ*P*_{E} (see above) and the weight of the magma, although these are negligible (Gudmundsson 1990), implying another contribution to floor depression, suggested by Lipman (1984) as ‘gravitationally-driven down-warping of rocks at lower crustal depths’. This idea has been proposed in various forms (Branch 1967; Bridgewater *et al.* 1974; Whitney & Stormer 1986; Cruden & McCaffrey 2001). In the model preferred here (Cruden 1998, 2006) it involves mechanical decoupling of sill underburden from sill overburden, and subsidence of the former into the source region by ductile shear of plastic lower crust and passive subsidence of brittle upper crust.

### Buoyancy pumping

Subsidence of the crustal region between source roof and sill floor (sill underburden) has numerous effects. It creates sill volume by floor depression, influencing its cross-sectional and plan geometry. By reducing direct support for the sill roof it suppresses laccolithic uplift. Foundering also provides a powerful mechanism for expelling source contents owing to the buoyancy of the latter and would simultaneously process non-melted protolith down through the melting zone. Subsidence would also tend fully to drain and compact source regions, eliminating relicts of their evolution. Finally, for a domical source region, subsiding sill underburden would experience a component of horizontal extension and reduction in σ_{H}, with implications for magma transport in dykes traversing the subsiding region.

Ductile subsidence thus explains shallow intrusion volumes and provides a major drive for draining source regions. The onset of subsidence defines a relaxation time for the crust determined by the geometry and properties of the subsiding region. However, since the source must already support the sill overburden hydrostatically, for magma within it to feel the incumbent load, underburden relaxation would probably always predate shallow sill intrusion. The effective viscosity of subsiding crust exceeds that of magma in the system by many orders. Thus, the former would always control the subsidence rate. If a conduit was kept open between source and pluton, magma would be expelled passively upwards by buoyancy. This process is designated *buoyancy pumping*. If, for some reason, the conduit were shut, the negative buoyancy of the subsiding underburden would be transmitted to magma within the compacting source as a hydraulic load, which would pressurize all magma within the isolated source and closed conduit system, tending to force the latter to reopen. However, the σ_{H}-drop affecting crust subsiding into a domical source would provide sufficient additional lateral space for conduits within the underburden as to make it unlikely that conduits would close during buoyancy pumping.

The subsided volume depends on source and sill geometry. It may be evaluated approximately using the expression for simple shear , where *Y* is the maximum vertical displacement and *a* the radius of subsidence. For *Y* = 3 km and *a*=60 km then ϵ is 0.025. Although strict knowledge of parameters is lacking for specific cases, it is noteworthy that at a typical ductile strain rate of 10^{−14} s^{−1} (Pfiffner & Ramsay 1982) such strain could occur in approximately 80 000 years. In this example, where the subsided volume is some 30 000 km^{3}, equivalent to the SJB, the downward volumetric underburden flux *Q*_{U} would be 0.4 km^{3} year^{−1} at 10^{−14} s^{−1} or 4000 km^{3} year^{−1} at 10^{−10} s^{−1}, consistent with estimates by Cruden (1998, 2006). Underburden subsidence at ductile strain rates can thus readily account for the short emplacement timescales of large intrusions.

A geometric consequence of subsidence into domical sources is that the greatest vertical movement is central, so the subsiding region undergoes subhorizontal elongation. This effect is exaggerated by the weakest underburden around the hot conduit relaxing first, consistent with funnel-shaped pluton feeders (Vigneresse 1995*a*) and centralized down-sagging of pluton floors (Bridgewater *et al.* 1974). Stretching of subsiding crust, in lieu of any tensile failure, thus yields a temporary reduction in the horizontal compressive elastic stress. This can be estimated for 2D simple shear by considering a region initially 100 km wide subsiding 3 km at one end, generating layer-parallel stretching ϵ_{H}=0.00045. Ignoring the slight non-coaxiality of this with respect to the original stress field, this gives a tensile subsidence-related elastic stress contribution −σ_{H}^{(subs)} of 22.5 MPa for *E*=5×10^{10} Pa, which would stabilize vertical cracks and/or create additional lateral space for dykes. The reduction of σ_{H} below lithostatic explains the significance of a major σ_{H}-drop for hypothetical magma ascent trajectories (Fig. 10b). In the example above, ϵ_{H}^{subs}=0.00045 equates to 45 m, so for our model case, subsidence provides 10–20 times more lateral space than the uplift responsible for rupture. Given efficient focusing of this space across the subsided region into existing conduits, sustained magma flow during buoyancy pumping is assured.

## System volume

The intrusive ‘space problem’ has long been recognized (e.g. Brown 1994). However, crustal decoupling at the intrusion level and subsidence into the source explain how, simultaneously, intrusion volume is created by floor depression, magma is transferred under buoyancy and the source is compacted (Cruden 1998, 2006). The key outstanding problem remains to explain how the fracture conduit system linking source to shallow intrusion initially develops given that Δ*P*_{V}, although critical for source rupture, also promotes reorientation of dykes to sills at some small height above source. Central to resolving this is that Δ*P*_{V}, as a hydraulic load, will decline to zero if the instantaneous *EMV** creating it declines to zero. This might occur if melting ceases (*EMV*→0) or if the cavity system enlarges volumetrically by non-recoverable strain by an amount equivalent to *EMV**. Under such circumstances Δ*P*_{V}→0, so the problem of magma wedging by dykes vanishes. We designate any intrusive system for which Δ*P*_{V}>0 as being in a state of *hydraulic inflation*.

### Hydraulic inflation

To consider how a vertically extensive conduit develops we disregard cessation of melting as an ad hoc function and because generation of the initial vertical crack system requires source rupture and ongoing dilative melting. The *EMV* of a developing source is accommodated by cavity system enlargement in three ways: (i) elastic compression of wall rocks, generating Δ*P*_{V} as stored elastic strain, the non-relaxed portion being *EMV**; (ii) inelastic source dilation, which for a low-aspect source involves essentially uplift; and (iii) volumetric growth of the crack system outside the source as dykes and sills. Of these sinks for *EMV*, we may write that *EMV**=(1 − [ϕ+η])×*EMV*, where ϕ and η are the proportions of the instantaneous *EMV* relaxed by (ii) and (iii), respectively. Hydraulic inflation may be understood by considering relations between ϕ, η and *EMV**, partitioning of which depends on their efficiency at converting *EMV* into work. Since each relaxes the same initial load, they may be qualitatively compared non-dimensionally. Figure 13 illustrates how, for an idealized magma system with constant volumetric melting rate, this partitioning of *EMV* relaxation defines four hydraulic subregimes.

#### Disequilibrium source dilation

At the onset of melting, prior to rupture, there is no crack system, thus η=0. *EMV* partitions between a non-relaxed portion (*EMV**) generating hydraulic overpressure (Δ*P*_{V}) of pore-bound magma and a portion relaxed by inelastic uplift ϕ (Fig. 13). Given that the onset of inelastic uplift requires a relaxation time for host crust, at the onset of melting ϕ=0. The *disequilibrium dilation* subregime occurs during this period of crustal relaxation and is thus characterized by an initially large *EMV** that declines abruptly as uplift (ϕ) begins. Once uplift has begun, the partition of uplift ϕ abruptly increases, while *EMV** concomitantly decreases owing to the crust's free surface, weakness in flexure and susceptibility to inelastic ϵ_{V}. Deformation during disequilibrium dilation thus involves spatially abrupt variations in Δ*P*_{E} and temporally abrupt variations in failure mode, geometry and deformation rate within the source, acting to smooth initially strong internal variations in σ_{d} caused by an uneven distribution of pore overpressure.

#### Equilibrium source dilation

Because ϵ_{V} readily occurs, uplift is favoured over elastic compression of wall rock to increase source volume and relax *EMV*. Thus, for a simple melting function, source swelling evolves to a partitioning of ϕ and *EMV** in which ϕ strongly dominates (Fig. 13). This subregime is characterized mainly by inelastic deformation and is designated *equilibrium dilation*. Expected rates of dilatant loading suggest that some dynamic overpressure (1−ϕ) always escapes relaxation during melting, however, so some Δ*P*_{V} is always maintained during this subregime. For a porous source in which *M*>*M*^{PT}, the hydraulic nucleus would therefore be maintained hydraulically inflated. This state would be near-isotropic except for a weak pressure gradient away from the most recently melted region near the source roof.

#### Disequilibrium cracking

Before rupture *EMV* is relaxed by uplift (ϕ), except for a residual portion *EMV** generating Δ*P*_{V}. Uplift induces a horizontal tensile stress contribution −σ_{H}^{(up)} that prepares the crust for vertical rupture initiated in the porous source. Crack growth, once initiated, then offers a new sink for *EMV* (η), defining a second transitional subregime *disequilibrium cracking*. Initially, the crack volume is small and limited to veins. This relieves little *EMV*, however, which continues to partition into ϕ and *EMV** (Fig. 13). However, the crack system may readily enlarge by extending outside the source through the conversion of stored Δ*P*_{V} into crack volume, as melt segregates into veins and these propagate as dykes. Under most conditions *EMV* is minimized much more efficiently by crack growth than by inelastic uplift, because the driving pressures for crack propagation are small owing to stress-concentration effects. The effective tensile stress −σ_{T}* at a crack tip is (Inglis 1913) −σ_{T}*=−σ_{T} (1+2*a*/*b*), where −σ_{T} is the effective driving pressure (*P*_{M}−σ_{3}) and *a* and *b* are the crack semi-axes. As restated by Pollard (1987), for cracks where *a*/*b*=1000 (e.g. 1 cm wide, 10 m high) −σ_{T}* is 2001 times −σ_{T}. Thus, given a small drop in σ_{H} through uplift, positive Δ*P*_{V} and a source height supporting veins up to kilometric half-height, crack growth is strongly favoured as the sink for new *EMV*. Disequilibrium cracking thus involves abrupt enlargement of the crack system and augmentation of η with a corresponding decrease in ϕ and *EMV** (Fig. 13). This surge of crack growth, powered by residual Δ*P*_{V}, would end once it balanced the instantaneous uplift-modified stress field.

#### Equilibrium cracking

Whereas inelastic deformation is too slow to relax new *EMV* instantaneously, crack growth is fast enough to compete with source loading. Thus, although melting be ongoing, *EMV** (and Δ*P*_{V}) would decrease to incipient levels during disequilibrium crack growth, terminating that subregime. This explains how we envisage the stalling effect of Δ*P*_{V} to be overcome. First, during disequilibrium cracking some *EMV* remains non-relaxed as elastic strain energy in the source walls, permitting an initial surge of disequilibrium cracking. During this stage the wedging effect of Δ*P*_{V} may induce initially vertical cracks to reorient to horizontal. Given high (although decreasing) Δ*P*_{V}, dilation of any deep horizontal crack segment so formed would locally restore the vertical load to a value above the horizontal, so alternating dyke and sill segments could be generated until instantaneous Δ*P*_{V}→0. Notwithstanding this, once *EMV** and Δ*P*_{V} had decreased to incipient levels, each new increment of *EMV** would partition directly into new crack growth. Strong Δ*P*_{V} would not become available so repeated stress-field switching by magma wedging could not reoccur. Moreover, given that some inelastic uplift would continue during ongoing melting, causing the crust to stretch and σ_{H} to decrease, vertical crack growth would continue to be favoured. Thus, once *EMV** reduced to incipient levels, Δ*P*_{V} would serve to keep the source and crack system hydraulically inflated while each new increment of *EMV** was converted directly to vertical crack growth through deflection of dyke walls and crack lengthening, to maintain the elastic crack geometry. This subregime, designated *equilibrium cracking* (Fig. 13), would facilitate the initial establishment of a vertically extensive fracture system.

### Transition to buoyancy pumping

For likely ϕ values of 0.9–0.999 during uplift, the absolute volumetric equivalence of *EMV** is at most a few tens of km^{3}, even for potent sources. This total *EMV** could be exhausted by just a few dykes tall enough to reach emplacement levels. For example, five triangular radial dykes, length and height 37.5 km, width 5 m, would accommodate 17.5 km^{3}. Hydraulically inflated dyke columns must therefore be thin if they are to be vertically extensive. Sill emplacement makes even greater demands on *EMV**. Deep, elastic sills have volume relations as for dykes, making alternating dykes and sills volumetrically inefficient at generating tall conduits. For shallow sills, because of the free surface, roof deflection yields even larger increases in system volume (Fig. 12) making laccolithic doming an even greater sink for *EMV**. Surface eruption, should it occur, represents an effectively infinite *EMV** sink. It is recalled that there is a critical depth for sill intrusion where the buoyancy of a magma column crosses the regional σ_{V} curve (Fig. 10). If the fracture system extended to this level during equilibrium cracking, sill emplacement would begin if the column then stalled, causing sill intrusion at depth *I*. If the column propagated initially to the surface, eruption would cause the freezing of magma, stalling the column, and explaining why hydraulically inflated systems do not erupt as floods of lava. Once shallow sill intrusion had initiated it would continue by buoyancy pumping. This is because depth *I* also marks the maximum overburden thickness supportable by the static buoyancy of the magma column alone. Only sills at this or shallower depth would permit crustal decoupling and buoyancy pumping to initiate.

Cessation of melting has not yet been considered. Once the causative thermal anomaly has stabilized the upper limit of the source then new *EMV* production ceases. For small sources, those with weakly dilative melting reactions or those in weak crust where inelastic deformation was unusually effective, *EMV* production may end before shallow sill injection occurred, so a transition to buoyancy pumping would not occur. Cartoons in Figure 14 illustrate the hydraulic inflation (a–c) and buoyancy pumping (d) regimes for small and large *EMV* reservoirs. Figure 14a shows hydraulic inflation before source rupture. Uplift (ϕ) is occurring, causing horizontal stretching. Figure 14b shows the end of the equilibrium cracking regime for a small *EMV* reservoir, where *EMV* is exhausted during vertical crack growth but before the crack system has reached depth *I*, so buoyancy pumping cannot initiate. Figure 14c shows the later stages of hydraulic inflation for a large *EMV* reservoir where the vertical crack system has reached the critical depth. Sill development has begun by elastic cracking and laccolithic uplift. Once sill intrusion has initiated, it is inconceivable that if melting then ceased 90–99% of the magma generated would remain in the source. Figure 14d shows how, once the crustal underburden has relaxed, magma remaining in the source would ascend by subsidence-powered buoyancy pumping. This is because the source supports the sill overburden hydrostatically, but if the pluton was at or above depth *I* and a conduit was open to the source then its floor is unsupported owing to the weakness of fluid magma in the source and conduit. Crustal decoupling and buoyancy pumping would commence such that the downward underburden flux, *Q*_{U}, controlled the upward flux of magma in the conduit, *Q*_{C}. Relaxation of the horizontal tensile stress contribution from subsidence −σ_{H}^{(subs)} widens the conduit significantly from the aperture generated during hydraulic inflation and ensures sustained magma flow.

## Magma flux rates during buoyancy pumping

Given buoyancy pumping, the proposed non-viability of felsic dykes (Rubin 1995) needs reassessment. Lateral space for vertical dykes is not only a function of Δ*P*_{E} and the crust's initial stress state, but includes contributions at the onset of both hydraulic inflation and buoyancy pumping from lateral stretching due to either uplift or subsidence. The initial hydraulic fracture system must be thin and focused to reach shallow depth, but is significantly widened during buoyancy pumping. Stretching of subsiding sill-underburden provides about 10 times more lateral space for dykes than the uplift responsible for cracking; several tens of metres across an extensive source region. Moreover, such space contributions would widen existing cracks, notably at intersections. For a domical source, rupture gives stellate vertical cracks (Fig. 9), focusing space contributions at the central nexus. In other stress-fields, orthogonal crack intersections would act in this way. Magma flow would localize at such intersections to generate pipe-like conduits by thermo-mechanical erosion of wall rock or by solidifying and freezing of magma at lateral dyke extremities.

For model purposes, flow up a fracture nexus approximates pipe flow. Potential values of magma flux up the conduit *Q*_{C} may thus be compared with *Q*_{U} during buoyancy pumping according to *Q*_{C}=π*r*^{4}/8µ_{M} (d*P*_{M}/d*z*), where d*P*/d*z* is the driving pressure gradient and µ_{M} is the magma viscosity (Turcotte & Schubert 1982). For d*P*/d*z*=400 Pa m^{−1} and µ_{M}=10^{6} Pa s, pipes of radii *r*=10, 20, 50 and 100 m support fluxes of 0.01, 0.15, 6 and 100 km^{3} year^{−1} respectively. At these rates a 30 000 km^{3} intrusion would require 3 Ma, 200 ka, 5 ka, or 300 years to fill, respectively. Given conduit widening by tens of metres during underburden collapse, pipe-like conduits at crack intersections permit fluxes adequate under reasonable pressure gradients to fill giant plutons in 10^{6} years or so. Note that resistance to magma flow does not limit source drainage rate. Rather, *Q*_{U} controls the process, while *Q*_{C} must balance *Q*_{U} for continuity. Subsidence also dampens fluctuations in *Q*_{C} provoked by vertically varying conduit radius.

### Complex pluton structure

Many intrusions evidence pulsed construction (Pitcher 1979) with nesting of phases, lobes (LG, SJB) or sheeting (LG). Multiple phases related to protolithic inhomogeneity or variation in melt proportion would be stacked in the vein system and drained sequentially. Lobate or sheeted structure is explicable by evolving stress conditions at the conduit neck or by multiple transitions from hydraulic to buoyant regimes. As underburden subsidence processes crust into the melting zone, then if a thermal anomaly persisted after subsidence initiated, renewed *EMV* production would reinflate the system during buoyancy pumping, generating an unpredictable interplay between magma production and source drainage. Possibly most plutons are generated in this way.

### Caldera formation and ‘super-volcanism’

Felsic volcano evolution is understandable in the model context. Hydraulically inflated conduits may intersect the surface before sill emplacement. However, *EMV* exhaustion by eruption strictly limits such volcanism except for the highest-*EMV* systems, e.g. CO_{2}-powered kimberlites. For others, once the initial conduit stalls it will intrude a sill owing to buoyancy. For powerful systems (Fig. 14c) volcanism may also accompany this hydraulic laccolithic stage, but this would be weak, as intrusive apophyses take up the horizontal slack (Parsons & Thompson 1991) preventing normal faulting and inducing mode 1 roof failure (Martí *et al.* 1994). Caldera collapse at this stage is prevented by hydraulic roof support, restorative hydraulic loading by subsiding blocks, roof binding by frozen magma and structural arching set at the inflation highstand.

After inflation, caldera-forming eruption becomes possible. The pluton may be isolated from source if the conduit freezes or closes beneath the sill. Dynamic support is withdrawn from the roof, which is supported only by magma and structural arching. If the intrusion depressurizes owing to, for example, eruption caused by volatile exsolution-driven internal overpressure (Burnham 1972; Jaupart & Tait 1990) then wholesale roof collapse is inevitable once it begins to shear into the magma layer. The structural pattern for roof failure is set during inflation (cf. Martí *et al.* 1994). Radial roof extension creates stress conditions for steeply inward dipping ring faults, a problem in caldera mechanics (Gudmundsson *et al.* 1997), which rotate to subvertical and activate when the roof sags, preventing mechanical wedging. Source and laccolithic doming at superimposed wavelengths (c.f. Gudmundsson *et al.* 1999) is proposed to explain this. Collapsing roofs come to rest on sill floors, so sill thickness limits caldera depth. Mismatches in caldera areas with parent batholiths by 2–4 times (e.g. Fig. 11a) is fixed by roof strength or the volatile-saturated roof zone.

Our model explains why, except in caldera-formation, only modest overpressures (<15 MPa) drive felsic volcanism (Jaupart & Tait 1990): during disequilibrium cracking Δ*P*_{V} is masked by viscous losses; during equilibrium cracking Δ*P*_{V} is reduced to incipient levels by volumetric crack growth; during buoyancy pumping Δ*P*_{E} is moderated by floor subsidence; while after pluton disconnection Δ*P*_{E} is limited by the roof strength.

## Acknowledgments

Professor R.S.J. Sparks and Dr S. Cruden are thanked for comments on an earlier version of the manuscript. S. Cruden and Dr E. Cañón-Tapia gave useful reviews.

## Appendix

For a mode 1 crack in an infinite medium under uniform internal Δ*P*_{E} the symmetric deflections (±*W*) of the roof and floor are (Sneddon 1951): *W*=4Δ*P*_{E} (1−*v*^{2}) (π*E*)×(*a*^{2}−*r*^{2}), where *r* is the radial co-ordinate. At the centre (*r*=0), *W*^{max} is: 4Δ*P*_{E} (1−*v*^{2})*a*/(π*E*), while *V*=16(1 −*v*^{2})Δ*P*_{E} *a*^{3}/(3*E*) and *a*={3*EV*/(16 [1−*v*^{2}]Δ*P*_{E})}.

The deflection of a thin, circular, simply supported plate of thickness *D* is (Ugural 1981): *W*=Δ*P*_{E} *a*^{4}/64ς[(*r*^{4}/*a*^{4})−2(3+v*r*^{2}/1+v*a*^{2})+(5+*v*/1+*v*)], where the flexural rigidity is ς=*ED*^{3}/12(1−*v*^{2}). For uniform Δ*P*_{E} then *W*^{max}=Δ*P*_{E} *a*^{4}/64ς[(5+*v*)/(1+*v*)]. An effective elastic thickness is substitutable for ς (Pollard & Johnson 1973). For a thick plate, *W*^{max} is: Δ*P*_{E} *a*^{4}/64ς [(5+*v*/1+*v*)+4*D*^{2}/(1−*v*)*a*^{2}], where the shear stresses (last term) make *W*/*W*^{max} greater.

- © The Geological Society of London 2008