Nomenclature

E

Young's modulus

E′

plane strain modulus

g

gravitational acceleration (9.81 m s⁻²)

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⁻¹. 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−ν²), 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Set-up for analogue experiments in PMMA and glass, where fracturing fluids were aqueous solutions of blue food dye and either glycerin or glucose, and the depth of the initial flaw varied 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).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Images of growing fractures for two tests: (a) showing the presence of significant fluid lag, after Bunger et al. (2005); and (b) illustrating elongation of an initially circular fracture to become egg-shaped.

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:

1. 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;

2. 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; and

3. The fracture size increases appreciably relative to its depth as t increases relative to t_mm̄.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Sketch of problem geometry, showing a cross-section of circular fracture.

Consider a sill growing at 2 km depth. Take the rock density as 2500 kg m⁻³ 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⁶ 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 𝒦.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Relation between elongation and the dimensionless toughness evolution parameter 𝒦. (a) Evolution of the elongation ratio R_max/R_min, with the mean radius to crack tip-depth ratio R/h for tests with differing mean 𝒦. (b) Correlation between R/h at onset of elongation and 𝒦. (c) Correlation between slope of R_max/R_min v. R/h after onset of elongation with 𝒦.

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Crack paths indicated by the crack-tip depth, h, and fracture radius, r, normalized by initial depth, H, for experiments performed in: (a) PMMA; (b) glass; and (c) a comparison of the results for the two materials.

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 β=ρ_ggH^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.