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₂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⁵ to 2.5×10⁵ Jkg⁻¹. 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₂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⁻² 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⁻² 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⁻² 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⁻² 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′/ρ′_Pc′, 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⁵ J kg⁻¹, ρ=2900 kg m⁻³, k=2.5 W m⁻¹ K⁻¹, and c=900 J kg⁻¹ K⁻¹. For solid protolith at T close to T_solidus, κ is approximately 8×10⁻⁷ m² s⁻¹, while in the melting zone, with average properties ρ′_P=2700 kg m⁻³ and k′=2.0 W m⁻¹ K⁻¹, then c′=2900 J kg⁻¹ and κ′=2.6×10⁻⁷ m² s⁻¹, 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⁻². 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⁻² 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²/κ, then heat focusing yields a modified timescale t′ for isotherm ascent through a melting region t′≈z²/κ′. For amphibolite this predicts a reduced isotherm advance rate of approximately 0.003 m year⁻¹. An idealized partial melting region would thus require approximately 8×10⁵ years to reach 2.5 km thickness. Given an appropriate thermal anomaly, the timescale for C-type source formation is thus about 10⁶ 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Melting controls on an idealized source internal structure. (a) Schematic cumulative melt productivity diagram for hypothetical protolith with a dominant hydrous mineral. M, melt fraction; T, temperature. (1) Solidus at temperature T = T₁; (2) percolation threshold (PT) here corresponding with the onset of the main melting pulse at T = T₂; (3) rheological transition (RT); (4) ‘pseudo-liquidus’ at T = T₃ representing the end of the main melting pulse with residual solid persisting; (5) liquidus at T = T₄. (b) Schematic section through an idealized C-type source showing thermo-rheological zones. Light shading, non-melted crust; dark shading, mantle; not shaded, crustal melting region. The base of the melting zone is the Moho, representing an abrupt increase in solidus temperature. An anatectic core is unlikely to develop because rupture would first permit segregation of melt from solid residue into veins. The hydraulic nucleus would then encompass the inner region of the source defined by PT. (c) The effect of pressurized magma in elongate cavities in the porous zone is to equalize the horizontal and vertical loads by magma wedging.

Figure 2a shows a hypothetical sigmoidal melt-productivity curve for a simple hydrous protolith. Melt forms initially at the solidus T₁, attains greatest productivity over the melting pulse T₂ to T₃, and ends production at the liquidus T₄. T₃ 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.

Crustal strength

Maximum differential stress (σ_d=[σ₁−σ₃]) 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⁻¹⁰ s⁻¹ (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⁻¹⁰ s⁻¹ 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Example depth–stress distributions for the model crustal system and representative long-term strength envelopes. Curve σ_V(z) gives the lithostat P=P_L(z). Curve σ_H1^(th) shows the small and temporary effect of the example thermal stresses (see the text). Curve σ_H3^(bo) is the instantaneous horizontal stress arising from the weight of the crust and is purely illustrative. Curve σ_H3^(bo+th) shows the additional effect of the example thermal stresses. Curve ΔP_B is the depth-integrated buoyancy overpressure of a column of magma of density 2300 kg m⁻³ in model crust (Table 1). Strength envelopes after Kohlstedt et al. (1995) are for a strain rate of 10⁻¹⁰ s⁻¹.

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 ρ_Rgz. 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.

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⁻³. 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_MM̄ 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¹⁰ 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¹¹ 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/dt∝dh/dt, where under stable (non-fractured) conditions dh/dt 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/dt for a developing source will always exceed dΔP_B^max/dt because it depends on the volumetric rather than linear melting rate. For a source with h=2500 m and Δρ=500 kg m⁻³ formed in 8×10⁵ years, then dΔP_B^max/dt is 16 Pa year⁻¹. If confined, with ΔV_M=0.1 and M=0.3, then if the whole overburden acted elastically dΔP_V^max/dt=125 Pa year⁻¹. If only a 1 km overburden layer absorbed this load then dΔP_V^max/dt would be 9.4 kPa year⁻¹. Instantaneous relaxation of such ΔP_V would require strain at 10⁻⁶ or 10⁻⁴ s⁻¹, 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⁻⁶ s⁻¹ (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) σ₃=τ_H; (ii) σ_d at failure=0.3–1τ_H; and (iii) σ₃=σ_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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Mohr diagram (shear stress, σ_S v. normal stress, σ_N) showing two stress fields (A, B) with, respectively, large and small differential stress (σ_d) to illustrate the effect of magma pore pressure. The effect is to reduce the normal stress, σ_N, on actual or potential fracture surfaces so that the effective stress, σ_N*, becomes σ_N−P_M. Stress field A cannot achieve tensile rupture once pore pressure P_M is accounted for, as σ_d is too great to allow intersection with the tensile part of the combined failure envelope (CFE). Stress field B can achieve tensile failure once σ₃ becomes equal to the tensile stress τ. Failure conditions are shown by the black dots. Note that the increase in principal stress values to generate a pseudo-lithostatic pressure spike around the source as a function of the magma pore overpressure has been omitted for clarity.

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 σ₃ (Anderson 1951), the effect of their internal pressure (P_M) being to augment σ₃ towards P_M, while the ambient contribution to P_M is σ₃ (Parsons & Thompson 1991). Thus, the effect of an internal excess magma pressure ΔP_E (=P_M−σ₃) in a crack is first to raise σ₃→σ₂ until σ₂≈σ₃, 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 σ₂ and σ₃, 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 (σ₂≈σ₃)→σ₁ until σ₁≈σ₂≈σ₃. 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=σ₃<σ₁) towards a magmatically elevated pseudo-lithostatic state (P_M=σ₃≈σ₁). However, as the reduction in σ_N by pressurized porosity would always reduce the effective σ₃ back towards 0 because P_M=σ₃+Δ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=σ₁, 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].

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Idealized rupture scenario for modestly σ_H-dominant conditions. (a) Depth–stress plot. (b) Mohr diagram. The plots show stresses at the vertical source axis. Models refer to an idealized, 5 km-thick, C-type source region in idealized lithosphere (Table 1). Mohr plots refer to stresses at the top of the porous zone.

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].

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

Idealized rupture scenario for modestly σ_H-dominant conditions with magmatically elevated lithostatic stress, modified by uplift. (a) Depth–stress plot illustrating rupture conditions after swelling and uplift. (b) Mohr diagram showing magmatic elevation of lithostat. (c) Mohr diagram showing rupture conditions. The model source is 5 km thick (Fig. 2).

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 7.

Idealized rupture scenario for σ_V-dominant conditions modified by source swelling. (a) Depth–stress plot illustrating prograde rupture by source swelling. (b) Mohr diagram.

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 (σ₁≠σ₂≠σ₃) 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=σ₃ 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 σ₃, 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=σ₁ 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 σ₃→σ₂ and subsequently σ₂→σ₁, 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 8.

Axisymmetric stress fields and preferred tensile crack-system geometries. Shading indicates magma-filled cracks. Stress fields (a) and (b) favour sills, (c) and (e) favour vertical radial dykes, and (d) and (f) favour ring dykes or cone sheets. Large arrows linking cases indicate principal stress-axis changes due to magma wedging (discussed in the text).

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 9.

Cartoon to illustrate preferred radial crack-system geometry propagated from a hydraulically inflated source region. Radial dykes meet at a nexus that efficiently focuses lateral dyke space across the projected area of the source into a central pipe-like conduit. The number and in-plane geometry of radial dykes is conjectural.

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).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 10.

Depth–stress plots illustrating factors influencing intrusive system architecture. (a) For rupture under modestly σ_H-dominant conditions. (b) Rupture under σ_V-dominant conditions accompanying a horizontal stress drop (see the text).

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 σ₃, 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.

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 12.

Maximum deflection W^max and radius, a, v. volume, V, for pressurized cracks in infinite host rock with v=0.2, E=5×10¹⁰ and ΔP_E=15 MPa. Thick crosses, parameter estimates for Lebowa Granite (LG) and San Juan Batholith (SJB). For crack radii and thicknesses equivalent to the Lebowa Granite or San Juan Batholith, W^max values are too low to explain them as single sills (although the approximate 100 m thickness of individual LG sills does coincide with predictions), while predicted radii for cracks with volumes of the LG and SJB are far greater than observed, unsurprising since the a/D ratios of these plutons violate elastic crack geometry owing to the free surface.

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⁻¹⁴ s⁻¹ (Pfiffner & Ramsay 1982) such strain could occur in approximately 80 000 years. In this example, where the subsided volume is some 30 000 km³, equivalent to the SJB, the downward volumetric underburden flux Q_U would be 0.4 km³ year⁻¹ at 10⁻¹⁴ s⁻¹ or 4000 km³ year⁻¹ at 10⁻¹⁰ s⁻¹, 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 1995a) 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¹⁰ 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.

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 13.

Partitioning of release mechanisms for potential energy represented by excess magma volume, EMV, between inelastic uplift, ϕ, elastic compression generating hydraulic overpressure, EMV*, and volumetric crack growth during hydraulic inflation, η. Inflation subregimes are labelled.

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³, 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³. 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 14.

Cartoons illustrating the proposed fracture-mediated intrusion model. (a) Pre-rupture evolution. (b) Ending of cracking for a source generating small EMV. I is the critical sill intrusion depth. (c) Cracking regimes for a source generating large EMV. (d) Buoyancy pumping regime. Note that lateral boundaries are vertical in a spherical Earth reference frame, to illustrate the mechanism by which uplift generates coeval horizontal stretching.

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₂-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.