Earth and Planetary Science Letters Igneous sills record far-ﬁeld and near-ﬁeld stress interactions during volcano construction: Isle of Mull, Scotland

Sill emplacement is typically associated with horizontally mechanically layered host rocks in a near- hydrostatic far-ﬁeld stress state, where contrasting mechanical properties across the layers promote transitions from dykes, or inclined sheets, to sills. We used detailed ﬁeld observations from the Loch Scridain Sill Complex (Isle of Mull, UK), and mechanical models to show that layering is not always the dominant control on sill emplacement. The studied sills have consistently shallow dips (1 ◦ –25 ◦ ) and cut vertically bedded and foliated metamorphic basement rocks, and horizontally bedded cover sedimentary rocks and lavas. Horizontal and shallowly-dipping fractures in the host rock were intruded with vertical opening in all cases, whilst steeply-dipping discontinuities within the sequence (i.e. vertical fractures and foliation in the basement, and vertical polygonal joints in the lavas) were not intruded during sill emplacement. Mechanical models of slip tendency, dilation tendency, and fracture susceptibility for local and overall sill geometry data, support a radial horizontal compression during sill emplacement. Our models show that dykes and sills across Mull were emplaced during NW–SE horizontal shortening, related to a far-ﬁeld tectonic stress state. The dykes generally accommodated phases of NE–SW horizontal tectonic extension, whereas the sills record the superposition of the far-ﬁeld stress with a near-ﬁeld stress state, imposed by emplacement of the Mull Central Volcano. We show that through detailed geometric characterisation coupled with mechanical modelling, sills may be used as an indication of ﬂuctuations in the paleostress state. © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Sheet intrusions represent an important and directional volumetric addition to the crust, with dykes accommodating horizontal extension, and sills commonly accommodating vertical thickening. Although vertical and lateral mafic magma transport through the crust is traditionally associated with dykes, recent studies have shown that sills can dominate the transport and storage network (e.g. Airoldi et al., 2011Airoldi et al., , 2016Eide et al., 2016;Magee et al., 2016). Dykes and sills are commonly observed in close spatial association with each other, which has led to the interpretation that regional sill complexes are fed by dykes, despite few supporting observations of this transition. Dykes are commonly used to infer crustal extension, with the idealised dyke plane developing perpendicular to the axis of minimum compressive stress (σ 3 ; here we reckon compressive stress positive, in which σ 1 > * Corresponding author. E-mail address: tls15@le.ac.uk (T.L. Stephens). σ 2 > σ 3 ), and parallel to the σ 1 -σ 2 plane. If sills are fed by dykes then a rotation of principal stress axes is required; from a horizontal σ 3 for dyke emplacement, to a vertical σ 3 for sill emplacement. A vertical σ 3 is typically inferred to be the result of a local stress perturbation at the interface between anisotropic mechanical layers (e.g. Gudmundsson, 2011), rather than a farfield stress state with a horizontal σ 1 -σ 2 plane. However, several studies of intrusive systems that include sills have identified the potential for a far-field (tectonic) stress state control on intrusion geometry (e.g. England, 1988;Chaussard and Amelung, 2012;Walker, 2016). This relationship between intrusion geometry and stress state is important for several key reasons: (1) far-field horizontal compression and shortening during intrusion may serve to inhibit vertical ascent of magma (via dyking) toward the surface (e.g. Chaussard and Amelung, 2012); (2) the dominant mechanism for sill emplacement will play an important role in the placement and distribution of intrusions (e.g. Maccaferri et al., 2011);and (3) sills may serve as a record for phases of horizontal shortening in regions that are otherwise considered tectonically inactive (e.g. Walker et al., 2017), or in cases where the scale of observa- tion precludes characterisation of minor strains (e.g. using remote sensing and geophysical imaging techniques).
To consider sill emplacement controls, we present a structural characterisation for sills that cut basement and cover sequences in the Loch Scridain area of western Isle of Mull, Scotland (Fig. 1a): The Loch Scridain Sill Complex. Notably, Precambrian basement rocks in the study area are vertically bedded (e.g. Figs. 2, 3c), whereas the Mesozoic sedimentary, and Paleogene volcanic, cover sequences are horizontally bedded (Fig. 3a, b), presenting a rare opportunity to investigate the role of mechanical layering in controlling sill geometry. For the first time, we use slip tendency, dilation tendency, and fracture susceptibility mechanical modelling, based on field data for intrusions and host rock fractures, to constrain the stress state during sill emplacement. Our field observations and mechanical models reveal that sill emplacement required a radial horizontal shortening, which we attribute to stress (a) NE-SW striking dykes are discordant to subvertical Moine foliation. Annotation gives strike, dip, and dip-direction, dykes are tinted red for clarity. (b) Dyke segments are locally parallel to foliation and abut subvertical fractures, 20 cm notebook for scale (dyke is tinted red). (c) Opposing dyke segments step against a pre-existing shear zone within the Moine host rock (dyke is tinted red). (d) Segments are discontinuous laterally along foliation planes. (e) Equal area lower hemisphere stereonet of dykes in cover sequences and basement, plotted with Moine foliation. Rose plots of all dyke strikes (f) and dips (g), divided into 10 • bins. Each concentric circle denotes one dyke measurement in (f), and two measurements in (g). (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.) state fluctuations arising from inflation of the Mull Central Complex.

Definitions for sheet intrusions
Igneous sheet intrusions are defined using a mixture of geometry, and their relationship to bedding (Walker, 1993): A sill is defined as a sheet intrusion that is horizontal (≤10 • ), but also concordant to bedding; a dyke is a sheet intrusion that is vertical (>70 • ), but also discordant to bedding. Inclined sheets are discordant to bedding with a dip range of 10 • -70 • (e.g. Muirhead et al., 2012;Marsh, 2015). Cone sheets can have sill-like and dyke-like geometries, tracking the imposed near-field stress axes related to a major intrusive body (e.g. Gudmundsson, 2011): Broadly they are conical and should dip toward their central source.
As most dykes and sills can be considered as magma-filled cracks, here we apply a geometric description for sheet intrusions based purely on the attitude of the crack semi-axes (where crack semi-axis lengths are a ≥ b c, with c oriented normal to the crack plane) at the time of emplacement. Thus, we define a sill as a sheet intrusion with an a-b plane that dips 0 • -30 • , and dykes as sheet intrusions that dip within a range of 60 • -90 • . Sheet intrusions that fall between these angular ranges (i.e. with a-b plane dips >30 • but <60 • ) we describe as an inclined sheet. These definitions are important here because the intrusions we present, and others like them, do not fit the current definition scheme, but also because the association with bedding has led to an inference that layering is the cause of sill emplacement in almost all instances. Intrusions may show local deflections that cross our angular ranges, or may deflect over kilometre length scales (e.g., as is the case for saucer-shaped sills). For the purposes of this paper, we use the overall intrusion geometry at the observation scale, to separate sills, inclined sheets, and dykes.

Geological background
The Isle of Mull hosts three major stratigraphic groups (Fig. 1a): (1) Precambrian Moine basement metasedimentary rocks; (2) Mesozoic sedimentary rocks; and (3) Paleogene lavas. Our study focuses on the Loch Scridain Sill Complex (LSSC) in western Mull, which cuts the entire outcropping sequence. Basement metasedimentary rocks outcrop in the southwest of the island and comprise quartzites, micaceous psammites, pelites, and gneissose pelites that have been folded to form a regional, tight and upright, SSWand NNE-plunging syncline (Holdsworth et al., 1987). Bedding in the basement strikes NNW-SSE to N-S, with subvertical dips. Folding has produced a pervasive axial planar cleavage which, in the study area, is transposed with bedding. The basement also hosts three fracture sets that each strike at ∼90 • to the subvertical Moine foliations, and which are separated here based on their dip: (1) vertical to subvertical (71 • -90 • dip); (2) inclined (21 • -70 • dip); and (3) horizontal to subhorizontal (0 • -20 • dip).
The cover sequences include Mesozoic shales and sandstones, overlain by Paleogene lavas and volcaniclastic sedimentary rocks  Paleogene volcanic activity on Mull has been dated to within a 2.52 ± 0.36 Myr range during Chron 26r (Chambers and Pringle, 2001), beginning with the extrusion of the lava pile onto the Moine basement and Mesozoic cover at 60.6 ± 0.3 Ma (Chambers and Pringle, 2001). The lavas have a present-day cumulative thickness of ∼1800 m (Bell and Williamson, 2002), comprising individual units that are typically less than 15 m thick (Skelhorn, 1969), and have an average dip of 3 • E. The Moine basement is juxtaposed against the lava cover across the ESE-WNW-striking Assapol fault (Fig. 1), which pre-dates lava extrusion (Williamson and Bell, 2012).
Magmatic activity on Mull produced a suite of large igneous intrusions -the central complex -and sets of dykes and sills (England, 1992). The central complex comprises three volcanic centres, which form a combined elliptical outcrop at surface with its long axis trending NW-SE (Fig. 1a). Emplacement of the central complex resulted in a radial shortening that affected the basement and cover sequences in central and south-eastern Mull, for distances of ∼4 km from the contact, forming a circumferential fold and thrust belt ( Fig. 1a; Skelhorn, 1969;Mathieu and van Wyk de Vries, 2009). This deformation has not been reported in the western Mull study area, which is ∼15 km from the central complex (Fig. 1). Dykes were emplaced continually throughout the period of magmatism on Mull and show a dominant NW-SE strike across the region (England, 1992). The LSSC is emplaced within the Moine basement, and the Mesozoic and Paleogene cover sequences. The sills have undergone fractionation and crustal contamination, but have primary compositions comparable to tholeiitic lavas associated with Centre 1 (Preston et al., 1998;Kerr et al., 1999). Those lavas, and the Loch Ba Ring Dyke (associated with Centre 3) are dated at 59.05 ± 0.27 Ma and 58.48 ± 0.18 Ma, respectively (Chambers and Pringle, 2001), therefore constraining the timing for intrusion of the sills to a 0.57 ± 0.45 Myr period.

Observations
Paleogene dykes cut the full stratigraphic sequence on Mull. Within the Moine basement, dykes strike NE-SW and NW-SE, and dip ∼80 • (Fig. 2); no preferred dip direction is noted, and no cross-cutting relationships were observed between the two dyke sets. NE-SW striking dykes cut the Moine foliation at a high angle (Fig. 2a), parallel to NE-SW subvertical fractures. NE-SW dykes do not intrude the Moine foliation. NW-SE dykes locally strike parallel to the Moine foliation (Fig. 2b), but generally cut the foliation at a shallow angle (Fig. 2c). NW-SE dyke tips are generally tapered .
Steps along the dyke contacts are parallel to inclined NE-SW striking fractures in the area (42 • -58 • dip) (Fig. 2b, c). Dykes within the Paleogene lava sequence are near-vertical (dipping SW) and strike NW-SE. Segments show deflections to produce minor variations in dip through the sequence (72 • -88 • ; Fig. 3a), with deflections in dip occurring at lava unit contacts, and within individual lava units. Where dykes and sills are observed together, sills always cut the dykes (Fig. 3a). We found no examples of dykes deflecting to subhorizontal intrusions, nor vice versa.

Summary & interpretations
The mean and modal strike of dykes across Mull is NW-SE (Jolly and Sanderson, 1995). We infer that this represents the axis of the far-field maximum horizontal compressive stress (σ H ) during dyke emplacement, with the axis of minimum horizontal compression (σ h ) oriented NE-SW. This is consistent with broader studies of dykes in Mull and on the Scottish mainland (see Jolly and Sanderson, 1995;England, 1988.) Optimally orientated pre-existing discontinuities (i.e., jointing, bedding contacts) allow for magmatic intrusion at lower magma pressures. Tapered dyke tips indicate that dilation and intrusion was possible along the Moine foliation, however, the overall strike of dykes in the Moine is discordant to the foliation (Fig. 2). Where dykes abut subvertical (80 • -89 • dip), NE-SW striking fractures (i.e. ∼90 • from σ H ), we infer that magma pressure did not exceed the normal stress acting on the fracture planes. Within the lava cover, local variations in dyke dip suggest that bedding may have influenced dyke propagation, but this was not sufficient to cause sill formation on a regional scale. Notably, in all observed cases, sills cut the dykes. Based on our field observations, we infer that the dykes observed at outcrop were not the feeders to sills in the same section.

Observations
The Loch Scridain Sill complex (LSSC) intrudes the entire stratigraphic succession on Mull. At the kilometre scale sills display subhorizontal dips, with mean dip values of 1 • -25 • (Fig. 3); the sills do not dip towards a central focus. At the sub-kilometre scale, the sills are consistently discordant to host layering, and cut the stratigraphy at a low angle as a series of steps, producing a 'flat' (4 • -17 • dip; mean is 8 • ) and 'ramp' geometry (11 • -33 • dip; mean is 22 • ) (Fig. 3a, b).

Sills in the Moine basement
The sills cut the Moine foliations at a high angle, with dips of 2 • -13 • SE to SSE and 3 • -4 • NW to NNW (Fig. 1b). In some locations, NW and SE dipping segments intersect to create an acute angle, with a horizontal bisector trending ENE-WSW (Fig. 4). Sill contacts show steps against pre-existing structures, including fractures (Fig. 5), and the Moine foliation ( Fig. 6), but igneous material does not occur within the subvertical structures above or below the sill contact step (Figs. 5c, 6b, c). These steps appear to represent local shear reactivation of pre-existing host rock structures that were not dilated during magmatic inflation. Sill tips are parallel to, and hosted within, subhorizontal to inclined (6 • -48 • dip) fractures (Fig. 5b, e), and in places, shear offsets of bedding and foliation are observed across the sill contacts (Fig. 6f).

Sills in the cover sequences
Sills cut through the Mesozoic sedimentary rocks and Paleogene lava sequence at a low angle (4 • -25 • : a greater range of dips than those emplaced within the Moine; Fig. 3), with a ramp-flat geometry (Fig. 3a, b). Sills at the boundary between the Mesozoic sedimentary rocks and the overlying lavas are roughly parallel to the formation boundary for a short distance (Fig. 3b). Ramp-flat transitions occur within units, and at unit interfaces ( Fig. 3a, b). Within the lavas, the flat sections are not consistently parallel to the host rock layering and cross-cut multiple unit interfaces (Fig. 3a). Where flat segments correspond to a single unit they occur mostly within the massive lava core, internal to the lava units, rather than following the upper, lower, or internal contacts. The lavas host vertical cooling joints, but we found no examples of sills intruded along joints.

Summary & interpretations
Sills that cut the vertically bedded Moine basement show similar attitudes to those that cut horizontally bedded Mesozoic sedimentary rocks and lava cover sequences (Fig. 3), indicating that host layering was not the dominant control on sill emplacement. Sills do not consistently dip toward a central focus, and the main dip directions (NW and SE) are not directed towards the central complex; hence we discount a cone sheet model in explaining their emplacement. Sills have consistent shallow dips to the NW and SE, and in some cases intersect, to form a continuous sill, with the contact forming an acute angle about the horizontal plane ( Fig. 4). We interpret this geometry as representing a conjugate structure, similar to extensional shear faults and veins (Sibson, 2003).
The stratigraphy on Mull was host to several discontinuities that were present at the time of intrusion, including the vertical Moine bedding and cleavage, horizontal bedding in the Mesozoic sedimentary rocks and Paleogene lavas, and vertical cooling joints within the lavas. Only subhorizontal fractures were intruded during sill emplacement. Sills cut cover sequence horizontal bedding at a shallow angle, and cut or abut pre-existing steeply dipping structures (>33 • ), including bedding, tectonic cleavage, fractures, and lava-hosted cooling joints. In all cases, host rock offset markers across the sills indicate near-vertical opening, regardless of the sill dip. We infer that this opening direction represents the minimum principal stress, σ 3 . For vertical structures to remain closed, we infer that the magma pressure during emplacement did not exceed the horizontal normal stress: σ 1 and σ 2 were horizontal and greater than the magma pressure. We suggest that steps in sill contacts record stages of propagation where, during sill emplacement, dilation of subhorizontal fractures was accompanied by inflation and shear on vertical structures to form a through-going sheet. Based on our observations, we propose a conceptual model for sill emplacement (Fig. 7) similar to models for sill emplacement during horizontal tectonic shortening presented by Walker (2016) for the Faroe Islands, European Atlantic margin, in which: (1) horizon- (2) those horizontal and subhorizontal fractures were intruded as magma propagated from source ( Fig. 7b, b ); and (3) inflation of sill segments resulted in fracture propagation and linkage to facilitate through-going magma flow and composite sill geometries (Fig. 7c, c ). Initial propagation as offset segments may have been impeded by low cohesion or cohesionless pre-existing structures (e.g. Gudmundsson, 2011), leading to blunt sill tips at foliation and fracture planes. As the sill inflated, pre-existing vertical structures were re-activated via shear, until inflation facilitated linkage with adjacent sill segments, producing a local step in the sill geometry (Fig. 7). Linked sill segments are observed in both dip-parallel (Fig. 5) and strike-parallel ( Fig. 6) orientations, suggesting that magma propagated in all directions within the horizontal plane. The formation of sills in this model requires that fractures became linked to a magmatic source; those that did not, remained closed.

Estimating the stress state during sill intrusion
The geometry and opening directions of sheet intrusions provide an indication of the stress state (i.e. the attitudes and relative magnitudes of the principal stresses) during emplacement, similar to constraining paleostress from fault and fracture systems (e.g. Baer et al., 1994;Jolly and Sanderson, 1997). Here, we use mechanical models of slip tendency (T s ), dilation tendency (T d ), and fracture susceptibility (S f ), to derive paleostress states from our field characterisation of sills, and the inferred reactivation state of pre-existing structures during sill emplacement in the Moine basement (Fig. 8).
Slip and dilation tendency show the attitude of planes susceptible to reactivation, via shear or dilation respectively. A plane is susceptible to slip when the ratio of shear stress (τ ) to the normal stress acting on the plane (σ n ) is large. Slip tendency is normalised here (relative to the maximum possible slip tendency, T s(max) ) to enable comparison between stress states (Eq. (1); Morris et al., 1996): A plane is susceptible to dilation when the difference between σ 1 and the normal stress acting on the plane approaches the magnitude of differential stress (σ D , where σ D = σ 1 -σ 3 ; Eq. (2) ;Ferrill et al., 1999): Fracture susceptibility (Eq. (3); e.g. Mildren et al., 2002) represents the magnitude of fluid pressure change ( P f : either magmatic pressure within a crack, or hydrous pore-fluid pressure within a crack or the host rock pore space) that is required to cause shear reactivation. Fluid-driven shear reactivation is dependent on the shear and normal stresses acting on the plane, as well as its cohesion (here assumed = 0), and static coefficient of friction (μ s ): (e) all exploited fractures and exploited tension fractures; (f) all stepped sill contacts against steep fractures and Moine foliation; and (g) fractures measured >10 m and <10 m from the sill contact (note that data number reflects accessibility and exposure).
Reactivation via fluid-driven dilation would require a fluid pressure greater than the fracture susceptibility.
Intrusion is favoured where T d is high and fracture susceptibility is low. Shear reactivation is favoured where T s is high and fracture susceptibility is low. The variation of T s , T d , and S f for all possible fracture attitudes in three dimensional (3D) space can be visualised for any given stress ratio (φ, where φ = (σ 2 − σ 3 )/(σ 1 − σ 3 )). Fig. 8a-c shows contoured pole-to-plane values, for T s , T d , and S f in three horizontally compressive fault systems. All planes were modelled as cohesionless surfaces and assigned an appropriate static coefficient of friction (μ s ) for reactivation (i.e. 0.6; Sibson, 2003). Due to this, the fracture susceptibility models represent the minimum change in fluid pressure required for reactivation.
Fig. 8a-c shows that at low stress ratios (φ = 0; Fig. 8a), planes susceptible to slip (T s → 1.0) have a radial distribution about a horizontal axis, that is elongate in the direction of σ 2 . The zone of high dilation tendency (T d → 1.0) forms a girdle centred about the σ 3 axis and is also elongate in the direction of σ 2 , suggesting dilation across a range of horizontal to vertical fractures (promoting sill and dyke emplacement). When φ = 0.5 (Fig. 8b), high T s forms bimodal ellipses. The zone of high T d is elliptical about the σ 3 axis, which limits dilation to subhorizontal and inclined planes; the ellipse long axis points toward σ 2 . When φ = 1 (Fig. 8c), planes susceptible to slip have a radial distribution about the σ 3 axis, and the zone of high dilation tendency is equidistant from σ 2 and σ 1 .
It should be noted that the regions for high T s partially overlaps with the area of high T d (Fig. 8), suggesting that fractures dipping at a low angle to the σ 1 -σ 2 plane could be reactivated as mixed-mode (hybrid) fractures (cf. Fig. 6c-f). A good fit between the models and the field observations requires a paleostress state that favours shear reactivation of vertical structures and dilation of subhorizontal structures.
Sills in the basement and cover sequences show consistent vertical uplift of the host rock. As previously mentioned, we found no examples of sills fed by dykes, and in all observed cases, dykes are cut by sills (Fig. 3a). Based on the dominant sill dip directions (NW and SE), and the observed intersection of these dipping sill  sets (Fig. 4), the sills are interpreted as conjugate structures, comparable to conjugate faults and veins. Conjugate fault geometries form when σ D is ≥4 times the tensile strength (T ) of the material. For the purposes of our illustrative mechanical model, an intact rock tensile strength of 3 MPa is assumed (typical ranges from 0.5 to 6 MPa; e.g. Gudmundsson, 2011); σ D is set at 4T (12 MPa), with the aim of creating extensional shear fractures. Holness and Watt (2002) estimate an emplacement depth for the sills of 2-3 km, based on mineral paragenesis in the metamorphic aureoles. Here, we used a fixed emplacement depth of 2 km. The overburden comprises basalt lavas, siliciclastic sedimentary rocks, and metamorphosed siliciclastic sedimentary rocks. Bott and Tantrigoda (1987) used a density of 2.5 g/cm 3 for the sedimentary rocks. Our density measurements of basalt lava and Moine samples, gave average densities of 2.81 g/cm 3 and 2.64 g/cm 3 , respectively (see Appendix A). For our models, we used the average overburden density (2.65 g/cm 3 ), which provided a vertical stress (σ z ) of 52 MPa. Jolly and Sanderson (1997) suggested that σ 3 will plot in the centre of a stereographic pole cluster of dilated fractures. Poles to the overall sill geometry and exploited fracture data, therefore, define the approximate extent of the high T d zone (0.8-1.0). The data fits a slightly elongate NE-SW ellipse (Fig. 8d), indicating a horizontal, radially-compressive stress state during emplacement (0.5 < φ < 1.0). The ellipse long axis gives the relative position of σ 2 (NE-SW) and therefore indicates that σ 1 trends NW-SE. This σ 1 axis is consistent with the Mull dyke swarm (see Jolly and Sanderson, 1995). Steps in the sill contact against vertical fractures (which strike ∼E-W) and subvertical foliation (which strikes ∼N-S), indicate that these planes were not intruded during sill emplacement (Figs. 5,6), requiring that the horizontal stresses were greater than the magma pressure. The stereographic distribution of poles to the steps is near-quadrimodal (Fig. 8d), which requires a near-radial zone of low T d (0.0-0.2), suggesting φ > 0.5. We find that φ = 0.75 provides a best fit for the field data ( Fig. 8d; also see Appendices B and C). The σ 3 axis is parallel to the central axis of the dilational ellipse which plunges 85 • towards 315 • , and σ 1 trends NW-SE, plunging 5 • towards 135 • . The measured dykes in the study area plot in the zones of low T s , low T d , and high S f , demonstrating that the best-fit stress state for sill emplacement was not favourable for simultaneous dyke emplacement.

The effect of fluid pressure
Intrusion of pre-existing structures is possible when the fluid pressure ( P f ) exceeds σ n on the fracture plane. We can further consider the relationship between magma pressure and the applied stress, as a ratio of driving pressure (R ; Baer et al., 1994): Equation (4) shows that when P f is less than, or equal to, σ 3 no pre-existing fracture will open (R ≤ 0) and when P f is greater than σ 3 (R > 0) pre-existing fractures will begin to open (Baer et al., 1994). Intrusion is therefore favoured when magma pressure is greater than σ 3 . Fluid pressure can also be considered in terms of the pore fluid factor (λ v ): a ratio of the fluid pressure to vertical stress (i.e. λ v = P f /σ z ) (Sibson, 2003). Mechanical models were run for 5 fluid pressures from R = 0 to R = 1 (Fig. 9) by imposing lithostatic (λ v = 1; P f = 52 MPa) and supralithostatic (λ v > 1; P f > 52 MPa) fluid pressure to the φ = 0.75 stress state shown in Applying a fluid pressure equal to σ 3 caused no dilation, but allowed shear reactivation of subhorizontal and inclined planes ( Fig. 9a). Fluid pressure less than σ 2 (i.e. σ 3 < P f < σ 2 ) only allowed opening of shallowly-dipping planes with poles approximately parallel to the σ 3 axis, and shear reactivation of inclined to subvertical planes (Fig. 9b). Fluid pressures greater than σ 2 (σ 2 ≤ P f < σ 1 ) allowed opening of low-angle and inclined planes, as well as steep planes with poles parallel to the σ 2 axis (Fig. 9c, d).
In each case, shear reactivation is restricted to subvertical and vertical planes. When fluid pressure was equal to, or exceeded, σ 1 ( P f ≥ σ 1 ) planes of all attitudes were dilated (Fig. 9e). Field observations show that vertical and inclined discontinuities were not intruded, neither in the lava cover, nor the Moine basement (e.g. Figs. 5, 6), indicating that the fluid pressures could not have been greater than σ 2 (i.e. 52 MPa < P f < 61 MPa). The best-fit to the field observations is therefore Fig. 9b, which allowed dilation of subhorizontal planes, and shear along inclined and subvertical planes.

Layering as a control on intrusion geometry
Host rock lithology, mechanical layering, and pre-existing discontinuities within the crust are commonly invoked as the primary controls on brittle sill emplacement (e.g. Kavanagh et al., 2006;Galland et al., 2009;Maccaferri et al., 2010;Galerne et al., 2011;Gudmundsson, 2011;Muirhead et al., 2012;Eide et al., 2016;Magee et al., 2016). These models include Cook-Gordon delamination (Cook et al., 1964), stress barriers (e.g. Gudmundsson, 2011), elastic mismatch (Dundurs, 1969) and material toughness (e.g. Gudmundsson, 2011). Each model provides a local mechanism for σ 3 to rotate from horizontal to vertical, causing a dyke to deflect into a sill. Each mechanism is typically modelled in a hydrostatic, or low deviatoric, stress state where the only loading in the system is magmatic overpressure (e.g. Kavanagh et al., 2006;Galland et al., 2009;Gudmundsson, 2011). Increasing the fluid pressure under hydrostatic conditions would eventually result in opening of all cohesionless discontinuities, hence we expect to see abrupt deflections in intrusion geometry. The LSSC occurs at the same stratigraphic levels as dykes, and yet where observed together, the sills consistently cut the dykes; there are no observed instances of a dyke adjoining a sill. The sills were emplaced within anisotropic host rocks, including the vertically bedded and foliated Moine basement, which is host to subhorizontal and subvertical fracture sets, and horizontally bedded sedimentary and volcanic cover sequences; the lavas in the sequence are vertically jointed. Only subhorizontal fractures have been intruded during sill emplacement (Figs. 5,6). It is important to note that even in cases where distinct mechanical discontinuities exist in the horizontal plane, such as the contact between Mesozoic sedimentary rocks and the overlying Paleogene lavas, sills appear to cut the interface at a shallow angle (Fig. 3). The juxtaposition of the vertically-bedded basement, and the horizontally bedded cover, preclude stress rotation due to layering, as has been invoked elsewhere (e.g. Gudmundsson, 2011;Magee et al., 2016). Host rock mechanical layering does not present the primary control on intrusion geometry within the study area. On a local scale, discontinuities and lithology play important roles in controlling the emplacement mechanisms and morphology of initial sill segments, however the composite geometry is likely controlled by the remote stress state.
Numerical and scaled analogue models suggest that when a homogeneous host is subjected to a remote horizontal compressive stress, rotation of a feeder dyke may occur over several hundreds of metres to kilometres (Menand et al., 2010;Maccaferri et al., 2011). It is not clear from our field observations whether the LSSC is fed by sills, or by dykes, and the vertical limitation in exposure may preclude this observation.

Causes of horizontal shortening during the emplacement of the Loch Scridain sill complex
Sills in the western Mull study area record a period of radial horizontal compression that is not related to layering-induced stress rotations. The LSSC was emplaced between ∼59.05 ± 0.29 Ma and 58.48 ± 0.18 Ma (Preston et al., 1998;Kerr et al., 1999;Chambers and Pringle, 2001), during phases of rifting and sea floor spreading associated with the formation of the North Atlantic. Extension during that time was generally NW-SE, leading to the initiation of sea floor spreading at ∼55 Ma (Bell and Williamson, 2002), suggesting a NW-SE oriented σ 3 during the emplacement of the Loch Scridain sills. However, the European Atlantic continental margin is host to a suite of structures that accommodated NW-SE contraction during this time, such as anticline-syncline pairs and reverse faults within the adjacent North Rockall Basin, dated to ∼59.2-28.1 Ma (Tuitt et al., 2010), and the NW-SE-striking Mull dykes (Jolly and Sanderson, 1995). Our mechanical models show that the best fit stress state for the sills involves a NW-SE σ 1 , but also favours a radial horizontal compression where φ = 0.75. Dykes were emplaced throughout the volcanic and magmatic period on Mull (England, 1992), however, in the western Mull study area, sills cut dykes in all instances.
Dyke emplacement requires a horizontal σ 3 axis, which in Mull  Fig. 8d). Planes to the left of the reactivation envelope (red line) can be reactivated in shear (yellow zone) or via dilation (grey zone), or a combination; planes to the right of the envelope (white zone) remain stable. The fluid pressure magnitude can be read from the graph where the contact between the zones of shear and dilation intersects the normal stress axis. Field data for stepped sill contacts, exploited fractures and overall sill geometries are plotted to show the range of structures that can be dilated at each fluid pressure. (For interpretation of the colours in this figure, the reader is referred to the web version of this article.) should be oriented NE-SW (Jolly and Sanderson, 1995;their S h ). Jolly and Sanderson (1995) inferred a maximum principal stress, associated with the dykes, oriented NW-SE and horizontal (their S H ), which is coaxial with our inferred σ 1 . We assume then, that σ 1 remained horizontal throughout dyke and sill emplacement. Therefore, a mechanism for rotation of σ 2 and σ 3 must be considered. We have demonstrated that layering is not a suitable mechanism for causing the necessary stress rotations in the study area; as such, we explore two additional possible causes for stress rotation: (1) dyke-induced, and (2) emplacement of the central complex.
Forceful injection of magma as dykes, in cases where the volume of added material is greater than the tectonically-induced extensional strain, can result in an induced, near-field, horizontal compression related to the directional addition of material, and lead to a switch in the principal stress axes (e.g. Anderson, 1951).  (Fig. 9b). Graphs show the change in magnitude of σ x , σ y , and σ z , with distance from Centre 1 to the SW (d) and NW (f), with a 6 MPa excess pressure. The shaded zones indicate a thrust fault regime (grey), strike-slip regime (white), normal fault regime (green), and the extent of sill exposure (red). Mohr diagrams show the calculated stress state at a distance of 3 radii, they indicate that conditions for sill emplacement occur in the SW (φ = 0.75, σ D = 4T ) (e), and for dyke emplacement in the NW (φ = 0.02, σ D ≈ 7T ) (g). (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.) The forceful injection of the NW-SE striking Paleogene dyke swarm on Mull could have induced horizontal shortening of the host rock, normal to the dyke margins, i.e. in a NE-SW direction. It is therefore possible that the added volume represented by the dyke swarm resulted in a rotation of the stress axes with a switch of the σ 2 and σ 3 axes; the σ 1 axis remained horizontal (England, 1988).
Based on dyke thicknesses, Jolly and Sanderson (1995) estimated between 2% and 20% crustal extension was accommodated by dyk-ing. Assuming the dykes were emplaced by forceful injection, and not as a response to crustal extension, the compressional effect of this additional volume should be at its greatest in areas of relatively dense dyke spacing (i.e. where dyke frequency is at its greatest), and less pronounced where the density of dyking decreases (i.e. where spacing increases). A large number of sills have been mapped in the west of Mull, where dyke spacing is very large. Very few sills have been mapped in the North of Mull, where dyke spacing is smaller (Fig. 1). This does not, therefore, support a dyke-induced stress rotation.
An alternative cause of horizontal shortening is presented by the central complex, which has been associated with a radial compression and shortening, forming a fold and thrust system circumferential to Centre 1, affecting a region up to ∼4 km away from its margins (Mathieu and van Wyk de Vries, 2009). Forceful emplacement of a large magmatic body into the crust can induce an outward radial (axial) compressive stress (σ rr ), with an induced tensile hoop (circumferential) stress (σ θθ ) (Fig. 10a). Here, we present a conceptual model to illustrate these effects superimposed onto an anisotropic far-field stress state; we consider the horizontal plane at the level of the Moine-hosted sills (at 2 km depth). The magma chamber is simplified as a completely molten, and instantly-emplaced cylinder. The basic formulation for the axial and hoop stresses at distance d from the central point of a fluid-filled body of radius r, can be derived from the Kirsch solutions for stress concentrations around a pressurised cylinder that is subject to a far-field anisotropic stress state, given by equations (5) and (6) (Jaeger et al., 2007): where P m is the total magma pressure (i.e. P m = excess pressure ( P e ) + lithostatic (overburden) pressure ( P l )), and θ is the clockwise angle from the x-axis to the line measuring distance (Fig. 10b). The vertical axis of the cylinder is parallel to the vertical stress (σ z ), and its x-and y-axes are aligned with σ x (NW-SE) and σ y (NE-SW), respectively. The radial stress is positive and at its maximum, at the intrusion wall (r = d). Once the intrusion has inflated, r may remain relatively constant, and since each part of Equations (5) and (6) requires a division by d, the induced perturbation stress will decay with distance from the magmatic body ( Fig. 10d, f). We consider a shield volcano above Centre 1 with a height of 1.5 km and radius of ∼11 km that extends to the fold and thrust region; beyond this region, we model a constant depth of 2 km (σ z = 52 MPa). This edifice would contribute an additional vertical load of 39 MPa, hence, we consider a conservative value for the magma pressure at 3.5 km depth, at 91 MPa (i.e. P e = P l ). Magma chambers are likely to rupture when the magma excess pressure equals or exceeds the tensile strength of the wall rock (Browning et al., 2015). For our conceptual model, we use a maximum excess pressure of 6 MPa (2T ). We find that increasing the excess pressure from a lithostatic value ( P e = 0 MPa, P m = 91 MPa) to an overpressure of 2T ( P m = 97 MPa) has a negligible effect on the principal stress magnitudes and stress perturbations around the chamber. As such, in Fig. 10 we only present the results for an excess pressure where P e = 2T . The sills are predominantly exposed beyond the radial fold and thrust region to a distance of d ≈ 24.5 km (Fig. 10a). The studied sills in the Moine and Lava are located about the 3 radii marker (d ≈ 21 km; Fig. 10a). Horizontal compression induced by Centre 1 would be directed radially and should therefore be evidenced in the north of Mull, as well as in the southwestern study area. To the SW of Centre 1 (θ = 90 • ), σ rr would act along a NE-SW axis (Fig. 10a, d, e), parallel to the inferred σ 3 axis for dykes in the region. The negative hoop stress would act in the vertical and along a NW-SE axis, parallel to the inferred σ 1 orientation. Using a far-field vertical radial extension (φ = 0.25), with a differential stress of 5T (where σ H (NW) = 63.25 MPa and σ h (NE) = 48.25 MPa), we find from Equations (5) and (6), that beyond the edifice, to the SW of Centre 1, a zone of horizontal shortening (thrusting) occurs that correlates with the zone of sill exposure (Fig. 10d). At a distance of 3 radii, this gives a stress ratio (φ) of 0.75 (Fig. 10d, e). With increasing distance, the phi value drops as the influence of the hoop and radial stress decrease and the stress ratio returns to the far-field value (φ = 0.25; Fig. 10d, e). Beyond the edifice, to the NW of Centre 1 (θ = 0 • ), the stress state is a consistent vertical radial extension with a NWdirected σ 1 (Fig. 10f). At the 3 radii marker, this produces a stress ratio of φ = 0.02 (Fig. 10f, g).
Based on our conceptual model, we envisage that emplacement of Centre 1 could have facilitated a switch in stress axes to the NE and SW of the centre, by locally increasing the horizontal NE-SW stress with addition of σ rr , and lowering the NW-SE and vertical stresses with addition of a negative σ θθ (Fig. 10d, e). The model indicates that in the SW, at 3 radii, the magnitude of σ 1 would be approximately equal to that of the far-field stress state.
At this distance, only σ 2 and σ 3 are switched due to the change in their relative magnitudes (Fig. 10d, e) promoting the emplacement of subhorizontal intrusions. In the north of Mull, σ rr would act along a NW-SE axis, parallel to, and adding to the magnitude of, σ 1 (Fig. 10f, g). The negative σ θθ in the vertical and NE-SWhorizontal axes would have the effect of reducing the magnitude of σ 2 and σ 3 such that their axes should remain constant, promoting dyke emplacement (Fig. 10f, g).
The induced near-field stress is sensitive to a number of independent parameters. Our illustrative model uses a low-angle conical edifice, a cylindrical magma chamber that was emplaced instantaneously, and which was subject to a far-field strike-slip stress state. The hosting material is treated as isotropic despite the presence of pre-existing structures and foliations. The model presents stress as a function of magma chamber radii, which would account for changes in the stress distribution owing to the size of the chamber. Material elastic properties can have an effect on the induced hoop stress, but this is minor compared to the effect of changing the shape of the chamber (e.g. Sartoris et al., 1990;Martí and Geyer, 2009;Davis et al., 2017), or the height of the volcano over the magma chamber; the latter has a marked effect on gravitational potential energy, and the induced compressive stress at the toe edge of the cone (e.g., Mathieu and van Wyk de Vries, 2009). Additionally, the model is also sensitive to the far-field stress state; switching between σ 2 and σ 3 is favoured when the differential stress is low (<6T ) and the φ-value is low (approaching 0), and becomes more unlikely as the φ-value and differential stress increase: σ 2 becomes much greater in magnitude relative to σ 3 .
Based on our illustrative model, we favour the inflation of Centre 1 as the cause of the stress rotation that enabled sill emplacement. Without the presence of a far-field NW-SE directed shortening (NE-SW extension), emplacement of the central complex would favour the formation of radial dykes, and cone sheets below the edifice. Our study demonstrates the importance of intrusion geometry in constraining the stress state during emplacement, and the relatively minor effect of pre-existing discontinuities, particularly the contacts between mechanical layers, in controlling the overall intrusion form.

Conclusion
The Loch Scridain Sill Complex records a previously unrecognised phase of radial horizontal compression, representing a perturbation in the far-field stress state. Cross-cutting relationships, and paleostress analysis of intrusion geometry, suggest that tectonically-controlled dyke emplacement, and NE-SW extension, was overprinted by emplacement of the Mull Central Complex, leading to a radial compression, and sill emplacement. The directional dependence imposed by superposition of the NW-SEdirected maximum tectonic stress, and the induced radial magmatic stress of the central complex, may have led to inhibited vertical ascent of magmas in the southwest of Mull, whilst promoting vertical ascent in the northwest and southeast. The results therefore have major implications to other systems in which volcano construction may overprint the far-field stress state (e.g., the East Africa Rift). This study shows how the combination of sill geometry and mechanical modelling can be used to interpret paleostress during intrusion and highlights the importance of field-based characterisation in constraining sill emplacement mechanisms.