## Abstract

The growth of shallow sills is studied in analogue experiments performed in polymethyl methacrylate (PMMA) and glass. The experimental fractures curve towards the surface to become saucer-shaped, which is consistent with many field observations of dolerite sills. The curvature of the saucer is shown to decrease as the *in situ* stress acting parallel to the surface increases relative to an estimate of the strength of the fracture-induced stress field. The initially circular fractures also elongate in plan view to become egg-shaped, a tendency that decreases with increasing importance of viscous dissipation in the growth process. Sill emplacement is further examined mathematically by considering a shallow, circular, fluid-driven fracture propagating in a homogeneous brittle elastic material. The fractures are shown to undergo three transitions related to the mechanics of sill growth. Each transition is associated with a characteristic time that is derived from analysis of the governing equations using scaling methods. These characteristic times provide an estimate of how long viscous flow is the dominant energy dissipation mechanism, how long significant lag between the fluid and fracture fronts is expected to persist, and how long the sill will take to attain an extent that is of the same order as its depth.

## Nomenclature

*E*- Young's modulus
*E*′- plane strain modulus
*g*- gravitational acceleration (9.81 m s
^{−2}) *H*- initial fracture depth
*h*- depth of the fracture tip
- ℋ, 𝒦, ℳ, 𝒮
- dimensionless evolution parameters
*K*_{Ic}(*K*′)- fracture toughness (alternate form)
*L*- scaling factor for fracture radius
*Q*_{o}- volumetric injection rate
*R*- fracture radius
*R*_{f}- radius of fluid-filled region
*R*_{max},*R*_{min}- maximum and minimum radial distance to fracture tip (experiments)
*W*- scaling factor for fracture width
*p*_{f}- fluid pressure
*p*_{v}- vapour pressure in lag region
*q*- volumetric fluid flux
*r*- radial co-ordinate
*t*- time
*t*_{mk},*t*_{om},*t*_{mm̄}- characteristic times
*w*- fracture opening
- χ
- dimensionless parameter controlling fracture curving
- Π
- dimensionless pressure
- Ω
- dimensionless fracture opening
- γ
- dimensionless fracture radius
- γ
_{f} - dimensionless fluid radius
- ϵ
- Small parameter relating
*p*_{f}and*E*′ - ϑ
- Scaling factor for ξ
_{f} - μ (μ′)
- dynamic fluid viscosity (alternate form)
- ν (ν′)
- Poisson's ratio (alternate form)
- ξ
_{f} - fluid fraction
- ρ
- normalized fracture co-ordinate
*r*/*R* - ρ̂
- normalized fluid co-ordinate
*r*/*R*_{f} - ρ
_{g} - rock density
- σ
_{o} *in situ*stress acting perpendicular to the specimen/Earth's surface- σ
_{r} *in situ*stress acting parallel to the specimen/Earth's surface

Magmatic sill intrusions are commonly saucer-shaped. This morphology is demonstrated by the great dolorite sills such as those in the South African Karoo Basin and numerous examples from seismic images (Francis 1982; Chevallier & Woodford 1999; Malthe-Sørenssen *et al.* 2004; Thomson & Hutton 2004; Hansen & Cartwright 2006). Sill intrusions are thought to have a profound impact on the mechanical properties, geomorphology and hydrocarbon development in their host basins (Malthe-Sørenssen *et al.* 2004) and may be source regions for hydrothermal vent complexes (Svensen *et al.* 2004). Saucer formation is also fundamental to industrial hydraulic fracturing in mining (Jeffrey & Mills 2000), environmental remediation (Murdoch & Slack 2002) and ion etching of semi-conductors (Gigùere *et al.* 2005, 2006). Previous studies have focused on describing the saucer shapes and discussing potential models (Francis 1982; Chevallier & Woodford 1999; Murdoch & Slack 2002; Thomson & Hutton 2004; Hansen & Cartwright 2006), solving models for flat, shallow fractures (Pollard & Holzhausen 1979; Murdoch 2002; Bunger & Detournay 2005), and obtaining numerical solutions that demonstrate the tendency of shallow fractures to become saucer-shaped (Fialko 2001; Malthe-Sørenssen *et al.* 2004). However, until now it has not been clear as to what geometric and *in situ* conditions are required for formation of the saucer shape and what information the saucer morphology may provide regarding conditions at the time of fracture growth.

This paper presents the results of analogue experiments performed using polymethyl methacrylate (PMMA) and borosilicate glass specimens. The focus of this paper is placed on understanding the processes that dictate two key aspects of the sill morphology; namely, the curvature of the sill that results in the saucer shape and the observed (Watson 1910; Murdoch & Slack 2002; Thomson & Hutton 2004; Bunger 2005; Hansen & Cartwright 2006) tendency of initially circular fractures to favour growth in one direction, thereby elongating to become egg-shaped. We make use of an approach to scaling of the mathematical model that simultaneously considers linear elastic fracture and one-dimensional laminar fluid flow (Savitski & Detournay 2002; Detournay 2004). By coupling this approach with the experimental data we are able to identify two dimensionless groups of parameters controlling the final geometry of the saucer-shaped fractures. The analysis also casts saucer-shaped growth into the context of the physical transitions a shallow, circular, fluid-driven fracture undergoes during its growth history, and thereby makes clear the conditions for which viscous dissipation is important and for which the fluid front is expect to lag significantly behind the fracture front.

## Experimental method

The experiments consist of propagating hydraulic fractures through PMMA and glass specimens using Newtonian solutions of blue food dye (FD & C Blue #1), water and either glycerin or glucose (Fig. 1). The fluids can be considered Newtonian in rheology and are characterized by a dynamic fluid viscosity, μ. By varying the fluid composition, μ ranged from 0.001 to 100 Pa s at room temperature. Fluid was injected by a positive-displacement stepping motor pump capable of producing a volumetric injection rate *Q*_{o} between 0 and 0.16 ml s^{−1}. The transparent, brittle elastic specimens consisted of either 150 mm-diameter cylinders or block-shaped specimens with sides ranging from 150 to 400 mm. A lateral compressive stress, σ_{r}, was applied to the cylindrical specimens using a Hoek-type triaxial compression cell (Hoek & Franklin 1968) and to prismatic specimens using water-filled flat jacks. The fracture toughness, *K*_{Ic}, for the PMMA and glass specimens is 1.3 and 0.6 MPa m^{1/2}, respectively. The plane strain modulus *E*′=*E*/(1−ν^{2}), where *E* is the Young's modulus and ν is Poisson's ratio, is 3.9 GPa for the PMMA and 65 GPa for the glass specimens. Fractures were initiated from an approximately 6 mm-radius initial flaw that was machined at the base of the injection tube, with the initial depth from the surface of the block, *H*, ranging from 12 to 30 mm.

The transparent specimens enabled direct video monitoring of the location of the fluid and fracture fronts, as shown in Figure 2a for a case where there was a significant lag region between the two. In addition, the transition from circular to elongating fracture growth could be readily observed (Fig. 2b). Furthermore, analysis of the intensity of light passing through the fluid-filled portion of the fractures can be used to measure the fracture opening, a method which has been shown to be accurate within about 10% (Bunger 2006).

In addition to video monitoring, the injection pressure and fluid temperature, which was used to obtain viscosity from calibration data, were measured using analogue transducers. The surface displacement was measured with a linear variable differential transformer (LVDT). The injection rate was nominally equal to the volumetric displacement rate of the pump; however, some variation was induced by the compressibility of the fluid and injection system, particularly very soon after initial fracture propagation. Hence, *Q*_{o} was determined either by computing the injected volume directly based on fracture opening measurements or by calibrating fluid flow to the pressure drop across the flow control valve (Bunger *et al.* 2004; Bunger 2005). The flow control valve was also used to moderate the compressibility effects on the injection rate by dissipating some of the elastic energy stored prior to initial fracture growth.

## Timescales and limiting regimes

In order to perform sound experimental design and to understand conditions under which analogue experiments are expected to bear physical similarity (Barenblatt 1996) to sills at the field scale, one must first carefully consider scaling. In this case the scaling takes advantage of a well-developed mathematical model for fluid-driven fractures that has evolved from Khristianovic & Zheltov (1955), and is in the same spirit as the magma ascent model of Lister & Kerr (1991). Deformation of the solid and fracture propagation are considered using linear elastic fracture mechanics (Rice 1968), which is coupled with modelling laminar flow of a Newtonian fluid in the fracture as a one-dimensional process according to Reynold's lubrication theory (Batchelor 1967).

A similar analysis is given for a plane strain fluid-driven fracture by Detournay & Bunger (2006); however, one must be clear that analysis of this sort is particular to the geometry under consideration. Here we consider a circular fluid-driven fracture situated relatively close to the Earth's surface, and this analysis is therefore unique to shallow fractures that grow with at least approximate radial symmetry. In addition, the existence of regimes of hydraulic fracture, dyke and sill propagation has been discussed by several previous authors (e.g. Barenblatt 1962; Spence & Sharp 1985; Lister & Kerr 1991), but only since Savitski & Detournay (2002) have they been cast in the context of small and large time regimes for an evolving fracture.

The analysis, detailed in the Appendix, illuminates three transitions that the shallow, circular fractures under consideration undergo during their growth: (1) from viscous- to toughness-dominated energy dissipation; (2) from the fluid occupying only a small portion of the fracture near the inlet to the fluid filling essentially the entire fracture; and (3) from a fracture that is small to one which is large relative to its depth. The nature of these transitions is determined by the presumed-constant fluid injection rate, *Q*_{o}, dynamic fluid viscosity, μ, fracture toughness, *K*_{Ic}, plane strain modulus, *E*′, initial fracture depth, *H*, overburden stress, σ_{o}, and vapour pressure in the region between the fluid and fracture fronts, *p*_{v} (Fig. 3). Note that for a non-curving fracture one may consider σ_{o} as a surface traction, as drawn, or equivalently as an overburden stress due to the weight of the overlying material (Zhang *et al.* 2002). However, when the fracture curves to approach the free surface, the variation of the stress with the depth may be important. For simplicity in clarifying some basic qualitative characteristics of these fractures, the scaling considers the mathematical model for a non-curving fracture propagating parallel to the surface. Also, equations are kept tidy by using the alternative parameters:
Each of the three transitions can be associated with a characteristic time that provides an estimate for how long each transition will take. These are given by:
1
One can then evaluate (see the Appendix) the role of each of these timescales in the scaled governing equations in order to better describe the physical transitions listed above. We find that:

The role of viscous dissipation diminishes as time (

*t*) increases relative to*t*_{mk}. Hence, for*t*≪*t*_{mk}viscous flow comprises the dominant energy-dissipation mechanism; that is, propagation is in the viscosity-dominated regime. In contrast, for*t*≫*t*_{mk}the fracture propagates in the toughness-dominated regime. In this case viscous dissipation can be neglected relative to the energy dissipation associated with the fracture toughness and so the fluid pressure, which still evolves in time, can be taken as uniform inside the fracture;The so-called fluid lag region, which is the region between the fluid and the fracture front, diminishes in size as

*t*increases relative to*t*_{mk},*t*_{om}or both; andThe fracture size increases appreciably relative to its depth as

*t*increases relative to*t*_{mm̄}.

Consider a sill growing at 2 km depth. Take the rock density as 2500 kg m^{−3} and the gas pressure in the lag region to be negligible relative to σ_{o}, and further let:
2
From equation (1) one finds that *t*_{mk} ≈ 10^{6} years, *t*_{om} ≈ 6 s, and *t*_{mm̄} ≈ 4 days. Note that this sill is predicted to require at least 4 years to grow to a radius of 20 km based on a lower-bound estimate provided by a zero-viscosity, elastically clamped circular plate model (Bunger & Detournay 2005). Hence, the observable life of this sill is expected to be characterized by viscosity-dominated fracture growth with negligible fluid lag (i.e. the magma penetrates nearly to the fracture tip) and considerable influence from the Earth's surface. The interested reader can easily verify that similar behaviour is expected for most realistic combinations of governing parameters for shallow sill growth. However, if gasses are released from the magma so that the net pressure in the lag region is comparable in magnitude to the overburden stress (Lister & Kerr 1991; Bunger & Detournay 2007), then the size of the lag region could remain significant for much of the fracture's life.

## Experimental propagation regimes

Because the overall goal in any analogue experiment is to attain physical similarity with the field-scale phenomena under examination, it would be ideal to replicate fracture growth in a regime typical of sill growth, as discussed above. However, in the experiments presented here, the overburden stress is due only to the atmospheric pressure acting on the block, so *t*_{om} is typically ≈75 days. An experiment typically lasts for only a few seconds or minutes. Hence, in the experiments one would expect either the fluid lag to be significant or the role of viscous dissipation to be negligible. Therefore, the propagation regime expected for field cases could only be reproduced if one were to modify the experimental set-up to simulate the presence of significant overburden stress but maintaining a condition of zero shear traction on the top surface of the specimen.

Because the focus of this experimental study is to examine the parameters controlling the geometry of saucer-shaped sills, the question of the effect of propagation regime should focus on its effect on the fracture path (curving) and the plan-view fracture shape (footprint). To this end, experiments were performed in the two regimes that were readily accessible in the laboratory, namely: (1) with *t*≈*t*_{mk}, *t*≈*t*_{om}, which is a regime with strong viscous dissipation and significant fluid lag (e.g. Fig. 2a); and (2) with *t*≪*t*_{mk}, *t*≪*t*_{om}, which is a regime with negligible viscous dissipation and hence negligible fluid lag (e.g. Fig. 2b). In both cases *t*≈*t*_{mm̄}. The results presented in the following section demonstrate that the regime of propagation, as defined here, determines the degree to which the fracture will persist in circular (plan-view) growth but has only a second-order effect on the fracture curving that leads to the saucer shape.

## Fracture curving and elongation

Analysis of the scaling of the governing equations illuminates the timescales associated with three physical transitions the fractures/sills undergo as they grow and provides a means for interpretation of the occurrence of significant fluid lag in some of the experiments. The analysis also sets the stage for understanding which parameters control two key aspects of the experimental saucer morphologies, namely fracture curving and elongation.

The analysis, presented in the Appendix, of the expected pressure drop along the fracture in terms of a dimensionless viscosity ℳ = (*t*/*t*_{mk})^{−2/5} (see equation 25) forms the basis for interpreting the tendency of the initially circular experimental fractures to eventually favour growth in one direction over the others, thereby producing an elongated egg-shaped geometry in plan view (Fig. 2b). But first it should be made clear that this discussion is pertinent to fracture elongation that can be interpreted to be caused by perturbations to the stress field or material properties in the case of a uniform *in situ* (horizontal) stress acting parallel to the surface. When the horizontal principal stresses differ appreciably, this stress difference may in fact control how the fracture will elongate.

Recall that the value of ℳ determines the importance of viscous dissipation in the growth process and formally sets the order of magnitude of the pressure gradient or the strength of the pressure drop arising from viscous flow. Qualitatively, the internal pressure loading is concentrated near the fracture centre when ℳ is large (*t*≪*t*_{mk}) and distributed uniformly along the fracture when ℳ is small (*t*≫*t*_{mk}). Alternately, we can consider the dimensionless toughness 𝒦= ℳ^{−5/18}, noting that 𝒦<0.5 corresponds to ℳ>10, and 𝒦>1.9 corresponds to ℳ<0.1. In general, 𝒦 varies from zero to infinity with increasing time, *t*. However, for a typical experiment useful data can only be obtained after 0.3 s. Similarly, the tests must be stopped when the fracture reaches the edge of the specimen, which occurs in a few seconds to a few minutes. Thus, 𝒦 covers only a limited range for a given experiment (i.e. 1.5<𝒦<1.8), so that a mean value of 𝒦 can be associated with each test.

Although further research is required to formalize the stability analysis for the shallow, curving fractures, our intuition may be guided by the stability analysis that Gao & Rice (1985) performed for perturbations to a circular fracture in an infinite medium. Their analysis shows that larger perturbations are required to induce unstable asymmetric growth when the internal pressure loading is more concentrated near the fracture centre. In other words, concentrated loads favour the stability of the circular shape.

Figure 4 demonstrates that the circular configuration was more stable in the experiments for smaller values of 𝒦; that is, when the pressure loading is expected to be more concentrated near the fracture centre. Figure 4a gives the evolution of the elongation ratio, *R*_{max}/*R*_{min}, where *R*_{max} and *R*_{min} are the maximum and minimum distance measured to the crack tip along radial lines emanating from the injection point. Note that *R*_{max}/*R*_{min}=1 corresponds to a perfectly circular fracture. Evolution is given in this case as a function of the mean fracture radius, *R*, to the depth of the crack tip, *h*. Three cases are compared. In all three, the tendency of *R*_{max}/*R*_{min} to remain the same or decrease through the first part of the experiment indicates a period during which the circular configuration is stable. However, at some point the elongation ratio begins to increase with a slope that corresponds to the degree to which the fracture is elongating. Hence, for the test with an average value of 𝒦=3.4 (ℳ=0.01), it is observed that the onset of elongation occurs when the fracture is smaller relative to its depth, and the elongation is more pronounced than when 𝒦=0.94 (ℳ=1.2). The case when 𝒦=2.4 (ℳ=0.4) is shown to be an intermediate case between the two. Figure 4b and c confirm the tendency throughout the experimental series of the fractures to become larger relative to their depths before onset of elongation and to elongate less dramatically for smaller values of 𝒦.

While the nature of the internal loading appears to play a key role in the degree to which an initially circular fracture will tend to elongate, the experiments show that it does not strongly change the curving of the saucer shape. Rather, fracture curving depends mainly on a dimensionless stress parameter:
3
which was first suggested by Zhang *et al.* (2002), but is here modified to clarify its dependence on the stress difference (σ_{r}−σ_{o}) rather than only σ_{r}, as defined in Figure 3. This parameter relies on the usual assumption in fracture mechanics that the crack path is determined by the near-tip stress field. Provided fracture propagation implies *K*_{I}=*K*_{Ic} and the fracture size is of the same order as the depth *H*, gives an estimate of the strength of the elastic near-tip stress field induced by the pressurized fracture. Because of the asymmetry introduced by the presence of the free surface, this induced stress field will be non-symmetric about the crack tip. The asymmetry can be considered as a consequence of the upper fracture face deforming more than the lower face (Zhang *et al.* 2002; Malthe-Sørenssen *et al.* 2004). It is this asymmetry in the induced stress field that drives fracture curving. However, the *in situ* stress difference (σ_{r}−σ_{o}) will oppose this tendency to curve because of the fracture's desire to remain oriented so that its opening is parallel with the least compressive *in situ* stress. Hence, χ can be interpreted as a parameter that compares the stress fields opposing and driving fracture curving.

Figure 5 shows the dependence of the saucer shape (i.e. the crack path) on the parameter χ. Figure 5a and b show the progressive flattening of the saucer for increasing values of χ for experiments performed in PMMA and glass, respectively. Figure 5c shows a comparison between experiments performed in these two materials for which χ was the same in spite of the fact that *K*_{Ic} is different for the two materials. It is also important to realize that these experiments represent cases with 𝒦 ranging from 0.5 to 4.5, so the details of the internal pressure loading are appreciably different from one test to another. Hence, the experimental results show that two fractures characterized by the same value of χ will have the same saucer shape up to a re-scaling by the initial depth *H* without respect to the values of the dimensional parameters or the details of the internal pressure loading.

## Discussion

In summary, the results presented here show that mechanical interaction of a sill with the Earth's surface is sufficient to generate saucer morphologies. Furthermore, the influence of naturally occurring perturbations to the stress field or material properties can lead to elongation of initially circular fractures. Hence, this simple model captures two of the obvious geometric features observed in nature. However, there remain some issues that this model does not address which are likely to be important to the growth of natural sills. We will discuss four issues that seem to be of the greatest importance.

First, we consider propagation in a homogeneous medium. In actuality, the sills often grow in stratified sedimentary basins, and the presence of these heterogeneities may control formation of stepping morphologies as the saucer transects the sedimentary layers (Francis 1982; Malthe-Sørenssen *et al.* 2004). Also, curving sills are not necessarily indicative of interaction of the sill with the surface. For example, a curving sill could result from the influence of a low-stiffness layer overlying the layer in which the sill forms (Selcuk *et al.* 1994).

Second, our experiments neglect the effect of the overburden; that is, the vertical stress component. This approach is reasonable for analysis of the experiments, but is not valid for field cases. Although the experiments consider only cases where the overburden stress, σ_{o}, is negligible, it can be shown from elastic fracture simulations that the parameter χ depends on the difference between the lateral and vertical stresses (σ_{r}−σ_{o}) (Martynyuk 2002; Xi Zhang pers. comm. 2006). However, the fact that these stresses increase with depth would have some effect on the crack shape over what is observed in the experiments and would depend on at least one additional parameter of the form β=ρ_{g}*gH*^{3/2}/*K*_{Ic}, where ρ_{g} is the overburden density and *g* is gravitational acceleration. Also, as demonstrated in the example calculation of the characteristic time *t*_{om}, the vertical stress for natural sills would usually be expected to cause the lag to be very small even when viscous dissipation is very large. In this case there is strong coupling between the fluid and solid in the tip region, and as a result the elastic stress singularity is changed from the usual linear elastic fracture case to a form that is unique to fluid-driven fractures (Spence & Sharp 1985; Lister & Kerr 1991; Desroches *et al.* 1994; Garagash & Detournay 2000; Bunger *et al.* 2005). Alteration of the near-tip stress field would probably alter the crack path; that is, the saucer shape. But by how much the path would change remains unknown. Finally, the influence of buoyancy forces associated with density differences between the magma and the rock remain unclear.

Third, it is important to acknowledge that important physical mechanisms associated with sill emplacement may be neglected both by the mathematical and analogue models considered here. For example, crystallization and rheological changes in the magma as it cools could have an important effect on the geometry observed in natural sills.

Fourth, we consider only the case where the lateral stress is the same in all directions. In this case the direction of eventual elongation of the experimental fractures is random and the curvature of the saucer in all directions is very nearly the same. However, in the Earth's crust one would expect directions of maximum and minimum horizontal stress to exist. Natural saucer morphologies are expected to reflect this fact. For example, the Golden Valley Saucer in South Africa's Karoo Basin (Malthe-Sørenssen *et al.* 2004) dips less steeply along its axis of elongation. Building on the relationship established by this research between the flatness of the saucer and χ, we expect and have observed in some preliminary experiments that the saucers will be elongated and flatter along an axis aligned with the maximum horizontal stress direction.

## Conclusions

Shallow sill growth has been modelled mathematically, and in analogue experiments as shallow, circular, fluid-driven fractures propagating in a brittle elastic, homogeneous material. The sills are shown to undergo three physical transitions as they grow. Each transition is associated with a characteristic time that is derived from analysis of the governing equations. These characteristic times provide an estimate of how long viscous flow is the dominant energy-dissipation mechanism, how long significant lag between the fluid and fracture fronts is expected to persist, and how long the sill will take to attain an extent that is of the same order as its depth.

The experimental fractures curve to become saucer-shaped and elongate to become egg-shaped in plan view. Each of these characteristics of the saucer morphology has been shown to correlate with a dimensionless group of parameters. Namely, the curvature of the saucer is shown to decrease as the *in situ* stress acting parallel to the surface increases relative to an estimate of the strength of the fracture-induced stress field. Further, the tendency of the fracture to elongate is shown to decrease with increasing importance of viscous dissipation.

Further work is underway in order to experimentally investigate near-surface fracture in biaxial stress states to develop better understanding of the asymmetric saucer shapes observed in sills. Such an investigation may lead to a method to constrain principal stress orientations at the time of sill emplacement from the sill geometry. In addition, this analysis may provide the way forward for constraining the magnitudes of the stresses. This proposed method would rely on estimating χ for a given sill by comparison with numerical or analogue model results and independently estimating the rock fracture toughness and initial emplacement depth.

## Acknowledgments

Funding has been provided by Australian Coal Association Research Programme (ACARP), Project C10010 with additional support from CSIRO Petroleum, the Theodore Bennett Chair and Schlumberger. Thank you to S. Polteau and Xi Zhang for helpful discussions and for bringing several references to our attention.

## Appendix

### Appendix: Analysis of intrinsic timescales

The problem consists of finding the fracture opening, *w*, and fluid net pressure, *p*, taken as the difference between the fluid pressure, *p*_{f}, and the overburden stress, σ_{o}. The location of the fluid and fracture fronts, *R*_{f}, and *R*, is also a part of the desired solution. Radial symmetry is assumed and, further, the fracture size is taken to be large relative to the dimension of the source so that injection is taken to be from a point source at the fracture centre. The solution is a function of time, *t*, the radial co-ordinate, *r*, the presumed-constant fluid injection rate, *Q*_{o}, the dynamic fluid viscosity, μ, the fracture toughness, *K*_{Ic}, and plane strain modulus, *E*′. In addition, we consider the solid to be prestressed with normal stress components σ_{r} and σ_{o}, as shown in Figure 3. Note that for a non-curving fracture one may consider σ_{o} as a surface traction, as drawn, or equivalently, as an overburden stress due to the weight of the overlying material (Zhang *et al.* 2002). However, when the fracture curves to approach the free surface, the variation of the stress with the depth may be important. For simplicity in clarifying some basic qualitative characteristics of these fractures, the scaling considers the mathematical model for a non-curving fracture propagating parallel to the surface. Also, equations are kept tidy by using the alternative parameters:

Letting the fluid flux, *q*, be described by the well-known Poiseuille equation:
4
and requiring local mass balance for the incompressible fluid
5
one arrives at Reynold's lubrication equation
6
For this problem, Reynold's equation is coupled with linear elasticity, which can be formulated as an integral equation embodying the superposition of appropriate fundamental solutions (Korsunsky 1994; Hills *et al.* 1996). For brevity it is expressed here in the general form:
7
The formulation is completed by the boundary conditions and moving boundary equations, given by (Savitski & Detournay 2002; Bunger & Detournay 2007):
8
9
10
11
where the overdot indicates differentiation with respect to time and *p*_{v} is the pressure in the lag region that is likely to be associated with gasses exolved from the magma or fluids that have entered from the host rock. The propagation condition (11) states the assumption that fracture propagation implies *K*_{I} = *K*_{Ic}, where *K*_{I} is the mode I stress intensity factor characterizing the near-tip stress field. In some cases *K*_{I} can be computed as a weight function integral (Rice 1968), while in many computations it is found from the asymptotic shape of the crack tip (Thomas & Pollard 1993). Here it is expressed generally in terms of the functional 𝒢{·}. It is also useful to note that the condition of the fluid front velocity (10) can be alternately accounted for using global volume balance:
12

The main points of the scaling of the governing equations are described below, leaving the details to Bunger (2005). The analysis begins with a substitution of the form:
13
Upon substitution one obtains non-dimensionalized governing equations that depend on seven dimensionless groups that contain the parameters {*t*, *Q*_{o}, μ′, *K*′, *E*′, σ_{o} − *p*_{v}} and the characteristic quantities {*W*, ε, *L*, ϑ}. Choosing the characteristic quantities involves setting four of the groups to unity. One choice (Savitski & Detournay 2002) leads to *W*=ε*L*, ϑ=1, ε=ε_{m} and *L*=*L*_{m} with:
14
The scaled governing equations and boundary conditions are then given by
15
16
17
18
19
20
The remaining groups, here named {ℳ, 𝒮, ℋ}, are evolution parameters that can be expressed as a power of the ratio between *t* and a characteristic time, namely:
21
where
22

Each of the evolution parameters can be associated with a physical transition that the fracture undergoes as it grows, and, as such, the characteristic times give an estimate of the ranges of time for which certain physical processes are expected to be important. Inspection of the scaled propagation condition (20) shows that 𝒦 can be considered as a dimensionless toughness. One can see, then, that the effect of the fracture toughness vanishes at early stage so that viscous flow is the dominant energy-dissipation mechanism. It is instructive to also consider a different choice of scaling (Savitski & Detournay 2002), similar to the one above but taking ε=ε_{k} and *L*=*L*_{k} with:
23
Reynolds equation is then written as:
24
where the subscript *k* is added to denote re-scaled quantities. Now a dimensionless viscosity:
25
appears and it becomes apparent that the pressure gradient ∂Π_{k}/∂ρ̂ vanishes when ℳ goes to zero, which is equivalent to 𝒦 going to infinity. Note that it is important to the validity of this analysis that the ‘M’-scaling (14) ensures that {Ω, Π, γ_{f}, γ} are order 1 quantities when 𝒦→0 and 𝒮, ℋ→∞. Similarly, the ‘K’-scaling (23) leads to {Ω_{k}, Π_{k}, γ_{fk}, γ_{k}} that are order 1 quantities when ℳ→ 0 and 𝒮, ℋ→∞. In fact, the reliance of this analysis on these properties of the scalings motivates the introduction of the dimensionless groups and the time-dependent scaling quantities in (14) and (23).

Some may find the association of the viscosity-dominated regime with small time counterintuitive because the length of the channel through which fluid is flowing increases with time. However, this analysis essentially shows that the increase in the opening *w* with time actually causes the effect of viscous dissipation to diminish in spite of the increase in radial extent. Equivalently, one can consider the transition from the viscosity-dominated to the toughness-dominated regime in radial fractures with a constant injection rate to result from the fact that the fracture-tip velocity decreases with time (Abé *et al*. 1976).

The remaining evolution parameters can be evaluated in a similar manner. The evolution parameter ℋ, which is like a dimensionless fracture depth, appears in the elasticity equation in such a way that the full-space problem is recovered as ℋ goes to infinity; that is, when time is small relative to *t*_{mm̄}. Conversely, ℋ becoming small corresponds to the case when the fracture is large relative to its depth. So the fracture changes from an effectively deep fracture to one which is effectively shallow with the characteristic time *t*_{mm̄}.

Finally, the physical transition associated with the parameter 𝒮 can be explained by considering a weight function integral representation of the propagation condition (20) for a deep fracture (ℋ→∞) (Rice 1968):
26
Employing the ‘M’-scaling (14) and taking the small-toughness limit 𝒦→0, we have (Bunger 2005)
27
where ξ_{f}=*R*_{f}/*R*. As previously discussed, the scaling ensures Π=*O*(1) when 𝒦≪1 and ℋ≫1, so it follows directly that satisfying (27) as 𝒮 goes from zero to infinity requires ξ_{f} to approach 1. In other words, the lag region between the fluid and fracture fronts diminishes as time increases relative to *t*_{om}. Physically, this transition corresponds to the net fluid pressure reducing in magnitude relative to the negative net pressure in the lag. Because the net pressure in the lag resists fracture extension, like a clamping in the tip region, the relative reduction of the net pressure will cause the fracture to stop growing pending sufficient increase of ξ_{f}. However, it is important to realize that the fluid lag not only diminishes as time increases relative to *t*_{om}, but also the lag diminishes along with the importance of viscous dissipation (Jeffrey 1989; Garagash 2005; Bunger & Detournay 2007). So the lag is expected to vanish when time is large relative to either *t*_{om} or *t*_{mk}.

- © The Geological Society of London 2008