The building blocks of igneous sheet intrusions: insights from 3D seismic reflection data

Dykes, sills, and inclined sheets may intrude as continuous, planar sheets, but can also form by the 15 amalgamation of discrete injections of magma. Describing the host rock deformation style around distinct geometric components within sheet intrusions can help us to determine emplacement mechanisms, how magma flow influences accumulation of economic resources, and how magma migration may translate into ground surface deformation prior to volcanic eruptions. Geometric components, such as ‘magma fingers’ and ‘segments’, form during sheet propagation by brittle 20 and/or non-brittle emplacement mechanisms (i.e., shear failure, host rock fluidization), and fracture segmentation (i.e., elastic-brittle fractures), respectively. Seismic reflection data provide unique opportunities to map the 3D geometry of igneous sills over large areas and in high spatial detail, allowing us to identify these distinct geometric components and constrain magma flow patterns. However, limitations in the resolution of seismic reflection data mean we cannot discern 25 the host rock deformation structures that allow the origin of these geometric components to be interpreted. Here, we use 3D seismic reflection data located offshore NW Australia to: i) introduce new terminology that defines the building blocks (i.e., elements) of sheet-like igneous intrusions and their connectors (i.e., H-type, S/Z-type), without tacitly linking their description to emplacement mechanisms; ii) quantify the 3D geometry of these elements and their connectors in 30 two sills; and iii) discuss how element and connector geometries can be related to emplacement processes. For example, our measurements of connector heights along connectors exclude synemplacement rotation of principal stress axes as a primary segmentation mechanism for the host Submitted for publication in Geosphere 2 sills. Based on seismic attribute analyses and the 3D geometry of elements, we conclude that potential brittle and/or non-brittle magma emplacement processes led to sill segmentation and to the formation of magma fingers. We show thickness varies across sills, and across distinct elements, which suggests that they grew and inflated independently with individual flow kinematics. Our data of element geometries and their connectors: i) permit comparison between 5 different subsurface and field datasets spanning a range of host rock types and tectonic settings; and ii) allow us to test predictions of numerical and physical analogue models, helping us better understand magma emplacement mechanisms. INTRODUCTION 10 Magma typically moves large distances laterally and vertically through the crust via networks of interconnected sub-vertical and sub-horizontal sheet intrusions (i.e., dykes and sills) (e.g., Anderson, 1937, 1951; Ernst et al., 2005; Muirhead et al., 2012; Magee et al., 2016; Cruden and Weinberg, 2018). For example, interconnected mafic sills and inclined sheets can transport magma laterally over >4100 km and vertically up to 12 km (Cartwright and Hansen, 2006; Leat, 2008; 15 Magee et al., 2016). Understanding the architecture of these magma plumbing systems and how magma migrates through them is essential because they influence: i) volcanic hazards and their warning signals (e.g., Sparks, 2003; Cashman et al., 2013); ii) carbon dioxide and sulfur dioxide degassing to the atmosphere (e.g., Hernández et al., 2001; Schmidt et al., 2015; Ilyinskaya et al., 2017); and iii) the accumulation of economic resources (e.g., Barnes et al., 2016; Spacapan et al., 2

resolution is often sufficient to allow us to identify these distinct elongate components and thereby map magma flow patterns over entire intrusion networks. Yet seismic resolution is limited so we typically cannot discern the centimeter-to-meter scale host rock deformation structures that would allow the origin of these components to be interpreted. Here, we introduce a new term that defines 25 the components (i.e., 'elements') of sheet-like igneous intrusions, without linking their description to emplacement mechanisms. Using 3D seismic reflection data from offshore NW Australia, we quantify the 3D geometry of these elements and their connectors within two sills and discuss how their shape may relate to emplacement processes. Based on seismic attribute analyses and our measurements of their 3D geometry, we conclude that the mapped elements likely formed by non-30 elastic brittle and/or non-brittle deformation ahead of the advancing sill tip, implying they are magma fingers. We show that thickness varies across sills, and across distinct elements, which we infer to represent flow localization and subsequent thickening of restricted areas. The quantification of element geometries is useful for comparisons between different subsurface and field-based datasets, spanning a range of host rock types and tectonic settings. This in turn facilitates the testing of magma emplacement mechanisms and predictions from numerical and physical analogue experiments.
Critically, the development and presence of elements can drive localization of magma flow within 25 individual intrusions, thereby controlling the distribution of magma in the subsurface (Guo et al., 2013;Magee et al., 2014). To understand the initiation and propagation mechanisms of elongate elements, as well as their flow kinematics and dynamics, it is crucial to link and compare numerical/laboratory experiments and analytical solutions with observations made in nature.
Despite the clear benefit and importance of an integrated, multidisciplinary approach to 30 understanding magma emplacement dynamics, we are currently faced with several challenges: i) we do not yet have fully coupled, fully 3D, numerical or laboratory models of magma emplacement that consider multiple emplacement mechanisms; ii) current numerical and laboratory models of sill emplacement mainly focus on the formation of saucer-shaped sills (e.g., Malthe-Sørenssen et al., 2004;Galland et al., 2009;Walker and Gill, 2020); iii) materials used as host rock analogues in laboratory models are limited in their rheological behaviors, resulting in prescribed magma 5 emplacement mechanisms; and iv) observations made in the field are often limited to relatively small outcrops compared to the whole intrusion, and in many cases only represent a twodimensional (2D) view of a 3D intrusion (e.g., Pollard et al., 1975). Seismic reflection data are a powerful tool to bridge the gap between observations made in analogue models and nature, as they allow mapping of sill morphologies in 3D down to tens-of-meters scale, which is at the upper limit 10 of elements that have been mapped in the field . However, in many cases sill thicknesses are below the vertical seismic resolution, meaning that sills are visualized in the seismic reflection data but their thickness and the size and style of surrounding host rock structures cannot be constrained (e.g., Thomson and Hutton, 2004;Hansen and Cartwright, 2006a;Thomson, 2007;Miles and Cartwright, 2010); this adds uncertainty to the interpretation of 3D intrusion 15 geometries. A lack of borehole data also means we often cannot directly determine intrusion thicknesses, and host rock lithologies and associated deformation mechanisms.
Here we use 3D seismic reflection data from the Exmouth Plateau, offshore NW Australia to study the geometry of igneous sills, focusing on their building blocks. More specifically we: i) introduce new descriptive terminology for the building blocks of sills (i.e., 'elements'); ii) quantify 20 the dimensions of 3D element geometries and their connectors; iii) identify intrusion geometries considered to be representative of flow kinematics; and iv) combine this information to produce a conceptual model of sill emplacement dynamics and related deformation mechanisms. By comparing sill thickness and reflection amplitude variations, we find that low-amplitudes within the Top-sill seismic reflection correspond to connectors between adjacent elements; these features 25 can be used to map elements both above and below the vertical seismic resolution. Distinct thickness variations within adjacent elements suggest that they grew and inflated independently.
Where adjacent elements are vertically offset, they can form connectors resulting in a step-like sill geometry. Following the assumption that vertical connector heights increase with increasing element length when fracture segmentation occurs due to a rotation of the principal stress axis 30 orientations (Fig. 1D) (e.g., Pollard et al., 1975;Rickwood, 1990;Hutton, 2009), we conclude that branching of the sills presented in this study was not caused by stress rotation. Our quantitative data of element geometries and their connectors can be used to compare different subsurface and field-based datasets with a range of host rock types, magma compositions, emplacement depths, and tectonic settings. The geometrical data presented here also allow us to test predictions from numerical and laboratory analogue experiments, contributing to a better understanding of sill 5 emplacement mechanisms at shallow depths in the Earth's crust.

ELEMENTS, MAGMA FINGERS, SEGMENTS, LOBES, AND THEIR CONNECTORS
Over the last few decades, several studies have shown that sills and dykes do not always propagate 10 as continuous planar sheets, but that they can branch into laterally and/or vertically offset magma fingers or segments at their propagating edges, depending on the syn-intrusive rheological behavior of the host rock (e.g., Pollard et al., 1975;Rickwood, 1990;Hutton, 2009;Schofield et al., 2010Schofield et al., , 2012aSpacapan et al., 2017;Galland et al., 2019;Magee et al., 2019b). As magma intrusion continues and these features inflate, they coalesce to form a continuous sheet intrusion 15 (e.g., Pollard et al., 1975;Baer and Reches, 1987;Hutton, 2009;Schofield et al., 2010Schofield et al., , 2012aGalland et al., 2019). We briefly summarize the formation and propagation of these features below in order to provide an overview of the building blocks of sheet intrusions, their geometries, and associated host rock deformation. Previously used structural terms for these features within sheet intrusions have acquired a genetic connotation and are therefore no longer descriptive. Because 20 the information required to demonstrate how these features form is commonly not resolved in seismic reflection data, we introduce new, descriptive terminology (i.e., 'element') that we propose should be used prior to the interpretation of emplacement processes. Mapped elements interpreted as magma fingers in a sill imaged in 3D seismic reflection data located in the Rockall Trough (modified after Thomson and Hutton, 2004;figure from Magee et 5 al., 2015), and in the Nequén Basin, Argentina (modified after Galland et al., 2019). (B) Segments; Opacity render of a sill located in the Flett Basin, NE Atlantic with large scale elements highlighted (dashed line). Corresponding seismic sections are interpreted to show segments and the growth of intrusive steps and bridge structures along their boundaries (A'-C) (modified after Schofield et al., 2012b). Field photograph shows vertically offset segments connected via steps emplaced into 10 sedimentary strata on Axel Heiberg Island, Canada (photo curtsey of Martin Jackson). (C) Schematic diagram showing multiple elements forming an element network. (D) Schematic diagram detailing how a rotation of principal stress axes orientation leads to fracture segmentation and the formation of vertically offset segments with parallel intrusion planes (after Hutton, 2009;Magee et al., 2019b;figure modified after Pollard et al., 1975). Note that interpreting elements as magma fingers or segments based on seismic reflection data can only be done if sufficient 5 information about according host rock deformation is available. To avoid misinterpretation, we suggest the use of the descriptive term element (refer to text for further discussion).
During sill emplacement, two different mechanisms can lead to branching of the propagating edge and the formation of discrete segments: i) a rotation of the principal stress axes 10 ahead of the propagating tip, resulting in en-échelon fractures with a consistent stepping direction (e.g., Pollard et al., 1982;Nicholson and Pollard, 1985;Takada, 1990;Schofield et al., 2012a;Magee et al., 2019b) (Fig. 1D); and ii) the exploitation of pre-existing, preferentially oriented weaknesses with a lower tensile strength and/or fracture toughness (e.g., bedding planes, fault planes), likely leading to an inconsistent stepping direction (e.g., Hutton, 2009;Schofield et al., 15 2012a;Magee et al., 2013b;Stephens et al., 2017;Magee et al., 2019b). Since adjacent segments propagate on vertically offset horizons, they either overlap or underlap and are originally separated from each other by 'bridges' of host rock (Fig. 1B) (e.g., Hutton, 2009;Schofield et al., 2012a;Magee et al., 2019b). Subsequent inflation may cause segments to connect due to brittle failure of the host rock bridge, forming a continuous sheet with a step-like geometry ( Fig. 1B) (e.g., 20 Rickwood, 1990;Hutton, 2009;Schofield et al., 2012aSchofield et al., , 2012b.
The thickness-to-width ratios of segments tend to be small compared to magma fingers. If the host rock deforms in a linear elastic manner, then: where W = characteristic thickness, P = magma overpressure relative to the minimum in situ stress, 25 L = characteristic crack length (radius for a circular sill; width (shorter dimension) for a segment), and the stiffness parameter E' = E/(1-v 2 ) where E = Young's modulus and v = Poisson's ratio (Olson, 2003). If crack propagation is governed by LEFM, then: where KC = fracture toughness and PC = critical overpressure (Olson, 2003). Under the assumption that the crack is propagating, $ = $ % , such that equation (1) and (2) can be combined. A range of reasonable values of KC and E' for host rocks such as sandstone (1.5 x 10^6 Pa m 1/2 and 20-55 GPa, respectively) or shale (1.4 x 10^6 Pa m 1/2 and 1-20 GPa, respectively) gives an upper W/L aspect ratio of ~0.00001-0.0004. 5

Elements and Their Connectors
Structural terms such as magma fingers and segments already imply potential emplacement mechanisms. However, magma fingers and segments cannot be distinguished in map-view, since information on host rock deformation mechanism and intrusion 3D geometry are missing. We 10 therefore propose the term element to generally describe the building blocks of sills (i.e., magma fingers, segments), without attempting to infer their 3D geometries or emplacement mechanisms.
Adjacent elements may be vertically and/or horizontally offset and separated by host rock, which we refer to as a host rock separator (Fig. 1A, B). When overlapping or underlapping elements inflate, brittle failure of the host rock separator may cause both elements to connect, forming 15 broken bridges and steps, respectively (e.g, Rickwood, 1990;Hutton, 2009;Schofield et al., 2012a). Broken bridges and steps have been used as evidence for magma propagation by tensile brittle fracture mechanisms (i.e., fracture segmentation) (e.g., Hutton, 2009;Schofield et al., 2012a). However, recent field studies show that magma fingers that propagated via brittle and/or non-brittle processes can also be vertically offset and connected, forming a step-like geometry 20 ( Fig. 1A) . We introduce the term H-connector, to describe connectors between overlapping elements (formerly known as a 'broken bridge') (e.g., Nicholson and Pollard, 1985;Schofield et al., 2012b) and S-, or Z-connector for underlapping elements (formerly known as 'steps') (e.g., Rickwood, 1990;Magee et al., 2019b), with a right-upwards and right-downwards oriented stepping direction, respectively, depending on the viewing direction (Fig. 1B). These 25 connectors form parallel to the long dimension of elements, such that their long axis can be used to determine the propagation axis of the elements, which in turn can help to identify potential magma feeder zones (e.g., Schofield et al., 2012b;Magee et al., 2014;Schofield et al., 2017;Galland et al., 2019;Magee et al., 2019b).
So far, we have described the smallest building blocks of sheet intrusions (i.e., magma 30 fingers and segments) and placed them within our new element-based terminology. In addition, numerous studies have shown that sheet intrusions can occasionally be sub-divided into larger sections, each created by the amalgamation of an element set (e.g., Thomson and Hutton, 2004;Hansen and Cartwright, 2006a;Miles and Cartwright, 2010;Schofield et al., 2012b;Magee et al., 2016a;Schofield et al., 2017;Fyfe et al., 2021). For example, where elements define a group of fanning magma fingers and/or segments, these larger sections have been termed lobes (e.g., 5 Thomson and Hutton, 2004;Hansen and Cartwright, 2006a;Schofield et al., 2012b). We designate larger sheet components, including lobes, also as elements and highlight that they can be classified into different orders. In this way, first-order elements comprise a group of smaller second-order elements, which potentially can be sub-divided into third-order elements, where magma fingers and segments represent the smallest end-members in the element hierarchy (Figs. 1A, C, and 4A). 10 A large, first-order element consisting of multiple orders of sub-elements regardless of their geometries is termed an element network (Fig. 1A, C).
The geometry of any element in map view has previously been shown to typically be lobate or finger-like (e.g., Thomson and Hutton, 2004;Hansen and Cartwright, 2006a;Miles and Cartwright, 2010;Schofield et al., 2012b). In this study, we classify element geometries into three 15 groups: i) elongate elements with parallel lateral edges (Type 1); ii) elongate elements that are widening in magma flow direction (Type 2; Fig. 4B); and iii) semi-radially fanning elements that rapidly increase in width (Type 3; Fig. 4B). All three types are scale-independent and are only used to describe the geometries, without implying their scale or hierarchy. Adjacent elements can form in response to: i) changing magma flow kinematics (e.g., Magee et al., 2016b); ii) discrete magma 20 injections (e.g., Thomson and Hutton, 2004;Horsman et al., 2016;Magee et al., 2016a); or iii) cycles of magma front cooling, inflation, and magma break-out (e.g., Hansen and Cartwright, 2006a;Miles and Cartwright, 2010;Currier and Marsh, 2015). Overall, the term 'element' is a purely geometrical description where the hierarchy refers to the scale of elements, thus excluding any interpretation of emplacement timing or deformation mechanisms. 25

GEOLOGICAL SETTING
The Exmouth Plateau is part of the North Carnarvon Basin, which is located on the rifted margin of northwest Australia and formed by multiple episodes of extension during the Late Paleozoic and Mesozoic ( Fig. 2A) (e.g., Symonds et al., 1998;Karner and Driscoll, 1999;Longley 30 et al., 2002;Stagg et al., 2004). The Exmouth Plateau is composed of a ~10 km thick sedimentary sequence atop a ~10 km thick basement of stretched continental crust (e.g., Exon and Buffler, 1992;Symonds et al., 1998). Formation of the Tethys Ocean during the Permian resulted in extensive stretching and thinning of the crust, followed by rapid subsidence and the deposition of fluvial deltaic siliciclastic pre-rift sediments of the Triassic Locker Shale and Mungaroo Formation, which form the host rocks of five sills in the study area (Figs. 2B-D) (e.g., 5 Exon and Buffler, 1992;Symonds et al., 1998;Stagg et al., 2004). Sedimentary rocks of this pre-rift sequence are offset by steeply dipping northeast-southwest striking normal faults that formed in the earliest rifting phase of the Exmouth Plateau during the Late Triassic-to-Early Jurassic break-up of Argoland from Greater India (Figs. 2C and D) (e.g., Exon and Buffler, 1992;Stagg et al., 2004;Black et al., 2017;Bilal et al., 2018). 10 Rift-related magmatism in the North Carnarvon Basin spanned the Late Jurassic-to-Early Cretaceous ( Fig. 2D) (e.g., Mihut and Müller, 1998;Symonds et al., 1998;Rey et al., 2008;Magee and Jackson, 2020). A large intrusive body, characterized by high seismic velocities that is inferred to comprise mafic-to-ultra mafic rocks, was emplaced into the lower crust ~165-136 Myrs ago, potentially acting as the source for the extensive network of sills and dykes observed across the 15 area ( Fig. 2D) (e.g., Mihut and Müller, 1998;Rey et al., 2008;Rohrman, 2013;Magee and Jackson, 2020). Dating of forced folds, hydrothermal vents, and fluid escape structures, all interpreted to be related to sill intrusion, suggest that sill emplacement across the North Carnarvon Basin started in the Kimmeridgian (~157 Ma; Fig. 2D) (e.g., Rey et al., 2008;Magee et al., 2013a;Rohrman, 2013;Magee et al., 2017;Mark et al., 2019;Norcliffe et al., 2021). These sills are occasionally cross-cut 20 by radiating dykes of the Tithonian (~152-147 Ma) Exmouth Dyke Swarm ( Fig. 2D) (Magee and Jackson, 2020). A final phase of magmatism was associated with continental break-up and the development of the continent-ocean transition zone outboard of the Exmouth Plateau ( Fig. 2D) (e.g., Symonds et al., 1998;Rey et al., 2008;Magee et al., 2013b). Our study focuses on Late Jurassic-to-Early Cretaceous sills that were emplaced into predominantly Triassic pre-rift 25 sedimentary rocks.

Seismic Reflection Data and Seismic Resolution
This study uses a zero-phase, time-migrated, 3D seismic reflection survey (Glencoe 3D) covering an area of approximately 4042 km 2 (Fig. 2). In-lines and cross-lines are oriented E-W and N-S, respectively and have a regular spacing of 25 m. These data image to a depth of 8 s two-way time 5 (TWT) and are displayed with SEG positive polarity; i.e., a downward increase in acoustic impedance corresponds to a positive reflection (red), whereas a decrease corresponds to a negative reflection (blue). Seafloor depths in the survey area range from 1335-1664 ms TWT and, assuming a water velocity of 1500 m s -1 , between 1001-1248 m, respectively.
Borehole data from wells near the studied sills were used to: i) constrain the ages of mapped 10 key horizons and their lithology; ii) perform a tuning wedge model; and iii) perform a time-depth conversion based on check-shot data, thereby allowing us to compare element geometries measured in seismic reflection data with field analogues and the predictions of both numerical and laboratory analogue experiments. We used a second order polynomial best-fit regression to the check-shot data to extrapolate the measurements to greater depths below the bottoms of the 15 boreholes (e.g., Reeves et al., 2018;Norcliffe et al., 2021). To most accurately constrain sill depths and thicknesses, we used different time-depth conversions for sills S1 and S2, using check-shot data from boreholes nearest to each sill (S1 = Chester1-ST1, Hijinx-1, Nimblefoot-1, Makybe-Diva-1, Lightfinger-1; S2 = Briseis-1, Glencoe-1, Toporoa-1; Table S2-1 and Jupyter notebook S2-1 in the Supplemental Material 1 ). The coefficient of determination (R 2 = 0.99) of the estimated 20 gradient between check-shots and depth measurements was calculated to quantify the accuracy of the simplified time-depth conversion, which we then used to determine the paleo-emplacement depth and vertical offsets of the sill elements in meters.
In addition to converting measurements in time to depth in meters, we use seismic velocity information to establish the resolution to the imaged sills. However, major sills (e.g., S1 and S2) 25 resolved in this survey mainly occur at depths of 3200-5150 ms TWT (~3000-6600 m) and are not penetrated by boreholes, meaning we have no direct constraints on their seismic velocities.
Strata-concordant, high-amplitude reflections intersected by Rimfire-1 at 3632 m depths do correspond to a ~20 m thick basaltic rock body, which we interpret as a sill. Based on reported interval velocities from mafic sills in other case studies (e.g., Skogly, 1998;Berndt et al., 2000;30 Smallwood and Maresh, 2002;Magee et al., 2019a), we assume a sill velocity of 5550 (± 10%) m s -1 to determine the seismic resolution in the vicinity of S1 and S2 (e.g., Skogly, 1998;Berndt et al., 2000;Smallwood and Maresh, 2002). The dominant seismic frequency (f) in the proximity of the sills is ~25 Hz, which in combination with the inferred sill velocities (v), suggests a dominant wavelength (λ) of ~224.5 m (λ = v / f ). Following Brown (2011), this leads to a limit of separability (λ/4) of ~56 m (±5.6) and a limit of visibility (λ/30) of ~7 m (±0.7); i.e., sills with a thickness >56 5 m (±5.6) are resolved with distinct top-and base-contact reflections, whereas sills that are <56 m (±5.6) and >7 m (±0.7) thick (limit of separability > intrusion thickness > limit of visibility) are not resolved but detected in the seismic data and expressed as tuned reflection packages ( Fig. 3) (e.g., Brown, 2011;Magee et al., 2015). In this case, discrete seismic reflections from the top and base intrusion contacts interfere on their return to the surface and cannot be distinguished, thus 10 making it difficult to quantify sill thicknesses (e.g., Jackson et al., 2013). Apparent thickness measurements of tuned reflection packages are slightly underestimated for thicknesses above maximum tuning, and overestimated for thicknesses below maximum tuning (Fig. 3). In both cases, we cannot identify a distinct Base-sill reflection and, instead, the lowermost, high-amplitude reflection with a negative polarity was mapped as the Base-sill reflector. This reflection may not 15 represent the true sill-host rock contact, resulting in both under-and overestimated apparent thicknesses (Fig. 3). When a seismic signal propagates through the Earth's crust, its dominant frequency decreases with increasing depth, resulting in an increasing dominant wavelength and thus in a decay of vertical resolution in depth (e.g., Bacon et al., 2007). The horizontal seismic resolution is approximated by λ/4 resulting in ~56 m (λ≈224.5 m) within the vicinity of the studied 20 sills (e.g., Lebedeva-Ivanova et al., 2018). shows potential sill pick of both the Top and Base-sill reflector, following the uppermost/lowermost, high-amplitude reflection with a positive/negative polarity, respectively. (C) Predicted normalized amplitude and apparent thickness plotted against true thickness (Jupyter notebook S2-2 in the Supplemental Material 1 ). 10

Mapping Sills and Primary Magma Flow Indicators in 3D Seismic Reflection Data
With the exception of a seismically unresolved sill intersected by Rimfire-1, no other boreholes in 15 the study area penetrate sills. To identify sills within the seismic reflection data, we therefore use the criteria defined by Planke et al. (2005) and classify reflections as sills or inclined sheet intrusions where they: i) have high amplitudes; ii) display a positive polarity (a downward increase in acoustic impedance); iii) locally transgress host rock strata; and iv) terminate abruptly.
However, igneous sills are often displayed as tuned reflection packages, which adds uncertainties to the interpretation of sill geometries and morphologies as geophysical artifacts are more likely to occur (e.g., Smallwood and Maresh, 2002;Magee et al., 2015). Tuned reflections allow sill 5 outlines to be relatively well constrained in map view, but only allow estimates of sill thicknesses to be between the limits of separability and visibility. Tuned reflections therefore add uncertainties to 3D geometry interpretations including those of elements and element connectors (e.g., Magee et al., 2015;Eide et al., 2018), hampering the ability to interpret magma emplacement mechanisms.
Synthetic seismic forward models have shown that image quality can also be influenced by the 10 host rock style (i.e., homogenous vs. interbedded) (Magee et al., 2015;Eide et al., 2018).
Reflections from the contacts of strata-discordant inclined sheets and those emanating from the bedding planes that cross-cut can interfere with each other on their return to the surface, resulting in geophysical artifacts expressed as an apparent step-like intrusion geometry (i.e., 'pseudosteps') and therefore may lead to misinterpretation of intrusion geometries (Magee et al., 2015;Eide et 15 al., 2018).
Elements and their connectors form parallel to the main magma flow direction and are primary magma flow indicators (e.g., Rickwood, 1990;Schofield et al., 2012b;Galland et al., 2019;Magee et al., 2019b). Identifying and mapping elements and H-and S-/Z-connectors in seismic data thus helps to: i) reconstruct magma flow pathways and potential sill emplacement 20 history; and ii) identify potential sill feeder zones, as well as larger, source magma reservoirs (e.g., Thomson and Hutton, 2004;Schofield et al., 2012bSchofield et al., , 2017Magee et al., 2014). Connectors between adjacent elements are typically expressed as low-amplitude stripes on top-sill amplitude maps (e.g., Schofield et al., 2012bSchofield et al., , 2017Magee et al., 2015), resulting from sill thickness variations (Magee et al., 2015). For example, where elements propagate along the same 25 stratigraphic level and coalesce, the sill is thinner at the point of coalescence (e.g., Fig. 4D) (e.g., Pollard et al., 1975). In case the sill thickness is below the maximum tuning thickness, amplitudes decrease with decreasing sill thickness and can therefore be used to infer relative thickness variation and to identify contact geometries. On the other hand, where elements are vertically offset and connected via H-or S-/Z-connectors, local sill thickening occurs due to the formation of sub-30 vertical connectors (Fig. 4D). This sill thickening causes a decrease in amplitude in case the sill thickness is above the maximum tuning thickness ( Fig. 3 and 4D) (Magee et al., 2015). These findings are based on synthetic seismic forward models (Magee et al., 2015) and suggest that lowamplitude stripes on Top-sill amplitude maps can indicate connector geometries (Figs. 3 and 4D).
In this study, we manually picked two igneous sills (S1 and S2; Fig. 2B), and key horizons in encasing host rock, on every 10 th inline and crossline resulting in a regular grid spacing of 250 5 m. These grids were carefully quality checked and used as seed points to propagate (i.e., autotrack) surfaces. We further identified and mapped the primary magma flow indicators of sill S1 and S2.

Quantitative Analyses 10
We used a combination of seismic horizons of sills (map-view) and related seismic sections (crosssection view) to map sill morphologies, including different orders of sill elements and their connectors. For each sill and its component elements, we manually measured maximum resolved lengths and their maximum resolved widths perpendicular to the length in map-view (Fig. 4A).
Errors in length and width measurements are estimated to be in range of the seismic horizontal 15 resolution (~56 m). It is important to note that sills and their elements could thin below the limit of visibility (~7 m (±0.7)); all values reported below should therefore be considered minima.
We classify the geometry of elements based on changes in the α angle, which spans between the left-and right-hand element margin looking along a reference line oriented parallel to either: i) the long axis of elongate elements with width-to-length aspect ratios < 1; or ii) the long 20 axis of a feeder conduit for elements with width-to-length aspect ratios > 1 (Fig. 4B). The angle α is not constant for a single element, but rather increases and decreases along the defined reference line, with α > 0 indicating widening (i.e., increase in element width) and α < 0 indicating narrowing (i.e., decrease in element width) (Fig. 4B). To ensure accuracy as well as to assess how α varies along the along the reference line, the element margins are partitioned into smaller sections (Fig. 25 4B). The angle between the related section and the defined reference line is calculated at each node. These angles are measured for the left-hand side (αL) and right-hand side (αR), respectively, and the sum of both angles at the same location along the reference line equals the total angle α at this location (Fig. 4B). Since the analyzed margin sections have irregular lengths (Fig. 4B), we multiply α for each section by a pre-factor (0-1) that represents the section's proportion of the total 30 element length. This results in an α angle that is weighted by the section length, which prevents the overinterpretation of small sections with relatively large or small α angles. The mean angles reported in this study are equal to the sum of all weighted α angles in each element (Jupyter notebook S2-4 in the Supplemental Material 1 ).
Where the sill S2 is expressed as two separate reflections, mapping the top and base reflector allows us to depth convert measured sill thicknesses in time to meters using the inferred 5 velocity (5550 m s -1 (± 10 %)). Thickness and amplitude measurements were exported every 25 m along three seismic sections oriented either perpendicular or parallel to the long-axis of the elements. Where sill S2 is imaged as tuned reflections and absolute thicknesses cannot be quantified, amplitude measurements were used to infer relative thickness changes as suggested by synthetic forward models (Figs. 3 and 4D) (e.g., Magee et al., 2015). However, assuming 10 thicknesses of tuned reflection packages at the outermost sill margin are below maximum tuning, increasing and decreasing amplitudes indicate relative sill thickening and thinning, respectively To quantify vertical offsets between adjacent elements, and the vertical height of their connectors, we measured the TWT at zero-crossings (i.e., null point of the amplitude) of lateral 15 element tips, which we then converted to depth (s TWT to m) (Fig. 4C, D). This data was collected in ~100 m spaced seismic sections, starting from the most inward part of the connector towards its tip, perpendicular to the strike of the mapped contact geometries (H-and S-Z-connector). The element outline is partitioned into small sections to quantify how α varies in long-axis direction of (i) the element itself where width-to-length aspect ratios < 1, and (ii) the feeder conduit where 5 element width-to-length aspect ratios > 1, also referred to as 'reference line'. Angles are measured between the section and the element long-axis at the node of each section. The sum of angles measured for the left-hand side (αL) and right-hand side (αR) at the same location on the reference line equals the total angle α. A positive angle α represents an element widening (i.e., increase in width), whereas a negative angle α represents an element narrowing (i.e., decrease in width). (C) 10 Schematic diagram showing two vertically offset elements that are connected via an S-Z-connector (modified after Schofield et al., 2012a). Gray-scale color bar highlights the variation of vertical offsets as well as the location of minimum and maximum offsets. Dotted lines (i) and (ii) indicate the location of schematic cross sections shown in (D). (D) Schematic cross section of (i) coalesced elements at the same stratigraphic level, and (ii) vertically offset elements that coalesced via an S-5 Z-connector, resulting in a strata-concordant and step-like sheet intrusion, respectively. Red and blue colors of schematic seismic responses indicate a positive and negative polarity, respectively, and zero-crossings represent the null point of the amplitude. Vertical offsets are defined as the height between the zero-crossings of adjacent elements. Schematic seismic amplitude plots are based on predictions of seismic forward models (Magee et al., 2015). TRP = Tuned reflection 10 package.

RESULTS
In order to investigate sill-morphologies and small-scale sill-geometries, we mapped two seismically resolved sills (S1 and S2), where the Top-sill contacts are expressed as laterally discontinuous high-amplitude reflectors with a positive polarity, highlighting the downward increase in acoustic impedance at host rock-sill contacts (e.g., Fig. 5). Base-sill contacts are 20 resolved as high-amplitude reflectors with a negative polarity (i.e., decrease in acoustic impedance at sill-host rock contacts). Where possible, we identified and mapped the Base-sill contacts, where the reflections are resolved as distinct seismic events. However, in large areas of sill S1 and parts of sill S2 distinct Base-sill reflections are not resolved, in which case the sills are imaged as tuned reflection packages. The mapped sills have different morphologies, including: i) strata-discordant 25 inclined sheets (S1); and ii) saucer-shaped sills, where the strata-concordant inner sill is flanked by strata-discordant transgressive limbs (S2) (Fig. 2B).

Sill S1
The S1 inclined sheet dips southeast and extends at least ~35 km down dip and ~25 km along strike 30 (Fig. 5). However, the complete geometry and size of S1 cannot be determined because the southern end is not imaged within the seismic survey. S1 transgresses from a maximum depth of 5.15 s TWT in the southeast up to 3.25 s TWT at the northern tip, which is located immediately below the Mungaroo Formation-Rhaetian Marl contact (Figs. 5A and D). S1 is expressed as distinct top and base reflections within its northern sector, where it is up to ~85 ms TWT (~236 m) thick (Figs. 5C-E). In contrast, the sill peripheries, as well as reflections at greater depths of > 4.2 s TWT, are below the limit of separability and thus imaged as tuned reflection packages (~7 m < sill thickness < ~56 m) (Figs. 5C-E and 6A-C). The Top-sill reflection is of relatively highamplitude at shallow levels (< 4 s TWT), decreasing with depth (Figs. 5B and D). It is important 5 to note that the general relationship between relatively high-amplitude reflections at shallow levels holds irrespective of whether the sill is characterized by distinct or tuned reflections (Figs. 5B and E). Top-sill horizon of S1, as well as the two-way time thickness (C) between the Top and Base-sill horizons. The Base-sill was mapped in areas where two distinct reflectors are resolved (black dashed line), whereas the lowermost, high-amplitude reflection with a negative polarity was mapped in tuned reflection packages. (D, E) Interpreted seismic sections detailing both sill 5 geometry and morphology as well as sill emplacement related host rock deformation. Seismic section (D) is oriented parallel to the dip-direction of S1 showing strata-concordant and stratadiscordant sill sections. Multiple orders of elements are shown in section (E), which is oriented NE-SW and covers an area of shallow-level sill emplacement. This section includes elements resolved as tuned reflections in the SW, and as distinct reflectors in the NE, both observed at a 10 depth of ~3.2-3.8 s TWT. 15

Elements and their Connectors
We classify large sections of S1 as large-scale, first-order elements, which are resolved as tuned reflection packages (Figs. 5B and 6A). Absolute peak amplitude maps highlight linear variation of relatively lower amplitudes at the lateral borders of each seismically resolved element, indicating contact geometries (H-and S-Z-connectors; Figs. 5E and 6C). At relatively shallow emplacement 20 depth (~3.8-3.5 s TWT depth), we observe four orders of sub-horizontal, bedding-parallel elements radiating from the inclined sheet, forming an arcuate element-network (Fig. 6A). This element-network has a westward termination and extends over ~8.5 km N-S and ~5.5 km E-W ( Fig. 6A). Elongated, fourth-order elements form the smallest seismically resolved member of this network and they emerge from the outer margin of third-order elements (Fig. 6B). Since elements 25 are vertically offset, adjacent elements either over-or underlap and often coalesce, forming H-and S-Z-connectors, respectively, with an inconsistent stepping direction (Fig. 6C). In order to measure their orientations and vertical heights, we mapped 72 connectors within the element network described above (Fig. 6A). Element-connectors have semi-radial orientations with strikes varying from SW to NW (Fig. 6D). Where adjacent elements are separated by a pre-5 existing fault plane, S-or Z-connectors coincide with the fault plane and strike NE-SW similar to the pre-existing structure (Fig. 6D). The maximum connector length ranges from 81-1808 m, with fault-related connectors tending to be longer (Fig. 6D). Where connectors coincide with pre-  of connectors that coincide with a pre-existing fault plane are highlighted in red (Table S2-3 and Jupyter notebook S2-3 in the Supplemental Material 1 ).

Sill S2
Sill S2 occurs at a depth of ~3.8-3.5 s TWT and extends at least 13.5 km N-S and 6 km E-W direction (Fig. 8A). Pre-existing, Jurassic faults define the sill borders, resulting in an overall NNEelongate geometry (Fig. 8B). Inwardly inclined, moderately dipping (~23-35º) transgressive limbs 15 are coincident with a major fault bounding the sill's eastern margin, and a more minor fault in the west. Both inclined limbs extend over vertical distances of up to 0.29 s TWT (~470 m), transgressing the Mungaroo Formation-Rhaetian Marl contact (Fig. 2C). Even though the bounding fault west of S2 is only resolved at the southern end of the sill, an inclined limb formed over a lateral distance of ~7.2 km towards the north-northeast (Figs. 8A and B). Where S2 is not bounded by an inclined limb in the west, two westward-diverging, sub-horizontal element-5 networks (EN1, EN2) form an arcuate geometry with a northwesterly termination (Figs. 8B and D). Sill S2 comprises two morphologies separated by a NE-SW striking normal fault: i) an elongate saucer-shaped sill in the south; and ii) a sub-horizontal, bedding-concordant sill in the north (Fig.   8B). Two distinct reflectors are observed along most of the elongate inner sill of S2, defining a maximum thickness of 98 ms TWT (~272 m (±27)). In contrast, the inclined limbs and the sub-10 horizontal outer margin in the south and the northwest (EN2) are imaged as tuned reflection packages (Figs. 8B and C). The elongate section of S2 thins from east to west (Fig. 8C). The TWTthickness-map (Fig. 8C) highlights these thickness variations but also shows that sill thickness can vary distinctively between neighboring elements at a range of scales (Figs. 8C and 9A).
Three different orders of elements were observed and mapped throughout S2. First-order 15 elements include the semi-radially fanning element-networks EN1 and EN2 in the northern section of S2, which we sub-divide into second-order elements (Fig. 8D). Elongate, third-order elements form the smallest seismically resolved member in the element-hierarchy of S2 and emerge from the outer margin of second-order elements of EN2 but also occur locally within EN1 (Fig. 8D).
Third-order elements are therefore comparable to the fourth-order elements described from S1 20 such that both represent the smallest seismically resolved order of elements (Figs. 6B and 8D).
Second-order elements, with occasional smaller, third-order elements, occur within the elongated, saucer-shaped section and the southern tip of S2.
Relatively low amplitude linear features define element connectors that have different orientations depending on the sill morphology (Fig. 8E). Where inclined limbs bound S2, element 25 connectors strike parallel to the NE-SW trending bounding faults and inclined limbs (Fig. 8E).
Connectors that coincide with pre-existing fault planes strike in the same NNE-SSW direction as the faults, whereas connectors that do not coincide with fault planes have strikes that vary from NNW-SSE to NE-SW (Fig. 8E). Element connectors within the two sub-horizontal, beddingconcordant element-networks (EN1, EN2) in the north have strikes that fan towards NW, ranging 30 from N to W (Fig. 8E).

Figure 8: (A-C) Maps showing the time-structure (A) and absolute-peak-amplitude (B) of the
Top-sill horizon of S2, as well as the two-way time thickness (C) between the Top-and Base-sill horizons. The Base-sill was mapped in areas where two distinct reflectors are resolved (black dashed line), whereas the lowermost, high-amplitude reflection with a negative polarity was mapped in tuned reflection packages. (D) Absolute-peak-amplitude map detailing two semi-5 radially fanning element networks (EN1, EN2) north of the potential feeder zone. Seismic sections as well as thickness and amplitude profiles along X-X', Y-Y', and Z-Z' are shown in Fig. 9. Rose diagrams (E) show the strike of mapped connectors between adjacent elements. Bins are colorcoded with respect to the connector's length. 10

Amplitude and Thickness Variations of Elongate Elements
Amplitudes of the Top-sill horizon and thickness measurements of S2 were collected along two 15 seismic lines oriented perpendicular to, and one parallel to the strike of the element connectors (location of lines shown in Fig. 8D). Profile X-X' (Fig. 9A) shows amplitude and thickness variations of elements in cross-section in the more inward part of EN1. In this area, the top and base of EN1 are resolved as two distinct seismic reflections. High amplitudes are mainly located in the center of elements, with minor fluctuations occurring along these high-amplitude plateaus 20 ( Fig. 9A; E1 and E4). Low amplitudes define segment connectors (i.e., S-Z-connectors) (Fig. 9A).
Sill thickness measurements based on depth-converted Top-and Base-sill horizons range from 123.5(±12.5) -185.5(±18.5) meters, with adjacent elements varying in thickness (Figs. 8C and 9A). Where single elements are resolved without overlapping neighboring elements, the peak thickness occurs in the center of the elements and thins out towards their margins ( Fig. 9A; E4, 25 E5). In some cases, peak thicknesses coincide with low amplitudes and the inferred location of element connectors ( Fig. 9A; contact between E2 and E3). Here, the top reflection of the stratigraphically higher element overlies the bottom reflection of the stratigraphically lower element, resulting in a measurement of the S-Z-connector thickness, rather than the single elements, and an apparent local sill thickening ( Fig. 9A; contacts between E2-E3 and E4-E5). 30 Thickness-to-width aspect ratios for the elements displayed in profile X-X' vary between 0.31 and 0.71.
Amplitudes periodically increase and decrease in a second-order element at the outer margin of EN2 (Y-Y'), where the sill is expressed as a tuned reflection package (Fig. 9B).
Relatively low amplitudes define element connectors (i.e., S-Z-connectors) (Fig. 9B). Low amplitudes within the second-order element are expressed as linear low-amplitude features in mapview and potentially highlight the contacts of strata-concordant, third-order elements, which within 5 their centers show relatively high amplitudes (Figs. 8D and 9B). Since EN2 is resolved as a tuned reflection package, thickness measurements do not represent the true sill thickness, and absolute thickness variations cannot be inferred from the data. Assuming a limit of visibility and separability of 7(±0.7) m and 56(±5.6) m, respectively, thickness-to-width aspect ratios of the elements in EN2 range from 0.02 to 0.35, which is smaller than those observed in EN1 (profile X-10 X').
Along seismic section Z-Z' parallel to the strike of element connectors, we show that thickness and seismic amplitude varies in both cases where the sill is resolved as separate reflectors (EN1) and tuned reflections (EN2; Fig. 9C). EN1 gradually thins towards its tip (132(±13) -99(±10) m), while the thickness drops below the limit of separability in EN2 (Fig. 9C). Minor 15 amplitude variations occur within EN1, and low amplitude values define the contact between EN1 and EN2. Amplitudes within EN2 fluctuate slightly, with a peak occurring at 2750 m along the profile, which decreases rapidly towards the sill tip (Fig. 9C).
A plot of normalized amplitudes versus the depth-converted thicknesses of all measurement points of the three seismic sections (Figs. 9A-C) indicates a high variance of 20 amplitudes for thicknesses greater than the limit of separability (Fig. 9D). With decreasing sill thickness, amplitudes increase below the limit of separability and are therefore expressed as tuned reflection packages (Fig. 9D). Amplitudes within these tuned reflection packages further increase until reaching the maximum tuning thickness, after which they decrease (Fig. 9D). The minimum apparent thickness observed is 38 m, which is only slightly below the maximum tuning thickness 25 of 51 m. However, we note that these apparent thickness measurements are likely to be overestimated (Fig. 3C). Predictions of a simplified tuning wedge model (Jupyter notebook S2-2 in the Supplemental Material 1 ) show similar amplitude trends, but the predicted maximum tuning thickness of 43 m is 8 m lower than observed in our data. Modeled amplitudes for thicknesses above the limit of separability are constant, contrary to data collected along the seismic sections 30 described above (Fig. 9D). host rock parameters are based on borehole data (Hijinx-1). Note that thicknesses below the limit of separability represent apparent thicknesses and are likely to be underestimated or overestimated (Fig. 3).

Element Geometries
We observe three types of element geometry in S1 and S2, including elongate elements with parallel lateral edges (Type 1), elongate and widening elements (Type 2), and fanning elements (Type 3) (Fig. 10A). Below, we describe and quantify the geometry of Type 2 and Type 3 to better constrain how elements widen and/or narrow (Figs. 4B and 10). 5 Type 2 -elongate, slightly to moderately widening elements. We subdivide these into Type 2A elements that steadily increase in width (α > 0); and Type 2B elements that first increase in width (α > 0), but then narrow along the long axis towards the tip (α < 0) (Fig. 4B, 10A). Even though the overall width of each Type 2 element consistently increases or decreases, local variation in α of up to 90º can occur on a small scale (Fig. 10A). We observed 38 Type 2 elements with 10 broad range of αmean values (-1-55º), where Type 2A elements (n = 29) tend to have larger αmean angles with a broader range (7-55º) compared to Type 2B elements (n = 9; αmean = -1-27º) (Fig.   10B). This occurs because Type 2B elements narrow at their tips with α < 0, resulting in smaller αmean angles (Figs. 3B and 10A). Since Type 2 elements have elongate geometries, their width-tolength aspect ratios are < 1 (0.27-0.93; average of 0.61), with element lengths ranging over two 15 orders-of-magnitude (181-5061 m) (Fig. 10C).
These elements are often fed via a relatively short and narrow conduit with a small to moderate α angle < ~50º, before a rapid increase in width with α up to 137º (Fig. 10A). We observed 14 Type 3 elements with a broad range of αmean angles (41-122º) (Fig. 10B). Their maximum width and 20 length range from 370-3940 m and 345-2348 m, respectively, resulting in width-to-length aspect ratios predominantly > 1 (0.77-2.37; average of 1.42), with aspect ratios < 1 (0.77 and 0.78) occurring in two elements with relatively long feeding conduits. Width-to-length aspect ratios of Type 3 geometries tend to increase with increasing element length (Fig. 10C).
Sills mapped in this study show multiple orders of elements where higher-order elements 25 contain at least one lower order of sub-elements. Long-axis orientations of sub-elements within Type 2 and Type 3 elements are predominantly sub-parallel to the margins of the higher-order element (Fig. 10A). Therefore, these long-axis orientations vary in the range of α, resulting in an approximately linear and a semi-radially fanning trend within high-order Type 2 and Type 3 elements, respectively (Figs. 10A and D).  Fig. 4B). (C) Cross-plot of 5 the maximum element length and width. Black trend lines showing the ratio between the width and the length (i.e., w/l = 0.25, w/l = 0.5, w/l = 1, w/l = 2), with aspect ratios < 1 indicating elongate geometries, and aspect ratios > 1 indicating wider geometries, highlighted by cartoons (solid black lines). (D) Cross-plot of width-to-length aspect ratios (C) and αmean angles. Plotted data shows a correlation to a second order polynomial fit (dashed line). Measurement errors in width and length are associated with the horizontal seismic resolution (i.e., 56 m). Error in width-to-length ratios indicate the range of ratios based on minimum and maximum width and length measurements. Since we cannot quantify the error of angle αmean measurements, the error is given as 5% of αmean. Different element geometries in (B-D) are color-coded (Table S2-4-1, S2-4-2, S2-4-3 and Jupyter notebook S2-4 in the Supplemental Material 1 ). 5

DISCUSSION 10
Because comprehensive analyses of the full 3D geometries of sub-horizontal sheet intrusions are limited to: i) geophysical techniques such as 3D seismic reflection data; and ii) numerical and laboratory analogue modelling, it is important to link observations of both approaches together with field observations to unravel magma emplacement processes. Our results reveal distinct element geometries, as well as thickness variations within different orders of elements. Below, we 15 discuss the implications of our results by considering how sill geometries and 3D seismic reflection data can help to infer mechanisms that may have led to sill branching. We then introduce a conceptual flow and emplacement model for the sills described in this study, based on sill and element geometries.

Can Seismic Reflection Data Be Used to Infer Sill Branching Mechanisms?
Field observations of host rock deformation adjacent to sills indicate the operation of: i) viscous emplacement mechanisms such as host rock fluidization (e.g., Kokelaar, 1982;Schofield et al., 2010); ii) elastic bending and uplift of the overburden strata potentially due to LEFM (e.g., Johnson and Pollard, 1973;Walker, 2016); and iii) viscous indentation and plastic shear failure (e.g., 25 Spacapan et al., 2017;Galland et al., 2019). All three deformation mechanisms described above may explain the breakdown of planar igneous sheets into elongate elements (e.g., Schofield et al., 2012a;Magee et al., 2019b and references therein). A key problem with seismic reflection data is that host rock deformation structures typically used to infer sill emplacement mechanisms are rarely imaged in sufficient detail to interpret their origin. Furthermore, boreholes rarely penetrate 30 sills and, where they do, provide only 1D information on rock properties. To interpret sill emplacement mechanics from seismic reflection data we therefore currently rely on measurement of overburden uplift (forced folds) above sills, generated by reverse faulting (e.g., Thomson and Schofield, 2008) and/or elastic bending (e.g., Hansen and Cartwright, 2006b;Galland, 2012;Jackson et al., 2013;Schmiedel et al., 2017;Reeves et al., 2018;Magee et al., 2019a;Norcliffe et al., 2021). For example, Jackson et al. (2013) suggested that discrepancies between sill thicknesses and forced fold amplitudes may be caused by compaction and/or host rock volume reduction due 5 to pore fluid expulsion, potentially indicating host rock fluidization and ductile flow of host rock.
In previous seismic reflection-based studies (e.g., Thomson and Hutton, 2004), contact geometries (i.e., H-and S-Z-connectors) between elements have been interpreted to form due to tensile brittle fracturing between two inflating, adjacent but vertically offset elements, suggesting magma emplacement through elastic-brittle fractures (e.g., Francis, 1982;Schofield et al., 2012a). 10 However, S-Z connectors have also been observed between vertically offset elements where host rock deformation suggests that magma emplacement processes were not dominated by elasticbrittle fracturing . For example, in the Neuquén Basin, Argentina, two slightly vertically offset magma fingers are connected via an S-Z connector and the dominant emplacement mechanism is interpreted to be viscous indentation (Fig. 1A, Finger 1 and 2)  . 15 Vertically offset elements therefore cannot be used to identify fracture segmentation processes without additional information about host rock deformation. This finding is especially important for seismic reflection data interpretation, where data and/or observations of host rock deformation is not available and highlights the importance of using the descriptive term element where mechanisms that accommodate both sill emplacement and branching cannot be identified. 20

Vertical Connector Heights Reveal Fracture Segmentation
Based on our observations of sills in the Exmouth Plateau, we suggest that detailed quantitative measurements of connectors can provide insights into magma emplacement mechanisms. For example, our data indicates that significant variations in vertical element offsets 25 and therefore in connector heights can occur along the length of elements (Fig. 7). Although recognition of vertical offsets cannot be used to determine whether emplacement involved elasticbrittle fracturing or brittle and/or non-brittle deformation mechanisms (i.e., viscous indentation, host rock fluidization), we can test whether sheet-segmentation was related to stress-field rotation.
In particular, if we assume host rock properties remain constant, we may expect that when 30 segmentation is due to stress-field rotation, the vertical offset of elements should increase continuously along the connector lengths, and the elements should step in a consistent direction ( Fig. 1D) (e.g., Pollard et al., 1975Pollard et al., , 1982Rickwood, 1990;Schofield et al., 2012a;Magee et al., 2019b). Conversely, when sill segmentation is caused by the exploitation of preferentially oriented, pre-existing structures such as bedding planes or faults, magma migrates along surfaces with relatively lower tensile strength and/or fracture toughness, which may result in an inconsistent 5 stepping direction, with vertical offsets varying irregularly along the connector length (e.g., Schofield et al., 2012a;Stephens et al., 2017;Magee et al., 2019b). Quantifying connector geometries and their vertical heights could therefore help to identify stress-field rotation as a fracture segmentation mechanism. Maximum vertical connector heights in S1 are evenly distributed along element lengths (Fig. 7A), indicating that syn-emplacement rotation of principal 10 stress axes orientations was not the dominant mechanism that led to sill branching of S1. This conclusion is supported by the inconsistent stepping direction observed within the S1 element network, indicating that exploitation of pre-existing weaknesses led to vertically offset elements (Figs. 5E and 6C) (Schofield et al., 2012a;Magee et al., 2019b).

Amplitude Variations May Imply Coalesced Magma Fingers
Since both elastic-brittle fracturing but also brittle and/or non-brittle emplacement mechanisms can lead to the exploitation of pre-existing weaknesses, we used seismic attributes to infer host rock properties to better assess the formation mechanisms of sill branches. For example, seismic reflections of the Top-sill horizon of S1 show significantly increased amplitudes at 20 shallower levels < 3.8 s TWT (Figs. 5A and B). Seismic reflections with relatively higher amplitudes are observed for both distinct and tuned reflections at similar stratigraphic levels (Figs. 5B and E), suggesting a relative change from a mechanically strong towards a weaker host rock.
Lithological changes plus increasing cementation and compaction with depth can increase the mechanical host rock strength (e.g., Kokelaar, 1982;Schofield et al., 2010Schofield et al., , 2012a. Such burial 25 effects are associated with a decrease in pore fluid volume which reduces the potential for host rock fluidization and thus non-brittle magma emplacement (e.g., Kokelaar, 1982;Schofield et al., 2010Schofield et al., , 2012a. We hypothesize that the relative change towards mechanically weaker host rocks at shallower emplacement depths of approximately 950-1150 m in our data resulted in non-brittle emplacement of S1 and, consequently, the formation of magma fingers. Our data further indicates that detailed analyses of seismic amplitudes can be used to identify coalesced magma fingers that propagated within the same seismically resolved stratigraphic interval (S2, Fig. 9B). Where the sill thickness is below the maximum tuning thickness, increasing and decreasing amplitudes indicate relative sill thickening and thinning, respectively, although the thickness cannot be measured directly. Since magma fingers are slightly 5 convex-upwards, they form cusp-shaped grooves at their point of coalescence, resulting in an undulating Top-sill geometry, where the sill is thinner at magma finger contacts (Fig. 1A, 4D).
Contact geometries of coalesced magma fingers are therefore expressed as linear features with relatively low amplitudes, oriented in the magma finger long-axis direction (Fig. 8D). These findings are consistent with predictions of synthetic seismic forward models (Magee et al., 2015). 10 Due to the limit of seismic resolution, small-scale elements may not be resolved and will remain undetected. We therefore note that the smallest-resolved element order in seismic reflection data may contain lower orders of unresolved sub-elements.

What Do Element Geometries Tell Us About Sheet Emplacement? 15
Below, we interpret element geometries described in this study and we discuss what these geometries reveal about the emplacement of igneous sheet intrusions. We then present a potential growth model for elements based on their geometries, which is discussed and compared with observations presented in existing field and seismic studies.

Variations Indicate Localized Magma Flow
Both magma emplacement and subsequent intrusion growth are often described as incremental processes, with numerous magma injections occurring over periods of months up to millions of years (e.g., Hurlbut Jr, 1939;Biggs et al., 2009;Annen, 2011;Annen et al., 2015;Magee et al., 2016a;Cruden et al., 2017;Cruden and Weinberg, 2018;Reeves et al., 2018). Distinct thickness 25 variations observed in our data suggest that S2 may be composed of three separate consecutive magma pulses, where the first pulse emplaced elements close to the feeder location (Fig. 8, 13).
When the second and third pulses of magma were injected it fed the propagating sill tip resulting in tip propagation, and led to sill inflation in areas that were previously injected (Fig. 8, 13). Based on the data available, we can neither determine the exact duration of sill emplacement nor the time 30 lag between the individual magma pulses. However, since the second and third pulses may have contributed to both sill tip propagation and vertical inflation of the previously injected magma body, these consecutive pulses were likely emplaced prior to cooling and solidification. We therefore suggest that S2 was probably emplaced over a relatively short period of months to years.
Our data indicates that similar patterns of thickness variation occur at smaller scales for the different orders of elements (Figs. 8C and 9A, B). These patterns can be accentuated by flow 5 localization, a process by which rapid solidification of thinner parts of a sheet intrusion (e.g., point of coalescence of elements) results in areas of focused magma flow and inflation (Fig. 11C) (Holness and Humpherys, 2003). Such internal flow localization could also produce thickness variations within a planar element or sill (Holness and Humpherys, 2003). However, these variations could be smooth and expressed as large-scale, undulating intrusion morphologies. We 10 therefore suggest that abrupt variations in thickness, as observed here, could be accentuated by but not produced due to internal flow localization (Fig. 11C). This also implies that once elements form, flow can be channelized within them (Fig. 11C, D) (Magee et al., 2013c). Identifying this channelized magma flow in sheet intrusions has important implications for: i) better understanding of the architecture of sill complexes and how they are built (e.g., Thomson and Hutton, 2004;15 Cartwright and Hansen, 2006;Hansen and Cartwright, 2006a;Guo et al., 2013;Magee et al., 2016a;Schofield et al., 2017); ii) assessing areas of potential fissure eruptions emanating from areas with localized magma flow (e.g., Pansino et al., 2019); and iii) exploration for Ni-Cu-PGE deposits, since sulfide liquids are commonly accumulated in elongate parts of intrusions, which may trap sulfide liquids (e.g., Barnes et al., 2016). 20 . When one element is fully solidified, magma flow can localize within the adjacent element, resulting in inflation (C, iii). Similar flow localization can occur when the duration between two magma pulses is long enough to fully solidify magma that was injected during the first pulse (D). 10 In this scenario, the second magma pulse propagates along the solidified first magma pulse, resulting in a similar final geometry as described in (C). (E) Cross-sectional field photographs of coalesced elements: i) coplanar magma fingers intruding sandstone at the outer margin of the Shonkin Sag laccolith, Montana; and ii) vertically offset segments intruding metasedimentary rocks on Ardnamurchan, NW Scotland (modified after Magee et al., 2019b). It is important to note 15 that the final geometry of the element network shown in (A) and (B, ii) is similar and that magma fingers or segments cannot be distinguished in map view.

The Growth of Elements
When a single element emerges from a planar sheet it grows in length and thins towards the propagating tip (Figs. 9C and 12C) (Horsman et al., 2005). With time and sustained magma supply, elements grow in width and inflate (Fig. 12C). We observe smaller thickness-to-width aspect ratios 10 closer to the sill tip than toward the center of S2, which further suggests that elements first grow in width before they inflate. This interpretation is in agreement with field observations made at the outer margin of the Shonkin Sag laccolith, where adjacent magma fingers rapidly increase in thickness while approaching their point of coalescence (Pollard et al., 1975). When an element grows in length and width and then inflates, the propagating tip can become unstable, eventually 15 breaking down into multiple smaller elements (Fig. 12B) (e.g., Saffman and Taylor, 1958;Pollard et al., 1982;Takada, 1990;Schofield et al., 2012a). We note that this can happen in response to both elastic and viscous instabilities (e.g., Pollard et al., 1975;Nicholson and Pollard, 1985;Takada, 1990;Schofield et al., 2010). The former element now acts as a feeder conduit for its newly formed sub-elements (Fig. 12A, B). Higher-order elements tend to be thicker close to their 20 feeder conduit and thin out towards their tips, analogous to lobes in pahoehoe lavas (e.g., Peterson et al., 1994), which are thicker close to their feeding vents, or flow inflation lobes (Hansen and Cartwright, 2006a;Miles and Cartwright, 2010). However, flow inflation lobes are thought to form where sills intrude unconsolidated sediments at shallow emplacement depths (<500 m) and therefore this mechanism may not account for the elements described here, which were emplaced 25 into consolidated sedimentary rocks at ~950-1150 m depth. Our described growth model is compatible with other sills mapped in 3D seismic reflection data.
For example, data from the North Rockall Trough images a sill with 'flow units' of varying orders, referred to as 'fingers' and 'en-échelon fingers' (Thomson and Hutton, 2004). The emplacement model proposed by Thomson and Hutton (2004) describes first-order flow units branching into 30 second-order flow units, which again may bud to form a third order. This succession implies a time-dependency between the varying orders, such that a higher order flow unit is older than the lower orders it contains. We emphasize that our use of the term 'element' is purely geometrical, implying neither a specific magma emplacement mechanism nor the relative timing of emplacement. The latter is important, as different sill emplacement histories can result in the same final sill geometry, which is recorded in 3D seismic reflection data (Fig. 11A, B). For example, consider a scenario where a first order element contains second and third order elements, with the 5 latter representing the smallest seismically-resolved order (Fig. 11). This element network may have formed in response to either: i) a single magma pulse (Fig. 11A); or ii) distinct magma pulses separated in time (e.g., Biggs et al., 2009;Annen, 2011;Magee et al., 2014;Reeves et al., 2018) ( Fig. 11B). In the second scenario, third order elements may have formed during magma pulse 1 (Fig. 11B, i), whereas an additional second order element formed during magma pulse 2 (Fig. 11B,  10 ii), such that third order elements can be older than second order elements within the same element network. However, neither hypothesis can be tested using seismic reflection data alone.
Our observations reveal that thicknesses can vary significantly between neighboring elements across multiple orders (Figs. 8C and 9). We suggest that: i) the emplacement of elements 15 is geometrically self-similar over multiple scales whereby elements with different orders may grow independently, with potentially independent flow kinematics; and ii) neighboring elements either inflate and/or deflate at different rates or different times, which implies that sill growth is incremental down to the smallest seismically resolved member in the element-hierarchy (Fig. 11,   12). 20 starting to emerge at the element tips. Contours are based on a numerical study on viscous fingering (modified from Ebrahimi et al., 2016). (C) Cross-section through a feeder conduit and emerging elongate elements, parallel (X-X') and perpendicular (Y-Y') to the element long axis. Schematic diagrams are not to scale. 10

Sill Emplacement History and Magma Flow Pathways in the Studied Sills
Seismically resolved sills in the Exmouth Plateau are not penetrated by wells and radiogenic isotope ages are therefore not available. However, paleo-seafloor surfaces mapped using sillinduced forced folds indicate that magmatism occurred during the Upper Jurassic-to-Early Cretaceous, which agrees with previous studies (e.g., Mihut and Müller, 1998;5 Symonds et al., 1998;Magee et al., 2013a;Mark et al., 2019;Magee and Jackson, 2020;Norcliffe et al., 2021). Interpretation of sill-emplacement related forced folds and onlapping sedimentary strata suggest a paleo-emplacement depth of ~ 4.7-0.4 km for the inclined sheet S1 (this study), and ~ 0.9 km for sill S2 (Norcliffe et al., 2021).
Sills in sedimentary basins often form extensive networks (i.e., sill complexes) (e.g., 10 Cartwright and Hansen, 2006;Magee et al., 2016a). Sill junctions form where two sills connect, and they can be used to reconstruct magma flow pathways and the overall plumbing system architecture (Hansen et al., 2004). Identifying the feeders of sills that are not fed by other sills is challenging since they often have relatively thin and steeply dipping geometries (e.g., dykes), and as such can only be recognized in high-resolution seismic reflection data as localized areas of low-15 amplitude reflections disrupting otherwise continuous host rock-related reflections (Magee and Jackson, 2020). In the field, extensive studies of magma flow indicators (e.g., magnetic fabrics, shape preferred orientation of minerals, slickensides, shear sense indicators) (e.g., Cruden and Launeau, 1994;Cruden et al., 1999;Walker et al., 2017;Magee et al., 2018;Galland et al., 2019;Martin et al., 2019) can help to reconstruct magma flow pathways and potential feeder locations. 20 However, intrusions are rarely fully exposed, which makes sample collection and a consistent spatial distribution of data challenging. 3D seismic reflection data allow for mapping of sill geometries, including element connectors (H-and S-Z-connectors), in high detail. Since element connectors form along the lateral margins of elements, it is widely accepted that their orientation is a magma flow pathway indicator. We highlight that contact geometries form between adjacent, 25 vertically offset elements (i.e., segments, magma fingers), but also between elements that propagate within the same stratigraphic interval (i.e., magma fingers), which allows us to interpret magma flow directions regardless of sill emplacement mechanisms. Here, we use quantitative measurements of these connector geometries to map magma flow pathways, infer potential feeder locations, and to reconstruct a potential emplacement history for S1 and S2.

Sill S1
Element-element contacts and thus magma flow directions in S1 can be precisely mapped within a fanning element network at shallower levels, indicating an overall westerly flow direction (Fig.   6). At greater depths, most likely due to a decrease in seismic resolution, only large-scale magma flow indicators are recognized in the form of potential first-order elements (Fig. 5B). However, 5 based on the available data, we infer that S1 is an inclined sheet that propagated upwards, crosscutting the host-rock strata, with a feeder(s) located at the deepest levels in the SE (Fig. 5). Magma flowed towards the N to NNE until it reached an area where relatively higher amplitude reflections may indicate a change in host rock lithology. We infer that this potential change in mechanical properties, possibly defined by a change to a weaker, less brittle host rock, may have allowed for 10 non-brittle deformation processes and potentially the formation of bedding-parallel magma fingers that fanned-out westwards forming a large-scale element network (Figs. 5 and 6). Magma fingers propagated along pre-existing weaknesses (e.g., bedding planes), resulting in vertical offsets and a step-like sill geometry with an inconsistent stepping direction (Fig. 6C). NE-SW striking Jurassic faults influenced magma flow pathways such that magma propagated along pre-existing fault 15 planes, resulting in a localized fault-parallel magma flow. With increasing element inflation, neighboring but vertically offset elements coalesced via H-and S-Z-connectors due to tensile brittle failure of the host rock separator, forming a broadly continuous sheet containing abrupt intra-sheet steps.

20
Sill S2 Radially spreading magma flow pathways suggest that even though clear evidence for magma feeder(s) in cross-section is missing, S2 was fed from a local source along the NE-SW striking fault that appears to splay off the major fault that bounds the eastern margin of S2 (t0; Figs. 8B and 13B-C). When the first pulse of vertically ascending magma interacts with the pre-existing fault, 25 local rotation of σ3 from the sub-horizontal towards an orientation perpendicular to the fault surface allows magma to propagate along the fault (Fig. 13A, i) (Magee et al., 2013b). Two potential models can explain the subsequent transition from a steeply dipping dyke to a sub-horizontal sheet intrusion (t1; Fig. 13A, iii): i) once the ascending magma reaches a horizon (e.g., bedding) with a tensile strength and/or fracture toughness that is lower than that of the intruded fault plane, magma 30 will propagate along the interface of host rock layers resulting in the transition to a sub-horizontal sill (t1; Fig. 13A, ii-iv) (e.g., Kavanagh et al., 2006Kavanagh et al., , 2015; or ii) increasing dyke thickening results in compressional stresses with σ1 orthogonal to the feeder wall and a sub-vertical σ3 at the feeder tip (e.g., Anderson, 1951;Parsons and Thompson, 1991), resulting in the emplacement of a subhorizontal sill. The obtuse angle between the feeding fault plane and the footwall resulted in a rotation of σ3 to a sub-vertical orientation, which forced the magma to primarily intrude the 5 footwall in a SE direction (Figs. 13A-C). Inflation and associated tensile failure may have caused sill propagation into the hanging wall towards the north (Fig. 13A, iii-iv). Due to a potential second magma pulse, S2 inflated and propagated SW along the major bounding fault in the east, forming an elongate, channel-like sill (t2; Fig. 13B). This highlights that elongate sills may not always be fed via an underlying dyke striking parallel to the intrusion long-axis (e.g., Galerne et al., 2011;10 Fyfe et al., 2021), and that pre-existing fault planes can have major influence on sill geometries.
Inclined limbs formed along the pre-existing fault plane in the east but also in the west of the sill, where only minor faults in the SW are imaged in the seismic reflection data (t3-t5; Figs. 13C and D). We interpret the straight geometry of the western inclined limb to be a consequence of magma following fractures and/or faults (e.g., Thomson and Schofield, 2008;Magee et al., 15 2013b) that are not resolved in the seismic reflection data, rather than being due to an asymmetric stress field at the sill tip (e.g., Malthe-Sørenssen et al., 2004). Three plausible scenarios may have formed the western inclined limb of S2: i) tensile failure induced by forced folding, whereby sill inflation lifts the overburden and magma follows fractures that formed at the point of highest flexure (e.g., Thomson and Schofield, 2008); ii) magma propagation along reverse faults that 20 formed due to overburden uplift during sill inflation (e.g., Thomson and Schofield, 2008;Haug et al., 2017); or iii) exploitation of pre-existing faults (e.g., Magee et al., 2013b) (Fig. 13D). Faults are below the seismic resolution in the latter two scenarios. Since Jurassic age faults on the Exmouth Plateau strike NNE-SSW, similar to the western inclined limb of S2, and given that forced folds are observed above S2, we cannot differentiate between the three mechanisms above 25 based on the available data. Detailed forced fold and fault analyses by Norcliffe et al., (2021) showed that sill inflation of S2 was likely accommodated by a combination of overburden strata uplift and, adjacent to the pre-existing bounding fault in the east, fault inversion, and compaction of overburden strata.
Bedding-parallel, in some cases slightly vertically offset, elements north of the potential 30 feeder zone formed the first-order element EN1 (t2; Fig. 13B). When elements within this network grew in width and thickness, brittle tensile failure of the host rock separators resulted in the formation of S-Z-connectors (t3; Fig. 13), as commonly suggested by field studies (e.g., Nicholson and Pollard, 1985;Rickwood, 1990;Hutton, 2009). During a potential third magma pulse, similar processes formed the first-order element EN2 (t4, t5; Fig. 13). However, in contrast to EN1, individual elements within second-order elements propagated along the same strata (Fig. 9B). 5 Based on this finding and in combination with thickness-to-width aspect ratios of 0.02-0.35, we suggest that EN2 comprises second-order elements of coalesced magma fingers, where minor vertical offsets between adjacent second-order elements resulted from the exploitation of preexisting structures.
A large variation of amplitudes is observed for thicknesses above the limit of separability 10 ( Fig. 9D), which we interpret as due to variations in host rock lithology, or sill composition. Since the latter scenario would likely result in accumulation of distinct amplitude clusters for portions of the sill characterized by a distinct igneous rock composition, we hypothesize that the broad scatter of amplitudes for thicknesses above the limit of separability is more likely to indicate varying host rock lithologies. These changes in host rock lithology could have promoted the formation of 15 elements and may have caused minor vertical offsets, which we observe between the second-order elements of S2. It is to note that color-coded areas in (B) indicate the outer margin of propagating magma at given timesteps; i.e., even though younger magma pulses seem to be accreted to the edges of the intrusion, they also contributed to sill inflation as highlighted in cross section (C), (D). Cross section locations of (A), (C), and (D) are indicated in (B). Schematic 10 cross-sections are not to scale. Please refer to 'Sill Emplacement History and Magma Flow Pathways in the Studied Sills' for a detailed discussion.

CONCLUSION
High-quality 3D seismic reflection data from the Exmouth Plateau, offshore NW Australia was used to: i) map basaltic sills and their elements; ii) document their geometries and seismic expressions; and iii) introduce new descriptive terminology to define different sill elements and connectors based on their geometries. Our key conclusions are: 5 1. Five Late Jurassic-to-Early Cretaceous basaltic sills intruded interbedded sandstones and claystones of the Triassic Mungaroo Fm. and the Rhaetian Marl at emplacement depths of ~4.7-0.4 km. The sills are seismically resolved in the Glencoe 3D survey, and they have inclined sheet, saucer-shape, and sub-horizontal sill morphologies.
2. Sub-horizontal, semi-radially fanning element networks occur in two sills where: i) an inclined 10 sheet intrudes into a potentially mechanically weaker host rock (S1); or ii) a saucer-shaped sill is bounded by only one inclined limb (S2).
3. Quantifying vertical heights of H-and S-Z-connectors in 3D seismic reflection data provides important insights into sill branching mechanisms. This approach helps to identify/exclude synemplacement rotation of principal stress axes orientations ahead of a propagating fracture as a 15 potential segmentation mechanism. 4. Elongate amplitude variations within elements can indicate coalesced sub-elements even where they are resolved as tuned reflection packages and emplaced at approximately the same stratigraphic level. These findings are in agreement with results of synthetic forward models (Magee et al., 2015) and permit the analysis of magma flow pathways in 3D seismic reflection 20 data at meter scale, and help to infer potential sill emplacement mechanisms.
5. Based on quantitative measurements of element geometries, we classified elements into: i) elongate elements with parallel lateral edges (Type 1); ii) elongate elements that consistently widen and have width-to-length aspect ratios < 1 (Type 2); and iii) semi-radially fanning elements often fed via a relatively short and narrow conduit, with width-to-length aspect ratios predominantly > 25 1 (Type 3). Quantifying element geometries permits comparisons between different subsurface and field-based datasets covering a range of host rock types and tectonic settings and further allows us to test the predictions of laboratory and numerical analogue experiments.
6. Sill thicknesses vary significantly in lateral directions. We characterized sill S2 into three main zones based on their thicknesses, where the sill is thicker close to its feeder. This suggests that 30 incremental emplacement of potentially three magma pulses led to the final sill geometry. 7. Thickness variations between adjacent elements indicate that these features grow and inflate independently and therefore may have had separate flow kinematics resulting in magma flow localization.
8. Thickness-to-width aspect ratios of individual elements at the outer margin of S2 are smaller than those measured towards the center of S2, suggesting that elements first grow in width before 5 they inflate.

ACKNOWLEDGMENTS
We are grateful to Geoscience Australia for making the seismic reflection and borehole data 10 publicly available from https://nopims.dmp.wa.gov.au/nopims. We also thank Down Under GeoSolution and Schlumberger for use of DUG Insight 4 and Petrel seismic interpretation software. We gratefully acknowledge helpful reviews by Eric Horsman and an anonymous reviewer, and additional commentary from Associate Editor Robert Miller. We thank Shanaka de Silva for his editorial handling of the manuscript. JK acknowledges a Monash Graduate 15 Scholarship. ARC, CALJ and CM acknowledge funding from ARC Discovery Grant DP190102422. We thank Andy Bunger for fruitful discussions on element geometries.