Seismic reflection data reveal the 3D structure of the newly discovered Exmouth Dyke Swarm, offshore NW Australia

Abstract. Dyke swarms are common on Earth and other planetary
bodies, comprising arrays of dykes that can extend laterally for tens to
thousands of kilometres. The vast extent of such dyke swarms, and their
presumed rapid emplacement, means they can significantly influence a variety
of planetary processes, including continental break-up, crustal extension,
resource accumulation, and volcanism. Determining the mechanisms driving
dyke swarm emplacement is thus critical to a range of Earth Science
disciplines. However, unravelling dyke swarm emplacement mechanics relies on
constraining their 3D structure, which is difficult given we typically
cannot access their subsurface geometry at a sufficiently high enough
resolution. Here we use high-quality seismic reflection data to identify and
examine the 3D geometry of the newly discovered Exmouth Dyke Swarm, and
associated structures (i.e. dyke-induced normal faults and pit craters).
Dykes are expressed in our seismic reflection data as ∼335–68 m wide, vertical zones of disruption (VZD), in which stratal
reflections are dimmed and/or deflected from sub-horizontal. Borehole data
reveal one ∼130 m wide VZD corresponds to an ∼18 m thick, mafic dyke, highlighting that the true geometry of the inferred
dykes may not be fully captured by their seismic expression. The Late
Jurassic dyke swarm is located on the Gascoyne Margin, offshore NW Australia,
and contains numerous dykes that extend laterally for > 170 km,
potentially up to > 500 km, with spacings typically < 10
km. Although limitations in data quality and resolution restrict mapping of
the dykes at depth, our data show that they likely have heights of at least
3.5 km. The mapped dykes are distributed radially across a
∼39∘ wide arc centred on the Cuvier Margin; we
infer that this focal area marks the source of the dyke swarm. We demonstrate
that seismic reflection data provide unique opportunities to map and quantify
dyke swarms in 3D. Because of this, we can now (i) recognise dyke swarms
across continental margins worldwide and incorporate them into models of
basin evolution and fluid flow, (ii) test previous models and hypotheses
concerning the 3D structure of dyke swarms, (iii) reveal how dyke-induced
normal faults and pit craters relate to dyking, and (iv) unravel how dyking
translates into surface deformation.


1a and c) (e.g., Buchan and Ernst, 2013;Baragar et al., 1996;Odé, 1957;Walker, 1986); and (iii) circumferential dyke swarms, which likely emanate from the lateral termination of a plume head, although the stress state controlling their emplacement remains poorly understood (e.g., Fig. 1d) (e.g., Buchan and Ernst, 2018a, b). Component dykes within dyke swarms can be up to 10's or 100's m thick and their emplacement is primarily accommodated by extending the host rock (e.g., Jolly and Sanderson, 1995;Rivalta et al., 2015;Gudmundsson, 1983;Paquet et al., 2007). Their geometry and scale 35 means dyke swarms can thus drive crustal extension, influencing plate tectonic processes on Earth and shaping other planetary bodies (e.g., Ernst et al., 2013;Ernst and Buchan, 1997;Wilson and Head, 2002;Ebinger and Casey, 2001;Halls, 1982;Paquet et al., 2007). Because they are typically emplaced over short timespans (≲5 Myr) and are sensitive to the prevailing stress field, dyke swarms also provide a record of syn-emplacement stress conditions and represent key spatial and temporal markers for palaeogeographic and palinspastic reconstruction (e.g., Hou et al., 2010;Bleeker and Ernst, 2006;Ju et 40 al., 2013;Peng, 2015;Halls, 1982). Furthermore, dyke swarms may be associated with the accumulation of critical economic resources (e.g., Jowitt et al., 2014;Ernst and Jowitt, 2013) and, if they feed extensive flood basalts, may contribute to climate change and related mass extinctions (e.g., Ernst and Youbi, 2017). Unravelling the emplacement history of dyke swarms and deciphering the processes controlling their intrusion and form, is therefore crucial to a wide range of pure and applied Earth Science disciplines. 45 Decoding dyke swarm emplacement requires knowledge of their 3D structure, which is typically inferred by quantifying and projecting downwards the plan-view morphology of dykes exposed at Earth's surface or identified in airborne/satellite imagery (e.g., Ernst and Youbi, 2017;Ernst, 2014;Bryan and Ernst, 2008;Coffin and Eldholm, 2005;Coffin and Eldholm, 1994;Bryan et al., 2010;Hou et al., 2010;Halls and Fahrig, 1987;Halls, 1982;Ernst and Baragar, 1992). Such inferences of 3D structure may be augmented by direct mapping of the local subsurface structure of dyke 50 swarms, or component dykes, intersected in mines or imaged in geophysical data (e.g., Kavanagh and Sparks, 2011;Keir et al., 2011;Wall et al., 2010). Whilst integrating these datasets has emphasised the lateral variability in dyke swarm architecture, a recent seismic reflection-based study has highlighted dyke swarm structure can also change with depth (Phillips et al., 2018). In particular, Phillips et al. (2018) demonstrated the width of a dyke swarm imaged offshore southern Norway increased with depth, implying the plan-view morphology of a dyke swarm may not be a proxy for its 3D geometry 55 (or total volume); i.e. the plan-view morphology of a dyke swarm is a function of its attitude relative to the present topography (see also Magee et al., 2019). We can use different physical, analytical, and numerical modelling approaches to evaluate the 3D geometry of dyke swarms, and to establish how their structure can be inferred from principally 2D, surfacebased analyses. However, these model predictions are difficult to validate without constraints on the true 3D form of natural dyke swarms (e.g., Paquet et al., 2007;Bunger et al., 2013;Jolly and Sanderson, 1995;Macdonald et al., 1988). Advancing 60 our understanding of dyke swarm emplacement thus requires a method for imaging their 3D structure in detail (e.g., Magee et al., 2019;Magee et al., 2018;Phillips et al., 2018).

Dataset and methods
Dykes are rarely imaged in seismic reflection data because their sub-vertical orientation preferentially reflects seismic energy deeper into the subsurface, rather than returning it to the surface to be recorded (e.g., Eide et al., 2018;Thomson, 2007).
Dykes identified in the field and/or in aeromagnetic data have been indirectly recognised in co-located seismic reflection 135 data where a localised reduction in returned seismic energy disrupts the continuity and strength (amplitude) of reflections associated with stratigraphic layering (e.g., Fig. 2) (e.g., Wall et al., 2010;Bosworth et al., 2015;Kirton and Donato, 1985;Ardakani et al., 2017); i.e. in these cases, dykes do not correspond to discrete reflections, but instead appear as 'vertical zones of disruption' (VZDs). Whilst dykes can thus be recognised in seismic reflection data, vertical strike-slip and normal faults, and non-magmatic fluid flow conduits (e.g., gas chimneys) may also be expressed as VZDs. To avoid interpretational 140 bias, we describe the features of interest in this study as VZDs, and collect additional data and make further observations to inform a critical discussion of their likely origin.
We use eight 3D and 63 2D, time-migrated seismic surveys to map 26 VZDs across ~40,000 km 2 of the North Carnarvon Basin (Figs 4a and b); the properties of each seismic survey are provided in Supplementary Table S1. Visual inspection of the data and extraction of variance volume attributes, which highlight trace-to-trace variations in seismic 145 wavelets to reveal structural (e.g., faults and VZDs) and stratigraphic (e.g. channel edges) discontinuities (Brown, 2011), allow us to identify VZDs in the 3D seismic volumes. These VZDs were mapped on sections oriented orthogonal to their strike every ~250-1200 m. In places, the VZDs were obscured by tectonic faults and could not be mapped at regular intervals. Along-strike projection of mapped VZDs outside of the 3D seismic volumes guided their interpretation on 2D seismic lines, where poorer data quality and/or lower resolution hindered their recognition; we were able to confidently 150 recognise VZDs in nine 2D seismic surveys (e.g., Figs 4c and d), although we cannot rule out their presence in other datasets.
In addition to mapping VZDs, we used biostratigraphic and well-log data from 24 wells to identify and interpret two key stratigraphic horizons across the study area: (i) the ~148 Myr near Base Cretaceous unconformity (BC); and (ii) the near Top Mungaroo Formation (TM), which is broadly equivalent to the Norian-Rhaetian boundary (i.e. intra-Upper Triassic) 155 ( Fig. 4a; Supplementary Fig. S1). Where we observed fluvial channels within the Triassic strata using variance time-slices (e.g., Fig. 5), we locally mapped intra-Mungaroo horizons to assess channel continuity across identified VZDs; this helped us assess VZD kinematics. We also interpreted key structures associated with the VZDs, including overlying normal fault systems, pipes, and sub-circular depressions.

Quantitative analysis 160
The mechanics and dynamics of dyke swarm emplacement controls and is reflected in the geometry of its component dykes (e.g., Bunger et al., 2013;Jolly and Sanderson, 1995;Mège and Korme, 2004;Gudmundsson, 1987). For example, dyke lengths within a swarm are expected to display a power-law distribution and may be used to differentiate feeder and nonfeeder dykes within a given population (Mège and Korme, 2004). The range of dyke thicknesses, which follows a Weibull distribution, and spacings can also provide insights into the strength of host rock and/or magma source conditions 165 (Krumbholz et al., 2014;Bunger et al., 2013). These predicted distributions for dyke properties within a swarm also allow us to test whether an observed dyke set comprises one or multiple generations of intrusion, perhaps originating from different sources (e.g., Krumbholz et al., 2014). We quantify VZD structure and compare our results to predicted distributions to help unravel the mechanics and dynamics of VZD formation.
We measured the plan-view, tip-to-tip length (L) and strike (S) of each VZD (Fig. 5). Many VZDs display minor 170 but abrupt changes in strike along their length (e.g., Fig. 5). These minor changes in strike sub-divide the VZDs into discrete planar segments, for which we measured strike (s) and length (l) (Fig. 5). Where coverage of 3D seismic volumes was sufficient, we also measured VZD width (w) and spacing (h) orthogonal to strike, on variance time-slices at 4.5 s two-way time (TWT) (Fig. 5); we specifically measured w and h, as well as the depth to VZD tips, along 35 ~E-W trending transects spaced ~4.7 km apart. Because data quality generally decreases with depth within individual seismic surveys, defining the 175 base of individual VZDs is problematic, making it difficult to ascertain whether most VZDs truly terminate downwards or if they extend below the survey limits. We therefore only qualitatively assess VZD height (H).

Seismic resolution
We used time-depth plots derived from the checkshot data available for the 24 wells to estimate seismic velocities (Supplementary Fig. S3 and Supplementary Table S2). Because the VZDs extend below the total depth of all wells, we 180 estimated seismic velocities (v) through the interval of interest by extrapolating a second-order polynomial trend-line through the cumulative checkshot data ( Supplementary Fig. S3). The dominant frequency (f) of the 2D and 3D seismic surveys broadly decrease with depth from a maximum of ~30-40 Hz at the top of the interval of interest (~2.8-2.9 s TWT; ~2.5-2.7 km) to a minimum of ~5-20 Hz at ~5.9-6.0 s TWT (~9.7-10.1 km). We calculated the average interval velocities for ~2.8-2.9 s TWT (~3.0 km s -1 ) and ~5.9-6.0 s TWT (~6.4 km s -1 ). Coupled with the dominant frequency data, these 185 average interval velocities allowed us to estimate the dominant wavelength (λ = v/f) of the data and constrain the limits of separability (~λ/4) and visibility (~λ/30) (Brown, 2011). The limit of separability corresponds to the minimum vertical distance between two interfaces required for them to produce distinct seismic reflections within a survey (Brown, 2011). If the vertical distance between two interfaces is between the limits of separability and visibility, their reflections will interfere and cannot be deconvolved; i.e. they produce tuned reflection packages (Brown, 2011). Interfaces separated by vertical 190 distances less than the limit of visibility will be indistinguishable from noise (Brown, 2011). Our calculations indicate the https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License.
limits of separability and visibility at the top of the interval of interest, within the Early Cretaceous Barrow Group, are ~19-25 m and ~2-3 m, respectively. Towards the base of the 3D seismic surveys at ~5.9-6.0 s TWT, the limits of separability and visibility decrease to ~80-320 m and ~11-43 m, respectively. The lateral resolution of the data similarly decreases with depth, from 10-12 m to ~31-62 m. 195

Errors
Here we carefully consider the errors associated with our quantitative analysis of VZD geometry. For example, synthetic seismic forward modelling indicates dyke-related VZD width is dependent on data quality and resolution, and thus likely does not equal dyke thickness (Eide et al., 2018). Data quality and resolution, in turn, is influenced by a range of geophysical (e.g., acquisition and processing parameters) and geological (e.g., faults may locally inhibit imaging) factors. The different 200 acquisition and processing histories of the seismic surveys we use, coupled with spatial variations in the geology of the study area, therefore makes it challenging to assess the likely errors associated with our measurements of VZD width and spacing; i.e. we cannot easily determine how closely the VZD geometry reflects the thickness and spacing of the structures they correspond to (e.g., Fig. 5). The local strike and dip of VZDs may also potentially differ from that of their corresponding structure(s) (e.g., Fig. 5), although we consider these variations to be negligible given their high length-to-width and height-205 to-width aspect ratios. Because we do not know how seismic velocity varies laterally away from areas of borehole control, we do not depth-convert the seismic reflection data, instead presenting measurements in time (milliseconds TWT) rather than depth (in metres). However, to help geoscientists working with field data and to provide an overall sense of scale, we provide approximate depth-converted value (in metres) for each measurement in time but we cannot ascertain the accuracy of these conversions. Overall, in an attempt to account for and visually communicate the potential errors described above, as 210 well as those introduced by human bias during interpretation and measurement, we conservatively consider that each quantitative parameter could have an arbitrary error of either: (i) ±10% if the property analysed is measured in time (e.g., VZD upper tip depth); or (ii) ±50 m if distances (e.g., VZD length, width, and spacing) are measured in plan-view.
the VZDs are subtle and typically only marked by a reduction in amplitude and/or minor geometrical distortion of the stratigraphic reflections they cross-cut (e.g., Figs 6d-f and 7). On some 2D and 3D seismic sections, particularly where data quality is poor and tectonic faults inhibit imaging, we could not recognise VZDs in locations where we predicted them to occur based on their along-strike projection (e.g., Fig. 7c). Conversely, we identified some additional VZDs on individual 225 2D seismic lines but could not map these on neighbouring sections located as little as 5 km along-strike (e.g., Fig. 6f); in these cases it was difficult to determine if the VZDs truly terminated along-strike, or whether they were simply not imaged on adjacent lines. Where VZDs cross-cut pre-existing fluvial channels or linear structures within the Mungaroo Formation, there is no resolved vertical or lateral offset of these potential host rock strain markers (e.g., Figs 5 and 8).

Geometry 240
In plan-view, the VZDs are linear, ranging in length (L) from ~4-171 km and with tip-to-tip strikes (S) between 353° and 021° (Figs 4 and 10b; Table 1). Overall, the VZDs have a mean S of 008° and broadly display a westwards progression from ~NNE-SSW striking to ~NNW-SEE striking (Figs 4 and 10b). Only the ~N-S striking (002°) VZD B intersects other VZD traces (i.e. VZDs C and D; Fig. 4); the resolution of the data is insufficient to determine whether the VZDs merge at these intersections, or if one cross-cuts and potentially offsets the other. Depending on their form between the Thebe and HEX03A 245 datasets, where tectonic faulting inhibits their imaging on the intervening 2D seismic lines, VZDs S-Y may also intersect or connect (Figs 4 and 7). Along most (94%) of the mapped VZDs, minor but abrupt changes in strike allow us to sub-divide them into numerous connected segments (Figs 4,5,and 9c). Across the mapped VZDs, we recognise 280 discrete segments (e.g., Dyke H.1 comprises 26 segments), which have strikes (s) between 350° and 044°, and lengths (l) of 0.4±0.05 to 33.1±0.05 km (Figs 4 and 10b; Supplementary Table S3). Both L and l display a relatively good-fit with log-normal and 250 negative exponential distributions, and poorer fits to normal and power-law distributions (Fig. 10c). within the Thebe and HEX03A 3D seismic surveys, which lie in the western part of the study area, occur at ~3±0.3 s TWT (~2.9 km) (e.g., Fig. 7). Regardless of their precise depth, VZD upper tips across the study area are consistently located ≳1 s TWT beneath the near Base Cretaceous unconformity (e.g., Figs 6 and 7). The expression of all VZDs, at some point along their length, continues below ~5 s TWT (~7 km), where they either appear to terminate or extend beneath the survey limit 260 (e.g., Figs 6 and 7). Although we cannot determine whether the observed lower tips of the VZDs truly mark the base of the structure they correspond to, our data suggests VZD heights are typically ≳1.5 s TWT (≳3.5 km) and potentially ≳3 s TWT (≳9 km) in places (e.g., Figs 6 and 7). Only on a few seismic sections, where data quality is high, do we observe undisturbed reflections directly beneath a VZD, allowing us to constrain its height (e.g., Fig. 6c). For example, the depth to the base of VZD E appears to decrease northwards from ≳5.8±0.55 s TWT to ~4.4±0.45 s TWT (~8.5-5.6 km) (e.g., Figs 6a-c). 265 The width (w) of the VZDs ranges from 68±50 m to 335±50 m ( Fig. 11b; Supplementary Table S4). In places, w could not be confidently measured because other structures (e.g., tectonic faults) locally inhibit VZD imaging. We note w varies between different 3D seismic datasets, each of which had different acquisition geometries, processing histories, and data quality (Fig. 11b). Regardless of these relatively short-wavelength changes in w, there is an apparent overall reduction in w northwards marked by a weakly negative trend-line for the combined dataset (Fig. 11b). 270 Spacing (h) between individual VZDs is variable across the measured transects but broadly increases northwards and is best either described by a log-normal or negative-exponential (Figs 11c and d; Supplementary Table S4). For example, h between VZDs D-E and G-H increases northwards from ~2.77±0.05 km to 4.90±0.05 km and ~6.17±0.05 km to 11.60±0.05 km, respectively (Fig. 11c). A prominent exception to this spatial trend in h, is the northwards reduction in h between VZDs C-D from ~6.80±0.05 to 3.29±0.05 km (Fig. 11C). For part of the lengths of VZDs B-D and B-E, h also 275 decreases northwards, although this is a function of the different orientation of VZD B relative to the other two VZDs (Figs 4 and 11c). Between physically unconnected VZD sections (e.g., G.1-G.6), h is ≲2.01 km, with a minimum of ~0.31±0.05 km ( Fig. 11c). Superimposed onto the large-scale variations in h are localised increases in h (Fig. 11c). The boundaries of these localised increases in h typically coincide with zones where physically unconnected VZD parts terminate, or where VZDs contain a short segment with a markedly different trend to its neighbouring segments (Figs 4 and 11c). There is a good-fit 280 between h and log-normal and negative exponential distributions, but the fit of h to normal and power-law distributions is poorer (Fig. 11d).

Structures associated with VZDs
Overlying and parallel to most VZDs are either one or two, large (up to ~170 km long) normal fault systems, which dip towards and typically converge on the uppers tips of the VZDs (e.g., Figs 6, 7, and 12). In plan-view, these broadly linear 285 normal fault systems extend along much of the length of the underlying VZD ( Sub-circular depressions occur within the graben and half-graben overlying the VZDs (Figs 12 and 13). These depressions are located at the near Base Cretaceous unconformity or at slightly deeper stratigraphic levels within the Dingo Claystone (e.g., Figs 6a, e, 7b, and 13). The depressions are up to ~0.5 km wide, ≲50 ms TWT (≲80 m) deep and infilled by overlying strata (e.g., Figs 6a, e, 7b, 12, and 13). Sub-vertical pipes, within which seismic reflections are displaced downwards relative to their regional trend, underlie each depression (e.g., Figs 6a, e, 7b, and 13). These pipes extend down 300 to the underlying VZD tip or terminate within the Mungaroo Formation above the corresponding VZD (e.g., Figs 6a, e, 7b, and 13).
We discount fluid escape as an origin for our VZDs because these events produce laterally restricted, pipe-like 315 conduits that are geometrically very different to the elongate planar features we observe here ( Fig. 4) (e.g., Cartwright and Santamarina, 2015;Jamtveit et al., 2004;Moss and Cartwright, 2010). We also demonstrate that fluvial channels and linear structures within the Mungaroo Formation are not vertically or laterally offset by cross-cutting VZDs (e.g., Figs 5 and 8), https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License.
indicating there is no evidence for strike-or dip-slip motion across the latter (cf. Harding, 1985). Plate reconstructions for the time of break-up between Greater India and Australia in the Late Jurassic-to-Early Cretaceous, informed by the orientation 320 of tectonic normal faults, seafloor spreading anomalies, and the Cape Range Fracture Zone, further suggest rifting was margin-parallel and thus unlikely to involve significant ~N-trending, strike-slip faulting (e.g., Heine and Müller, 2005). We therefore consider it unlikely that the VZDs are faults.
We interpret the VZDs as igneous dykes because: (i) their seismic expression appears similar to dykes in other real and synthetic seismic datasets (cf. Figs 2, 6, and 7) (e.g., Wall et al., 2010;Ardakani et al., 2017;Kirton and Donato, 1985;325 Holford et al., 2017;Minakov et al., 2018;Plazibat et al., 2019;Eide et al., 2018); and (ii) the geometry of individual VZDs, as well as that of the array they comprise, are akin to the morphology of dyke swarms exposed at Earth's surface (cf. Figs   1A-B and 4) (e.g., Jowitt et al., 2014;Ernst et al., 2001;Halls, 1982). The ~48 m thick basalt interval intersected by the Chester-1 ST1 well, which occurs within VZD H.1, may further support our interpretation that the VZDs correspond to igneous dykes ( Fig. 9). However, to attribute the recovered basalt cuttings to a dyke, we first need to assess whether the well 330 could instead have penetrated a lava flow or sill. Based on an interval velocity of ~4.7±0.5 km s -1 and a dominant frequency of ~20 Hz around the intersected basalt, we calculate that the limits of separability and visibility are locally ~59±6 m and ~8±1 m, respectively. Given these limits of separability and visibility, coupled with the higher density and seismic velocity of the basalt compared to the surrounding sedimentary rocks (Fig. 9d), a ~48 m thick lava flow or sill should be seismically expressed as a high-amplitude, positive polarity, tuned reflection package (e.g., Eide et al., 2018;Rabbel et al., 2018). Yet 335 the intra-Mungaroo seismic reflection coincident with the basalt in Chester-1 ST1 has a negative polarity and is of moderate amplitude (Figs 9A and B). These observations suggest the basalt intersected by Chester-1 ST1 does not come from a lava flow or sill, but instead supports our interpretation that the coincident VZD H.1, and likely other VZDs, are igneous dykes.
We concede that limitations in seismic and well data resolution mean we still cannot determine whether individual VZDs correspond to a single dyke or multiple, closely spaced dykes (e.g., Fig. 5). 340 Our interpretation that the VZDs correspond to igneous dykes raises the question as to whether the observed overlying normal fault systems and pipes, which converge on the inferred dykes, were genetically related to magmatism (e.g., Figs 4, 6, 7, 12, and 13). For example, normal fault systems and sub-circular depressions similar to those we describe have been observed above dykes on Earth, other planetary bodies, and in physical and numerical models (e.g., Wyrick et al., 2004;Wilson and Head, 2002;Trippanera et al., 2015b;Trippanera et al., 2015a;Rubin and Pollard, 1988;Pollard et al., 345 1983;Rubin, 1992;Hardy, 2016;Wyrick and Smart, 2009;Okubo and Martel, 1998). Numerical and analytical models suggest normal faulting above intruding and widening dykes is driven by the concentration of tensile stress at the dykes upper tip and at the contemporaneous surface (e.g., Rubin and Pollard, 1988;Pollard et al., 1983;Rubin, 1992). Shear failure within this local dyke-induced stress field produces graben-or half graben-bounding, dyke-parallel normal faults that dip towards and converge on the dykes upper tip (e.g., Rubin and Pollard, 1988;Pollard et al., 1983;Rubin, 1992;Trippanera et 350 al., 2015b); these faults are termed 'dyke-induced normal faults'. Dyke intrusion and widening can also locally produce cavities through the accumulation and release of magmatic volatiles at its upper tip, or the heating and escape of pore fluids https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License.

Timing of dyke emplacement
Radiometric dates are unavailable to constrain the emplacement age of the studied dykes, so we have to apply seismicstratigraphic techniques. Each dyke intrudes and terminates within the Mungaroo Formation, indicating their emplacement occurred during or after the Triassic (e.g., Figs 6, and 7). The dykes also cross-cut and thus post-date sills intruded within the Triassic Mungaroo Formation (e.g., Figs 6b, c, and f). Although we have no constraints on the age of these sills cross-cut by 365 the dykes, it is likely they were emplaced during a regional phase of Late Jurassic-to-Early Cretaceous magmatism (e.g., Symonds et al., 1998;Rohrman, 2013;Magee et al., 2017;Magee et al., 2013b;Magee et al., 2013a). Onlap of overlying strata onto intrusion-induced forced folds suggest sill emplacement elsewhere in the North Carnarvon Basin may have begun in the Kimmeridgian (Magee et al., 2017;Magee et al., 2013a).
The near Base Cretaceous unconformity (~148 Ma) is the youngest stratigraphic horizon deformed by most of the 370 interpreted dyke-induced normal fault systems and pit craters (e.g., Figs 6 and 7). Where dyke-induced normal fault systems and pit craters are observed elsewhere on Earth or other planetary bodies, they deform the surface, contemporaneous with dyke intrusion (e.g., Wyrick et al., 2004;Wilson and Head, 2002;Trippanera et al., 2015b;Trippanera et al., 2015a;Rubin and Pollard, 1988;Pollard et al., 1983;Rubin, 1992;Hardy, 2016;Wyrick and Smart, 2009;Okubo and Martel, 1998). Our seismic-stratigraphic observations therefore suggest the near Base Cretaceous unconformity (~148 Ma) likely marked the 375 palaeosurface during dyking, indicating emplacement principally occurred during or after its development, but ceased before the overlying Barrow Group was deposited. Some pit craters terminate within rather than at the top of the Dingo Claystone (e.g., Fig. 13), suggesting dyking may have initiated in the Late Jurassic before the near Base Cretaceous unconformity formed at ~148 Myr. The apparent extension of some dyke-induced normal faults into the ~146.7-138.2 Ma Barrow Group, which is located above the near Base Cretaceous unconformity, may be indicative of renewed, post-Valanginian dyking 380 (Figs 6d, e, and 7a-c). An alternative suggestion is that the upward extension of the dyke-induced normal faults into the Barrow Group simply reflects fault reactivation and/or dip-linkage during later polygonal fault formation (i.e. these fault extensions are unrelated to dyking). Such reactivation or dip-linkage of the dyke-induced normal faults is supported by the: (i) reduced dip of many dyke-induced faults segments above the near Base Cretaceous unconformity (e.g., Figs 6d, e, and 7a-https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License. b); and (ii) similar extension of some tectonic normal faults above the near Base Cretaceous unconformity, occasionally to 385 just below the seabed. Overall, we propose all dykes were likely intruded during a short period in the latest Jurassic, probably during the Tithonian (~152-147 Ma), before the onset of Barrow Group deposition at ~146.7 Ma (Reeve et al., 2016); we name this newly discovered suite of igneous dykes the Exmouth Dyke Swarm.

Dyke swarm structure
To understand the kinematics and mechanics governing dyke swarm emplacement, we typically rely on measuring the 390 geometrical properties (e.g., length, width, and spacing) of dykes exposed at the Earth's surface (e.g., Paquet et al., 2007;Jolly and Sanderson, 1995;Gudmundsson, 1983). A potential problem with these analyses is that we can only measure the surface, principally 2D expression of dykes and dyke swarms, which may not equal their true 3D geometry. For example, seismic reflection data from offshore southern Norway reveal the width of an imaged dyke swarm increases with depth, implying the dimensions of dyke swarms measured at the surface depend partly on erosion level and may therefore not 395 capture the true swarm geometry (Phillips et al., 2018). Seismic reflection data thus provide a unique opportunity to examine and quantify the 3D structure of a dyke swarm independent of the potential bias introduced by the processes (e.g., erosion) controlling how dyke swarms intersect the surface. Here, we specifically discuss how our measurements of dyke length, thickness, and spacing compare to predicted distributions of these geometrical properties derived from surface-and physical, numerical, and analytical modelling-based studies. 400

Dyke length
Lengthening of fractures is commonly facilitated by linkage between individual segments (e.g., Mège and Korme, 2004;Schultz, 2000;Cladouhos and Marrett, 1996;Gudmundsson, 1987). The evolution of a fracture population can be unravelled from its length distribution if we can ascertain whether linked or closely spaced fractures should be treated as one or several structures (e.g., Schultz, 2000;Mège and Korme, 2004); i.e. does the length-frequency distribution of a fracture population 405 change through time in response to linkage modifying the behaviour of the system, or is it scale invariant? Dykes are magma-filled fractures and can broadly be considered to intrude instantaneously and independently formed fractures (i.e. they do not interact), implying the length-frequency distribution of a dyke swarm should preserve the initial configuration of the fracture population (Mège and Korme, 2004). Comparing data from fracture and dyke populations reveal their lengthfrequency distributions are both broadly power-law, suggesting mechanical linkage of fractures does not modify system 410 behaviour (e.g., Mège and Korme, 2004;Schultz, 2000;Cladouhos and Marrett, 1996;Gudmundsson, 1987;Paquet et al., 2007). Here we use our data, assuming the dykes are Mode I fractures, to examine whether: (i) measurement of dyke-surface intersections introduces bias to length-frequency distributions; and (ii) dyke segmentation, which may be indicative of noninstantaneous and non-independent fracture growth, also display a power-law length-frequency distribution (cf. Mège and Korme, 2004). 415 https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License.
Cumulative length-frequency plots for all measured dyke lengths (L), which comprise connected and/or closely spaced but physically unconnected segments, initially appear to fit a log-normal or negative exponential, rather than a power-law distribution (Fig. 10c) (cf. Paquet et al., 2007;Mège and Korme, 2004). Dyke segment length (l) data display similar log-normal and negative exponential distribution characteristics (Fig. 10c). However, power-law distributions can be fit to L values between 20-160 km and l values of 5-20 km (Fig. 10c). The population exponents (C) for the L and l datasets 420 are 1.29 and 2.85, respectively, consistent with values derived from the analysis of other fracture and dyke populations (Fig.   10c) (see Mège and Korme, 2004 and references therein). The observed departure of our measured L and l values from a power-law distribution at small and large length-scales could indicate bias in the data. For example, restrictions in dyke imaging and 2D seismic line spacing may mean: (i) the dykes are likely longer than mapped; (ii) some dykes (e.g., VZDs X and Z) may be connected along-strike, thereby increasing their length (Figs 4b-d); and (iii) small dykes and/or dyke 425 segments are difficult to recognise or may not be imaged because they occur between 2D seismic lines outside of areas imaged by the 3D surveys. We contend that our data could thus be considered consistent with previous studies in describing dyke length distributions as power-law, indicating processes controlling dyke length (e.g., segmentation) are scale invariant (Mège and Korme, 2004). Furthermore, our results suggest the free-surface intersection of fractures or dykes is, at least typically, representative of a population's length distribution. 430

Dyke thickness
The thickness of a dyke, or cumulative thickness of a dyke swarm, influences a variety of processes, including eruption rates and crustal extension (e.g., Krumbholz et al., 2014). Statistical analyses of dyke thickness distributions derived from surfacebased measurements inform dynamic models of dyke emplacement, shedding light on the processes controlling dyke thickness (e.g., Krumbholz et al., 2014;Jolly and Sanderson, 1995;Klausen, 2006;Klausen, 2004). Resolving the 3D 435 structure of dyke swarms in seismic reflection data provides a potential opportunity to examine both lateral and vertical variations in dyke thickness distribution. However, synthetic seismic forward models suggest the width of VZDs corresponding to sub-vertical dykes is greater than the true thickness of the dyke itself (Eide et al., 2018). Furthermore, because VZD width is partly controlled by the acquisition and processing properties of the data in which they are imaged in (e.g., frequency; Eide et al., 2018), evidenced by the marked differences in VZD width between different seismic surveys 440 (Fig. 11c), it is difficult to determine how VZD width and true dyke thickness are related. We show VZD widths measured across multiple 3D seismic surveys gradually decrease northwards (Fig. 11c). Because the northwards decrease in VZD width is consistent across multiple seismic surveys, which each have different acquisition and processing parameters, we suggest this trend marks a similar northwards decrease in true dyke thickness (Fig. 11c). From the Chester-1 ST1 well, which likely intersects a 48 m long section of basalt, we calculate the dyke has a true thickness of ~18 m, assuming its orientation is 445 parallel to that of the ~130±50 m wide VZD it relates to (Fig. 14). This well data confirms synthetic seismic forward model predictions that dyke-related VZD width is much greater than true dyke thickness (Eide et al., 2018). Further work in https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License.
understanding how dykes are expressed in seismic reflection data is required before these data can be used to quantify dyke thickness distributions.

Dyke spacing 450
Plan-view sections through dyke swarms reveal individual dykes are typically regularly spaced, with the spacing (h) of radiating swarms increasing away from their focal area (e.g., Bunger et al., 2013;Jolly and Sanderson, 1995;Ernst et al., 1995). Identifying controls on h is fundamental to understanding why dykes occur in swarms and, thus, how they interact with and/or drive crustal extension on Earth and other planetary bodies (Bunger et al., 2013). Analytical predictions suggest first-generation, laterally propagating dykes will have energetically optimal spacings that are related to dyke height (H) and 455 magma source conditions (Bunger et al., 2013). For dykes emanating from a constant pressure magma source (i.e. an infinitely large, compressible reservoir), h/H is expected to be ≈ 1, whilst those from a constant influx magma source (i.e. a small, incompressible reservoir) will have either a h/H of ≈ 2.5 or ≈ 0.3 (Bunger et al., 2013). Constraining the relative age of dykes is critical to testing these analytical predictions because second-generation or younger may preferentially intrude between first-generation dykes, thereby reducing the apparent spacing (Bunger et al., 2013). 460 Dyke spacing within the Exmouth Dyke Swarm ranges from ~22.4±0.05 km to 0.3±0.05 km, with a geometric mean of ~4.1 km, and broadly increases northwards (Figs 4 and 11d). This northward increase in h, coupled with apparent northwards reductions in dyke thickness and abundance, implies extension accommodated by the Exmouth Dyke Swarm similarly decreased northwards (Figs 4 and 11d). To test analytical predictions using our measured h values, it is first important to recognise key limitations in our dataset: (i) not all dykes within the swarm may be imaged by the seismic 465 reflection data, suggesting our h measurements are likely only maximum values; (ii) H is difficult to quantify because a reduction in data quality with depth likely means we cannot accurately pick the lower tips of dykes, some of which may extend beneath the seismic surveys (e.g., Figs 6 and 7); and (iii) it is challenging to ascertain whether all dykes were emplaced simultaneously or not during the latest Jurassic dyking event. Because the dykes are typically >1.5 s TWT tall (e.g., Figs 6 and 7), we use extrapolated checkshot data to estimate the average H is at least ~3.5 km ( Supplementary Fig.  470 S3). As a maximum estimate for average H, we consider the dykes could extend upwards from the base of the crust, which across the Exmouth Plateau is likely ~20-28 km beneath the present day seabed (e.g., Reeve et al., 2016;Tindale et al., 1998;Stagg et al., 2004;Stagg and Colwell, 1994;Mutter and Larson, 1989). Given the upper dyke tips broadly occur at ~3.7±0.37 s TWT, equivalent to a depth of ~4.1 km, we therefore suggest the maximum average H could be up to ~24 km.
Assuming dyking was instantaneous and using the geometric mean for h (~4.1 km), we calculate h/H ≈ 1.17-0.17. 475 The calculated h/H values of 1.17-0.17 are broadly consistent with and cannot be used to discriminate between the constant pressure (h/H ≈ 1) and constant influx (h/H ≈ 0.3) end-member source conditions (Bunger et al., 2013). However, if seismically unresolved dykes are present in the study area, we may expect h to be less than that measured and thus more consistent with h/H ≈ 0.3, implying the dykes were fed from a constant influx magma source (Bunger et al., 2013).
we would expect h≳4.1 km for the first-generation dykes; this would imply the original maximum h/H ratio could be ≈ 1.
Potential evidence for incremental emplacement of the Exmouth Dyke Swarm includes: (i) the relatively good fit of h to a negative-exponential distribution (Fig. 11d), which suggests h is random and likely results from incorporation of different dyke sets into the data; and (ii) the observation that some pit craters occur within (rather than at the top of) the Dingo Claystone (e.g., above Dyke F; Fig. 13), suggesting their associated dykes were emplaced before the formation of the near 485 Base Cretaceous unconformity (~148 Ma). For example, if we hypothetically consider VZDs C, F, H, and I were emplaced first, their geometric mean h of 12.4 km implies h/H ≈ 3.54-0.52, which again could be considered consistent with a constant pressure (h/H ≈ 1) or constant influx (h/H ≈ 2.5) source (Bunger et al., 2013). Mapping the occurrence and distribution of pit craters formed before the near Base Cretaceous unconformity may allow us to identify first-generation dykes and thereby constrain dyke source conditions. 490

Emplacement of the Exmouth Dyke Swarm
We mapped the Exmouth Dyke Swarm, as well as associated dyke-induced normal faults and pit craters, across a ~40,000 km2 part of the North Carnarvon Basin (Figs 4-7 and 12). Long, linear graben, containing sub-circular depressions, similar to the dyke-induced normal faults and pit craters we identified, occur at the near Base Cretaceous unconformity elsewhere in the North Carnarvon Basin (e.g., Fig. 15) (Velayatham et al., 2019;Velayatham et al., 2018). The formation of some of these 495 other depressions has been linked to fluid escape following faulting of overpressure strata, and not dyking (Velayatham et al., 2018). However, their geometrical similarity to and occurrence at the same structural level as the dyke-induced normal fault systems and pit craters described here suggests they could be the palaeosurface expression of the Exmouth Dyke Swarm (cf. Figs 12 and 15) (see also Velayatham et al., 2019). This potential distribution of dykes (except for VZD K), dyke-induced normal fault systems, and pit craters across the North Carnarvon Basin appears to describe a giant radial dyke 500 swarm (cf. Figs 1c and 15c) (cf. Ernst, 2014;Ernst et al., 2001;Ernst et al., 1995;Halls and Fahrig, 1987). Projecting the inferred dykes to a common focal area, which is located on the Cuvier Margin, suggests the Exmouth Dyke Swarm could be >500 km long and distributed around a ~039° (perhaps up to ~054°) arc (Fig. 15c). To unravel the origin of the Exmouth Dyke Swarm, we first discuss evidence for magma propagation direction and syn-emplacement stress conditions.

Dyke propagation direction 505
The radiating form of the Exmouth Dyke Swarm suggests individual dykes may have been sourced and thus flowed laterally northwards from the northern sector of the Cuvier Margin (Fig. 15c) (see also Velayatham et al., 2019). Lateral propagation of the dykes to the north is supported by the: (i) maintenance of dyke upper tip depths (Figs 6, 7, and 11a), consistent with the expectation that horizontally emplaced dykes have fixed upper and lower tip positions (e.g., Townsend et al., 2017); (ii) subtle northwards decrease in VZD width (Fig. 11b), which we suggest reflects thinning of dykes, perhaps towards their 510 laterally propagating tip (e.g., Healy et al., 2018); and (iii) minor but abrupt changes in the strike of connected dyke https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License. segments (Figs 4 and 5), which are reminiscent of the kinked geometry attained by the Bárðarbunga-Holuhraun dyke during its incremental, lateral propagation (Woods et al., 2019;Sigmundsson et al., 2015).

Palaeostress conditions during dyke emplacement
The orientation and structure of dykes and dyke swarms is commonly used to reconstruct syn-emplacement stress and 515 magma conditions (e.g., Odé, 1957;Hou et al., 2010;Grosfils and Head, 1994;Lahiri et al., 2019;Jolly and Sanderson, 1997;Jolly and Sanderson, 1995). Deriving these overarching controls on dyke emplacement assumes that dykes preferentially develop orthogonal to σ3 within the σ1-σ2 plane (e.g., Anderson, 1951). Although the orientation of dykes and dyke segments studied here is variable, they are broadly N-to NE-trending and sub-vertical (~80-90°), suggesting an average syn-emplacement σ3 currently oriented 100/00° (Fig. 16). Mutually orthogonal to the calculated σ3 on a lower-520 hemisphere, equal area stereographic projection are two axes, at 010/00° and 280/90° respectively, which can be ascribed to σ1 or σ2 depending on their proximity to the cluster of measured dyke poles (e.g., Jolly and Sanderson, 1997; Lahiri et al., 2019). Specifically, the angle measured along the σ1-σ3 plane between the cluster of dykes and σ1 (i.e. θ2) will be greater than that measured along the σ2-σ3 plane between the data and σ2 (i.e. θ1; Fig. 16) (e.g., Jolly and Sanderson, 1997; Lahiri et al., 2019). Our data thus suggests that during dyking, the overarching stress field in the study area was extensional with a 525 vertical σ1 (000/90°) and horizontal, N-trending σ2 (010/00°) (Fig. 16). The syn-emplacement, ~W-trending, horizontal σ3 axis we define is comparable to suggested W-to NW-trending extension directions, estimated from tectonic fault orientations and seafloor spreading patterns, for Late Jurassic-to-Early Cretaceous rifting and break-up offshore NW Australia (e.g., Hopper et al., 1992;Heine and Müller, 2005;Driscoll and Karner, 1998). Where NW-trending dykes may dominate to the west of the study area (Fig. 15c) (Velayatham et al., 2018), we anticipate the horizontal principal stress axes 530 (σ2 and σ3) were oriented NW-SE and NE-SW, respectively, whilst σ1 remained vertical.

Tectono-magmatic setting and source of the Exmouth Dyke Swarm
Magmatism across the North Carnarvon Basin has been attributed to decompression melting during rifting (Karner and Driscoll, 1999), coupled rifting and small-scale convective partial melting (e.g., Hopper et al., 1992;Mutter et al., 1988;Mihut and Müller, 1998), and/or mantle plume activity (e.g., Rohrman, 2015Rohrman, , 2013Müller et al., 2002;Mihut and Müller, 535 1998). We show emplacement of the Exmouth Dyke Swarm occurred during the latest Jurassic (~152-147 Ma), after intrusion of extensive sill-complexes (e.g. , Figs 6 and 7). Individual dykes likely propagated laterally away from a source focal area, which we infer was located on the Cuvier Margin, SSE of the study area (Fig. 15c). Dyking and earlier sill emplacement thus predated the main phase of igneous activity recorded across the North Carnarvon Basin, which was associated with formation of the ~136-130 Ma continent-ocean transition zones bordering the Gascoyne and Cuvier margins, 540 and ultimately continental break-up in the Hauterivian (e.g., Rey et al., 2008;Symonds et al., 1998;Mihut and Müller, 1998;Reeve et al., 2019;Direen et al., 2007;Robb et al., 2005). Seismic reflection data also reveal there was little upper crustal normal faulting or rifting across the Exmouth Plateau in the Late Jurassic, immediately prior to and during dyking (e.g., https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License. Driscoll and Karner, 1998;Bilal et al., 2018). These age relationships suggest the Exmouth Dyke Swarm and earlier sills were likely not associated with rift-related melting, which appears to have initiated in the Early Cretaceous (cf. Mihut and 545 Müller, 1998;Karner and Driscoll, 1999;Hopper et al., 1992;Mutter et al., 1988). Instead, the large extent and radial disposition of the Exmouth Dyke Swarm suggests it may have been sourced from either a regional, thermal mantle anomaly (e.g., a plume or small-scale convection cell) or a large volcanic system (e.g., Ernst and Buchan, 1997;Ernst et al., 1995;Odé, 1957;Speight et al., 1982).
Any process invoked to explain the origin of a thermal mantle anomaly in the Late Jurassic, and potentially into the 550 Early Cretaceous, needs to account for (e.g., Rohrman, 2015Rohrman, , 2013Müller et al., 2002;Mihut and Müller, 1998;Hopper et al., 1992;Mutter et al., 1988): (i) the latest Jurassic distribution of magmatism across the Gascoyne and Cuvier margins; and (ii) denudation patterns and formation of contemporaneous regional unconformities (e.g., the near Base Cretaceous unconformity). Two possible mantle plume sites on the Cuvier Margin have previously been proposed, with one located on the Bernier Platform, initiating at ~136 Ma, and the other active on the conjugate to the Cuvier Margin near the current Cape 555 Range Fracture Zone between ~165-136 Ma (e.g., Fig. 17a) (cf. Rohrman, 2015;Müller et al., 2002;Mihut and Müller, 1998). Mantle plume activity has previously been discounted as a viable source for Late Jurassic-to-Early Cretaceous magmatism because no clear hotspot tracks have been identified (e.g., Müller et al., 2002), although Rohrman (2015) argued the Quokka Rise and Zenith Plateau are part of such a track (Fig. 17a). An alternative interpretation to a mantle plume source is that melting reflects small-scale mantle convection instigated by juxtaposition of thick and thin lithosphere across a 560 transform margin (e.g., the Cape Range Fracture Zone) (e.g., Müller et al., 2002;Mutter et al., 1988). Because the formation of transform margins along the NW Australian Shelf occurred during break-up of Greater India and Australia in the Early Cretaceous (~136-130 Ma), coincident with the age of the proposed Bernier Platform mantle plume, it seems unlikely these processes could have generated the latest Jurassic Exmouth Dyke Swarm (cf. Müller et al., 2002;Mihut and Müller, 1998).
The interpreted age and distribution of the Exmouth Dyke Swarm thus fits best with the mantle plume model proposed by 565 Rohrman (2015).
Within the framework of the mantle plume model proposed by Rohrman (2015), melting is expected to have initiated ~165 Myr ago, leading to emplacement of a mafic-to-ultramafic, high-velocity magmatic body near the Moho and formation of the Callovian unconformity during associated uplift (i.e. vertical σ1; Fig. 18a). This high-velocity magmatic body likely fed the Late Jurassic sill-complex prior to emplacement of the Exmouth Dyke Swarm (Fig. 18a) (e.g., Rohrman, 570 2013;Magee et al., 2017;Magee et al., 2013a;Symonds et al., 1998). We suggest that emplacement of this sill-complex occurred as plume activity waned and uplift ceased, causing the regional stress to relax such that the vertical principal stress axis became σ3 and basin subsidence initiated (Fig. 18b); this change in stress orientation could explain why the ascent of buoyant magma from the high velocity body formed a sill-complex rather than a vertical dyke swarm. Layering in the sedimentary basins may also have favoured sill emplacement (Fig. 18b) (see Magee et al., 2016b and references therein). 575 The apparent transition from sill-complex formation to intrusion of the Exmouth Dyke Swarm in the latest Jurassic marks an abrupt change in emplacement conditions. To generate the Exmouth Dyke Swarm, which broadly coincided with a phase of https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License. uplift and denudation (i.e. formation of the Base Cretaceous unconformity), we show σ1 had become vertical and σ3 was circumferential to the swarms focal area (Figs 16a, 18c, and d). We suggest these conditions, which favoured dyking rather than sill-complex emplacement, could have been instigated by a renewed influx of plume material, with the swarm fed 580 either: (i) directly from a thermal mantle anomaly (Fig. 18c); or (ii) via a large intrusive centre located at the southern boundary of the Exmouth Sub-basin, which manifests as a sub-circular (~20 km diameter), positive magnetic anomaly and a zone of disturbance in seismic reflection data (e.g., Figs 17 and 18d) (Müller et al., 2002). Cessation of plume activity immediately after dyking, following removal or reduction of the thermal anomaly, may explain the rapid subsidence (i.e. <0.24 mm yr -1 ) required to accommodate the Late Jurassic-to-Early Cretaceous Barrow Group (cf. Reeve et al., 2016). 585 Overall, our data seemingly support the presence of a mantle plume offshore NW Australia during the Late Jurassic-to-Early Cretaceous (e.g., Rohrman, 2015Rohrman, , 2013Müller et al., 2002). However, it remains uncertain whether igneous activity coincident with Hauterivian break-up was also related to the presence of a mantle plume or not.

Implications and future studies
Giant dyke swarms are recognised worldwide onshore (e.g., Ernst and Youbi, 2017;Ernst, 2014;Bryan and Ernst, 2008;590 Coffin and Eldholm, 2005;Coffin and Eldholm, 1994;Bryan et al., 2010;Hou et al., 2010;Halls and Fahrig, 1987;Halls, 1982;Ernst and Baragar, 1992). Projection of these onshore dyke swarms and the known importance of dyking to break-up and formation of magma-rich margins suggests dyke swarms should also be prevalent on offshore continental shelves (see Magee et al., 2019 and references therein). Our work extends a growing consensus that vertical dykes can be recognised in seismic reflection data imaging continental margins (e.g., Wall et al., 2010;Bosworth et al., 2015;Malehmir et al., 2018;595 Holford et al., 2017;Kirton and Donato, 1985;Ardakani et al., 2017;Jaunich, 1983;Plazibat et al., 2019). Key criteria for defining vertical dykes in seismic reflection data include: (i) identification of thin, long, tall, typically sub-vertical zones of disturbance within otherwise sub-parallel reflections defining the host rock (e.g., Figs 6 and 7) (e.g., Eide et al., 2018;Minakov et al., 2018;Wall et al., 2010); (ii) lack of lateral or vertical offset of host rock strata, best revealed by mapping piercing points (e.g., fluvial channels, pre-existing structures) across inferred dyke-like features (e.g., Figs 5 and 8), which 600 suggests the features are not strike-slip or steeply dipping normal faults; and (iii) potential association with overlying pit craters or dyke-induced normal faults, which are likely easier to resolve and map in seismic reflection data compared to dykes (e.g., Figs 6, 7, 12 and 13). By increasing our collective awareness of how these criteria can be used to identify dykes in seismic reflection data, we expect more dyke swarms will be revealed across continental margins worldwide. Recognition of dyke swarms within seismic reflection data will help us produce better physical models of the subsurface, aiding our 605 understanding of a margins thermal history, and fluid and/or gas plumbing systems of sedimentary basins.
We also demonstrate that mapping dykes, dyke-induced normal faults, and pit craters across vast areas using seismic reflection data provides unprecedented opportunities to resolve and quantify their natural structure in 3D (e.g., Figs 4-13). Future work should focus on: (i) unravelling the geophysical expression of dykes, such that additional and more accurate quantitative data (e.g., dyke thickness) can be recovered; (ii) deciphering the kinematic history of dyke-induced 610 https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License.
normal faults, which we may expect should relate to and thus inform dyke structure and emplacement dynamics; (iii) quantifying the geometrical relationship between pit craters and the dyke intrusions driving their formation; and (iv) determining whether dyke-induced normal faults and pit craters can be used to constrain the temporal evolution of a dyke swarm. These four initiatives will provide new insights into and allow us to test hypotheses concerning the 3D structure and growth of dyke swarms, and their associated structures. We envisage that these findings will improve how we can invert the 615 surface expression of active or ancient dyke swarms, i.e. dyke-induced normal faults and pit craters exposed at the surface of Earth or other planetary bodies, to recover more information on their otherwise inaccessible subsurface structure and the processes that formed them.

Conclusions
Dyke swarms are ubiquitous on Earth and other planetary bodies. Yet we know little of the 3D structure of dyke swarms 620 because the pseudo-2D nature of planetary surfaces means we can typically only access their plan-view morphology, and then only at the given erosion level. Here we use a suite of seismic reflection datasets from the Exmouth Plateau offshore NW Australia, to map 26, latest Jurassic (~152-147 Ma) dykes in 3D across ~40,000 km 2 ; we name this the Exmouth Dyke Swarm. The mapped dykes correspond to ~N-to NE-trending, vertical zones of disturbance within the seismic reflection data that are can be up to 171 km long, ≲355 m wide, likely ≳9 km high, and can be sub-divided into smaller segments with 625 subtly different orientations. Directly above the dykes are a series of graben-bounding normal fault systems, which dip towards and converge upper dyke tips, and sub-vertical pipe-like features; we interpret these structures as dyke-induced normal faults and pit craters. Our quantitative analyses reveal dyke length broadly follows a power-law distribution consistent with previous studies, whilst dyke spacing conforms to a negative-exponential distribution, which we attribute to sampling of different dyke generations. Across the study area, dyke orientations are consistent with an ENE-trending, 630 horizontal and a vertical minimum and maximum principal stress axes, respectively. However, recognition of possible dykeinduced normal faults and pit craters elsewhere on the Exmouth Plateau suggest dykes are distributed radially across a 39° arc, implying the minimum principal stress axis was circumferential, centred on the Cuvier Margin to the south. This focal area on the Cuvier Margin likely marks the dyke swarm source, which is consistent with evidence the dykes propagated laterally northwards. Overall, we suggest emplacement of the Exmouth Dyke Swarm related to renewed activity of a mantle 635 plume located on the Cuvier Margin between ~165-136 Ma. Our work demonstrates seismic reflection data can be used to identify vertical dykes across vast areas on continental margins, whilst providing unprecedented into the 3D structure of these natural systems. By defining a series of criteria that can be used to interpret dykes in seismic reflection data, we anticipate future studies will: (i) recognise dyke swarms across continental margins worldwide, providing new insights into basin evolution (e.g., thermal histories) and controls on fluid flow; (ii) provide more robust constraints on dyke swarm 640 geometry, allowing previous models and hypotheses of their 3D structure to be tested; (iii) reveal how dyke-induced normal faults and pit craters are kinematically linked to dyking; and (iv) demonstrate how dyke swarms may be expressed at the syn-https://doi.org/10.5194/se-2019-201 Preprint. Discussion started: 13 January 2020 c Author(s) 2020. CC BY 4.0 License. emplacement surface, meaning we can improve inversions of such surficial features observed on Earth and other planetary bodies to better predict underlying dyke structures.

Data availability 645
The seismic reflection and well data used (see text for details) are publically available through Geoscience Australia at http://www.ga.gov.au/nopims. All measurements presented are within, or based on data within, the compiled Table and Supplementary Tables.                     Continent-ocean transition zone (COTZ) Direen et al., 2008) Tectonic element boundary Sill-complexes (modified from Symonds et al., 1998)