Mechanical models to estimate the paleostress state from igneous intrusions

Dikes and sills represent an important component of the deformation history in volcanic systems, but unlike dikes, sills are typically omitted from traditional paleostress analyses in tectonic studies. The emplacement of sheet intrusions is commonly associated with mode I 10 fracturing in a low deviatoric stress state, where dilation is perpendicular to the fracture plane. Many natural examples of sills and dikes, however, are observed to accommodate minor shear offsets, in addition to a component of dilation. Here we present mechanical models for sills in the San Rafael Subvolcanic Field, Utah, which use field-derived measurements of intrusion attitude and opening angles to constrain the tectonic stress axes during emplacement, and the 15 relative magma pressure for that stress state. The sills display bimodal dips to the NE and SW and consistent vertical opening directions, despite variable sill dips. Based on sill attitude and opening angles, we find that the sills were emplaced during a phase of NE-SW horizontal shortening. Calculated principal stress axes are consistent (within ~4  ) with paleostress results for penecontemporaneous thrust faults in the area. The models presented here can be applied 20 to any set of dilational structures, including dikes, sills, or hydrous veins, and represent a robust method for characterising the paleostress state in areas where other brittle deformation structures (e.g. faults), are not present.


Introduction
Sills and dikes are traditionally treated as extension fractures with a dilation vector normal to the fracture wall, i.e. they are extension fractures (Mode I, e.g.Anderson, 1951).This assumption has important implications for the use of sheet intrusions in constraining tectonic stress states, because the minimum compressive stress (σ 3 ) is perpendicular to the extension fracture walls; local deflections of the intrusion attitude are commonly inferred to represent local rotations of the stress axes.This is most commonly attributed to mechanical layering, and the presence of pre-existing structures (e.g.Rubin, 1995;Gudmundsson, 2002;2011a;Magee et al., 2016).This model implies that intrusions can locally propagate out of the regional σ 1 -σ 2 plane, via Mode I failure of intact rock, or through Mode I dilation of pre-existing structures, producing intrusions that display variable dilation vectors along a single intrusion.Notably, many field examples of sills and dikes exhibit near-parallel dilation vectors, regardless of the intrusion attitude (e.g.Hoek, 1991;Walker, 1993;Airoldi et al., 2011;Martinez-Poza et al., 2014;Muirhead et al., 2014;Eide et al., 2016;Walker et al., 2017).Intrusions that demonstrate shear-offset of markers across their margins indicate that during emplacement the dilation vector was inclined from plane-normal (Muirhead et al., 2014;Stephens et al., 2017;Walker et al., 2017); this obliquity of opening can be characterised by the opening angle (Fig. 1).Delaney et al. (1986), Baer et al. (1994) and Jolly and Sanderson (1997), who applied mechanical methods to estimate paleostress states using sheet intrusion attitudes.These mechanical models have been adopted for statistical constraints on paleostress axes, and the paleostress state, from vein or dike data (e.g.Sato et al., 2013;Yamaji, 2016).Although several methods exist to determine paleostress axes, and a state of paleostress via fault, fracture, or dike data, subhorizontal sheet intrusions (sills) are typically omitted from such analyses.This is because sill emplacement is commonly associated with a local rotation of the principal stress axes, as a result of bedding or other layer heterogeneities, during regional extension or a low-deviatoric stress state (i.e.where σ 1 ≈ σ 2 ≈ σ 3 ; e.g.Gudmundsson, 2011a;Magee et al., 2016).Several studies have demonstrated that sills may also form during tectonic horizontal shortening, where the minimum principal stress is vertical (Walker, 2016;Walker et al., 2017;Stephens et al., 2017); hence sills can be used as paleostress indicators, provided they are accompanied by detailed kinematic characterisation.

Dilated structures have been studied in detail by
Here, we present mechanical models, based on those of Jolly and Sanderson (1997), to determine paleostress state using the attitude of dilated fracture sets; we verify this using the measured opening angle of the intrusions.The method is applicable to any dilated fracture; here we focus on the stress state associated with sill emplacement in the San Rafael Subvolcanic Field, Utah, and compare these results to fault data in the same area, to demonstrate the particular importance of subhorizontal igneous intrusions as records of paleostress.

Dilation of pre-existing fractures
Dilation of a cohesionless pre-existing planar structure occurs when the fluid pressure (Pf) exceeds the normal stress on the plane (σ n ; Fig. 1) (Delaney et al., 1986;Jolly and Sanderson, 1997).Normal stress is related to the plane attitude (θ), and the maximum (SH) and minimum (Sh) principal stresses acting on the plane (Jolly and Sanderson, 1997): with the range of possible dilatant fracture attitudes controlled by the stress ratio (), and the driving pressure ratio (R') (Baer et al., 1994): The stress ratio, ϕ, is a non-directional value that describes the relative magnitudes of the principal stresses: σ 3 ≤ σ 2 ≤ σ 1 (here, compressive stresses are positive).The driving pressure ratio R' describes the relative magnitudes of the fluid pressure and the remote stress state: When R' < 0 (i.e., Pf ≤ σ 3 ) there is no dilation, and when R'  1 (i.e., Pf ≥ σ 1 ) all fracture attitudes could dilate.This assumes that an intrusion is emplaced during a deviatoric stress state (i.e. where σ 1  σ 2  σ 3 ); if Pf > σ 3 during a non-deviatoric stress state (i.e.where σ 1 = σ 2 = σ 3 ), all pre-existing fracture attitudes could dilate.
Pre-existing fractures, with poles parallel to the σ 3 axis, will show plane-normal dilation (i.e.dilation parallel to the normal stress vector).Fractures inclined to the plane of σ 3 resolve a shear stress on their surface (Delaney et al., 1986).During dilation, shear stress is reduced to zero through a plane-oblique dilation vector (extensional-shear).The angle between the dilation vector and normal stress is defined as the opening angle (µ; Fig. 1a-c) and represents the ratio of shear to dilation (Delaney et al., 1986).Extensional-shear therefore acts to reduce the amount of dilation.For an intrusion comprising inclined and subhorizontal sections, the inclined sections will therefore be thinner than subhorizontal sections oriented perpendicular to σ 3 (e.g.Fig. 1c, d) for a given fluid pressure, providing the stress state is deviatoric and 0 < R' < 1 (England, 1988).The dilation vector can be measured using traditional compass techniques in the field based on offset piercing points or unique contact geometries, such as recognisable corners in the intrusion walls (e.g. Figs 1c,2a,b).True offsets result from extensional-shear opening (Fig. 2c, d), however apparent shear offset can be produced by marker units that are oriented at an oblique angle to the intrusion plane (Fig. 2e), or via dilation of pre-existing faults, which produce a larger than expected opening angle (Fig. 2f).To determine whether offset and the opening angle is true, or apparent, multiple manual measurements of opening angle and intrusion thickness must be made along strike of an intrusion (Figs 1c,2g).The opening angle (µ) of a dilated fracture can be calculated as the inverse cosine of the true thickness (t, normal to the plane) and the vertical thickness (tv, parallel to the dilation vector; Fig. 1c): Alternatively, µ can be measured as a rake of the obtuse angle between the intrusion contact and the dilation vector, minus 90°.Importantly, this measurement only accounts for the opening angle in two dimensions, movement in or out of the face will not be recorded; therefore dip-parallel sections should ideally be used.
The opening angle is related to the shear stress (τ), normal stress (σ n ), and fluid pressure (Pf) acting on that plane at the time of intrusion (Delaney et al., 1986;Jolly and Sanderson, 1997): ). (5) Equation ( 5) shows that if the fluid overpressure (Pf -σ n ) is equal to the shear stress, the opening angle is 45º, and the shear to dilation ratio is unity.If the overpressure is greater than the shear stress, the opening angle is less than 45º, and the fracture will show a greater component of dilation to shear.When µ is negative, the fracture will remain closed as the fluid pressure did not exceed the normal stress.An intrusive segment, however, may inflate against a closed fracture (where Pf < σ n ), causing a local contractional shear, and a blunt intrusion tip (e.g.Fig. 2g).

Mechanical models for fracture dilation and opening angle
Equation ( 5) can be represented in three-dimensional (3D) space, for any given stress state, stress ratio and fluid pressure: Figure 3 shows the opening angle (µ) of all possible fracture planes in 3D space, plotted as colour-contoured pole to plane values of µ, on equal area lower hemisphere stereographic projections and on 3D Mohr Diagrams.Fluid pressure was calculated at five equal intervals relative to the ambient stress state, which are expressed as R' values (Fig. 3a). Figure 3 shows plots of opening angles for three tectonic regimes where ϕ is 0.5: (1) a thrust-fault regime (σ v = σ 3 ) (Fig. 3b); (2) a strike-slip regime (σ v = σ 2 ) (Fig. 3c); and (3) a normal-fault regime (σ v = σ 1 ) (Fig. 3d).The principal stress attitudes are constant in all models, with the azimuth for the maximum horizontal stress (σ H ) trending E-W.All planes are modelled as cohesionless surfaces.
The models complement the results of Jolly and Sanderson (1997), and demonstrate that within a given stress state, increasing the fluid pressure increases the range of pre-existing fracture attitudes that can dilate, provided the dilated fracture is linked to the magmatic source.
The models also show that for any given fracture, increasing Pf decreases µ.Where σ 3 < Pf ≤ σ 2 (Fig. 3b-d, i-ii), the dilation zone, delineating the poles to fractures that are predicted to dilate, forms an ellipse about the σ 3 axis and is elongate in the direction of the intermediate stress (σ 2 ); only fractures with poles parallel to the σ 3 axis show Mode I opening.If σ 2 < Pf < σ 1 (Fig. 3biii-diii) the dilation zone forms a girdle parallel to the axis of σ 2 , with two defined zones of near Mode I opening (0-10º) surrounding the σ 2 and σ 3 poles.Fractures of all attitudes can dilate if Pf = σ 1 (Fig. 3biv-div).The models suggest that Mode I opening of pre-existing fractures should only be common if the fluid pressure exceeds σ 1 (e.g., .For dilation of misoriented cohesionless fractures, with planes oriented at a high angle to σ 3 , a higher fluid pressure is required to overcome the normal stress acting on their surface (e.g.Gaffney et al., 2007).Importantly, if misoriented planes do not become linked to the fluid network, high fluid overpressures will be relaxed by dilation of preferentially oriented pre-existing structures, or via failure of intact rock (e.g.Sibson, 2012).Dilational fracture data should therefore plot stereographically with the greatest density of poles to planes about the σ 3 axis (Jolly and Sanderson, 1997;Yamaji et al., 2010).Propagating intrusions may cut across or terminate against misoriented structures, which plot in the zone of no dilation (white areas on the Mohr circles and stereonets; Fig. 3), as the normal stress acting on these planes exceeds the fluid pressure, preventing dilation (e.g.Gaffney et al., 2007;Gudmundsson, 2011a).Terminated intrusions may inflate against closed fractures if magma supply continues; dilation will be accommodated by local shear along the closed fracture producing a blunt intrusion tip (Baer, 1991;Stephens et al., 2017).Through measurement of the opening angles of veins or intrusions, and determining the distribution of fracture attitudes that are in contact with intrusions, but not dilated, it is therefore possible to constrain the fluid pressure relative to the ambient stress state.
In addition to changes in fluid pressure, fracture dilation is also sensitive to changes in remote stress, particularly at lower Pf.For instance, where σ 3 < Pf ≤ σ 2 (Fig. 3b-d, i-ii) each stress regime is defined by a particular range of fracture attitudes, with a unique opening angle pattern.Where σ 2 < Pf ≤ σ 1 (Fig. 3b-d, iii-iv), the thrust-fault regime (Fig. 3biii-iv) and strikeslip regime (Fig. 3ciii-iv) have similar opening angle patterns, whereas the distribution of attitudes and opening angles for the normal-fault regime (Fig. 3diii-iv) is unique.Applied carefully, using detailed kinematic field analyses of sills, and their associated deformation, the method has the potential to discriminate between intrusions related to different remote stress states.Here, we have applied this method to igneous sills in the San Rafael Subvolcanic Field, Utah (SRSVF).
Dilated fractures, created during the same dilational event, will produce a unique distribution of opening angles when plotted stereographically as poles to planes (e.g.Jolly and Sanderson, 1997;Yamaji et al., 2010;Sato et al., 2013).Whether this is a tight cluster or a distributed set of poles, the pattern is usually interpreted to represent a single event, governed by one fluid pressure (e.g.Jolly and Sanderson, 1997); however, it is possible that the data represents multiple events caused by fluid pressure pulses of varying magnitude (Yamaji et al., 2010).As the opening angle represents the ratio of shear to dilation (Delaney et al., 1986), magmatic events of varying magnitude, within the same governing stress state, could produce cross-cutting intrusions with similar geometries, but different opening angles (e.g.Fig. 3).The attitude and opening angle of dilated fractures, along with cross-cutting relationships, can therefore be used to identify whether a network of dilated fractures represents one, or multiple, fluid pressure pulses and the stress state during each pulse.

Stress state model for sill emplacement: San Rafael Subvolcanic Field
The San Rafael Sub Volcanic Field (SRSVF) is located in the north-western Colorado Plateau, central Utah (Fig. 4a), and comprises numerous, sills, dikes and volcanic breccia bodies (Fig. 4b, c).The close spatial association between these different intrusive bodies has been suggested to reflect emplacement during a low deviatoric stress state (Delaney and Gartner, 1997;Richardson et al. 2015).No crustal magma chambers have been identified as the source for the intrusive complex (Delaney and Gartner, 1997).The exposed intrusive system was emplaced at ~1 km depth, within the Mid Jurassic sedimentary rocks of the San Rafael Group, which comprises four formations; the Carmel Formation (siltstone and shale), Entrada Sandstone (interbedded sandstone, siltstone and claystone), Curtis Formation (interbedded glauconitic sandstone, and siltstone), and the Summerville Formation (siltstone, sandstone and shale) (Gartner, 1986;Richardson et al., 2015).The sequence represents a paralic environment, with near-shore and shallow marine deposits (Gartner, 1986;Peterson, 1988).The sills are mainly emplaced within the Entrada Sandstone, however they also cut across formation boundaries to intrude the Carmel and Summerville Formations (e.g.Fig. 4d; Walker et al., 2017).Richardson et al., (2015) highlighted an example of a dyke-to-sill transition for the Cedar Mountain sills (Fig. 4b).Our study focuses on sills in the southern SRSVF (Fig. 4b), and we cannot confirm that relationship here.In our study area, sills and dikes display mutual cross-cutting relationships (e.g.Fig. 4e; Gartner, 1986;Walker et al., 2017), and we have observed no dykesill transitions for the major sills (e.g., Fig. 4c).Delaney and Gartner (1997) determined ages of dike emplacement between 4.6 and 3.7 Ma; the field relationships therefore constrain sill emplacement to a ~1 Myr interval.Walker et al. (2017) interpret the sills in the SRSVF as low-angle conjugate intrusions, and concluded that the sills record a phase of horizontal tectonic shortening, rather than relating to local deflections due to material layering.The sills range from <10 cm to ~30 m thick, and are discordant to bedding, with dips of 1 -25º (Fig. 4d, f).Based on sill attitude measurements over km-scale outcrops, Walker et al. (2017) showed that sill poles cluster about a near-vertical axis with two defined clusters: a NW and SE dipping set in the northern SRSVF, and a NE and SW dipping set in the southern SRSVF.This study focuses on the NE and SW dipping sills in the southern SRSVF (Fig. 4b, c).
Subvertical to vertical fractures (>70º dip) were not intruded (Fig. 7a) suggesting that fluid pressure was less than the normal stress on those planes during intrusion, and the opening angle was < 0°.
Since the range of dilated sill attitudes is known, we can determine the relative stress and fluid pressure, using the parameters ϕ and R'.We calculate these parameters stereographically, using an adapted method from Jolly and Sanderson (1997).The local sill contact data are plotted as poles to planes and coloured by their µ-value; the pole cluster showing Mode I opening (µ=0-10º) is fitted with an ellipse (Fig. 7b).For data with a clustered distribution, this ellipse geometry provides a guide for the total range of dilated fractures, where blunt tips define the limit of fracture attitudes able to dilate; beyond this is the zone of no opening (where µ<0º; Fig. 7b).The minimum compressive stress (σ 3 ) plots in the centre of the Mode I ellipse (which likely coincides with the centre of the data cluster); σ 1 and σ 2 are mutually orthogonal to this (Fig. 7c).If the data has a girdled distribution, the 0-10° ellipse should be used as an approximation for the geometry of the zone of no opening.In these cases, there will be two clusters of poles where µ=0-10° (see Fig. 3biii-diii), the larger ellipse signifies the location of σ 3 .A distributed set of dilated fractures may contain several zones of Mode I opening (µ=0-10°; see Fig. 3): σ 3 will plot in the larger of these ellipses and σ 1 in the zone of highest opening angle, or the smallest Mode I ellipse when Pf exceeds σ 1 (see Fig. 3bv-bv).
Field data for sills in the SRSVF can be fitted to an elliptical region on the stereonet with a NW-SE long axis (Fig. 7b), giving horizontal NE-SW maximum compression where: σ 1 plunges 3° towards 068°; σ 3 is vertical, plunging 87° towards 265°; and σ 2 is horizontal plunging 1° towards 158° (Fig. 7c).Using the stereonet data in Figure 7c we derive the angles θ1 and θ2 which are used to calculate the stress ratio (ϕ) and driving pressure ratio (R') by Jolly and Sanderson (1997): and where θ1 is the angle between the σ 2 axis and the perimeter of the dilational ellipse and θ2 is the angle between the σ 1 axis and the dilational ellipse (Jolly and Sanderson, 1997).When an ellipse can only be fitted to the zone of no dilation, Jolly and Sanderson (1997) provide an alternative method for calculating ϕ.The R'-value is two-dimensional (2D); it does not take into account the magnitude of σ 2 .The three-dimensional relationship between the fluid pressure and all principal stresses is illustrated by the ellipse geometry on the stereonet (Fig. 7c, d), and can also be visualised through construction of a 3D Mohr circle (Fig. 7e).
The calculated ϕ (0.77) and R' (0.68) values define the stress ratio and driving stress used to create the opening angle mechanical model (Fig. 7d).For the model, we assigned a minimum, vertical stress of 25 MPa to simulate a ~1 km emplacement depth (Gartner, 1986).
Development of extensional-shear fractures requires a low differential stress (σ D : σ 1 -σ 3 ): 4T < σ D < 5.66T; compressional shear faults require that σ D >5.66T (T is the host rock tensile strength; Sibson, 2003).Due to the bimodal (conjugate) dip distribution of the sills, their consistent near-vertical opening direction, and the mutual cross-cutting and intrusive relationship between thrust faults and sills we estimated that σ D = 6T.We estimate the tensile strength of the host rock at 1 km depth to be 3 MPa, giving σ D = 18 MPa.It should be noted, however, that due to the nature of the model, providing all parameters are scaled relative to T, the value of T does not change the resulting opening angle pattern, only the magnitudes of the principal stresses and fluid pressure.ϕ (0.77) and R' (0.68) derived in the model indicate a mild horizontal radial compression during emplacement, with fluid pressure less than σ 2 (Fig. 7d,   e).
The measured opening angles of the studied sills fit with the contouring of a single ellipse (Fig. 7d), suggesting that during the intrusion of sills in the study area, the fluid pressure and stress state remained relatively constant.However, we acknowledge that a broader study across the SRSVF may reveal spatio-temporal fluctuations in either the stress regime or magma pressure.

Sill geometry as a record of far-field stress
Sills in the SRSVF are mildly transgressive, in that although they have intruded parallel to bedding at a local scale (e.g., Fig. 5b, c, f), at the field-scale they cut across the stratigraphy at a low angle (Figs 4d and 5a-b;Walker et al., 2017).Transgressive sill geometries are typically inferred to represent either the peripheral section of a saucer-shaped sill (e.g., Malthe-Sørenssen et al., 2004), or a low-angle cone sheet emanating from a magma chamber (e.g.Gudmundsson, 2006;Martí and Geyer, 2004).In either case, the transgressive segments should have an overall centripetal dip distribution.No magma chamber has been identified in the SRSVF, and the sills do not dip towards a central source, nor do they have a saucer-shaped geometry.Sills in the SRSVF show consistent bimodal dip patterns, with mutual cross-cutting relationships, at the outcrop to field scale (i.e. for sills that are centimetres-to tens-of-metres thick).We infer, as noted by Walker et al. (2017), that the bimodal sill geometry is representative of low-angle conjugate faults, which required a deviatoric stress state (σ D > 5.66T; Sibson, 2003) during emplacement.
At the local (cm -m) scale, the SRSVF sills display a range of attitudes with consistent vertical opening directions, and attitude-dependent thickness variations.In a low-deviatoric stress state, local variations in intrusion attitude are traditionally inferred to represent rotations of the principal stress axes, which may be caused by the effects of mechanical layering (e.g.Gudmundsson, 2011a), or pre-existing faults and fractures (e.g.Fossen, 2010;Magee et al., 2013).Local stress rotations would therefore cause intrusions of any attitude to dilate in a Mode I sense, producing a polymodal distribution of opening directions across an intrusive suite.In that case, the aperture (true thickness) of the segment would be dominantly controlled by the host rock elastic properties (mainly Young's Modulus, E) (e.g.Brenner and Gudmundsson, 2004).Units with low E are less stiff than those with high E; i.e., a low E material will accrue more strain for a given stress.Accordingly, when a fluid-filled fracture, with constant fluid pressure, cross-cuts units with varying elastic properties, it will have a larger aperture in units with low E, and a smaller aperture in units with higher E (Brenner & Gudmundsson, 2004;Gudmundsson, 2011b).Local attitude variations of the SRSVF sills do not correspond solely to bedding interfaces.Although we infer that the inclined sill sections are dilated pre-existing fractures, the consistent vertical opening directions suggests that inclined fractures were subject to extensional-shear, rather than Mode I, dilation.This is supported by the attitude-dependant (rather than lithology-dependent) thickness variations; horizontal sill segments are consistently thicker than adjoining inclined segments, regardless of the hosting lithology.
Areas that host several intrusion sets, with different attitudes, are commonly interpreted to represent discrete and separate intrusive events (e.g.Delaney and Gartner, 1997;Walker, 1993).Our study, however, suggests that the observed range of sill attitudes is indicative of magmatism during a deviatoric far-field stress state (e.g., Walker, 2016;Walker et al., 2017).
The sills may, therefore, be equivalent to reverse faults in the area, accommodating horizontal shortening, and vertical thickening.

Rafael Sub-Volcanic Field
The SRSVF is host to dikes and sills; the sills cut, and are cut by, thrust (and reverse) faults that dip to the NE and SW (Figs 5f, 8a;Walker et al., 2017).Thrusts, deformation bands, and gypsum veins in the region form conjugate sets (Fig. 5d-i), which record a coaxial horizontal shortening and vertical thickening (Walker et al., 2017).Right Dihedra paleostress analysis (using the method of Delvaux and Sperner, 2003) of the thrust faults and deformation bands in the San Rafael show a close correlation with the calculated stress state results for sills, with an angular mismatch between the principal stress axes of ~4º; ϕ-values calculated for the faults are lower than for the sills, at 0.53-0.56(Fig. 5h, i).For comparison, Bingham analyses (Yamaji, 2016) of the local sill contact data derives an average ϕ-value of 0.63 and principal stress axes creating a 6º angular mismatch with the opening angle and Right Dihedra models (Fig. 8b).To further constrain the stress ratio, and the relationship between intrusions and contractional shear structures, we input the thrust fault, deformation band, gypsum vein, and overall sill geometry pole data into mechanical models of normalized slip tendency (Ts; Fig. 8c) and dilation tendency (Td; Fig. 8d).Although typically used to assess the reactivation potential of pre-existing structures in a present-day stress state (e.g.Ferrill et al., 1999), the models can also be used to fit a paleostress state to field data (e.g.Stephens et al., 2017).
Normalized slip tendency (Ts = (τ / σ n )/ Tsmax; Morris et al., 1996) and dilation tendency (Td = (σ 1 -σ n ) / (σ 1 -σ 3 ); Ferrill et al., 1999) are calculated from the stresses acting on a plane, or potential plane.This can be used to predict the attitude of both potential failure planes, and preexisting fractures that are susceptible to reactivation in a given stress state.Notably, the zone of high Ts (0.8-1.0) overlaps with the zone of high Td (0.8-1.0), suggesting a zone of potential extensional-shear.To enable conjugate shear failure, we used the same differential stress of 6T, where T = 3 MPa, as used previously.By fitting the zones of high slip tendency to the combined thrust fault and deformation band data, and zones of high dilation tendency to the overall sill geometry and gypsum vein data, we were able to derive a best-fit ϕ-value of 0.65, with a horizontal NE-SW σ 1 that plots within the 6° angular mismatch.The ϕ-value is consistent with the Right Dihedra (ϕ=0.53;Fig. 5h, i), Bingham analysis (ϕ=0.63;Fig. 8a) and our opening angle models (ϕ=0.77;Fig. 7).The field observations combined with the paleostress analyses suggest that sill emplacement took place during a state of mild horizontal radial compression.
A key difference between the strains recorded by sills and thrusts is the dominance of While the opening angle of dilated fractures has been used here to characterise the paleostress state during dilation of pre-existing cohesionless fractures, we envisage that it could also be applied to extension or extensional-shear fractures, formed via failure of intact rock.In these cases, we would expect the fracture network to resemble a mesh structure (e.g.Sibson, 1996), which comprises predominantly parallel (Mode I), and bimodal (conjugate extensionalshear) fractures, which dip at a low angle to σ 1 , with opening angles of <20° (e.g.Hancock, 1985;Ramsey & Chester, 2004).Importantly, the models presented in this study can be used to determine whether intrusive suites record changes in the tectonic stress state, or to identify fluid pressure pulses of varying magnitude in a single governing far-field stress state.These attributes have significant implications for improving our understanding of the development of past and present-day magmatic systems.

Conclusions
Our mechanical models build upon the work of Delaney et al. (1986) and Jolly and Sanderson (1997) to use fracture geometry and opening angles to derive the principal stress axes during sill emplacement, and provide crucial new constraints on the stress state and fluid pressure, applicable to dikes, sills, and veins.The geometry of sills in the SRSVF record a continuous horizontal shortening that is otherwise accommodated by contractional faults in the area.
Contoured regions on stereonets for opening angles suggest that in a high deviatoric stress state, it should be relatively rare for intrusions to be purely Mode I structures.In such settings, it should be commonplace for intrusions to accommodate a component of shear during their emplacement.Our opening angle model is particularly useful in determining paleostress states for regions where there is little brittle deformation (i.e.faulting), other than intrusions, and it may therefore present a useful and important tool in tectonic and magmatic studies.(Delvaux and Sperner, 2003).Calculated stress axes given as plunge/ trend measurements, ϕ-value is the stress ratio (σ 2 -σ 3 / σ 1 -σ 3 ).
dilation during sill emplacement, and compressional shear during thrusting.The range of local dips, and consistent vertical opening, suggests that the sills were emplaced by a combination of brittle failure and dilation of pre-existing structures (including thrust faults); where the fluid overpressure accommodates both shear and dilation.Consistent with the interpretation ofWalker et al. (2017), our field observations and mechanical model results suggest that the sills represent conjugate intrusions, which record the continuity of horizontal shortening during periods of elevated magmatic pressure.

Figure 1 .
Figure 1.Schematic diagrams to show the dilation direction (yellow arrow) and opening angle (µ) of fractures in different stress regimes: (a) normal-fault regime (where SH = σ v ); (b) thrust-fault regime (where Sh = σ v ); the plane attitude is described by θ, the angle between the normal stress and the maximum stress.(c) Example of crosscutting sills; the middle sill (outlined in blue) shows a consistent vertical opening direction: Mode I opening on horizontal planes, and extensional-shear opening on inclined planes.(d) Graph showing changes in vertical thickness, tv, and true thickness, t, with changes of apparent dip, for the sills shown in (c).

Figure 2 .
Figure 2.Figure 2. Schematic diagrams to show the relationship between opening mode and apparent shear.(ab) Initial pre-intrusion geometries with a marker unit oblique to the fracture/ fault plane.(c-d) Extensional-shear dilation of a pre-existing fracture and fault.(e-f) Mode I dilation of a pre-existing fracture and fault, creating an apparent shear offset.(g) Dilation directions and opening angles along linked segments can be used to infer the true dilation direction and fluid pressure during intrusion.μ(ap) = apparent opening angle.

Figure 3 .
Figure 3. Mechanical models showing contoured fracture-opening angles for various stress and fluid pressure conditions, using the equation from Delaney et al. (1986), projected onto lower hemisphere, equal area stereonets and Mohr Circles (colour scheme from Thyng et al., 2016).The dashed line on the Mohr Circles indicates the fluid pressure magnitude.Five increments of fluid pressure (ai-v) have been modelled with three stress states: (b) thrust-fault regime (σ v = σ 3 ); (c) strike-slip regime (σ v = σ 2 ); (d) normal-fault regime (σ v = σ 1 ).Colour contours correspond to pole-to planes.Poles that plot in the yellow zone (where μ = 0-10°) will display approximately plane-normal opening, poles with larger opening angles will display extensional-shear dilation and have a smaller true thickness than adjoining segments that plot in the yellow zone.

Figure 4 .
Figure 4. Location maps for the San Rafael Sub-Volcanic Field in Utah amended from Walker et al. (2017).(a) Digital elevation Model for Utah, showing major structural and depositional areas of the Colorado Plateau.Solid black line shows province boundaries.Dashed black line is a region of lower-crustal delamination and crustal thinning detailed in Levander et al. (2011); dashed white line is their outline of a downwelling body at 200 km depth, estimated from body wave tomography.(b) Aerial imagery for the San Rafael Sub-Volcanic Field (SRSVF)highlighting location and distribution of intrusive bodies; the white star marks the Cedar Mountain sill studied byRichardson et al. (2015).(c) Hillshaded digital elevation models for thick sills in the southern SRSVF, coloured to show extrapolated elevation data for sill top contacts (modified fromWalker et al., 2017).(d) Field photograph of the Last Chance Sills, which transgress through the stratigraphy; the upper sill cuts across the Carmel-Entrada Formation boundary.(e) Sills cut dike and cut, and abut against, subvertical fractures (sills outlined in black, dike in white).(f) Lower hemisphere stereographic projections show sill top contact polygon attitudes (extrapolated from c) as great circles, and poles to planes for each sill system (data fromWalker et al., 2017).

Figure 5 .
Figure 5. Sill geometry and paleostress analyses of deformation structures in the San Rafael.(a) km-scale segmented sills show en-echelon stepping consistent with conjugate faults (sills shaded red).Inset shows schematic interpretation.(b) 100 m scale sills showing NE and SW dips (sills outlined in black).(c) <100 m scale sills with bimodal dips (sills outlined in black).(d) Deformation bands and low-angle fractures parallel gypsum veins.(e) Gypsum veins with conjugate geometries.(f) Sill intruding and cut by a thrust fault and low angle fractures (sills outlined in black).(g-i)Lower hemisphere equal-area projections, showing: (g) gypsum veins and fractures, (h) thrust faults, and (i) deformation band data; principal stress axes calculated for (h) and (i) using the Right Dihedron method(Delvaux and Sperner, 2003).Calculated stress axes given as plunge/ trend measurements,

Figure 6 .
Figure 6.Sill-sill cross-cutting relationships.(a) Younger horizontal sill maintains a constant thickness when cutting older sill.(b-c) Inclined sill sections are thinner than horizontal sections.All sills show near-vertical opening directions, regardless of attitude.

Figure 7 .
Figure 7.The opening-angle mechanical model.(a) Field example of thin sills (<1m thick) displaying a range of local contact dips, and a consistent vertical opening.(b) Lower hemisphere, equal area stereonet with local sill contacts plotted as poles to planes, coloured relative to their opening angle, with Mode I and total dilational ellipses constructed.(c) Determination of principal stress axes and θ angles.(d) Mechanical model results using the stress ratio and driving pressure (ϕ=0.77,R'=0.68; see text for details) derived from (c).(e) 3D Mohr circle showing how θ1 and θ2 can be used to calculate relative fluid pressure; Mohr Circle is colour contoured for values of μ, local sill contacts are plotted as poles and coloured by their opening angle.

Figure 8 .
Figure 8. Sill geometry and paleostress analyses of deformation structures in the San Rafael.(a) Field example for comparison between sill and thrust fault geometry.(b-d) Lower hemisphere equal-area projections with calculated stress axes given as plunge/ trend measurements: (b) Local sill-host rock contact data, Bingham analysis (Yamaji et al., 2016) used to calculate principal stress axes, shown with 95% confidence regions and density contours (shaded region).(c) and (d) are mechanical models contoured for slip tendency (Ts) and dilation tendency (Td), respectively.