Confidential manuscript submitted to JGR Solid Earth The roles of elastic properties, magmatic pressure, and tectonic stress in saucer- shaped sill growth

Near-surface igneous sills commonly exhibit saucer-like shapes, formed due to interaction with the Earth’s surface once some critical sill length is reached relative to its depth. Sill geometry has been strongly linked to the host material conditions, particularly in terms of the elastic properties and shear cohesion of the host rock, operating as primary controls on sill geometry. Here we conduct dynamic numerical simulations for sill growth in the near surface, in which we vary the host rock properties, magma pressure profile internal to the sill ∆PP, and the externally applied tectonic stress σσrr, to consider their contributions to sill geometry. We find that elastic properties alone have little impact on sill geometry. Saucer shapes reflect the additive stress components of the magma overpressure within the sill, and the tectonic stress, controlled by the locus of the maximum energy release rate GGmmmmmm. Initially GGmmmmmm is in-plane with the sill, but deflects to ~25° at a critical base length rrcc relative to depth, due to interaction between the sill and the free surface. Increasing σσrr decreases this angle and increases rrcc of the sill. ∆PP controls sill growth rate, but has little effect on overall geometry. Host rock cohesion and elastic properties control the absolute magnitudes of σσrr required to affect a change in sill geometry. Plain Language Summary Horizontal magma pathways – sills – are a crucial part of volcanic plumbing systems, acting as potential feeder conduits to volcanic eruptions, and as magma storage systems. Saucer-shaped sills, which exhibit a flat inner region and inclined outer region, are a common type of magma pathway in the shallow crust. Models for saucer-shaped sills, both as physical analogue models, or numerical simulations, typically under-predict the length of the inner flat region, and overpredict the outer inclined region; models are typically too short and too steep. Here we use numerical simulation to investigate parameters that may control sill shape. We find that the dominant controls on sill shape are the competing effects of: (1) bending of the rocks above the sill, which promotes a transition to inclined growth, typically at ~25°; and (2) plate tectonic shortening, which serves to decrease the angle of incline, towards 0° when the horizontal force is high. Increasing the applied horizontal tectonic force can produce sills that are up to five times longer in the inner region, before growing as inclined sills at ~5°. This matches very closely with observations of natural sills, indicating that tectonic forces are an important consideration in the growth of sills.


Introduction
Igneous sills have played a primary role in the formation of Earth's oceanic crust (Maclennan, 2018) and the growth of evolved continental crust (Jackson et al., 2019). Sills have been shown to represent an important part of the plumbing system in volcanically active regions (Amelung et al., 2000;Sigmundsson et al., 2010) and in ancient systems (Leat, 2008;Muirhead et al., 2014). Sills are common in association with horizontally-layered sedimentary rocks (e.g., in basin settings), and typically considered as initially flat-lying intrusions that are parallel to bedding, representative of extension (mode I) cracks. This association between sills and horizontal layering has led to numerous models that invoke layering as the primary cause of sills; typically involving a vertical feeder dyke that instantaneously transitions to a layer-parallel sill at some layer interface (Kavanagh et al., 2006;Galland et al., 2009;Gudmundsson, 2011). Gudmundsson (2011) summarises three mechanisms for plane-deflection in such cases: (1) Cook-Gordon delamination, in which tensile stress ahead of the ascending dyke tip causes the horizontal interface to open, and eventually being intruded; (2) stress barriers, in which the principal stresses are rotated by layers of varying mechanical properties; and (3) elastic mismatch, whereby dykes ascending through materials with a low Young's Modulus may meet a material of high Young's Modulus, causing deflection (e.g., Kavanagh et al., 2006). Notably, in all such models, it is assumed that the ambient stress state for sills is a modified form of that required for dyke emplacement; either the maximum compressive stress, 1 , is vertical, with the intermediate and minimum principal stresses ( 2 and 3 respectively) being horizontal, or the principal stresses are equal. This inference is important as it requires the maximum and minimum principal stresses to switch, to facilitate the transition from vertical dyke emplacement to horizontal sill emplacement. As an alternative, horizontal extension and dyke emplacement could occur if 1 and 3 are horizontal; forming sills would, therefore, require the maximum and intermediate principal stresses to switch (e.g., Walker, 2016;Stephens et al., 2017). Although the effect of an ambient tectonic stress state has been investigated for dyke growth for several decades (e.g., King Hubbert and Willis, 1957;Maccaferri et al., 2010), few studies have focused on the tectonic stress state control on sill propagation and geometry (see Bunger et al., 2008).
Although sills are typically considered as layer-parallel and horizontal, sills that cut bedding are observed, at a range of scales (metre to kilometre), including: (1) gently transgressive sills inferred to relate to horizontal shortening (tectonic) strains, which accommodate mode I-II extensional-shear opening (Fig. 1;Walker, 2016;Walker et al., 2017;Stephens et al., 2017Stephens et al., , 2018; and (2) saucer-shaped sills, defined by a broadly convex profile that consists of a flat inner region, with transgressive flanks ( Fig. 1; Malthe-Sorenssen et al., 2004;Planke et al., 2005;Polteau et al., 2008;Hansen et al., 2011). In the latter case, the characteristic saucer shape, with transition from a flat to inclined sill, is typically related to asymmetry in silltip stress during propagation, due to inflation-induced bending of the overburden (Malthe-Sorenssen et al., 2004;Galland et al., 2009). Other proposed mechanisms for saucer-shaped sill formation include the depth-dependent increase in Young's Modulus (Hansen, 2015), or shear failure of the overburden (Haug et al., 2017(Haug et al., ,2018. Analogue models (Galland et al., 2009;Galland and Scheibert, 2013) and numerical simulations of sill growth (Chen et al., 2017) using isotropic host materials, have demonstrated that the saucer shape is a fundamental geometry.
Transition from the flat inner region to inclined outer region of natural saucer-shaped sills (e.g., Figs. 1 and 2) typically occurs at some critical radius length scale ( ) relative to the initial depth of emplacement ( ), typically where > 3 (Fig. 3). The incline section of natural sills typically rises at an angle < 30° (Fig. 1), which initiates at a low angle and increases in dip towards the extremities, defining a convex profile . By contrast, and in general, most modelled saucer-shaped sills differ from natural sills in three critical ways: (1) ≪ 3 (Galland et al., 2009); (2) ≥ 30° (Haug et al., 2017); and (3) sill incline profiles are concave, with initially steep, before decreasing at distance from the sill centre (Malthe-Sorenssen et al., 2004;Galland et al., 2009;Haug et al., 2017Haug et al., ,2018. Although recent models for saucer-shaped sills have investigated the role of layering in controlling sill geometry (e.g., Chen et al., 2017), few studies have investigated the role of tectonic stress as a geometric control (see e.g., Bunger et al., 2008). Here we use axisymmetric finite element (FE) simulations to consider the effects of the magma pressure profile within sills, and of horizontal stress on sill growth rate, dimensions, and geometry. To isolate the effect of Confidential manuscript submitted to JGR Solid Earth these parameters from those of layer-controls, we use an isotropic host material. We show that deviatoric stress (here, where = ≠ ) is important in controlling intrusion geometry. Figure 1. Natural sill geometries. (A) Flat to gently-transgressive sills in the San Rafael Sub-Volcanic Field, Utah USA, cuts the formation boundary separating Carmel Formation mudstones and siltstones below, and Entrada Formation sandstones and siltstones above (Stephens et al., 2018). (B) The Golden Valley saucer-shaped sill, South Africa, emplaced at ~2.4 km depth into Permian-Triassic age, sub-horizontal sedimentary rocks (Polteau et al., 2008). The distance between the white stars is ~7 km (GoogleEarth image; map data: © CNES/Airbus, © Landsat/Copernicus, © 2019 Digital Globe). (C) Saucer-shaped sills in the Faroe-Shetland Basin, imaged in seismic datasets (after Moy and Imber, 2009;Moy, 2010). Left-hand image is a 3D two-way-time (TWT) map of the sill top reflector. Middle and right-hand images are 2D topsill time structure and dip-azimuth maps respectively.  Comparison of natural and modelled sill profiles. Axes of sill radius and depth ℎ are normalized to the depth of emplacement of the initial sill, for comparison (see Bunger et al., 2008). Profiles shown are for analogue experiments of (A) Galland et al. (2009) and (B) Bunger et al. (2008), and numerical result of (C) Hansen (2015) and Haug et al., (2017). Panel (D) shows all datasets together. Solid black profiles are for the Golden Valley Sill and Tulipan Sill, as labelled in D.

Previous analytical and numerical models for sills
There is a rich history of study into the emplacement of dykes and the dyke to sill transition, via analogue experiments Takada, 1990;Menand and Tait, 2002;Menand et al., 2010;Kavanagh et al., 2015Kavanagh et al., , 2018Pollard and Townsend, 2018), and numerical models (Pollard, 1973;Gudmundsson, 2011;Barnett and Gudmundsson, 2014;Townsend et al., 2017).It has long been proposed that the geometry of dykes and hydraulic fractures is strongly controlled by the ambient regional stress field (e.g., Anderson, 1936;King Hubbert and Willis, 1957;Sibson, 2000;Maccaferri et al., 2010). Dykes are commonly inferred to represent magma-filled mode I (extension) fractures that accommodate crustal extension. The dyke plane in such a model lies parallel to the vertical plane containing 1 and 2 , and therefore normal to 3 . In contrast, if sills represent extension fractures, the plane containing 1 and 2 should be horizontal, and with a vertical 3 axis. Dykes and sills are commonly found in close spatial association, and it is assumed in many cases that dykes feed sills, even if the dyke is not observed.
Several mechanisms have been proposed for dyke to sill transition, such as the level of neutral buoyancy (Francis 1982), the effect of dyke-related volume increase on local stresses (Anderson, 1936;Wyrick et al., 2015) and several layer-related controls such as weak interfaces or contrasting layer properties (Kavanagh et al., 2006;Gudmundsson, 2011), and weak host rocks (Schofield et al., 2012). No studies have reproduced sills from dykes at a level of neutral buoyancy; rather, on reaching such a level, dykes propagate laterally but remain vertical (Lister and Kerr, 1991;Takada, 1999). Layering is typically considered the most important control on dyke to sill transition. In most cases, models for sills are conducted in initially isotropic conditions (i.e., no layering, and no applied stress), with the exception that the sill is typically seeded at some initial depth by a material interface (Galland et al., 2009) or existing crack (Bunger et al., 2008), or the sill already exists in the model (which is the case for our simulations). Initial stresses in the models are typically hydrostatic (i.e. = = , where is the vertical stress, and are the horizontal stresses) rather than lithostatic ( = = =

1−
, where is Poisson's ratio, the ratio of lateral to axial strain; = , where is the rock density, is the gravitational constant 9.81 / 2 , and is the depth). A hydrostatic case is valid if either = 0.5, or the deviatoric stress has been removed by creep. Considering the association of dykes as feeders to sills, and the general association of dykes with horizontal extension, this lack of an applied tectonic stress within models is notable. Figure 3 shows sill profile results from a number of analogue studies (Bunger et al., 2008;Galland et al., 2009) and numerical studies (Hansen, 2015;Haug et al., 2018), with comparison to examples of well-constrained natural sill geometries. Galland et al., 2009 used silica powder and vegetable oil as the host rock and magma analogues respectively. Bunger et al., (2008) used glass or polymethyl methacrylate (PMMA) as the rock analogue, and glycerine-water as the magma analogue. Their tests involved a deviatoric stress, applied as a crack-parallel radial stress during sill growth, and introduced a dimensionless parameter χ to consider scaling, where = √ / ; is the horizontal stress in excess of , is the initial depth of the sill, and is the fracture toughness of the material. The fracture toughness was 1.3 and 0.6 MPa m 1/2 for their PMMA and glass respectively. Bunger et al. (2008) showed that the application of a radial compressive stress during analogue sill models, in which = > reduced the effect of anisotropic stress ahead of the propagating tip, producing saucer-shaped sills with inclines 30°> > 5°, more closely resembling natural saucer-shaped sills (Fig. 3). Haug et al., (2017) simulated static cracks in a homogenous Mohr-Coulomb material (Fig. 3C); individual plots show the damage associated with cracks of different starting lengths, each with a starting depth = 2 and sill thickness = 50 . Hansen (2015) calculated inflation sequences ( ) for elastic displacements above ( ) and below ( ) a sill intrusion centre line for isotropic host rocks as a function of the depth-dependent increase in Young's Modulus , where = + 2 ; is simply a multiplication factor used to vary . Results are shown for his models as a function of , with a magma overpressure = 10 ( Fig. 3C). However, in general we note that for modelled sills (1) the critical sill radius is too short, (2) is too high, and (3) the sill profiles are mostly concave (Fig. 3).

The Model
We employ a dynamic model that allows a sill to determine its own path through the crust as it extends over time. Three mutually dependent physical processes are involved in the sill intrusion process: (1) elastic deformation of the host rock, which varies as a function of the rock elastic properties, the overburden tectonic stresses, and the magma overpressure distribution within the sill; (2) fracture of the host rock; and (3) viscous flow of magma through the fracture. These are represented by partial differential equations which we solve here via the finite element (FE) method, using the commercial software package COMSOL Multiphysics v5.3. An advantage of this method is that it is particularly well suited for solving complex, strongly-coupled, globally-connected problems. In addition, and unlike several previous models for sills, our models represent dynamic simulations of growing sills, rather than static-sill-related stress or strain.
Our sill models are simulated in axisymmetric space (Fig. 4). An initial sill intrusion is introduced at a depth = ( Fig. 4; nomenclature is summarized in Table 1). The lithostatic pressure ( ) is assumed to be known and hydrostatic ( = = ). Therefore, only the unknown deviation from the lithostatic stress state, Δ , is calculated in the simulation. The actual stress is the sum of the two components: where the indices , = , , and is the identity tensor (1 if = and 0 otherwise). Compressive stress here is reckoned positive; tensile stress is negative. This separation of the stress state into two components allows better representation of the conditions at the sill tip, which is critical for modelling the evolution of sills. In general, the magma pressure in the sill and the applied tectonic stress will enhance the stress at the sill tip, whereas the lithostatic pressure will not.

Failure Criteria
The sill intrusion process is analogous to the growth of a crack filled with a pressurised fluid (Anderson, 1936;King Hubbert, 1951;Bunger et al, 2008). Here it is assumed that the stress state near the sill (crack) tip cannot be resolved accurately, partly because the appropriate lower length scale cannot be resolved in a computationally tractable model without adaptive meshing schemes (e.g., Shen and Lew, 2014), which would be more difficult to employ for an evolving crack along which a fluid flows, but also because the exact nature of the local crack tip conditions is not well understood. For example, heating due to the magma may lead to localised non-linearity (e.g. viscoplasticity; Souche et al., 2019), thermoporoelasticity (Zoback, 2007, p. 83), and water vaporisation in wet rocks (e.g., Pollard et al., 1975;Schofield et al., 2010;Galland et al., 2019), and the true dimensions of a tip lag region cannot be resolved. The stress at the sill tip is expected to be considerably larger than the stress resolved in a discrete continuum model with a fixed mesh size. Sill propagation depends on the maximum stress, and the finite resolution of the model only averages the stress over a region comparable in size to the sill thickness.
Here it is assumed that the actual stress at the sill tip is proportional to the stresses at the tip resolved in the model through a constant > 1, such that the actual locally enhanced stresses at the sill tip are given by: The constant therefore represents a calibration between the FE model and a real material. Two failure criteria for the rock at the sill tip are adopted. For brittle tensile fracture, we assume that failure occurs when: where 3 is the minimum compressive stress (where principal stresses 3 ≤ 2 ≤ 1 ) and < 0 is the tensile failure strength of the material. Using Equation 1 gives the equivalent condition in the model as: To induce sill growth in the model, a fixed source-overpressure is applied to the root of the initial intrusion at = 0 ( Fig. 4). In general, the overpressure is not known, as the value required to initiate sill intrusion depends on the depth at which this occurs, the precise physical conditions at the sill tip (for instance tip geometry, and temperature) and the fracture toughness of the rock. In this model, the choice of can be made arbitrarily without affecting the final outcome; it is simply a reference stress against which the relative magnitude of the other stresses in the model are defined. When no tectonic stress is applied, the observed stress required to initiate tensile failure in the model will therefore be proportional to this overpressure, i.e. Δ 3 = − , where is a constant.
The sill starts at a depth so that at first failure Equation 4 gives a sill tip stress enhancement factor of = ( − ) , and the tensile failure condition becomes: which reduces as the sill climbs towards the free surface (at = 0). The FE superscript is used to indicate properties that are calibrated for the FE simulation. They are not to be confused with the actual physical parameter itself, although these are used to determine the parameter in the FE simulation relative to the chosen magma overpressure .
Although the sill generally extends as a mode I crack, the rock at the sill tip may also fail by micro-cracking along lines of shear according to the Mohr-Coulomb relationship where the failure shear stress depends on the cohesive strength , the angle of internal friction , and the normal stress acting on the plane. The optimal slip condition is: where the maximum shear stress max is Eq. 8 and the hydrostatic stress max is Therefore, the Mohr-Coulomb condition is written as: As in the tensile failure case, this is also scaled relative to the overpressure. For an assumed friction angle of = 30 , this gives The critical failure stresses for determining the initial growth and onset of climbing for a sill are those determined at the starting depth . Given that typically ≪ , the ratio of the shear failure stress to the tensile failure stress, , can be expressed for a general friction angle as Eq. 12 This factor is a function of depth and determines whether the shear (small ) or tensile (large ) failure criteria is dominant.

The damage model
Rock fracture is simulated by a crack-band damage model (e.g., Bazant and Oh, 1983), whereby the damage parameter 0 ≤ ≤ 1 characterises the level of fracture at a point, with = 0 representing fully intact rock and = 1 representing completely failed rock. This method allows the sill to determine its own path, without placing prescriptions on its growth. Damage typically only occurs in elements adjacent to the sill tip where the stress is enhanced (by factor ). The damage parameter evolves according to: where the rate constant is = 1000 −1 if either failure condition (Equation 5 or 11) is satisfied in the tip region, and zero otherwise. The value of is capped so it cannot exceed one. The choice of the rate constant is not important, as long as it is large enough to ensure that fluid flow dictates the sill growth rate, not the rate of damage. The damage affects the linear elastic analysis through adaptation of the Young's modulus (e.g., Heap et al., 2010) Eq. 14 where is the modulus of the unfractured rock. This means that the fractured rock has no mechanical stiffness once it has failed, i.e. it represents the interior of a crack. The failure conditions, Equations 5 and 11, are determined using stresses calculated directly from the strains using the undamaged modulus, . This is necessary as the reduced modulus of Equation 14 causes the stresses in a damaged region to decrease causing failure to halt, whereas the strains increase (as the crack opens) allowing failure to continue.

Magma pressure and fluid flow
We assume that the magma exhibits laminar flow down a naturally-evolving pressure gradient in accordance with standard lubrication theory (e.g., Rubin, 1995) such that the fluid flux is given by where is the width of the magma channel and is the magma viscosity (Witherspoon et al., 1980;Chen et al., 2017) and the only non-zero contribution from the acceleration is that due to gravity such that = − . We define the total pressure as = + Δ . Mass conservation requires an incompressible fluid satisfies: This stiff system is difficult to implement numerically and it is more practical to consider the magma as a slightly compressible fluid such that where is the bulk modulus of the fluid and the parameter is zero until an element is fully damaged, at which point it is set to unity. The introduction of therefore restricts the fluid flow to the damaged (intrusion) zone. The magma overpressure Δ is added to the stress state of the simulation, where here ensures that the overpressure Δ is only applied to the magmatic region. A positive (upward) body force of = Δ is also added directly to the magmatic body, to model the buoyancy force of the magma, where Δ = − is the difference between the density of the rock ( ) and that of the magma ( ). A positive buoyancy acts to widen the channel and increase flow to the tip. Rubin (1995) has shown that this makes only a small contribution relative to the magma overpressure for intrusions such as those considered in this paper which only rise a distance of the order of a kilometre or so through the crust.
The growth of the sill is strongly coupled to the magma pressure and host rock elasticity (cf. Pinel et al., 2017) through the dependence of the channel width ( 3 ) on the elastic opening of the channel, where 3 is the maximum (elongation) principal strain in the damaged (sill) region. The magma flows to accommodate this strain and fill the enlarged channel (note then that the strain refers to the channel opening, and not the magma). The exact conditions that drive magma through intrusions are not fully understood. Two common approaches in modelling are: (i) constant flux of magma from the source and (ii) constant magma pressure at the source. These two bounding cases are both explored here. Rubin (1995) derived the profile for the pressure along a 2D vertical dike for constant source pressure. For short intrusions (<1km) they ignored buoyancy, and as such we consider the analysis is equally valid for 2D horizontal sills. That study postulated the presence of a cavity at the intrusion tip, and investigated the effect of the pressure in that region. It was shown that such an intrusion propagates with a velocity proportional to its length; hence constant source pressure is consistent with an exponentially increasing intrusion length that can be attributed to an ever-widening channel. This is considered to be unrealistic for longer dikes (>1km in height) and it is assumed that the availability of magma is limited once the vertical extent of the intrusion exceeds this length scale: i.e. the model proposes a switch from constant source pressure to constant flux.
For constant magma flux, the intrusion velocity will decrease over time. This is consistent with the concept of an intrusion eventually halting due to cooling and solidification, whereas the accelerating growth of a constant source pressure intrusion would be unlimited in this respect. This suggests that these two cases are extreme examples which define the limits of what is likely to be physically correct. Here, we consider both cases by describing the channel width as = 0 ( 0 + 3 ), where 0 is a characteristic channel width; the choice of this parameter is not important geometrically, and simply determines the time scale. For constant flux it is assumed that the channel opening is fairly insensitive to strain, such that 0 = 1. This effectively leads to a constant source flux model and is consistent with a sill that has a blunt tip of width 0 (even if there were to be a sharp tipped cavity ahead of the magma). Constant source pressure is modelled with 0 = 0, which is equivalent to the source pressure model of Rubin (1995). The difference in the pressure profiles generated by the two models is illustrated and discussed in section 4.2, but sensitivity analysis (see supplementary files) shows that both models generate very similar sill profiles; i.e. the sill path is insensitive to the magma pressure distribution. For this reason, only the constant flux model is used in the calculation of subsequent results. This model is selected as it produces a sill that propagates slower as it gets longer, as opposed to the constant source pressure model, which as predicted (Rubin,1995), has an exponentially increasing growth rate, which is not considered to be physically probable.

Simulation parameters
The model has been implemented using the commercial FE software package COMSOL Multiphysics v5.3. Linear 4-noded square elements of side length = 1 50 (Fig. 4) are employed for all three physical processes, where is the initial depth of the intrusion. The simulations are axisymmetric and conducted within a large domain 8 deep by 16 radius to avoid boundary effects ( Fig. 4; see supplementary files for sensitivity tests of the domain size). An initial horizontal rectangular intrusion of width and radius 5 is introduced at the depth of . The Young's modulus ( ) of the rock is taken as 1 GPa, with a Poisson's ratio ( ) of 0.4, to simulate a weak mudstone (Hobbs et al., 2002), which is a common host for basaltic sills (Mark et al., 2019). Sensitivity tests for these parameters, in the range of 1-100 GPa, and of 0.1-0.5, show that only the upper extreme values (i.e., = 100 and/or = 0.5) result in any appreciable difference in sill geometry (see supplementary files). Magma viscosity ( ) is constant at 10 6 Pa⋅s and the bulk modulus of the magma ( ) is 0.1 GPa, in accordance with calculations for shallow (2 km) magmas with volatiles (Huppert and Woods, 2002). The rock tensile strength ( ) is -3 MPa, and the Mohr-Coulomb parameters are a shear cohesion ( ) and an angle of internal friction ( ) of 30°. The rock density ( ) is 2.5 g/cm 3 and the magma density ( ) is 2.3 g/cm 3 . This value for magma density would be reasonable for high silica-content magmas, but low for basaltic melts (Stolper and Walker, 1980;Bottinga et al., 1982); notably the 'magma' in our model has not traversed the column of denser crust and mantle lithosphere at depth, and we consider that a relatively low density is justified here. To simulate tectonic stress in the model, a horizontal strain = (1− ) is applied to the outer radius of the simulation, where is the applied horizontal tectonic stress excess to the lithostatic pressure. In our simulations, is applied in the range of 0-5 MPa via a horizontal radial contraction.
We take a reference magmatic overpressure ( ) of 5 MPa, with the two failure conditions defined by Equations 5 and 11. The failure stress required to induce tensile failure at the tip of the initial intrusion in the FE model ( ( ) = − ) is -6 MPa with this overpressure, giving = 1.2. At a depth of = 2 , is 50 MPa giving a stress enhancement factor at the crack tip of = 50+3 6 ≈ 9. This implies that the actual stresses at the sill tip required to activate tensile fracture are roughly nine times higher than those resolved in the simulation. This is to be expected given the actual tip shape may be much sharper than we can simulate here, and shows how this approach mitigates the approximation incurred by not modelling the sill tip in detail. The applied source pressure is constant in our models, with constant flux or pressure conditions enforced as stated in section 2.5, though it should be noted that neither the magma pressure, nor the magma flux should be expected to remain constant in nature (Rivalta, 2010;Schopa and Annen, 2013;Kavanagh et al., 2015); here, we generate variable pressure decay profiles, which are strongly coupled to the modelled sill growth.
The critical stress in the FE model required to instigate Mohr-Coulomb failure at the chosen overpressure in the absence of an applied tectonic stress is found to be 3.8 MPa, compared with -6 MPa required to induce tensile failure. The critical value of , defined by Equation 15, at which both shear and tensile failure are equally likely in the absence of tectonic stress is therefore = 3.8 6 = 0.63. The value of this parameter depends on the ratio of . We therefore consider two cases which cover the expected range of this parameter (Haug et al, 2018): (1) with a lower shear strength of = 8 MPa, Confidential manuscript submitted to JGR Solid Earth and = 2km, giving = 0.16 and = 0.63 (Fig. 4); and (2) with a higher shear strength of = 20 MPa, and = 1km, giving = 0.8 and = 1.19 (Fig. 5).

Magma pressure distribution along the sill
As noted above, sill growth in the models is determined by the distribution of maximum stress at the sill tip. The model is fully coupled in that the damage depends on the stress; the stress depends on the pressure distribution in the sill; the pressure distribution in the sill depends on the elastic deformation of the host material; and the elastic deformation depends on the stress, ad infinitum. In section 2.5, two magma flow models were introduced: (i) constant source flux and (ii) constant source pressure. The two produce similar sill profiles ( Fig. 5C) but different pressure distributions (Fig. 5A) and radically different propagation rate characteristics (Fig. 5D). Figure 5A shows the two magma overpressure distributions (Δ ) within a flat intrusion subjected to zero tectonic stress. The constant source pressure at the root of the sill is fixed at = 5 ; Δ therefore refers to pressure within the conduit. The pressure drop along the sill causes magma to flow towards the tip. The pressure profile is convex for case (i) constant source flux and concave for case (ii) constant source pressure and comparable with Rubin (1995) (i.e., their figure 6, although this is for 2D planar sills not the axisymmetric sills modelled here). It can be seen that a small negative pressure has developed at the tip of the intrusion for case (i) and a much larger negative pressure for case (ii). This indicates that the sill crack grows at the tip faster than magma can readily accommodate the increase in volume, causing the magma to be sucked in to the tip region (cf. Rubin, 1995). The negative pressure will act to close the crack tip, but it is held open by the much more significant action of the positive pressure in the sill behind it. Alternatively, this could be considered a condition for separation of the fracture front and the magma front, i.e. the crack grows faster than the sill producing a tip cavity, which is not possible in this model but may occur in scaled analogue models and nature (e.g., Rubin, 1993;Liss et al., 2002;Kavanagh et al., 2006;Bunger et al., 2008). Figure 5C shows that the sill profiles are not strongly influenced by the pressure distribution (a fact that is independently verified by the energy release rate model in section 4.2), with only a slightly shorter critical radius for the constant source pressure model. The time scale in Figure 5D is normalised so that both sill models reach about 5 km in length at the same time. In practice the constant source pressure model is three orders of magnitude faster than the constant flux model. The constant source pressure model exhibits exponential growth, as expected (Rubin, 1995). The constant flux model exhibits a constant volumetric growth rate, such that the magma volume increases linearly with time, commensurate with the sill radius increasing as 1 3 for self-similar sill shapes. Given the small difference between the sill profiles resulting from the two models, the constant flux model is adopted for the rest of the paper, as a slowing sill growth rate seems more probable than an accelerating one.

Material properties
To calibrate the simulation results against reality it is necessary to determine the effective fracture toughness of the rock in the FE model. This can be calculated from the critical energy release rate, . This method avoids the need to consider the details of the stress field around the sill tip. The critical energy release is defined as = , where is the total change in elastic energy associated with an increase in the crack area . The energy release rate in the simulations was constant, indicated by the constant slope in the versus plot in Fig. 5B. This shows that the magma pressure distribution evolved to keep the stress intensity at the crack tip constant, i.e. = . This indicates that the sill growth is controlled by the toughness of the rock rather than the viscosity of the magma. The assumption that intrusions grow under constant is supported by the crack length versus thickness calculations of Scholz (2010) for dykes and veins (from Schultz et al., 2008, and references therein). The effective (mode I or II) fracture toughness can be inferred from the standard elastic relationship (Zehnder, 2012) = 2 Eq. 18 yielding = 26 √ , where the FE superscript has again been employed to demonstrate that this is a simulation parameter not a physical one. This is verified via a second method. We also use the stress intensity factor for a penny-shaped crack in an infinite solid subjected to a varying internal pressure distribution = 4 √ ∫ ( ) √ 2 − 2 0 (Sneddon, 1946) to generate an estimate for the stress intensity at the sill tip, where is the current crack length. The constant flux pressure profile in Fig. 5A yields a constant value for the stress intensity factor during the extension of = = 25 √ hence the two methods for calculating are in good agreement. Figure 6 shows profile plots which summarise the effect of a range of material properties for sills emplaced at 2 km, in which rock and magma properties, including , , , , are set as individual variables. In each case the material parameters are as described in Section 2.6, with the exception of the displayed parameter. In all cases, the sills propagated as a flat sheet to some critical radius ( ), before climbing upwards towards the surface. The variation of each material parameter generally imposes a minor control on the geometry in terms of (~1.25H +/-0.1H), and sill incline (~23° +/-2°). Figure 6A shows that increasing the host rock Young's Modulus can result in a flattening of the sill incline on the order of ~5° from = 1 GPa and 50 GPa; is also seen to increase by 0.25H across this range in . Young's modulus also has an effect on the sill thickness and roof uplift (Fig. 7), with increases in decreasing the thickness of the sill and amount of uplift as it propagates (Fig. 7C versus 7E). Sill thickness and surface uplift are not significantly affected by magma buoyancy. Numerous numerical tests were also conducted to explore the influence of other parameters, including the magma bulk modulus and mesh size ( = � , , � ), but none of these parameters significantly affected the overall sill profile geometry.   Figures 8 and 9 show simulation results for sills emplaced at 2 km and 1 km respectively; in both cases the key variable is . In both experiments, sills intruded under lithostatic conditions resulted in the initial growth of a flat sill until some critical radius (1.81-1.25H: Table 2, Figs. 8-11), after which the sill climbed at a consistent angle of ~23-24°. Increasing in both experiments resulted in an increase in , to a maximum of ~5 , and a decrease in , to ~1° (Table 2,  . In the lower shear strength case ( = 0.16 and = 0.63), increasing above 4 MPa resulted in the onset of a decrease in and an increase in (Figs. 8, 10A, and 11A,B). Increasing and appears to delay the onset of this transition beyond the range of magnitudes tested here (Figs. 9, 10B, and 11C,D). Figure 9 demonstrates a transition in the overall sill profile as a function of . For experiments with relatively low host cohesion (Fig. 9A), sills emplaced where < 4

Tectonic stress
show a concave profile in which decreases as the length of climb increases; the curve of the profile is a 4 th -5 th order polynomial. Where ≥ 4 , sills show a convex profile, in which increases with increasing length ; the curve simplifies to a 3 rd order polynomial. For experiments with a relatively higher cohesion host, and at shallower depth, the transition from convex to concave occurs later, with only the = 5 sill showing a concave form.      Bunger et al., (2008) conducted analogue experiments for sill intrusion, monitoring the axisymmetric crack radius and height over time as pressurised fluid was pumped into the crack at a fixed rate. Two materials were selected for the analogue rock (Glass and PMMA) with glycerine or glucose solution for the analogue magma. The sill was initiated from a flaw of 6 mm radius machined at the end of the injection tube, with sample depths in the range of 12 to 30 mm, i.e. = 0.2 to 0.5 initially. The effect of a horizontal stress, equivalent to , on sill profiles was investigated. They introduced a dimensionless parameter as a means to define scale-independent factors affecting sill profiles

Comparison with analogue models
where in the case of Bunger et al., (2008) is specifically the mode I fracture toughness of the host material.
To compare the simulation results with the analogue model of Bunger et al (2008) it is necessary to determine the effective fracture toughness of the rock in the FE model, . This functional form of the failure stress is given by = √ , where the initial depth of the sill, , is the characteristic length scale for the problem, is a constant geometric factor, and is the tensile failure stress in the FE model. It was shown in section 3.2 that = 26 √ , yielding a value of = 0.097 for our chosen FE mesh. Thus, given ( ) = − , we can write the effective value determined from the FE model as To generate results comparable with the analogue model of Bunger et al., (2008), the effects of lithostatic pressure and buoyancy are removed from the simulation, and failure of the host material is limited to mode I tensile fracture. The results are shown for PMMA and glass in Figure 12. It can be seen that the slopes of the sills are in very good agreement between the two models, although the difference in horizontal section length prior to this is considerable. There are a number of potential explanations for this difference: (1) the horizontal section in the analogue results is typically of the order of the initial flaw size (0.2 < < 0.5) suggesting that there is little to no horizontal growth in their models; (2) the analogue cracks initiate from a machined notch of finite thickness which may not induce an initial horizontal crack; (3) the analogue cracks are not completely axisymmetric; (4) the fracture is seen to extend ahead of the fluid in the analogue case, possibly due to the action of capillary forces which are not significant at the km scale; and (5) the FE model requires a finite perturbation in the stress to destabilise horizontal sill growth. The FE model of Equation 20 provides scaling behaviour for a blunt sill and predicts that the value of =8.6 corresponds most closely with the Tulipan Sill (see Fig  10A) for which it is expected that the tectonic stress has the same value as the magma source overpressure ( = ). The overpressure is not known accurately, but values of 0.5 to 6 MPa have been estimated elsewhere (Gudmundsson, 2012) with an average of 3 MPa and a maximum of 9 MPa. The higher end of this scale (≥ 3 MPa) would appear to be a reasonable differential relative to a lithostatic pressure of 50 MPa at a depth of 2km, as one would expect it to be larger than natural variations in the hydrostatic stress state (Suppe, 1985).
Using the expression of Bunger et al., (2008) from Equation 18, with a typical fracture toughness for sandstones at about 2 km depth of roughly 3 √ (Stoeckhert et al, 2016) predicts an applied tectonic stress of around 0.6 MPa for a value of = 8.6. This is at the lower end of the expected range given above, and is a very small stress in comparison to the additionally imposed lithostatic pressure of 50 MPa. This analysis is based on a crack tip that is sharp, whereas in reality the tip of the sill would not be as sharp as a crack. The FE model does not assume a sharp sill tip, leading to the more realistic value of 3-9 MPa predicted by Equation 20.

Energy release rates for different sill geometries
The intrusion simulations have shown that the model is quite insensitive to the material input parameters, with the sill incline demonstrating a strong preference for a particular orientation at a given tectonic stress. This suggests the elastic driving force due to the presence of the free surface and the tectonic stress are the primary controls on saucer-type sill shape. Here we employ a different method to separately explore and verify this response. The energy release rate of pressurised intrusion-like cracks as a function of their geometry and the pressure distribution within them is investigated for a fixed range of shapes. The axisymmetric sill geometry for a sill at depth is now described in terms of just three parameters, illustrated in Figure 2: the radius of the flat section of the sill, , the length of the inclined section of the sill, , and the climb angle of this section, . The normalised energy release rate is determined for different geometries and loading conditions for host rock with a Poisson ratio = 0.3. A sill is expected to propagate in the direction which maximises the energy release rate.

Effect of internal pressure distribution
It has been shown in the FE model of section 3.1 that different boundary conditions at the magma source results in different internal pressure distributions in the magma (e.g., Fig.  5A). However, the resulting sill profile is only weakly affected. This is investigated by proposing three different pressure distributions, each with smooth gradients; we do not consider abrupt pressure changes along the length of the sill, which would reflect significant viscosity variation within the conduit, which has been shown to have a marked effect on intrusion growth and geometry . We use the parameter 0 ≤ ≤ 1 to denote the location along the sill, with = 0 at the root (at = 0), increasing linearly over the total sill length of + , to = 1 at the tip. The pressure distributions are: (1) constant pressure, ( ) = 0 ; (2) linearly decaying pressure, ( ) = 3 0 (1 − ); and (3) quadratic pressure decay, ( ) = 6 0 (1 − ) 2 . The pre-factors have been chosen to keep the total force roughly similar between cases. These are expected to cover the extreme range of expected pressure distributions, with the quadratic pressure distribution being closest to that observed in the constant flux model (see Figure 5A), and a constant pressure being closest to that postulated for the constant source pressure model.
A sill is expected to grow at an inclined angle once the energy release rate is no longer maximum along the horizontal ( = 0 ). Figure 13A shows the normalised energy release rate, ̅ = 2 , as a function of for = 1.5 , a length at which the modelled sills are seen to initiate an incline. The energy release rate peaks at about 25-30 o suggesting this is the angle at which the sill would climb if there were no other forces acting on it (e.g., buoyancy forces). This is also typically the angle found in the sill simulations in the case of zero applied tectonic stress (Figs. 10 and 11). The angle of inclination is largely insensitive to the pressure distribution, which explains the insensitivity of the simulation to the various parameters that might affect this, such as buoyancy forces, magma bulk modulus, and the magma permeability tensor. Figure 13B shows the change in the normalised energy release rate, ̅ , when a flat sill begins to climb at 25° as a function of its flat section length; a positive energy change indicates that the switch is energetically favourable. It can be seen that the energetics for a switch are favourable for = 0.8 for the quadratic pressure distribution, rising to = 1.15 for the constant pressure distribution. In the FE model a sill length of = 1.5 is typically seen for zero tectonic stress (Fig. 11). This suggests that there is a finite energy barrier to be overcome to switch to inclined growth in the FE model, probably due to the sill having a blunt tip. In practice, there will also be localised stabilisation of horizontal natural sills due to host rock mechanical layering (e.g., Mudge, 1968;Johnson and Pollard, 1973;Burchardt, 2008).
In summary, the energy release rate calculations show that, in a linear elastic body, an axisymmetric sill will grow horizontally until its length exceeds 0.8-1.15 times its depth, which is dependent on the magma pressure distribution within the sill. In the absence of a tectonic stress, the sill will then climb at roughly 25-30 o to the horizontal (Figs. 8-11). This is consistent with the naturally evolving sill profiles resulting from the case (i) and case (ii) magma flow models, which are largely insensitive to the form of the pressure distribution (Fig. 5).

Effect of tectonic stress
A constant horizontal stress was applied to the remote vertical boundary, and the internal pressure on the crack removed. The applied stress will lead to a stress concentration at the crack tip, which if compressive will hinder the growth of the crack, and if tensile will assist its growth. The net differential stress (above lithostatic) is assumed to remain tensile, i.e. the pressure in the sill always acts to open the crack. Figure 10C shows the effect of the applied stress on the normalised energy release rate, ̅ = 2 , for a crack of length = 1.5 for different inclination angles =15 o , 30 o , 45 o and 60 o as a function of the inclined sill length . The lines are almost linear such that Equation 18 suggests that the normalised stress intensity at the sill tip due to the tectonic stress is where ℎ = is the height of the sill tip above its starting depth. This means that the applied stress has little effect on a horizontal sill, but that a compressive stress can stabilise horizontal sill growth by reducing the stress intensity at the sill tip as it starts to rise. This equation implies that the applied stress has greater effect as the sill grows towards the surface. This analysis does not identify how the mix of tensile and shear loading changes at the tip, with the contribution from the shear component expected to increase as the angle of incline gets larger. Figure 13D shows the dimensionless pre-factor � , � for two different sill base lengths. This shows the expected result that a tensile differential tectonic stress ( < 0) will act to drive an intrusion to climb at or near 90 o , as this assists the stress intensity at the sill tip most at this angle. In compression ( > 0), the reverse is the case, with the applied stress decreasing the stress intensity at the sill tip most effectively at 90 o , hence driving the crack towards 0 o as this reduces the stress intensity the least, allowing sills to continue to grow at these low angles.
These findings are summarised in the schematic illustration of Figure 14, which demonstrates that the differential stress (beyond lithostatic) at the tip of an inclined sill comes from two main contributions: (1) The internal magma pressure ∆ generates the majority of the stress that fractures the host rock, and drives the sill forward. The stress intensity at the sill tip is proportional to the magma overpressure times the root length, in which the characteristic length is the starting depth. It has been shown that in a lithostatic stress state, the asymmetric stress field at the tip drives the sill to climb at an incline of about 25 to the horizontal (Fig. 14A). This is the typical sill angle observed in models for sills conducted in the absence of a tectonic stress (e.g., Galland et al., 2009).
(2) The imposition of a horizontal tectonic stress creates an additional stress field that is superimposed on the tip (Fig. 14B), with the stress intensity proportional to the applied stress times the root of the height above the initial sill. If the differential tectonic stress is compressive ( > 0) then it has been shown that the stresses at the sill tip are reduced, with the maximum reduction applying to vertical intrusions and the minimum reduction for horizontal intrusions. A compressive applied tectonic stress therefore drives the sill to a shallower angle of incline, and a tensile tectonic stress drives the intrusion towards higher angles of incline.

Fig. 14.
Schematic illustration of the competing sill tip stress contributions: (A) the internal magma pressure has been shown to drive the sill towards the surface at an angle of about 25 to the horizontal once it has exceeded a length of ≈ . This is the dominant stress in driving the sill forwards. (B) A compressive applied differential tectonic stress (above lithostatic) acts to reduce the tensile stress at the sill tip caused by the magma pressure. The reduction is minimal in the horizontal plane (as shown) and hence acts to lower the angle of inclination of the sill.

Conclusions
Our numerical simulations of saucer-shaped sills show that sill geometry is generally insensitive to the host rock elastic properties and shear cohesion. Sill geometry is instead controlled by the competing contributions of the magma pressure profile within the crack, and the magnitude of tectonic stress. We applied two end-member magma flow models -constant source flux and constant source pressure -which produce very different magma pressure profiles within the conduit. Pressure distribution in the sill has little effect on the growth path of a sill, but results in markedly different propagation rates as a function of the sill length: sill growth slows as length increases in the case of constant flux, whereas constant source pressure models show an exponential increase in growth rate. In both cases, magma pressure within the sill drives growth, initially as a flat plane. On reaching a critical base length radius, relative to the emplacement depth, growth transitions to an inclined sheet, following the angle of maximum energy release rate at 25-30° to the horizontal. Compressive tectonic stresses subdue the effect of asymmetric tip stress generated by interaction with the free surface, allowing the sill to grow as a flat plane for greater lateral distances, and reducing the angle of incline once it exceeds its critical base length. Our results are in good agreement with analogue models of fluid-filled cracks presented by Bunger et al., (2008), indicating this is a scalable phenomenon. In our simulations, sills are most closely matched to natural examples, when the magnitude of tectonic stress above lithostatic is equal to the magma overpressure within the sill; applied tectonic stresses above or below this value result in steeper sills.