Segment tip geometry of sheet intrusions, I: Theory and numerical models for the role of tip shape in controlling propagation pathways

Inferences about sheet intrusion emplacement mechanisms have been built largely on field observations of intrusion tip zones: magmatic systems that did not grow beyond their observed state. Here we use finite element simulation of elliptical to superelliptical crack tips, representing observed natural sill segments, to show the effect of sill tip shape in controlling local stress concentrations, and the potential propagation pathways. Stress concentration magnitude and distribution is strongly affected by the position and magnitude of maximum tip curvature κmax. Elliptical tips concentrate stress in-plane with the sill, promoting coplanar growth. Superelliptical tips concentrate maximum tensile stress pσmaxq and shear stress out-of-plane of the sill, which may promote non-coplanar growth, vertical thickening, or coplanar viscous indentation. We find that σmax “ Pep1` 2 ? aκmaxq, where Pe is magma excess pressure and a is sill half length. At short length-scales, blunted tips locally generate large tensile stresses; at longer length-scales, elliptical-tipped sills become more efficient at concentrating stress than blunt sills. Non-technical summary Magma sheets underground grow in stages: small isolated sheets grow, link, and inflate to become larger/longer sheets. Preserved tips at the ends of sheets are inferred to represent these growth stages. However, a number of processes may alter the shape of sheet tips, which operate after growth. Here we model idealised shapes representing natural sheet tips, to determine the efficiency of those shapes—how likely they are to facilitate growth—and the likely pathway for growth based on the greatest efficiency. We find that growth efficiency and direction is linked to the magnitude and position of the maximum tip curvature. Elliptical tips have one maximum curvature at the tip, and will promote continued growth along the long axis of the sheet: small sheets will grow to become longer sheets. Squared-ellipses can be efficient, but a single tip displays two maximum curvature zones, which may change the


Introduction
Igneous sheet intrusions are a fundamental component of volcanic plumbing systems, representing conduits that supply volcanic eruptions, nascent magma storage systems in the crust, and fundamental building blocks of oceanic [Gass 1968;Maclennan 2019] and continental crust [Jackson et al. 2018]. Sheet intrusions that grow close to Earth's surface can induce large earthquakes and surface ruptures, in some cases leading to volcanic eruption [Wright et al. 2012;Sigmundsson et al. 2015]. The distribution of earthquakes during intru-Corresponding author: rich.walker@le.ac.uk sion, and the nature of surface deformations (e.g. doming, graben formation, and fissuring) is proposed to be related to the shape of the intrusion and the stages of its growth: the style and magnitude of surface deformations relate to the length, thickness, and tip geometry of the magma sheet [Mastin and Pollard 1988;Rubin and Pollard 1988].
Magma sheet growth models and emplacement mechanisms are critical to contemporary studies of active intrusions, since they provide key insights into the signatures of seismicity during active magma sheet growth [Pallister et al. 2010], the role of magma sheets in volcano flank instability [McGuire 1996], and surface topographic effects of intrusion in planetary rift systems [Rubin 1992;Wilson and Head 2002;Wyrick and Smart 2009]. A number of models have been proposed for sheet intrusion propagation [Rivalta et al. 2015;Galland et al. 2018;Stephens et al. 2021], including: (1) the splitting model [Lister 1990;Rubin 1993], related to tensile stress concentration at the tip according to linear elastic fracture mechanics (LEFM) (e.g. Figure 1); (2) the LEFM-Barenblatt cohesive zone model [Rubin 1993]; (3) host rock shear-failure models, via brittle and/or ductile tip zone faulting [Pollard 1973], or viscous indentation [Merle and Donnadieu 2000;Spacapan et al. 2016] (e.g. Figure 2E-H); and (4) host rock fluidisation [Schofield et al. 2010] (e.g. Figure 2C). These models are founded largely in field observations of natural intrusions; however, it is not possible to directly observe the growth of sheet intrusions.
Several recent laboratory and numerical models, coupled to field observations of intrusions in relatively weak host rocks [e.g. Haug et al. 2018;Poppe et al. 2019;Souche et al. 2019] have highlighted the role of host rock material properties in controlling the intrusion emplacement mechanism. Sill emplacement in mudstones can be accommodated through host rock shear failure [Haug et al. 2017;Galland et al. 2019]. Shear failure typically associates with blunted fingerlike intrusion geometries (e.g. Figure 2), rather than the classical bladed shapes (e.g. Figure 1) that have formed the basis for propagating models using a LEFM framework for the past 50 years [e.g. Pollard 1973;Mastin and Pollard 1988;Gudmundsson 2002;Magee et al. 2019]. Tip blunting in these cases is probably a consequence of the low shear strength of the host rock, with elevated tip stress resulting in ductile deformation as distributed intergranular fracturing [Magee et al. 2019]. Alternatively, blunted tips may arise from: (1) magma emplacement into hot and/or deep, ductile host rocks [e.g. Bons et al. 2004]; (2) emplacement of viscous magma within the whole intrusion or only within the tip region [e.g. Figure 2B; Wilson et al. 2016]; and (3) abutment of the tip against a pre-existing structure such as bedding interfaces, fractures, or faults [Figure 2I;Baer and Beyth 1990].
The purpose of the present two-part contribution (see Stephens et al. [2021]) is: (1) to demonstrate the effects of sill tip geometry on propagation potential, using Finite Element (FE) models, for a range of idealised non-bladed tip shapes (e.g. Figure 2A); (2) to critically assess the inference that intrusion segments represent the early stages of the intrusion process; and (3) to propose an evolutionary model for intrusion emplacement [Stephens et al. 2021].
2 Background: Relating crack shape to induced stress in the host rock Emplacement of sheet intrusions has been widely considered to occur in stages, with (1) initial fracture prop-agation, followed by (2) magma flow through the fracture, and (3) fracture inflation by magma pressure inside the fracture [Baer 1995]. Initial fractures are segmented, and may exhibit minor offsets in the position of fracture planes (e.g. Figure 1B, C,F and Figure 2I); their propagation leads to interaction and linkage, which preserves these offsets, as steps in the sheet wall, or remnant and abandoned tips [e.g. Figure 1B,H; Pollard et al. 1975;Francis 1982;Rickwood 1990]. Preserved tips and isolated segments are therefore treated as being representative of the growth stages of sheet development [Magee et al. 2019]. This classical model for sheet emplacement has been generally considered as an elastic-brittle process, in which magma pressure in the intrusion performs work on the intrusion walls in a hydraulic fracture [Rubin 1995;Gill and Walker 2020]. Deformation ahead of the tip occurs as elastic bending of the host rock to accommodate the gradual gradient of the intrusion thickness ( Figure 1A, F-H). The intrusion fracture will grow when the tip stress intensity factor K equals the material fracture toughness K c . For a nongrowing crack subjected to a uniform internal pressure ∆P , where a is the major semi-axis of the crack [Rubin 1995]. If the crack is growing, Equation 1 must be expanded to account for a spatially variable excess pressure [Lister and Kerr 1991]. The K " K c fracture criterion is commonly applied in the sheet intrusion literature because dykes and sills have generally been treated as a Linear Elastic Fracture Mechanics (LEFM) problem. The only geometric term considered is the crack length, since the assumption is that the crack is sharp; hence K does not directly apply to non-bladed (i.e. blunt) cracks. Crack shape has long been considered in terms of stress concentration with respect to tensile loading of an infinite elastic medium containing an elliptical hole [Kirsch 1898;Inglis 1913;Eshelby 1957], where the local maximum tensile stress at the tip due to an applied remote stress is where a and b are respectively the major and minor semi-axes of the hole (e.g. Figure 3B), and σ 8 is the applied remote tensile stress at distance from the hole, applied parallel to b [Inglis 1913]. This reduces to σ max " 3σ 8 in the case of a circular hole (a " b), whereas for cracks, b Ñ 0, hence σ max should theoretically increase towards infinity at the crack tip. Since materials cannot withstand infinite stress, the stress concentration is capped by non-elastic deformation in the tip region, which is commonly referred to as the Note that, here we reckon tensile stress positive and compressive stress negative, such that the maximum tensile stress is σ 1 and the least tensile stress is σ 3 [Pollard and Fletcher 2005]. process zone. We can define tip curvature for an elliptical crack (e.g. Figure 2B) as κ " a{b 2 (where κ " 1{ρ, the radius of curvature; see e.g. Figure 3B). For a sharp crack tip, the stress intensity factor due to an external load is K " σ 8 ? πa, and for an internally applied pressure it is K " P ? πa, which suggests that the internal load can simply be substituted for the external load in Equation 2. Hence, we propose that the tip zone stress concentration related to internally-applied loads for elliptical cracks can be written as where P is the pressure internal to the crack, which in the case of the sills considered here is the excess pressure P e relative to the hydrostatic pressure (see Section 3 below). The aim here is to test this hypothesis for generalised tip shapes, as well as elliptical cracks. Note that Equation 3 is not expected to perform well for circular features because-as with Equation 2-it predicts σ max " 3P for a circle, whereas in practice σ max " 2P in this case. This is because there is a significant horizontal load component acting on the circular surface, whereas the analogous external load of Equation 2 is purely vertical. For high-aspect ratio horizontal features such as sills, the internal pressure will act almost Segment tip geometry of sheet intrusions, I  entirely vertically, and in the sharp crack case it is entirely vertical, which is more applicable to the scenario envisaged for the application of Equation 3.
Several studies of intrusion tip segment shapes, particularly on the lateral periphery of igneous sills, have shown that segments can exhibit relatively blunted forms (e.g. Figure 2), with superelliptical [i.e. Lamé 1818] cross sections in the tip region (compare Figure 1 and Figure 2). The superellipse curve form can be defined by´x a¯n`´y b¯n " 1 (4) where n " 2 will give a circle if a " b, or an ellipse if a ą b (Figure 2A and 3B). Superellipses have geometries ranging from a cross (n Ñ 0) to a square/rectangle (n Ñ 8, e.g. Figure 3B). Superellipse geometries where n ą 2 do not exhibit a single tip maximum curvature κ max (where κ max " 1{ρ min ; see Figure 3B), nor a sim-ple relationship between the axis length-scales and tip curvature. Superellipses therefore pose issues for calculation of stress intensity K, and there are few analytical solutions for the stress around superelliptical tips [e.g. Obert and Duvall 1967]. Solutions for stress concentration around superelliptical cracks subjected to an external load [e.g. Greenspan 1944] show that tip form affects the magnitude and distribution of σ max around the tip, with increasingly superelliptical tips (i.e. higher n) moving σ max increasingly out of plane of the crack. Here we apply this principle to consider the effect of superelliptical tips in controlling the magnitude and distribution of σ max for fluid pressure internal to a crack, and therefore the potential pathways of intrusion growth.

Methods
To consider the effects of superelliptical tip geometries on the stress induced in the surrounding host rock, we present 2-D linear-elastic models for a range of sill tip shapes, including elliptical and superelliptical (Figure 2A and 3B) curves. The models can be applied to dykes by rotating the reference frame by 90°, but here we will consider horizontal intrusions. The models are constructed using the finite element (FE) software COMSOL Multiphysics (v5.5). Intrusion tip geometry variation is a consequence of some combination of material and environmental conditions and the resulting dominant deformation processes, and can be modified by non-elastic processes [Rubin 1995] such as: (1) mechanical and/or thermal erosion [Delaney and Pollard 1981;Bruce and Huppert 1990]; (2) distributed faulting [Rubin 1992]; and (3) syn-and post-intrusion viscous deformations in the host rock [Rubin 1993;Schmiedel et al. 2017]. Here, we prescribe a static sill geometry within a linear-elastic medium, and apply a constant magma pressure P m " 1 MPa internal to the sill. To isolate the effect of sill tip geometry, we do not apply any other stress state to the material volume (i.e. the host rock lithostatic pressure P L " 0 MPa); the stress perturbation in the host rock is solely the result of excess pressure P e (" P m in this case) within the crack.
The 1ˆ1 m 2-D model space comprises an isotropic and homogenous elastic medium. The model space is sufficiently large compared to the sill size to avoid stress interaction with the boundary: in all models, the sill major semi-axis a is 15 % of the x dimension of the model space, and stress is zero at the model extremities. Figure 3 shows the model symmetry axes. Material and model parameters are chosen to avoid significant strain within the material for the applied magnitude of P m . The host rock has a Poisson's ratio ν " 0.25, and Young's modulus E " 40 GPa, though since these are linear elastic models, these properties only affect the scale of distribution of stress amplification: larger ν and/or smaller E would amplify the effect of an applied stress. As noted above, the pressure applied within the sill represents an excess pressure P e " P m´PL . The magma pressure serves as a reference pressure to which all other stress is relative, and changes in the applied stress would result in a linear change in resolved stresses within the host rock. Stresses in FE simulation are dependent on the model size and the mesh size. As a consequence of the finite model size, and the mesh size, the infinite stress predicted for a bladed crack is reduced to a finite value. As noted above, this limit would be imposed anyway by non-elastic deformation within the process zone ahead of the tip.
Modelled sill dimensions are varied to consider (a) the effect of tip shape on the magnitude and distribution of stress in the host rock, and (b) the relative effect of tip shapes for sills of different lengths as a function of sill aspect ratio R " a{b. Tip shape is varied within the models by changing the n exponent in Equation 4, in the range 2 to 100, and we consider each sill tip shape for a range of aspect ratios R " t2, 5, 10, 50, 100, 1000u. For an ellipse-a superellipse where n " 2, which separates cross-like superellipses (n ă 2) from squared-appearance superellipses (n ą 2)-the curvature is given by: The minimum radius of curvature for an ellipse (Figure 2A) is given by ρ min " b 2 {a, and κ max " 1{ρ min . The relationship for a superellipse where n ą 2 is more complex, as κ max splits to a position above and below the centre-line of the crack (Figure 2A). In Cartesian form, the curvature of a superellipse is The parametric form for a superellipse is x " a cos m θ and y " b sin m θ, where m " 2{n, and the curvature κ is given by where x 1 " dx{dθ etc. Hence, we can write κpθq " " ηp2´mq am where η " b a " 1{R and f pθ, ηq " sin m´4 θ cos m´4 θ pcos 2m´4 θ`η 2 sin 2m´4 θq This is minimised when η 2 " pm´4q`2pm´1q tan 2 θ tan 2m´4 r2pm´1q`pm´4q tan 2 θs .
This is a non-linear function in tan θ so there is no analytical solution. We solve by calculating η 2 pθq and then given the value of η 2 interpolate the function to find the critical value of θ (see Supplementary files). We validate calculated κ max by plotting the superellipse curve and taking the curvature as a step function, which we perform in MATLAB. This method, and resulting value for κ max , is limited as it is dependent on the step size, and for superellipses with large n, the step size must become vanishingly small. Comparison of results for κ max based on Equation 6 in the range of n and R where stepsize is not a limitation, and results from Equations 7-10 are in agreement to three decimal places, which is more than sufficient for our purposes. Figure 3C shows the results of calculated maximum curvature as a function of n for a range of aspect ratios (1-1000).

Results and Interpretation
Propagation of a sill through our modelled isotropic host medium would be most likely where the tensile stress and/or the shear stress concentration is at its maximum (i.e. σ max and/or τ max ), and is therefore most likely to exceed the host rock strength. The mode of failure (extension, extensional-shear, or compressionalshear) is controlled by the conditions of stress at the tip, though it should be noted that our models describe the approach to failure, and not the failure itself. Increasing the aspect ratio of modelled elliptical sills (n " 2) from R " 1 to R " 1000 causes an increase in the stress concentration at the position of maximum curvature κ max , at θ " 0˝, as is expected from Equation 1 and 3 (Figure 4 and 5): from Figure 5E we find that ∆σ max " σ max´Pe " A ? κ max , where A " 2P e ? a in accordance with Equation 3. Inspection of Figure 5D suggests that deviation from the power law trend at the upper extremity of the regression is a mesh size-related underestimation of stress for the very sharp tip region of the R " 1000 ellipse.
Changing the tip region shape to superelliptical geometries, as a function of the n exponent, results in two key changes to stress concentrations in the host rock: (1) changes to the maximum stress generated; and (2) changes to the position of that maximum stress with respect to the horizontal centre line of the sill (Figure 4 and 5). For superellipses with aspect ratio R " 2, increasing n from 2 (elliptical) to 3 (a slightly-squared ellipse; e.g. Figure 3B) results in a decrease in resolved tensile stress in the host rock ( Figure 5A-B). Increasing n from 3 to 100 (a rectangle at the scale of observation; e.g. Figure 3B) results in a significant increase in the maximum stress generated in the host rock, such that where n ě 5, σ max is larger than that recorded for an elliptical intrusion of the same aspect ratio (Figure 5A-B); σ max for an ellipse is~19 % of the magnitude recorded for a superellipse where n " 100 ( Figure 5B).
Increasing superellipse aspect ratio results in a relatively diminished effect of increasing n, such that where R ě 50, an ellipse (n " 2) records a larger stress concentration than superellipses in the study range and at the same aspect ratio ( Figure 5A-B). The relative effect of sill aspect ratio on stress concentration as a function of n is highlighted by normalising stress concentrations relative to those recorded for R " 2 (Figure 5D). We attribute this to the change in maximum curvature, which is a function of the sill aspect ratio and the tip shape ( Figure 5E). Increasing the aspect ratio of an ellipse results in an increase in the maximum curvature by about six orders of magnitude from R " 2 to R " 1000 ( Figure 3C and 5E). Superellipses where R " 2 and n Ñ 100 already have a large κ max at the apparent-squared corners of the tip. At any aspect ratio, a superellipse where n " 100 appears indistinguishable from a rectangle at the scale of observation: Increasing the aspect ratio for large n segments (from  Figure . The colour scale data range is uniform for all plots of ∆σ 1 shown in [A], and uniform for all plots of ∆τ xy shown in [B]. Note that maximum ∆σ 1 , i.e. σ max , is along the sill centre axis for n " 2, and splits to create two zones of maxima that are increasingly off-axis with increasing values for n. ∆τ xy shows a similar splitting with increasing n; stress maxima (σ max and τ max ) in the modelled host rock correspond to the position of maximum curvature κ max .
[C] Potential shear (slip) planes [c.f. Pollard ] plot onto ∆τ xy for n " 2 and n " 100. Light-grey shear-couple arrows show potential slip sense for growth by vertical inflation of the sill, whereas darkgrey arrows show slip sense for horizontal growth by viscous indentation, in both cases accommodated by faulting in the host rock.

Figure :
The relationship between modelled stress change and tip zone curvature.
[A] Maximum modelled tensile stress change for sills of aspect ratio ranging from to , with curve forms varied as a function of n.
[B] σ max normalised to the value of the maximum σ max for each given aspect ratio.
[C] The y-coordinate position of σ max for each given aspect ratio, normalised to the minor semi-axis b. A value of . corresponds to the sill centre axis, and a value of . is the full thickness of the intrusion. [D] σ max for aspect ratios of five and above, normalised to the σ max values for R " 2.
[E] σ max versus calculated maximum curvature κ max for the modelled sill dimensions. Note that increases in σ max show a consistent power law relationship with increasing κ max , ranging from κ 0.42 max to κ 0.55 max (i.e. " ? κ max ).
R " 2 to R " 1000) increases the maximum curvature only by about three orders of magnitude ( Figure 3C). As with ellipses, from Figure 5E we find that for superellipses 2 ă n ď 100, ∆σ max " A ? Figure 5E shows that when fitting a set of results over a small range of curvatures, the exponent for each fit deviates from a square-root scaling by´0.08 and`0.05 (i.e. from κ 0.42 max to κ 0.55 max ), but overall the fit to Equation 3, across the entire set of results, is excellent. Figure 4A and 4B show that for superellipses where n ą 2, σ max and τ max move increasingly out-of-plane as n increases. Figure 5C shows the position of σ max in the y coordinate normalised to the half-thickness (b) of the intrusion; i.e. where y{b " 0, the maximum stress is along the sill centre axis. For elliptical sills (n " 2) we expect coplanar growth initiated at the site of σ max , producing a longer, flat sill (y{b " 0). Increasing n moves σ max progressively away from the sill centre axis (0 ă y{b ă 1). The non-coplanar distribution of maximum stress for superellipses may result in a number of growth scenarios, which we discuss below.
Our linear elastic models for stress distributions around static sill segments show that tip shape-here simplified to elliptical to superelliptical curve formsis important in controlling the magnitude of stress concentration at the tip, and the position of that maximum stress. Blunted tip shapes concentrate stress out of plane of the sill, with the position of maximum stress corresponding to the position of the maximum tip curvature (Figure 4). For small aspect ratio segments, blunted tips may be more efficient at concentrating stress than ellipses (Figure 5B), as a consequence of the larger curvature associated locally with the squared superellipse corners. Increasing the aspect ratio of an ellipse or superellipse characterised by a small nexponent results in a significant increase in the maxi-mum stress, as a consequence of the increased tip zone curvature. In contrast, increasing the aspect ratio of superellipses characterised by a large n-exponent, results in a minor change to the tip curvature relative to short segments.

Tip shape and stress concentration in the host rock
Segment tip-zone curvature imposes a sill-lengthdependent control on the maximum stress generated in the host rock ( Figure 5), and controls the position of that maximum stress change within the host rock relative to the centre line of the sill major axis (Figure 4 and 5). These findings have implications for (1)  [2018] found that short sills (in their case, a " 200 m) of constant thickness (b " 25 m), with a semi-circular tip region (ρ min " 25 m and κ max " 0.04), required an excess pressure of >900 MPa to cause shear failure of the host rock; longer sill segments (a>4 km) required much less magmatic overpressure (i.e. < 60 MPa). This decrease in the required excess pressure is to be expected given the relationship between stress intensity and crack length (e.g. Equation 1) and σ max (e.g. Equation 3), since the key length scale controlling tip stress is the fracture length (i.e. half-length a in those equations). Their sills were modelled at~2 km depth, where a natural host rock in situ tensile strength should be expected to be on the order of 3-5 MPa, and the shear stress at failure may be <50 MPa. As noted by Haug et al. [2017], it is unreasonable to expect such large magma overpressures within a static sill, at such a depth on Earth; estimates for the likely range of possible magmatic overpressure in sills is 1-20 MPa [Rubin 1995]. Although not directly comparable, the models in Haug et al. [2017] and Haug et al. [2018] are similar to the relatively blunted (superelliptical) tips modelled here since the tip curvature is less than it would be for an elliptical sill of the same dimensions. Using the example above (a " 200 m and b " 25 m), the equivalent elliptical crack would have a maximum tip curvature κ max " 0.32 as opposed to their κ max " 0.04. The maximum stress is proportional to the square root of the maximum curvature, hence an elliptical sill would require significantly lower magma overpressure to grow.
The chosen crack dimensions and geometry of the crack tip imposes a control on the magnitude and distribution of stress in the host rock, and this must be taken into account when interpreting model results.
Our results indicate that blunted tip shapes (n ą 2), which are tip-forms that require non-elastic processes within the host rock, will produce stress maxima out of plane from a sill. There are four main possibilities for continued sill growth: (1) intrusion by tensile or mixed mode failure in the host rock, at the site of κ max and σ max ( Figure 4A); (2) horizontal growth by viscous indentation of the host rock, activating low-angle shear planes (e.g. as shown in Figure 4C for n " 100); (3) roof uplift and vertical inflation of the sheet, activating steep shear planes (e.g. as shown in Figure 4C for n " 100; c.f. Figure 2D for natural sills and Figure 2I for natural dykes); or (4) arrest of the blunted sill tip, leading to growth at some other position on the intrusion, or formation of a new sill plane or planes (e.g. Figure 2D). This final option could result in a large number of short intrusions that may be distributed at different levels within the host rock, and may therefore account for: (a) observed inconsistent stepping of sill segments [e.g. Schofield et al. 2012;Magee et al. 2019]; (b) a large number of short sills particularly in host rocks that regularly show blunted tips (e.g. mudstones: [Mark et al. 2018;Galland et al. 2019]); and (c) long, finger-like segments with blunted lateral tip zones [Pollard et al. 1975;Wilson et al. 2016].
For finger-like intrusions, the terminal (growing) end of the intrusion may not necessarily match the geometry of the lateral tips. Previous constant-pressure static sill models [Haug et al. 2017] indicated that sills emplaced by shear-failure would not grow as flat coplanar sheets; paradoxically, each increment of sill growth in their models stepped out of plane, and it is therefore difficult to envisage how flat intrusive sheets might grow for long distances via a macro-scale sheardominated failure mechanism. Large tensile and shear stresses at the blunted tip corners (Figure 4) may facilitate opening or mixed mode host rock failure, in which case-should the magma pressure exceed the normal stress on the newly-formed plane-the intrusion would grow as an inclined sheet; at least initially.
Alternatively, viscous indentation and/or fluidisation of the host rock ahead of a blunted tip, may facilitate coplanar growth, and field observations [e.g. Figure 2C, 2G;Spacapan et al. 2016] and numerical simulation [e.g. Souche et al. 2019] confirm this is a viable sill-lengthening emplacement process, at least for short distances. Growth may be promoted by slip planes emanating from the large curvature zones at the corners of blunted, superelliptical tips ( Figure 4C), however, it is not an efficient growth mechanism. Unlike elastic splitting, newly-formed fault surfaces would remain in contact, and slip would be subject to residual friction on the interface. Again, if magma enters the slip planes to part the surfaces, growth may become more efficient. Most natural examples of sills do not preserve magmatic injection into such slip planes; they are compressional-shear faults, and again it is difficult to envisage this emplacement mechanism operating as a long-lived process for all but very weak rocks (e.g. mudstones: Souche et al. [2019]).
We suggest that evolving conditions at the tip, including changes to the host rock and magma properties, would lead to an evolution of the dominant emplacement mechanism as the intrusion grows and arrests [Stephens et al. 2021]. This evolution may be important to the preserved geometry of sheet intrusions, and the features of the tip region, with the potential for late-stage processes and features to overprint those from earlier stages.

The field record of igneous sheet intrusions
Intrusion segments are generally treated as representing the staged growth of through-going sheets, with isolated segments inferred to represent the earliest stages of the process [Delaney and Pollard 1981;Baer and Reches 1987;Rickwood 1990;Schofield et al. 2012], before linkage to create the sheet-like form of sills, dykes, or inclined sheets. Field-based observations of intrusive segment tips tend to be the lateral edges of the segments as viewed in the axis of bulk magma transport [e.g. Galland et al. 2019], rather than the terminal end tip region at the distal extremities of the segment [e.g. Poppe et al. 2020]. In either case, observations of segment tips are of a system that did not grow further than the observed state, and therefore represent the arresting form of segments, rather than the growing form. This is important for interpretation, since the features associated with the preserved segments are biased towards processes that led to an end of growth for that particular observed segment; field-observationbased interpretations of growth processes are subject to preservation bias. We explore this point in more detail in part II, using field-based observations [Stephens et al. 2021].
To consider the growth stages of intrusions, it is critical to identify and separate early-stage processes from features associated with late stage modification, in the conduit [Wadsworth et al. 2015;Walker et al. 2017], and in the host rock [Delaney and Pollard 1981;Bruce and Huppert 1990;Rubin 1993;Bons et al. 2004;Wilson et al. 2016]. As magma solidifies it undergoes physical property alteration, particularly in terms of an increase in viscosity. Reduction in temperature, increases in crystal proportion, and/or exsolution of volatiles may each increase viscosity [Lejeune and Richet 1995;Hess and Dingwell 1996;Ishibashi and Sato 2007], increasing viscous dissipation and potentially slowing magma flow to the fracture tip. At the same time, intrusion growth may lead to changes in host rock cohesive strength. For instance, a reduction in host rock cohesion may arise from pore fluid boiling, causing dis-tributed rock fracture [Schofield et al. 2010]. Fractured rock cannot withstand elevated tensile stresses, and may instead deform by distributed slip and/or granular flow, leading to a relative blunting of the tip zone. These processes may be more prevalent in areas of the host rock that have sufficient contact time with the magma so as to allow conductive heat transfer, such as at the lateral edges of a sill, and in rocks with sufficiently low permeability as to allow pore fluid pressure to increase (e.g. mudstones). Lateral-edge tip zones, rather than frontal tip zones, are probably more commonly found in natural dyke and sill networks, and may have experienced long-lived heat transfer as magma flowed along the conduit to a frontal propagating tip region. In each of these cases, the final stages of intrusion could contribute to a change in the material properties relative to the growth stages of intrusion; those early stages of growth would be overprinted by the final stage of emplacement, leading to a bias in preserved emplacement mechanisms.

Conclusions
Crack tip geometry imposes a critical control on the magnitude and distribution of stress concentration ahead of a pressurised crack; a principle that we apply here to igneous sheet intrusions. Elliptical tip zones concentrate stress in-plane with the sheet and therefore growth should continue as a flat plane. Superelliptical ("blunt") tip zones concentrate stress out of plane, which may lead to (a) inclined growth with respect to the sheet plane, (b) coplanar growth via shear in the host rock, or (c) thickening of an arrested sheet. These findings are important to the results of numerical models for sheet intrusions, in which the tip zone geometry is prescribed, and to the potential for preservation of tip features indicative of sheet growth in the rock record.

Acknowledgements
This study was conducted without allocated funding. Catherine Greenfield is a Daphne Jackson Trust Research Fellow funded by NERC and the University of Leicester. Simon Gill is supported by a Royal Society Apex Award. Sam Poppe was supported by a frs-FNRS postdoctoral fellowship while at UlB and currently by a ULAM scholarship (NAWA Poland) and NCN Poland grant 2020/37/K/ST10/02447 at SRC PAS. The authors thank Tim Davis and Steffi Burchardt for discussions that benefitted the formulation of this paper, and thank the three anonymous reviewers and Editor-Jamie Farquharson-for review commentary that helped to improve the clarity of this paper.

Author contributions
RW, TL, and DH formulated the research concept and contributed field data and examples. CG, SG, and RW developed, tested, and processed numerical models and results. SP contributed to discussion topics and field examples. All authors have contributed to the theoretical background, interpretation, and discussion.