## Abstract

Crustal magma transport is typically described using a complex, non-linear model associated with fluid-driven fracturing, and therefore fundamentally sound modelling forms the basis for interpretation of magmatic intrusions. One of the most basic considerations is that magma-driven sills can be broadly categorized based on the energy dissipation mechanism that is predominant during intrusion growth. In cases where either viscous flow or overcoming fracture toughness strongly dominates fracture behaviour, it is typical to speak of viscosity-dominated or toughness-dominated regimes, each of which defines a class of fracture propagation with significant implications for modelling. This paper presents a straightforward and geometry-independent means for local determination of the expected propagation regime based on an experimentally verified mathematical analysis of the multi-scale, coupled mechanics that govern the near-tip region. The propagation regime is then related directly to the ratio between a characteristic length associated with the near-tip physics compared with the size of the fracture/sill. Sill growth is shown to be expected in or near the viscosity-dominated regime and hence modelling generally must take into account the complexity of the near-tip region rather than relying solely on the tip behaviour implied by linear elastic fracture, although toughness-dominated mafic intrusions can also be anticipated if fracture toughness increases sufficiently rapidly with the intrusion size.

Crustal magma transport is typically modelled as a fluid-driven fracturing process (e.g. Spence *et al.* 1987; Lister 1990; Lister & Kerr 1991). These systems are complex owing to the presence of multiple, coupled physical processes such as rock fracturing and viscous fluid flow in the growing fracture. They are also non-linear owing both to a non-linear relationship between the fluid flow and the fracture opening and to the presence of moving boundaries associated with the fluid front and propagating fracture tip. Because of the complexity and non-linearity of these systems, the behaviour predicted by a given mechanical model can vary widely. One must therefore take great care in developing solution methods or approximations because, for example, a solution derived under the assumption that the viscous dissipation is negligible can be shown to give vastly different predictions of intrusion geometry or growth rate than would be derived under the assumption that the energy associated with fracture growth is negligible (e.g. Lister & Kerr 1991; Savitski & Detournay 2002).

In devising a modelling approach upon which geological interpretation can be based, one of the basic questions is whether energy is primarily dissipated due to viscous fluid flow in the fracture or due to the process of breaking material bonds in the crack tip region, a process that is often expressed via a critical energy release rate or fracture toughness. Hence it is useful to speak of viscosity- and toughness-dominated regimes of fracture/sill growth, where viscosity-dominated refers to the regime in which most of the energy dissipation is due to viscous fluid flow, and toughness-dominated refers to the regime in which most of the energy is dissipated due to the process of rock fracturing. These regimes are useful, because they define classes of fracture behaviour where the contribution of phenomena related to fracture toughness or fluid viscosity, respectively, can be neglected, thereby providing significant simplification for fracture modelling and analysis. In many cases it is these simplifications that make the underlying equations tractable (e.g. Spence & Sharp 1985; Detournay 2004), but it then remains vital that the propagation regime associated with a given approximate solution is well-understood, if it is to be applied in the interpretation of geological features.

A number of past contributions present analyses of propagation regime based on relative values of quantities derived from a coupled elastic fracture/lubrication theory mathematical model (after Khristianovic & Zheltov 1955) for a variety of idealized geometries (e.g. Spence & Sharp 1985; Lister & Kerr 1991; Adachi 2001; Savitski & Detournay 2002; Detournay 2004; Garagash & Detournay 2005). Yet the regime of propagation is still generally determined only from analysis of groups of parameters that are specific to particular idealized geometries.

A peculiar aspect of fluid-driven fractures is the multi-scale nature of the tip region. For example, the aperture *w* in the tip region has been shown to be given by linear elastic fracture mechanics (LEFM) in the toughness-dominated regime, so that *w*∝*s*^{1/2}, where *s* is the distance from the tip. Conversely, for the viscosity-dominated regime, *w*∝*s*^{2/3} when viewed at the scale of the fracture (e.g. Detournay 2004). Each of the aforementioned contributions, at varying levels, discuss the dependence of the behaviour of the near-tip region on the regime of propagation. However, the details of the complex multi-scaled near-tip region require careful consideration, and it is the solution resulting from this tip analysis that leads to the theoretical foundation for the present work.

The purpose, then, of the present work is to utilize the previously derived (Garagash & Detournay 2000, 2005) and recently experimentally verified (Bunger & Detournay 2008) near-tip solution in order to present a rigorous and straightforward method for determining fracture propagation regime in a manner that can be applied locally to a growing fluid-driven fracture, provided that the fracture's leading edge is sufficiently smooth (i.e. there are no sharp corners in the overall fracture shape). After reviewing the relevant theoretical and experimental contributions and describing the new method for determining propagation regime, the paper concludes with a discussion of the implications of this analysis for numerical modelling of fluid-driven fractures, such as magma-driven sills, with the elastic fracture/lubrication theory model.

## The sill growth model

One of the most common approaches to modelling magma transport in Earth's crust is to use a model that considers elastic fracture mechanics coupled with laminar fluid flow within the fracture as described by lubrication theory. This model can be traced back to the seminal work of Khristianovic & Zheltov (1955), and has been proposed for modelling magma-driven fractures in a number of contributions (e.g. Spence *et al.* 1987; Lister & Kerr 1991). A form of this model will be presented below. The viability of this model, in light of its simplifying assumptions, will also be discussed.

At the outset the model will neglect thermal effects on the host rock, solidification processes in the magma, and the permeability of the host rock. This is akin to assuming sufficiently rapid propagation so that: (1) the fracture tip remains ahead of the thermal and fluid diffusion fronts; (2) nearly all of the fluid (i.e. magma) supplied, for example from a feeder dyke, is stored in the growing sill rather than diffusing into the surrounding rock; and (3) insufficient time is allowed for the heat transfer required to bring about a phase transition in the magma or host rock during the time of sill extension. Of course these assumptions must be critically evaluated on a case-by-case basis. Further, it is assumed that the effects of gravity are negligible and therefore consideration does not cover the growth of features such as buoyancy-driven dykes.

The mathematical model then begins with Reynold's lubrication equation, which relates the fracture opening *w* to the fluid pressure *p*_{f} by combining the mass balance equation for an incompressible fluid with the fluid flux expression given by the Poiseuille equation. Hence,
1
where μ′=12 μ with dynamic fluid viscosity μ, δ(*x,y*) denotes the Dirac delta function with the origin of the system of coordinates (*x,y*) taken to coincide with the injection point, *Q*_{o} is the mean volumetric injection rate and ψ(*t*) expresses the variation of the injection rate with time (see Table 1). Although the Newtonian lubrication model is an accepted approach for a wide range of applications, it should be noted that magmas containing a significant crystal fraction could be expected to deviate significantly from Newtonian behaviour. However, non-Newtonian versions of this model can be developed and the basic elements of the approach remain unchanged (e.g. Desroches *et al.* 1994). Additionally, compressibility is sometimes considered to be important to certain occurrences of magma-driven fracture behaviour (e.g. Woods *et al.* 2006); however, this is expected to be mainly relevant to the problem of magma ascent in dykes, which has already been omitted from consideration under the present assumptions.

The lubrication equation is coupled to a description of the fracture opening based on linear elastic fracture mechanics (LEFM), which includes both a propagation condition (to follow) and a relation coupling the fracture opening with the internal pressure and far-field stress. Taking a mean far-field stress σ_{o} acting perpendicular to a fracture with the a shape in plan view (i.e. the fracture footprint) denoted by *S*(*t*), the elasticity equation can be expressed as
2
where *E*′ is the so-called plane strain elastic modulus of the solid, *x*′ and *y*′ are variables of integration, and ϕ(*x*,*y*) embodies spatial variation in the far-field stress. Closed-form solutions have been derived for some very simple geometries and pressure/stress distributions; however, this relation usually requires numerical solution.

For opening mode fractures, well-known asymptotic analysis of the elasticity equations imply that the near-tip opening is given by
3
where *K*_{I} is the mode I stress intensity factor, *s* the distance from the fracture tip (with the *s*-axis directed inwards), and *L* is the fracture length. Fracture propagation is then taken to require the availability of sufficient energy to break the rock just a head of the fracture tip. Within LEFM this is conveniently expressed by *K*_{I}=*K*_{Ic}, where *K*_{Ic} is the fracture toughness of the material. The basic assumption for this approach to fracture propagation is that the inelastic damaged zone near the fracture tip is sufficiently small so that there exists an intermediate scale, outside of the damaged zone, but small relative to the fracture size at which the energy release rate associated with the fracture can be computed with a desired accuracy using only the leading order term of the asymptotic expansion of the elasticity equation, i.e. Equation 3 and its corresponding *s*^{−1/2} singular term for the near-tip stresses.

A couple of additional points should be made regarding the application of LEFM in this context. Firstly, laboratory experiments and field observations give strong evidence that *K*_{Ic} is not a constant rock parameter, but instead it can increase significantly with increasing fracture size and *in situ* stress (e.g. Rubin 1993; Carpinteri 1994; Fialko & Rubin 1997). In spite of a number of important contributions, it remains difficult to confidently select, based on the available data, an appropriate value for *K*_{Ic} that would govern a large scale/large *in situ* stress process. One may hope that careful ongoing experimentation and modelling will eventually provide some insight; however until then, we move ahead using Equation 3 with *K*_{I}=*K*_{Ic}. Future work may lead to important refinements to this approach.

While this first issue is primarily concerned with the physics of rock fracturing, a second issue is related to the mathematical properties of the elasticity equations leading to Equation 3. For illustration, consider a two-dimensional fracture of length 2*L* under plane strain conditions, which is subjected to an internal uniform fluid pressure *p*_{f} and an isotropic *in situ* compressive stress σ_{o} (Fig. 1). The Williams (1957) type expansion of the near tip stresses is given by
4
where *f*_{ij}(θ) is a known function of the angular coordinate θ and *O* is used to indicate terms which are of this order and higher in *s*/*L*. This representation of the near-tip stresses has been used to argue that the region that is dominated by leading term, which is the one retained in LEFM analysis, requires that *s*/*L*≪(*p*_{f}−σ_{o})^{2}/*p*_{f}^{2}, which, for typical values of the *in situ* stress and fluid (magma) pressure, would render the zone governed by the LEFM term so small that it would typically be irrelevant (Rubin 1993, 1995). However, it is easily shown using the well-known *J*-integral method Rice (1968, 1974) that the uniform stress term (second term on the right-hand side of Equation 4) does not contribute to the computation of the energy release rate. Therefore, provided that one is considering a problem for which LEFM considerations are used solely as a tool for determining whether a fracture has sufficient energy for propagation (i.e. to enforce the condition *K*_{I}=*K*_{Ic}), the applicability of this theory is not undermined by the presence of large *in situ* stresses, provided that the inelastic zone near the tip is sufficiently small and that it is possible to appropriately account for (e.g. Khazan & Fialko 1995) the expected dependence of the fracture toughness on the *in situ* stress.”

Returning to the mathematical description, it can be seen immediately that this problem involves at least one moving boundary. In fact, generally the fluid and fracture fronts do not coincide with one another, thus requiring the tracking of two moving boundaries (e.g. Bunger 2005; Lecampion & Detournay 2007). However, it can be demonstrated (Garagash & Detournay 2000) that satisfying the condition
5
with *V* the fracture tip velocity, implies that the lag between the fluid and fracture fronts can be neglected. In this case the mathematical model is completed by the addition of the boundary conditions at the moving front
6
7
where an alternative form of the fracture toughness *K*′=*K*_{lc}(32/π)^{1/2} has been introduced to keep the equations uncluttered. Here the first equation expresses the LEFM propagation condition discussed above, while the second simply expresses a no-flux boundary condition at the fracture tip.

## The near-tip region

In spite of the geometric specificity required in order to solve the above system of equations, it can be demonstrated that the region close to the fracture tip (Fig. 2) has a universal form (Garagash & Detournay 2000, 2005). This universal form will eventually be used to obtain a geometry-independent method for determining the propagation regime, but its basic results must first be examined. Non-dimensional forms of the distance from the tip , the fracture opening , and the internal net pressure may be introduced according to 8 with 9 Noting that variation in the far-field stress can typically be neglected when examining only the near-tip region, it can be shown rigorously (Bunger & Detournay 2008) that the above governing equations for fluid-driven fracture growth degenerate near the fracture tip to 10 11 12

Obviously, the fracture opening nearest the tip is described by the limit in Equation 12. However, simultaneous solution of Equations 10 and 11 implies another asymptotic form given by
13
which in dimensional form is given by
14
This 2/3 asymptote of the fracture tip opening is recognized as an asymptotic form associated with fluid-driven fracture propagation when *K*_{Ic} is equal to zero (Spence & Sharp 1985; Lister & Kerr 1991; Desroches *et al.* 1994). Its reconciliation with the LEFM asymptote Equation 12 under conditions of ‘small toughness’ relies on the realization that the near-tip region has the form of a boundary layer solution when the fracture is propagating in the viscosity-dominated regime, with the LEFM asymptote providing the inner behaviour ( → 0, closest to the tip) and the 2/3 asymptote providing the outer behaviour ( → ∞, but still near the tip relative to the rest of the fracture) (Garagash & Detournay 2000). Hence, Garagash & Detournay (2005) present a numerical solution connecting the inner and outer asymptotics so that the universal form of the fluid-driven fracture tip is known (Fig. 3).

The validity of the autonomous generalized tip asymptote () has recently been confirmed with a series of laboratory experiments that consisted of driving circular hydraulic fractures through polymethyl methacrylate (PMMA) or glass specimens using fluids which were solutions of water, blue food dye and either glycerine or glucose (see Bunger & Detournay 2008 for details of the design and results of these experiments). Experimental data pertaining to 10 different tests are presented in Figure 4 together with the tip solution (). The experimental results exhibit some scatter, mainly due to the fact that the fracture opening becomes very small in the tip region, which can be to the detriment of the signal to noise ratio for the measurements. Nonetheless, the close agreement between the experimental and analytical results for the fracture tip opening uphold the multiscale tip asymptote that has been developed to describe the tip region of hydraulic fractures.

Physical understanding of the near-tip behaviour is aided by recognition of the fact that the fracture tip region is governed by two physical processes, which in this case happen to be energy dissipation mechanisms. Here energy is dissipated both in viscous flow of the fluid and in breaking of bonds just ahead of the fracture tip. As is typically the case when multiple physical processes are coupled, each process carries with it its own characteristic length scale (or scales), so that a multi-scale problem naturally arises. For the present problem, LEFM processes occupy the region closest to the tip, with coupling between the fluid and solid governing the outer tip behaviour. Furthermore, the tip length scale ℓ gives a characteristic length associated with the transition between the region that is LEFM dominated and the one that is dominated by fluid–solid coupling (and hence the 2/3 asymptote).

## Determining propagation regime

As will be demonstrated in the sections to follow, the regime of fluid-driven fracture propagation has a profound impact on approaches to modelling fracture growth numerically. Furthermore, for simplified geometries it may be desirable to use an analytical solution based on some simplified form of the governing equations; however, the development and application of these analytical solutions depends crucially on the propagation regime. Hence the motivation is to devise a method for estimating the propagation regime that is both straightforward and that does not rely on quantities that are associated with a particular idealized geometry.

A few details of the criterion proposed here depend on how one defines the tip region. Here a definition of the fracture tip is considered that is relevant to numerical modelling. In particular, the tip region is taken as a single point located at *s*_{t}=0.01 *L*. For example, from a numerical modelling perspective this could correspond to the centre of the discretization element nearest the fracture tip, if the fracture is meshed using 50 equal-sized elements.

The question in determining propagation regime then becomes which, if either, of the above tip asymptote Equations 6 or 14 give an accurate estimate of the fracture opening at *s*_{t}. The propagation regimes can then be defined in terms of the thresholds
15
for the toughness-dominated regime and
16
for the viscosity-dominated regime. Choosing the desired accuracy α=0.02, making use of the numerical tip solution (Fig. 3) and the tip scaling Equation 9, one arrives directly at conditions in terms of the tip length scale ℓ and an estimate of the overall fracture extent *L*, namely

toughness-dominated, ℓ>3800

*L*;viscosity-dominated, ℓ<0.03

*L*.

A graphical representation of these criteria is portrayed in Figure 5. Determining the propagation regime then requires estimating ℓ and *L*, either *a priori* from physical reasoning and parameter value estimates, or *a posteriori* from the solution to a numerical model or analysis of data from analogue experiments. It is important to note that this method relies on a local estimate of the propagation regime based on a local value of the fracture tip velocity *V*, which may be different at different locations around the fracture front. However, for reasonably simple geometries it may be valid to estimate the overall propagation regime based on an average expected tip velocity. Alternatively, bounds could be placed on the regime of propagation based on upper and lower bound estimates on the local tip velocity.

It is apparent that the definitions of the regimes of propagation are specific to the definition of the tip region *s*_{t} and the desired accuracy α. Hence these regimes are not defined in absolute terms, but instead require one to decide what is meant by ℓ/*L* being either large enough or small enough to warrant a given classification. This issue is not unique to the method of determining regime presented here. It also arises in methods that are based on parameters derived from the global geometry (e.g. Lister & Kerr 1991; Detournay 2004), in that one must always decide how small or large the value a group of parameters must be in order for the fracture to be in the viscosity or toughness dominated regime. In some situations a different choice of *s*_{t} or α may be appropriate. Tables 2 and 3 gives the coefficients *c*_{k} and *c*_{m} for different combinations of *s*_{t} and α, where ℓ>*c*_{k}*L* and ℓ<*c*_{m}*L* define the toughness- and viscosity-dominated regimes, respectively.

## Expected regime for sill growth

The challenge remains to acquire appropriate estimates of the necessary parameters. At this point such choices, and even the degree to which the feature of interest may satisfy the necessary assumptions for application of the theory, are specific to each field case. However, one can attempt to make a few general, albeit somewhat cautious, statements regarding sill growth based on expected ranges of parameter values.

First a range of magma viscosity values is taken as 100≤μ≤2000 Pa s based on the discussion of Lister & Kerr (1991). Here the selection is essentially limited to relatively low viscosity basaltic magmas. In the case of granitic magmas, which are expected to have considerably higher viscosity values than the range considered here, significant questions remain regarding how to interpret the mechanism of emplacement based on field data. Whether one is considering diapirism or episodic tabular injection (e.g. Petford *et al.* 2000), it leads to the same conclusion that the emplacement mechanism for granitic bodies may in many cases lie outside of what can be captured by the mathematical model under consideration here.

The host rock is then characterized by a plane strain modulus taken from the range 10≤*E*′≤70 GPa and a fracture toughness taken from the range 1≤*K*_{Ic}≤120 MPa m^{1/2}. When discussing the mechanical properties of rocks it is impossible to escape the issue of size dependance. As previously discussed, there remains large uncertainty regarding how to quantify the increase of *K*_{Ic} with fracture size and *in situ* stress as one pushes the theories to geological scales. For the basic estimates desired, the considered range accounts for a possible factor of 1000 increase in the critical fracture energy *G*_{c} (noting the well known relation *K*_{Ic}∝*G*_{c}^{1/2}) at the scale of sill growth relative to the laboratory-scale range, which is taken up to *c*. 4 MPa m^{1/2}.

Finally, the tip velocity *V* and fracture extent *L* must be considered. In field settings the fracture extent *L* is typically known from observation. However, estimating the tip velocity requires one to assume a model. For example one could make a rough estimate based on a known *L* and bounds (e.g. Petford *et al.* 2000) on the total emplacement time. Alternately, in this case *V* will be estimated based on a solution for a viscosity-dominated, penny-shaped hydraulic fracture presented by Savitski & Detournay (2002) in which the fracture extent (here the radius) is given by
17
Here the magma supply rate *Q*_{o} is considered to be constant and in the range 1≤*Q*_{o}≤10^{5} m^{3} s^{−1} (after Lister & Kerr 1991). The tip velocity is, of course, given by simple differentiation of Equation 17 with respect to time. The time of injection *t* is then selected so as to give sills which fall within the range 10 m≤*L*≤10 000 m based on Equation 17. Upon subsequent computation of the tip length scale ℓ for 40 000 randomly selected combinations of parameters from within the ranges listed above, one can identify the area in the propagation regime diagram in which natural sills are expected to lie. Tabulation of these test cases is presented in Figure 6. Here we see that approximately 90% of the realizations fall in the viscosity-dominated regime; a fraction of 1% fall in the toughness-dominated regime, with the remainder lying in the transitional regime.

Note that a solution for a viscosity-dominated fracture was chosen at the outset for estimation of the fracture size *L* and tip velocity *V*. On the other hand, one could have chosen a toughness-dominated solution presented by Savitski & Detournay (2002). In this case one would have predicted, for a given injection rate and time, a larger fracture size and tip velocity than with the viscosity-dominated model. This would drive the propagation regime estimates for the parameter ranges given here even further to the viscosity-dominated regime so that, on examination of these results in a diagram such as Figure 6, one would immediately see that the viscosity-dominated fracture model is more appropriate for the majority of the realizations.

It is important to reiterate that the current results have been obtained based on a somewhat arbitrary assumption about how fracture energy increases with fracture size. For example, if *K*_{Ic} is taken in a range up to 400 MPa m^{1/2}, then only 50% of the realizations lie in the viscosity dominated regime, 1% lie in the toughness-dominated regime, and the rest are in the transitional regime. It should be further noted that it is possible, then, that the lower end of the fracture toughness range may be irrelevant at the scale considered, in which case toughness-dominated behaviour could in fact be more important than is indicated by this analysis. Nonetheless, it is clear that analysis requires careful consideration of regime, and that in most cases it will be required to consider the effects of viscous dissipation. This work has also made it clear that development of models for analysis of geological data is reliant on careful experimentation and analysis in order to better understand how fracture energy/toughness scales fracture size.

## Implications for numerical modelling

In order to track the moving fracture front, fluid-driven fracture models almost invariably involve imposing a condition of the fracture tip opening (or near tip stress field) based on an expected asymptotic behaviour. These methods are well accepted in fracture mechanics modelling, where models based on the LEFM asymptote (Equation 6) form an important and very commonly employed class of propagation criteria. One temptation that follows from this history is to bring an approach from LEFM to the problem of fluid-driven fractures without accounting for the consequences of fluid–solid coupling in the fracture tip region. Of course, if the fluid can be assumed inviscid, this approach is valid. However, it has been demonstrated in detail above that, for many, if not the majority of the cases which are of interest to the geosciences (and indeed industry also), fluid–solid coupling is important, and as a result the region that is governed by the LEFM asymptote may be very, very small relative to the size of the fracture that is being modelled.

One further temptation is to attempt to address the fact that the LEFM asymptote may only govern a very small region near the fracture tip by simply using a large number of elements. Based on the analysis presented above, one can find directly that maintaining 2% accuracy at the tip element with a uniform meshing of the fracture requires
18
elements. Hence it can be seen that accurate modelling of the magma-driven sills using LEFM considerations only would typically (based on 90% of the realizations) require *N*<10^{8} elements *in each modelling dimension*, assuming a uniform mesh. Obviously, grading the element size from the tip towards the inlet would significantly decrease the number of elements needed, although at the expense of substantial algorithmic complexities. Also, the number of elements required would, of course, decrease upon relaxation of the accuracy requirement, but only by about an order of magnitude if one were to require only 10% accuracy. Furthermore, it is important to realize that a strong increase in the fracture toughness with fracture size could also drive the sills towards the toughness-dominated regime and thus reduce the number of elements required for accurate modelling based on LEFM considerations.

Nonetheless, it is clear that in most cases modelling fluid-driven fractures using a tip condition based on LEFM is expected to be, at best, computationally impractical and at worst a brute forcing of the wrong physics at the wrong scale. In contrast, it has been demonstrated that accurate numerical solutions (compared with known analytical solutions for simple geometries) can be obtained for viscosity-dominated fractures using very coarse meshing, i.e. *N*≈10 elements in each dimension for the two-dimensional model, provided that the 2/3 (Equation 14) asymptote is imposed at the tip elements rather than the LEFM solution (Bunger *et al.* 2007).

## Implications for predicting the thickness of sills

Although modelling the aspect ratios observed in particular natural sills is beyond the scope of this paper, it is useful to determine whether the present model is able to produce aspect ratios which are within the same orders of magnitude as are observed in the field. According to graphical catalogs (McCaffrey & Cruden 2002; Breitkreuz & Petford 2004; Cruden & McCaffrey 2006) of natural intrusions spanning the classifications of minor intrusions, sills, laccoliths, and plutons, the characteristic thickness *w*_{o} to length *L* ratios (or just ‘thickness ratios’) are expected to lie in the range 0.0005<*w*_{o}/*L*<0.5. Recall that the present analysis neglects: (1) formation of tabular bodies by periodic injection of magma; (2) the potential for complex, non-elastic deformation of the host rock; and (3) the large viscosity range which would typically be associated with the magmas that form granitic intrusions. Hence the present analysis is expected to provide a lower bound to the observed range of thickness ratios.

As discussed above, provided that *K*_{Ic} does not become too large as the fracture size and *in situ* stress increase, one would expect the majority of sills to grow in the viscosity-dominated regime. In this case, using the scaling relations derived by Savitski & Detournay (2002), it is easily shown that the relationship between *w*_{o} and *L* is estimated by
19
It is shown here that fracture growth is expected to be characterized by decreasing thickness ratios or lateral spreading, as its length *L* increases. Taking the parameter ranges above and *L*=10 m, one finds that 0.001<*w*_{o}/*L*<0.01, which is consistent with the available data for very small intrusions. However, if one considers cases for which *L*=10 km, one finds that *w*_{o}/*L*<0.0004 so that the predicted range of thickness ratios is smaller than what has been measured on natural sills.

Interestingly, the departure of the predictions based on viscosity-dominated fractures from the available data corresponds approximately to the point at which intrusions can safely be expected to have grown to a lateral extent that is significantly greater than the depth of emplacement. Again, making use of previously derived scaling relationships, in this case for a shallow (i.e. *L*>5*H*, where *H* is the emplacement depth) circular fracture that is governed by the elastic fracture/lubrication model and is propagating in the viscosity-dominated regime, one finds that (Bunger 2005)
20
For the parameter ranges used above one can demonstrate thickness ratios up to *w*_{o}/*L*≈0.01. Recall, however, that the ranges above do not account for the large viscosities associated with granitic magmas. Taking, for example, μ=10^{6} Pa s one obtains a thickness ratio of *w*_{o}/*L*≈0.07. So consideration of the effect of the Earth's surface on fracture growth would be expected to give sill thickness ratios that are within the range of the observed data. Furthermore, one finds that these shallow, viscosity-dominated fractures are expected to be characterized by growth with a constant thickness ratio *w*_{o}∝*L*, which is reasonably consistent with the scaling observed for large intrusions (McCaffrey & Petford 1997).

Promising comparisons such as this should drive a methodical and rigorous improvement of the mathematical model for the emplacement of sills and other intrusions. For example, one may hope that as we come to better understand large-scale rock fracture and the formation of tabular/sheeted intrusions, this model may be extended to encompass the range of intrusion geometries observed in Nature.

Others may suggest that elastic fracture mechanics should be replaced by approaches that account for the details of the plastic deformation in the fracture tip region, but it is important to recognize that such an approach requires plasticity at a sufficiently large scale relative to the fracture size. Hence one must use caution owing to the fact that it is possible to have inelasticity large relative to the observational (i.e. outcrop) scale, but still small relative to a kilometre-sized intrusion, so that elastic fracture mechanics may in fact still be applicable. Specifically, the present theory is built on the premise that there exists a lag region between the fluid and fracture fronts which is very small compared with the fracture length, but much larger than the zone of near-tip plastic deformation. Although the interplay among length scales in the near-tip region is a complex matter, a reasonable starting point is to require that the intrusion length be about 1000 times or more than the observable zone of plastic deformation.

## Conclusions

The lubrication/elastic fracture model, in various forms, is perhaps the most common approach to modelling magma transfer in the Earth's crust. While solving this coupled, non-linear system of equations for non-ideal geometries is in general an arduous task, the near-tip region possesses a universal form that is leveraged in the present work in order to devise a geometry-independent method for determining the regime of propagation for magma-driven sills. In so doing, it has been demonstrated that viscous dissipation is expected to be an important consideration in modelling sill formation, owing to the fact that most of the combinations of relevant governing parameters lead to sills that are viscosity-dominated or in the transition regime. As a result, numerical solution methods must be devised so as to account for the effect of fluid–solid coupling on the near-tip behaviour, for example by imposing the 2/3 asymptote of Equation 14 rather than the LEFM asymptote (Equation 6) at the fracture tip in the case of viscosity-dominated fractures.

A number of implications for evaluation of geological data based on mechanical modelling have been discussed. Firstly, it has been demonstrated that a full understanding of the propagation regime relevant to field cases requires a clearer understanding of how the fracture energy/toughness scales with fracture size. Secondly, it has been shown that naive application of the LEFM tip asymptote in modelling intrusion growth can lead to application of the wrong physical process governing fracture extension, and therefore to erroneous results. Finally, promising results have been presented, which indicate that the linear elasticity/lubrication model considered here is likely to be capable of producing realistic predictions of the thickness to length ratios for magmatic intrusions, provided that viscous dissipation and near-surface effects are considered.

## Acknowledgments

The present work arose out of a series of discussions with Emmanuel Detournay, for which I am sincerely grateful. Funding for this work has been provided by the Petroleum Research Fund administered by the American Chemical Society (grant no. ACS-PRF 43081-AC8), with additional support from Schlumberger. This support is gratefully acknowledged.

- © The Geological Society of London 2008