The propagating margins of igneous sills (and other sheet intrusions) may divide into laterally and/or vertically separated sections, which later inflate and coalesce. These components elongate parallel to and thus record the magma flow direction, and they can form either due to fracture segmentation (i.e., “segments”) or brittle and/or non-brittle deformation of the host rock (i.e., “magma fingers”). Seismic reflection data can image entire sills or sill-complexes in 3-D, and their resolution is often sufficient to allow us to identify these distinct elongate components and thereby map magma flow patterns over entire intrusion networks. However, 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 the components (i.e., “elements”) of sheet-like igneous intrusions without linking their description to emplacement mechanisms. Using 3-D seismic reflection data from offshore NW Australia, we quantify the 3-D 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 3-D geometry, we conclude that the mapped elements likely formed through non-elastic-brittle and/or non-brittle deformation ahead of the advancing sill tip, which implies 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 data sets that span 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.

Magma moves large distances laterally and vertically through the crust via networks of interconnected sheet intrusions (i.e., dykes and sills) (e.g., Anderson, 1937, 1951; Ernst et al., 2005; Muirhead et al., 2012; Magee et al., 2016a; 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 (Elliot and Fleming, 2004; Cartwright and Hansen, 2006; Elliot and Fleming, 2008; Leat, 2008; Magee et al., 2016a). Understanding the architecture of these magma plumbing systems and how magma migrates through them is essential because they influence: (1) volcanic hazards and their warning signals (e.g., Sparks, 2003; Cashman and Sparks, 2013); (2) 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 (3) the accumulation of economic resources (e.g., Barnes et al., 2016; Spacapan et al., 2018). There is a growing consensus that sheet intrusions within magma plumbing systems do not always propagate as continuous sheets but rather as separate, elongate sections that inflate and coalesce (e.g., Pollard et al., 1975; Schofield et al., 2012a; Magee et al., 2016a; Galland et al., 2019). These components, which we refer to here as “elements,” have been studied using field observations (e.g., Pollard et al., 1975; Schofield et al., 2010; Spacapan et al., 2017; Galland et al., 2019), laboratory and numerical analogue models (e.g., Chanceaux and Menand, 2014, 2016; Bertelsen et al., 2018; Souche et al., 2019), and three-dimensional (3-D) seismic reflection data (e.g., Thomson and Hutton, 2004; Guo et al., 2013; Magee et al., 2014). Critically, the development and presence of elements can drive the localization of magma flow within individual intrusions and thereby control 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 understanding magma emplacement dynamics, we are currently faced with several challenges: (1) we do not yet have fully coupled, fully 3-D, numerical or laboratory models of magma emplacement that consider multiple emplacement mechanisms; (2) 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); (3) materials used as host rock analogues in laboratory models are limited in their rheological behaviors, which results in prescribed magma emplacement mechanisms; and (4) observations made in the field are often limited to relatively small outcrops compared to the whole intrusion and in many cases only represent a two-dimensional (2-D) view of a 3-D intrusion (e.g., Pollard et al., 1975). Seismic reflection data are a powerful tool for bridging the gap between observations made in analogue models and nature, as they allow mapping of sill morphologies in 3-D down to the tens-of-meters scale, which is at the upper limit of elements that have been mapped in the field (Galland et al., 2019). However, in many cases sill thicknesses are below the vertical seismic resolution, which means 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 3-D intrusion 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 3-D seismic reflection data from the Exmouth Plateau offshore NW Australia to study the geometry of igneous sills with a focus on their building blocks. More specifically, we: (1) introduce new descriptive terminology for the building blocks of sills (i.e., “elements”); (2) quantify the dimensions of 3-D element geometries and their connectors; (3) identify intrusion geometries considered to be representative of flow kinematics; and (4) 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 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 that result 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 orientations (Fig. 1D) (e.g., Pollard et al., 1975; Rickwood, 1990; Hutton, 2009), we conclude that the 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 data sets 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, which contributes to a better understanding of sill emplacement mechanisms at shallow depths in the Earth's crust.

Over the last few decades, several studies have shown that sills and dykes do not always propagate as continuous planar sheets; 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., 2010, 2012a; Spacapan 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 (e.g., Pollard et al., 1975; Baer and Reches, 1987; Hutton, 2009; Schofield et al., 2010, 2012a; Galland et al., 2019). We briefly summarize the formation and propagation of these features below 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 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.

Magma Fingers

Magma fingers are relatively narrow features that are elongated in the magma flow direction with blunt and/or bulbous tips (e.g., Pollard et al., 1975; Schofield et al., 2010; Spacapan et al., 2017; Galland et al., 2019). Fingers mostly propagate along the same stratigraphic level, but small vertical offsets can occur if magma exploits preferentially oriented, pre-existing weaknesses (Fig. 1A) (e.g., Schofield et al., 2010; Galland et al., 2019). Thickness-to-width aspect ratios of magma fingers measured in the field vary between 0.09 and 0.67 (e.g., Pollard et al., 1975; Horsman et al., 2005; Spacapan et al., 2017; Galland et al., 2019). Host rock deformation related to finger emplacement is observed to mainly involve: (1) brittle faulting and cataclastic flow in the compressive regime between inflating magma fingers (e.g., Pollard et al., 1975; Wilson et al., 2016, 2019; Spacapan et al., 2017; Galland et al., 2019); and/or (2) non-brittle processes (e.g., host rock fluidization, magma-sediment mingling) closer to the finger margins (e.g., Schofield et al., 2010, 2012a; Eide et al., 2017). The growth of magma fingers is therefore not dominated by a linear elastic fracture mechanics (LEFM) process. Where two fingers horizontally coalesce to form a sheet, cusp-shaped grooves containing deformed host rock form at the top and bottom sill-host rock contacts (Fig. 1A) (e.g., Pollard et al., 1975; Schofield et al., 2010, 2012a). Magma finger initiation mechanisms are poorly understood and only generally described as being due to instabilities between a propagating magma front and the host rock (e.g., Saffman-Taylor instability) (e.g., Pollard et al., 1975; Schofield et al., 2010; Spacapan et al., 2017).

Segments

Similar to magma fingers, segments are also elongated in the magma flow direction but are commonly vertically offset (Fig. 1B) (e.g., Francis, 1982; Rickwood, 1990; Schofield et al., 2012a). Unlike magma fingers, segment growth is thought to be a LEFM process, where the intrusion is treated as a fluid-filled fracture (e.g., Pollard, 1973; Delaney and Pollard, 1982; Pollard et al., 1982; Rickwood, 1990; Schofield et al., 2012a). LEFM has been used to explain the emplacement of sheet intrusions such as dykes, sills, and laccoliths (e.g., Pollard and Johnson, 1973; Nicholson and Pollard, 1985; Takada, 1990; Lister and Kerr, 1991; Kavanagh et al., 2006; Bunger and Cruden, 2011). In this scenario, the host rock is assumed to deform in a linear elastic manner, and the fluid-filled fracture propagates by mode I failure (i.e., tensile opening) of the host rock that is driven by the magma overpressure (e.g., Pollard, 1973).

During sill emplacement, two different mechanisms can lead to branching of the propagating edge and the formation of discrete segments: (1) a rotation of the principal stress axes ahead of the propagating tip that results 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 (2) the exploitation of pre-existing, preferentially oriented weaknesses with a lower tensile strength and/or fracture toughness (e.g., bedding planes, fault planes) that likely lead to an inconsistent stepping direction (e.g., Hutton, 2009; Schofield et al., 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, which forms a continuous sheet with a step-like geometry (Fig. 1B) (e.g., Rickwood, 1990; Hutton, 2009; Schofield et al., 2012a, 2012b).

The thickness-to-width ratios of segments tend to be small compared to those of magma fingers. If the host rock deforms in a linear elastic manner, then:
formula
where W = characteristic thickness, P = magma overpressure relative to the minimum in situ stress, L = characteristic crack length (radius for a circular sill; width [shorter dimension] for a segment), and the stiffness parameter E′ = E/(1–v2), where E = Young's modulus and v = Poisson's ratio (Olson, 2003). If crack propagation is governed by LEFM, then:
formula

where KC = fracture toughness and PC = critical overpressure (Olson, 2003). Under the assumption that the crack is propagating, P = PC, 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 × 106 Pa m1/2 and 20–55 GPa, respectively) or shale (1.4 × 106 Pa m1/2 and 1–20 GPa, respectively) gives an upper W/L aspect ratio of ~0.00001–0.0004.

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 3-D geometry is missing. We therefore propose the term element to generally describe the building blocks of sills (i.e., magma fingers and segments) without attempting to infer their 3-D 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 (Figs. 1A–1B). When overlapping or underlapping elements inflate, brittle failure of the host rock separator may cause both elements to connect, which forms broken bridges and steps, respectively (e.g., Rickwood, 1990; Hutton, 2009; Schofield et al., 2012a). Broken bridges and steps have been used as evidence of 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, which forms a step-like geometry (Fig. 1A) (Galland et al., 2019). 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-connector 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 connectors form parallel to the long dimension of elements such that their long axes can be used to determine the propagation axes 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 fingers and segments) and placed them within our new element-based terminology. In addition, numerous studies have shown that sheet intrusions can occasionally be subdivided into larger sections, each of which is 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., Thomson and Hutton, 2004; Hansen and Cartwright, 2006a; Schofield et al., 2012b). We also designate larger sheet components, including lobes, 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 in which magma fingers and segments represent the smallest end-members in the element hierarchy (Figs. 1A and 1C). A large, first-order element consisting of multiple orders of sub-elements regardless of their geometries is termed an element network (Figs. 1A and 1C).

The geometry of any element in map view has 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 groups: (1) elongate elements with parallel lateral edges (Type 1); (2) elongate elements that are widening in magma flow direction (Type 2); and (3) semi-radially fanning elements that rapidly increase in width (Type 3). 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: (1) changing magma flow kinematics (e.g., Magee et al., 2016b); (2) discrete magma injections (e.g., Thomson and Hutton, 2004; Horsman et al., 2016; Magee et al., 2016a); or (3) 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 and thus excludes any interpretation of emplacement timing or deformation mechanisms.

The Exmouth Plateau is part of the North Carnarvon Basin, which is located on the rifted margin of northwest Australia and formed following multiple episodes of extension during the late Paleozoic and Mesozoic (Fig. 2A) (e.g., Symonds et al., 1998; Karner and Driscoll, 1999; Longley 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; Exon et al., 1992; Symonds et al., 1998). Formation of the Tethys Ocean during the Permian resulted in extensive stretching and thinning of the crust, which was 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–2D) (e.g., Exon and Buffler, 1992; Exon et al., 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–2D) (e.g., Exon and Buffler, 1992; Stagg et al., 2004; Black et al., 2017; Bilal et al., 2018).

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 ultramafic rocks was emplaced into the lower crust ~165–136 m.y. ago and potentially acted as the source of the extensive network of sills and dykes observed across the 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 of which are interpreted to be related to sill intrusion, suggests that sill emplacement across the North Carnarvon Basin started in the Kimmeridgian (ca. 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 by radiating dykes of the Tithonian (ca. 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 sedimentary rocks.

Seismic Reflection Data and Seismic Resolution

This study uses a zero-phase, time-migrated, 3-D seismic reflection survey (Glencoe 3D) covering an area of ~4042 km2 (Fig. 2). In-lines and cross-lines are oriented E-W and N-S, respectively, and have regular spacing of 25 m. These data image to a depth of 8 s two-way time (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 ms TWT to 1664 ms TWT and, assuming a water velocity of 1500 m s–1, between 1001 m and 1248 m, respectively.

Borehole data from wells near the sills studied were used to: (1) constrain the ages of mapped key horizons and their lithology; (2) perform a tuning wedge model; and (3) 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 boreholes (e.g., Reeves et al., 2018; Norcliffe et al., 2021). To accurately constrain sill depths and thicknesses, we used different time-depth conversions for sills S1 and S2 with 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 Material1). The coefficient of determination (R2 = 0.99) of the estimated 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) resolved in this survey mainly occur at depths of 3200–5150 ms TWT (~3000–6600 m) and are not penetrated by boreholes, which means 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 that 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; 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 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 are 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, which thus makes 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 represent the true sill-host rock contact, which results 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, which results in ~56 m (λ≈224.5 m) within the vicinity of the sills studied (e.g., Lebedeva-Ivanova et al., 2018).

Mapping Sills and Primary Magma Flow Indicators in 3-D Seismic Reflection Data

With the exception of a seismically unresolved sill intersected by Rimfire-1, no other boreholes in the study area penetrate sills. Therefore, to identify sills within the seismic reflection data, we use the criteria defined by Planke et al. (2005) and classify reflections as sills or inclined sheet intrusions where they: (1) have high amplitudes; (2) display a positive polarity (a downward increase in acoustic impedance); (3) locally transgress host rock strata; and (4) 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 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 3-D geometric interpretations including those of elements and element connectors (e.g., Magee et al., 2015; Eide et al., 2018), which hampers the ability to interpret magma emplacement mechanisms. Synthetic seismic forward models have shown that image quality can also be influenced by the host rock style (i.e., homogenous versus 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, which results 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 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: (1) reconstruct magma flow pathways and potential sill emplacement history; and (2) identify potential sill feeder zones as well as larger source-magma reservoirs (e.g., Thomson and Hutton, 2004; Schofield et al., 2012b, 2017; Magee et al., 2014). Connectors between adjacent elements are typically expressed as low-amplitude stripes on top-sill amplitude maps (e.g., Schofield et al., 2012b, 2017; Magee et al., 2015) and result from sill thickness variations (Magee et al., 2015). For example, where elements propagate along the same 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-vertical connectors (Fig. 4D). This sill thickening causes a decrease in amplitude if the sill thickness is above the maximum tuning thickness (Figs. 3 and 4D) (Magee et al., 2015). These findings are based on synthetic seismic forward models (Magee et al., 2015) and suggest that low-amplitude 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 tenth inline and crossline, which resulted in a regular grid spacing of 250 m. These grids were carefully quality checked and used as seed points to propagate (i.e., auto-track) surfaces. We further identified and mapped the primary magma flow indicators of sills S1 and S2.

Quantitative Analyses

We used a combination of seismic horizons of sills (map-view) and related seismic sections (cross-section view) to map sill morphologies and included 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 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: (1) the long axis of elongate elements with width-to-length aspect ratios of < 1; or (2) the long axis of a feeder conduit for elements with width-to-length aspect ratios that are > 1 (Fig. 4B). The angle α is not constant for a single element but rather increases and decreases along the defined reference line, and α > 0 indicates widening (i.e., increase in element width) and α < 0 indicates narrowing (i.e., decrease in element width) (Fig. 4B). To ensure accuracy and to assess how α varies along the reference line, the element margins are partitioned into smaller sections (Fig. 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 margin sections analyzed 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 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 [see footnote 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 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 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 (Fig. 3C).

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 element tips, which we then converted to depth (s TWT to m) (Figs. 4C–4D). These data were collected in ~100 m spaced seismic sections starting from the most inward part of the connector toward its tip and perpendicular to the strike of the mapped contact geometries (H- and S-/Z-connector).

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, which highlights the downward increase in acoustic impedance at host rock-sill contacts (e.g., Fig. 5). Base-sill contacts are 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 that include: (1) strata-discordant inclined sheets (S1); and (2) 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 (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 5D). 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–5E). 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–5E and 6A–6C). The top-sill reflection is of relatively high-amplitude at shallow levels (< 4 s TWT) and decreases with depth (Figs. 5B and 5D). It is important to note that the general relationship of relatively high-amplitude reflections at shallow levels holds irrespective of whether the sill is characterized by distinct or tuned reflections (Figs. 5B and 5E).

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, which indicates contact geometries (H- and S-/Z-connectors; Figs. 5E and 6C). At relatively shallow emplacement depth (~3.8–3.5 s TWT), we observe four orders of sub-horizontal, bedding-parallel elements radiating from the inclined sheet, which form 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 are vertically offset, adjacent elements either overlap or underlap and often coalesce, which forms H- and S-/Z-connectors, respectively, with an inconsistent stepping direction (Fig. 6C).

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-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 m to 1808 m, and fault-related connectors tend to be longer (Fig. 6D). Where connectors coincide with pre-existing fault planes, their maximum vertical height ranges from 22 m (±17) to 257 m (±17), whereas connectors that are not coincident with fault planes have maximum vertical heights of 4–121 m (±17) (Fig. 7A). Minimum vertical offsets mostly occur at the most inward part of the elements (55%) or at their propagating tips (21%) (Fig. 7B). Most minimum vertical height measurements (78%) range from 0 m to 25 m; however, greater heights of up to 154 m (±17) locally occur along the connector's length (Fig. 7B). Maximum vertical heights are evenly distributed along connector lengths with a slight increase at the element tips (18%) (Fig. 7A). Maximum connector heights range up to 257 m (±17) with the majority (94%) being <100 m.

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 (Fig. 8A). Pre-existing Jurassic faults define the sill borders, which results in an overall NNE-elongate geometry (Fig. 8B). Inwardly inclined, moderately dipping (~23–35°) transgressive limbs 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) and transgress 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 toward the north-northeast (Figs. 8A–8B). Where S2 is not bounded by an inclined limb in the west, two westward-diverging, sub-horizontal element networks (EN1, EN2) form an arcuate geometry with a northwesterly termination (Figs. 8B and 8D). Sill S2 comprises two morphologies separated by a NE-SW–striking normal fault: (1) an elongate saucer-shaped sill in the south; and (2) 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 and define a maximum thickness of 98 ms TWT (~272 m (±27)). In contrast, the inclined limbs and the sub-horizontal outer margin in the south and the northwest (EN2) are imaged as tuned reflection packages (Figs. 8B–8C). The elongate section of S2 thins from east to west (Fig. 8C). The TWT thickness 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 elements include the semi-radially fanning element networks EN1 and EN2 in the northern section of S2, which we subdivide 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 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 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, bedding-concordant element networks (EN1, EN2) in the north have strikes that fan toward NW and range from N to W (Fig. 8E).

Amplitude and Thickness Variations of Elongate Elements

Amplitudes of the top-sill horizon and thickness measurements of S2 were collected along two 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, and minor fluctuations occur along these high-amplitude plateaus (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 m (±12.5) to 185.5 m (±18.5) 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 toward their margins (Fig. 9A; E4, 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, which results in a measurement of the S-/Z-connector thickness rather than measurements of the single elements and an apparent local sill thickening (Fig. 9A; contacts between E2–E3 and E4–E5). 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 map view and potentially highlight the contacts of strata-concordant, third-order elements, which within 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 m (±0.7) and 56 m (±5.6), 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-X′).

Along seismic section Z-Z′ parallel to the strike of element connectors, we show that thickness and seismic amplitude vary in both cases where the sill is resolved as separate reflectors (EN1) and tuned reflections (EN2; Fig. 9C). EN1 gradually thins toward its tip (132 m (±13) – 99 m (±10)), while the thickness drops below the limit of separability in EN2 (Fig. 9C). Minor amplitude variations occur within EN1, and low-amplitude values define the contact between EN1 and EN2. Amplitudes within EN2 fluctuate slightly, and a peak occurs at 2750 m along the profile that decreases rapidly toward 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–9C) indicates a high variance of 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 the maximum tuning thickness is reached, after which they decrease (Fig. 9D). The minimum apparent thickness observed is 38 m, which is only slightly below the maximum tuning thickness 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 [footnote 1]) show similar amplitude trends, but the predicted maximum tuning thickness of 43 m is 8 m lower than is observed in our data. Contrary to data collected along the seismic sections described above, modeled amplitudes for thicknesses above the limit of separability are constant (Fig. 9D).

Element Geometries

We observe three types of element geometry in S1 and S2 that include 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).

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 toward the tip (α < 0) (Figs. 4B and 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 a 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°) than Type 2B elements (n = 9; αmean =–1–27°) (Fig. 10B). This occurs because Type 2B elements narrow at their tips with α < 0, which results in smaller αmean angles (Figs. 3B and 10A). Since Type 2 elements have elongate geometries, their width-to-length aspect ratios are < 1 (0.27–0.93; average of 0.61), and their element lengths range over two orders of magnitude (181–5061 m) (Fig. 10C).

Type 3—semi-radially fanning elements with a moderate to large α angle (~30–137°). These elements are often fed via a relatively short and narrow conduit with a small to moderate α angle of < ~50° before there is 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 widths and lengths range from 370 m to 3940 m and 345–2348 m, respectively, which results in width-to-length aspect ratios predominantly > 1 (0.77–2.37; average of 1.42), and aspect ratios < 1 (0.77 and 0.78) occur 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 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 α, which results in an approximately linear and a semi-radially fanning trend within high-order Type 2 and Type 3 elements, respectively (Figs. 10A and 10D).

Because comprehensive analyses of the full 3-D geometries of sub-horizontal sheet intrusions are limited to: (1) geophysical techniques such as 3-D seismic reflection data, and (2) numerical and laboratory analogue modeling, it is important to link observations of both approaches 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 discuss the implications of our results by considering how sill geometries and 3-D seismic reflection data can help to infer the 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: (1) viscous emplacement mechanisms such as host rock fluidization (e.g., Kokelaar, 1982; Schofield et al., 2010); (2) elastic bending and uplift of the overburden strata potentially due to LEFM (e.g., Johnson and Pollard, 1973; Walker, 2016); and (3) viscous indentation and plastic shear failure (e.g., 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 sills and, where they do, they provide only 1-D information about rock properties. To interpret sill emplacement mechanics from seismic reflection data, we therefore currently rely on measurement of overburden uplift (forced folds) above sills, which is 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 to pore fluid expulsion, which potentially indicates host rock fluidization and the 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, which suggests magma emplacement through elastic-brittle fractures (e.g., Francis, 1982; Schofield et al., 2012a). 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 elastic-brittle fracturing (Galland et al., 2019). 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, Fingers 1 and 2) (Galland et al., 2019). 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 are not available, and highlights the importance of using the descriptive term element where mechanisms that accommodate both sill emplacement and branching cannot be identified.

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 indicate that significant variations in vertical element offsets 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 elastic-brittle fracturing or brittle and/or non-brittle deformation mechanisms (i.e., viscous indentation and 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 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., 1975, 1982; Rickwood, 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 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), which indicates that syn-emplacement rotation of principal 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, which indicates 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 and 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 shallower levels < 3.8 s TWT (Figs. 5A–5B). Seismic reflections with relatively higher amplitudes are observed for both distinct and tuned reflections at similar stratigraphic levels (Figs. 5B and 5E), which suggests a relative change from a mechanically strong host rock toward a weaker rock. Lithological changes and increasing cementation and compaction with depth can increase the mechanical host rock strength (e.g., Kokelaar, 1982; Schofield et al., 2010, 2012a). Such burial 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., 2010, 2012a). We hypothesize that the relative change toward mechanically weaker host rocks at shallower emplacement depths of ~950–1150 m in our data resulted in non-brittle emplacement of S1 and, consequently, the formation of magma fingers.

Our data further indicate 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 convex upwards, they form cusp-shaped grooves at their points of coalescence, which results in an undulating top-sill geometry, where the sill is thinner at magma finger contacts (Figs. 1A and 4D). Contact geometries of coalesced magma fingers are therefore expressed as linear features with relatively low amplitudes; they are 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). 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?

Below, we interpret element geometries described in this study and 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.

Thickness 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, 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 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). 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). Based on the data available, we can neither determine the exact duration of sill emplacement nor the time 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 indicate that similar patterns of thickness variation occur at smaller scales for the different orders of elements (Figs. 8C and 9A–9B). These patterns can be accentuated by flow 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 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 (Figs. 11C–11D) (Magee et al., 2013c). Identifying this channelized magma flow in sheet intrusions has important implications for: (1) better understanding of the architecture of sill complexes and how they are built (e.g., Thomson and Hutton, 2004; Cartwright and Hansen, 2006; Hansen and Cartwright, 2006a; Guo et al., 2013; Magee et al., 2016a; Schofield et al., 2017); (2) assessing areas of potential fissure eruptions emanating from areas with localized magma flow (e.g., Pansino et al., 2019); and (3) 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).

The Growth of Elements

When a single element emerges from a planar sheet, it grows in length and thins toward 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 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 agrees 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 and eventually break 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 (Figs. 12A–12B). Higher order elements tend to be thicker close to their feeder conduits and thin out toward 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 into consolidated sedimentary rocks at ~950–1150 m depth.

Our growth model is compatible with other sills mapped in 3-D seismic reflection data. For example, data from the North Rockall Trough images a sill with “flow units” of varying orders that are 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 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 and implies 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 3-D seismic reflection data (Figs. 11A–11B). For example, consider a scenario where a first-order element contains second- and third-order elements, with the latter representing the smallest seismically resolved order (Fig. 11). This element network may have formed in response to either: (1) a single magma pulse (Fig. 11A) or (2) 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, 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: (1) the emplacement of elements is geometrically self-similar over multiple scales whereby elements with different orders may grow independently with potentially independent flow kinematics; and (2) 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 (Figs. 1112).

Emplacement History and Magma Flow Pathways in the Sills Studied

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 sill-induced forced folds indicate that magmatism occurred during the Upper Jurassic to Early Cretaceous, which agrees with previous studies (e.g., Exon et al., 1992; Mihut and Müller, 1998; 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., 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 they can only be recognized in high-resolution seismic reflection data as localized areas of low-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. However, intrusions are rarely fully exposed, which makes sample collection and a consistent spatial distribution of data challenging. Three-dimensional 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 orientations are a magma flow pathway indicator. We highlight that contact geometries form between adjacent, vertically offset elements (i.e., segments and 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, which indicates an overall westerly flow direction (Fig. 6). At greater depths, only large-scale magma flow indicators are recognized in the form of potential first-order elements, which is most likely due to a decrease in seismic resolution (Fig. 5B). However, based on the available data, we infer that S1 is an inclined sheet that propagated upwards and cross-cut the host-rock strata with at least one feeder located at the deepest levels in the SE (Fig. 5). Magma flowed toward 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, which was possibly defined by a change to a weaker, less brittle host rock, may have allowed for non-brittle deformation processes and potentially the formation of bedding-parallel magma fingers that fanned out westward and formed a large-scale element network (Figs. 56). Magma fingers propagated along pre-existing weaknesses (e.g., bedding planes), which resulted 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 planes, which resulted 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, which formed a broadly continuous sheet containing abrupt intra-sheet steps.

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–13C). When the first pulse of vertically ascending magma interacts with the pre-existing fault, local rotation of σ3 from the sub-horizontal toward an orientation perpendicular to the fault surface allows magma to propagate along the fault (Figs. 13A and 13I) (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): (1) 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 will propagate along the interface of host rock layers and result in the transition to a sub-horizontal sill (t1; Fig. 13A, ii–iv) (e.g., Kavanagh et al., 2006, 2015); or (2) 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), which results in the emplacement of a sub-horizontal 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 footwall in a SE direction (Figs. 13A–13C). Inflation and associated tensile failure may have caused sill propagation into the hanging wall toward 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, which formed 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's long axis (e.g., Galerne et al., 2011; Fyfe et al., 2021) and that pre-existing fault planes can greatly influence 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–13D). 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., 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: (1) 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); (2) magma propagation along reverse faults that formed due to overburden uplift during sill inflation (e.g., Thomson and Schofield, 2008; Haug et al., 2017); or (3) 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 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 and in some cases slightly vertically offset elements north of the potential 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 is 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). 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 pre-existing structures.

A large variation in amplitudes is observed for thicknesses above the limit of separability (Fig. 9D), which we interpret as due to variations in host rock lithology or sill composition. Since the latter scenario would likely result in the 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 elements and may have caused minor vertical offsets, which we observe between the second-order elements of S2.

High-quality 3-D seismic reflection data from the Exmouth Plateau, offshore NW Australia, were used to: (1) map basaltic sills and their elements; (2) document their geometries and seismic expressions; and (3) introduce new descriptive terminology to define different sill elements and connectors based on their geometries. Our key conclusions are:

  1. Five Late Jurassic to Early Cretaceous basaltic sills intruded interbedded sandstones and claystones of the Triassic Mungaroo Formation 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-shaped, and sub-horizontal sill morphologies.

  2. Sub-horizontal, semi-radially fanning element networks occur in two sills where: (1) an inclined sheet intrudes into a potentially mechanically weaker host rock (S1), or (2) a saucer-shaped sill is bounded by only one inclined limb (S2).

  3. Quantifying vertical heights of H- and S-/Z-connectors in 3-D seismic reflection data provides important insights into sill branching mechanisms. This approach helps to identify/exclude the syn-emplacement rotation of the orientations of principal stress axes ahead of a propagating fracture as a 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 agree with the results of synthetic forward models (Magee et al., 2015) and permit the analysis of magma flow pathways in 3-D seismic reflection data at the meter scale and help to infer potential sill emplacement mechanisms.

  5. Based on quantitative measurements of element geometries, we classified elements into the following groups: (1) elongate elements with parallel lateral edges (Type 1); (2) elongate elements that consistently widen and have width-to-length aspect ratios of < 1 (Type 2); and (3) semi-radially fanning elements that are often fed via a relatively short and narrow conduit with width-to-length aspect ratios predominantly > 1 (Type 3). Quantifying element geometries permits comparisons between different subsurface and field-based data sets 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 divided sill S2 into three main zones based on their thicknesses, where the sill is thicker close to its feeder. This suggests that 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 that resulted in magma flow localization.

  8. Thickness-to-width aspect ratios of individual elements at the outer margin of S2 are smaller than those measured toward the center of S2, which suggests that elements first grow in width before they inflate.

We are grateful to Geoscience Australia for making the seismic reflection and borehole data publicly available at https://nopims.dmp.wa.gov.au/nopims. We also thank Down Under GeoSolution and Schlumberger for the 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. J. Köpping acknowledges a Monash Graduate Scholarship. A.R. Cruden, C.A.-L. Jackson, and C. Magee acknowledge funding from Australian Research Council Discovery Grant DP190102422. We thank Andy Bunger for fruitful discussions on element geometries.

1Supplemental Material. Uninterpreted seismic reflection images (Figs. 5, 6, and 8); Jupyter notebooks and associated input- and output-files to: (1) describe the time-depth conversion of data, (2) create tuning wedge models, (3) quantify connector heights, and (4) quantify element geometries; and a QGIS project including all shapefiles and georeferenced maps used for this study. Please visit https://doi.org/10.1130/GEOS.S.16959634 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: Shanaka de Silva
Associate Editor: Robert B. Miller
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.