Evolution of normal fault displacement and length as continental lithosphere stretches

Continental rifting is accommodated by the development of normal fault networks. Fault growth patterns control their related seismic hazards, and the tectonostratigraphic evolution and resource and CO2 storage potential of rifts. Our understanding of fault evolution is largely derived by observing the final geometry and displacement (D)‐length (L) characteristics of active and inactive fault arrays, and by subsequently inferring their kinematics. We can rarely determine how these geometric properties change through time, and how the growth of individual fault arrays relate to the temporal evolution of their host networks. Here we use 3D seismic reflection and borehole data from the Exmouth Plateau, NW Shelf, Australia to determine the growth of rift‐related, crustal‐scale fault arrays and networks over geological timescales (>106 Ma). The excellent‐quality seismic data allows us to reconstruct the entire Jurassic‐to‐Early Cretaceous fault network over a relatively large area (ca. 1,200 km2). We find that fault trace lengths were established early, within the first ca. 7.2 Myr of rifting, and that along‐strike migration of throw maxima towards the centre of individual fault arrays occurred after ca. 28.5 Myr of rifting. Faults located in stress shadows become inactive and appear under‐displaced relative to adjacent larger faults, onto which strain localises as rifting proceeds. This implies that the scatter frequently observed in D‐L plots can simply reflect fault growth and network maturity. We show that by studying complete rift‐related normal networks, rather than just individual fault arrays, we can better understand how faults grow and more generally how continental lithosphere deforms as it stretches.


| INTRODUCTION
The formation of extensional basins is controlled by the development of normal faults. Large, basin-bounding fault arrays typically grow by the initiation, propagation, interaction and linkage of smaller fault segments (e.g. Anders & Schlische, 1994;Dawers & Anders, 1995;McLeod et al., 2000;Peacock & Sanderson, 1991;Trudgill & Cartwright, 1994). A kinematically linked network of fault segments and arrays comprise a fault network. Fault networks evolve in response to co-and interseismic stress feedbacks between the constituent segments and arrays; this can lead to the temporal progression from numerous, short, low-displacement segments to a few, long, highdisplacement arrays ('strain localisation'; Cowie, 1998;Gawthorpe & Leeder, 2000) (Figure 1).
As strain localisation is often associated with an increase in basin subsidence rate, the temporal and spatial evolution of fault networks are reflected in changes in rift physiography, and the location and rate of generation of accommodation (e.g. McLeod et al., 2000). The way in which fault networks grow therefore strongly influences the size, shape and distribution of sedimentary depocentres and drainage catchments, which in turn control the location of various energy resources, groundwater reservoirs and potential CO 2 storage sites (e.g. Gawthorpe & Leeder, 2000).
Mechanical interaction between adjacent faults profoundly affects the stress-dependent nature of earthquakes in potentially hazardous, seismically active regions (e.g. Allmendinger & Shaw, 2000;Dolan et al., 2007;Harris, 1998;Nicol et al., 2010;Mildon et al., 2019;Stein, 1999). An improved understanding of both, long-term fault displacement patterns from ancient fault networks imaged in seismic reflection data, and with short-term earthquake slip rates typically derived from regions of active extension, would allow us to characterise fault growth and interactions over multiple timescales, from a few seconds up to tens-to-hundreds of millions of years (Allmendinger & Shaw, 2000;Nicol et al., 2006Nicol et al., , 2010Nicol et al., , 2019Scholz & Gupta, 2010). However, our current understanding of fault growth largely stems from detailed studies of individual fault arrays (e.g. Childs et al., 2017;Fossen and Rotevatn, 2016;Jackson & Rotevatn, 2013;Jackson et al., 2017;Morley, 1999Morley, , 2002Nicol et al., 2005;Nicol et al., 2016), rather than rift-scale fault networks, thereby often overlooking the important role fault interactions and stress feedbacks may play in controlling the early stages of fault and rift development and, ultimately, continental breakup.
Historically, fault growth models have been proposed based on field and subsurface observations of the final geometry of ancient faults, often in the absence of kinematic constraints on how these faults evolved. The key F I G U R E 1 Fault network evolution, modified after Gawthorpe and Leeder (2000). Faults grow from an initiation stage (a) with isolated, short (20-10 km) and low-displacement fault segments into a period of interaction and linkage (b), and finally strain localises onto a few long (>10 km), high displacement 'through-going' fault zones (c). Typically, studies investigating fault growth measure the geometry of fault arrays when they have reached (c) geometric relationship that has been used is the relationship between fault length (L) and maximum fault displacement (D): where c represents the scaling D/L ratio and n represents an exponent value that may vary from a value of 1 indicating a linear scaling law (i.e. self-similarity) to 2 (Cowie & Scholz, 1992;Dawers et al., 1993;Marrett & Allmendinger, 1991;Schlische et al., 1996;Walsh & Watterson, 1988;Watterson, 1986). The exponent n remains unresolved due to the significant amount of scatter in global D-L profiles ( Figure 2) (Cartwright et al., 1995;Cladouhos & Marrett, 1996;Cowie & Scholz, 1992;Dawers et al., 1993;Nicol et al., 1996;Schlische et al., 1996Scholz et al., 1993, which may reflect natural geological variations between study areas (e.g. host rock mechanical properties), and uncertainties in data quality and sampling approaches (Gillespie et al., 1993;Kim & Sanderson, 2005;Walsh & Watterson, 1988;Walsh et al., 2017). It is, however, possible that variations in n reflect fundamentally differing styles of fault growth. For example current conceptual models suggest that individual faults may grow by either the: (1) 'propagating model', in which fault displacement and length simultaneously increase, producing a linear growth trajectory (Figure 2; e.g. Cartwright et al., 1995;Walsh & Watterson, 1988;Walsh et al., 2003); or (2) the 'constant length model', in which faults reach their near-full lengths rapidly before significant displacement accumulation, producing an initially shallow gradient trajectory on D-L plots as the fault lengthens, followed by a near-vertical trajectory as the fault accumulates displacement (Figure 2; e.g. Childs et al., 2017;Jackson et al., 2017;Meyer et al., 2002;Nicol et al., 2005;Nicol et al., 2016;Rotevatn et al., 2019;Walsh et al., 2002). The two fault models likely represent end-member styles of fault growth (Childs et al., 2017;Fossen 2016;Nicol et al., 2020), with an intermediary, hybrid growth model stating that faults initially lengthen via propagation and linkage for 20%-30% of the fault's lifespan, before later accumulating displacement for the final D max = cL n , F I G U R E 2 Global displacement-length (D-L) data of normal faults. D-L relationships are presented in log-log space, therefore the line of best fit results in the scaling exponent, n. Static geometric observations are indicated as a scatter point, and kinematically constrained data are depicted as fault growth lines. Faults may grow via the 'isolated' model where increases in displacement and length occur sympathetically (e.g. Walsh & Watterson, 1988), producing the blue growth trajectory; or the constant-length model where the final length is established early in slip history, prior to accumulation of significant displacement (e.g. Meyer et al., 2002;Walsh et al., 2002), producing the green trajectory. The 'hybrid' model (Rotevatn et al., 2019) in orange characterise a rapid stage of lengthening followed by a stage of displacement accrual 70%-80% (Rotevatn et al., 2019). However, it remains unclear how the hybrid model is represented in current global D-L data.
Our current understanding of fault segment, array and network growth is largely geometric and based on observations from modern extensional basins or ancient basins imaged in seismic reflection data. There are very few studies that provide a temporal framework for the growth of normal fault networks using their D-L geometries (e.g. Meyer et al., 2002;Tvedt et al., 2016). This is mainly due to the lack of age-constrained growth strata to record the timing of fault activity (e.g. Dawers et al., 1993;Schlische et al., 1996;Peacock & Sanderson, 1991). There are even fewer attempts to determine the evolution of entire fault networks that likely dissect the entire crust and which bound tens-of-kilometres-wide, graben and half-graben basins (e.g. Claringbould et al., 2017). We suspect this may reflect the time-consuming effort required to: (1) manually interpret many seismically imaged normal faults; (2) subsequently extract geometric properties (i.e. displacement, length) from the rift-related fault network; and (3) map spatially related growth strata and use these to assess fault kinematics at a range of scales. As a result, most studies focus on the growth of large, basin-bounding structures (i.e. major fault arrays; e.g. McLeod et al., 2000), whilst ignoring smaller, intra-basin ones (i.e. fault segments; e.g. Tvedt et al., 2016).
Here we use 3D seismic reflection and borehole data from the Exmouth Plateau, NW Shelf, Australia to determine the growth of rift-related, crustal-scale fault arrays and networks over geological timescales (>10 6 Ma). This study area was subject to Late Triassic-to-Early Cretaceous extension, forming extremely well-defined and seismically well-imaged normal faults, half-grabens and grabens (Exon & Willcox., 1978;Mutter & Larson, 1989). The relatively large extent of the study area (1,200 km 2 ) is comparable to major sub-basins in areas of active continental extension, which contain large, seismogenic fault networks (e.g. East Africa; Albaric et al., 2010;Baker et al., 1972;Poggi et al., 2017;Central Italy (Collettini & Barchi, 2002;Cowie et al., 2013;Luccio et al., 2010;Roberts and Michetti, 2004), Central Greece (Bell et al., 2009;King et al., 1985;Taylor et al., 2011), SE Russia (Logatchev & Zorin, 1992;Mats, 1993), western Turkey (Seyitoglu, 1996 and the western USA (Hamilton, 1987;Thatcher et al., 1999). This makes our study area an excellent ancient analogue with which to better understand the geological evolution of areas of active continental extension and areas that ultimately proceed to full lithospheric breakup (e.g. many continental passive margins).
To the best of our knowledge, this study is the first to quantify the temporal evolution of such a large, rift-scale fault network using geometric observations of all seismically resolvable faults. We first define the present-day structure of the fault network (i.e. the final product of rifting). Second, we use depocentre mapping (via the use of time-thickness maps known as isochrons) and fault displacement backstripping to reconstruct the geometry of the entire fault network at an earlier stage of rifting. This allows us to determine the growth trajectory of individual fault segments and arrays, which we then compare to global D-L datasets ( Figure 2). Based on our results, we propose a generic model for the growth of fault networks, which although based on our observation from NW Australia, might be more broadly applicable. Our model links our current knowledge of the styles of segment-and array-scale growth (i.e. propagating vs. constant length models) to the development of the larger (host) fault network. Our model also provides an explanation for the long-recognised scatter observed in global D-L datasets and sheds light on how continental lithosphere stretches during early rifting (e.g. Cowie & Scholz, 1992;Cartwright et al., 1995;Kim & Sanderson, 2005;Rotevatn et al., 2019).

| GEOLOGICAL HISTORY
The study area is located on the Exmouth Plateau (Exon and Willcox 1978) portion of the Northern Carnarvon Basin, NW Shelf of Australia (Geoscience Australia 2014) ( Figure 3a). Formation of the Northern Carnarvon Basin was initiated during the Carboniferous to Permian in response to the breakup of Pangea, with episodic rifting continuing until the Late Cretaceous (Driscoll & Karner, 1998;Colwell & Stagg, 1994). Late Triassic rifting produced NE-trending sub-basins (Figure 3a), until a phase of regional uplift in the Callovian (Barber, 1988;McCormack & McClay, 2013;Metcalfe, 2013;Pryer et al., 2002). Extension renewed in the Late Jurassic as the Argo block separated from Gondwana, culminating with seafloor spreading and the opening of the Argo Abyssal Plain during the Hauterivian (Heine & Müller, 2005;Gibbons et al., 2012;Metcalfe, 2013;Tindale, 1998). The period of Late Jurassic-to-Early Cretaceous rifting was accommodated by the formation of NE-SW-striking faults, locally reactivating Late Triassic faults (e.g. Jitmahantakul and McClay, 2013;Magee et al., 2016).
The sedimentary succession of the Exmouth Plateau comprises thick, Upper Triassic pre-rift succession of the fluvio-deltaic and marginal marine deposits of the Mungaroo Formation (Figure 4) (Driscoll and Karner, 1988;Longley et al., 2002Tindale et al., 1998. During the Late Triassic to Late Jurassic, a syn-rift succession comprised of shallow marine claystone and limestone of the Brigadier Formation (Late Triassic), and marine claystones of the Athol and Murat formations (Early to Middle Jurassic), was deposited within normal fault-bound depocentres (e.g. Longley et al., 2002;Colwell & Stagg, 1994;Stagg et al., 2004;Tindale et al., 1998). During the Late Jurassic to Early Cretaceous, regional uplift resulted in a decrease in water depth and an influx of interbedded shale and fluvial deltaic sands of the Barrow Group ( Figure 4) (Boyd et al., 1992). A major transgression during a subsequent period of subsidence resulted in the deposition of a thick succession of post-rift marine deposits of the Muderong Shale, Windalia Radiolarite and Gearle Siltstone, which are collectively known as the Winning Group (Driscoll & Karner, 1998;Hocking, 1987) (Figure 4).

| Seismic reflection and well data
Our database comprises a pre-stack time-migrated 3D seismic reflection survey (HEX07B) located in the WA-346-P permit area, and two wells (Thebe-1 and Thebe-2) ( Figure 3b). The seismic data image to a depth of ca. 4.5 s TWT (ca. 3.7 km) and covers 45 km × 39 km, providing a total areal extent of 1,200 km 2 . Crosslines are spaced at 25 m and inlines at 12.5 m. The dominant seismic frequency (extracted at the position of Thebe-1) is 45 Hz at the depth of the studied fault network. From time-depth plots derived from well checkshot data, we estimate the average velocity to be 2,480 m/s within the Jurassic syn-rift units; when combined with frequency data, this results in an estimated spatial vertical resolution of ca. 14 m (λ/4). Overall, the data quality is excellent throughout and enables a highresolution analysis of fault structure and evolution.
The Thebe-1 and Thebe-2 wells are drilled in the footwalls of tilted fault blocks ( Figure 3b) and contain a standard well-log suite (e.g. gamma ray, density, sonic, checkshot and biostratigraphic data). We present our measurements of fault throw in milliseconds two-way time (ms TWT) rather than depth, given we are primarily interested in the relative rather than absolute changes of throw in a fault network that is buried at a fairly consistent depth across the study area (cf. Jackson et al., 2017;Tvedt et al., 2016). To estimate depth in metres from observations in ms TWT, we use checkshot data from wells that define a polynomial relationship between time and depth. We have not attempted to decompact our throw values as the loss of throw is likely low (<15%) in the mixed sand-mudstone growth sequences that characterise much F I G U R E 3 The location of our study area: (a) location of the Thebe field seismic dataset within the Exmouth Plateau, North Carnarvon Basin, North West Shelf of Australia. Sub-basin boundaries and the main structural NE-trends are taken from Geoscience Australia (2011); and (b) major faults F1, F2 and F3 mapped within the Thebe dataset. The segments of major faults are labelled in colour, correlating to their respective throw-profiles in Figure 9. Minor faults are coloured in light grey. Wells Thebe-1 and Thebe-2 are shown of the syn-to post-rift succession of the Exmouth Plateau (cf. Taylor et al., 2008). A lack of information on sediment compaction parameters from our two wells also hampers decompaction of throw values. For better comparison with global D-L compilations, we converted our throw data to displacement data by assuming an average fault dip of 60°.

| Structural framework
We mapped ten horizons throughout the seismic reflection volume (labelled H1-H10; see ages in Figures 4 and 5) to delineate the present basin structure and then unravel the fault evolution. Fault interpretation was conducted throughout the seismic volume, aided by overlaying a dip attribute (i.e. the deviation of a seismic reflection from a horizontal plane; Chopra and Marfurt, 2007) onto mapped horizons. We measured the tip-totip lengths of all seismically resolvable, individual fault segments at the base of the syn-rift, given this structural level captured most of the rift-related strain (see Figure 1 for definitions). In D-L plots, we show the scaling properties of individual fault segments and arrays (F1-F3). We classify faults arrays as major (F1-F3) if their throw F I G U R E 4 Horizons and their respective formation ages. H1 and H10 are pre-rift and post-rift reflectors, whereas H2-H9 define the syn-kinematic unit of the study area. H4, H5 and H6 do not contain biostratigraphic ages, therefore age is inferred through the assumption of constant sedimentation. Our different methodologies to understand rift evolution are resolved to different stages in time, which has been demonstrated on the right is >150 ms (TWT) and minor if their throw is <150 ms TWT (labelled on Figure 3b).
Fault throw was constrained by measuring the vertical distance between horizon cut-offs in the footwall and hanging wall on seismic profiles trending perpendicular to local fault strike (cf. McLeod et al., 2000;Tvedt et al., 2013). Values are taken every 100-400 m along the fault depending on the structural complexity and throw variability along the fault. In the presence of erosional footwall fault scarp degradation and fault-related ductile (continuous) deformation ( Figure 5), the horizon cut-offs are projected towards the fault plane along the regional structural trend (e.g. Barrett et al., 2020;

| Constraining rift evolution
To constrain the fault network evolution, we use isochron (time-thickness) mapping and throw backstripping techniques. The temporal framework is provided by six seismic horizons that are directly age-constrained by biostratigraphic data from wells (H2, H3, H4, H8, H9 and H10; Figure 4). The ages of a further three syn-rift horizons are not constrained by the wells (H5, H6 and H7); their ages are therefore inferred based on an assumption of constant sedimentation rates between the bounding, age-constrained horizons ( Figure 4). We define H2 (Top Mungaroo; 208.5 Ma) as the base syn-rift (i.e. the start of rifting), consistent with other publications (e.g. Stagg  , 2004;Tindale et al., 1998), given we observe characteristic systematic upward decreases in throw from H2 to H9 ( Figure 6) and no across-fault sediment thickness changes between H1 and H2 ( Figure 5). Fault activity ceases at H9 (Top Muderong; 113 Ma), based on observations of no across-fault changes in thickness above H9, which we thus define as top syn-rift (i.e. the end of rifting; Figure 6). Although minor (ca. 10 m; Figure 6) reactivation of some large, basement-involved faults (F1-3) may have occurred at the same time as widespread polygonal faulting in strata above H10, we consider the throw addition associated with this to be negligible (ca. 2%) compared to the faults' total displacement (650 ms TWT). Our ages for top and base syn-rift thus give a total rift duration of ca. 85.5 Myr (Figure 4). We produced six isochron maps within the syn-rift interval (SU1-6; Figure 10), which we use to constrain patterns of fault-controlled subsidence and thus fault growth (Jackson & Rotevatn, 2013;McLeod et al., 2000;Thorsen, 1963).
We use fault throw backstripping to constrain the style and patterns of fault growth (Childs et al., 1993;Dutton & Trudgill, 2009;Jackson et al., 2017;Rowan et al., 1998). Throw backstripping involves the subtraction of fault throw at multiple stratigraphic levels to determine the throw accumulated during deposition of the backstripped unit (see Chapman & Meneilly, 1991). Faults can be backstripped using the 'original method' or the 'modified method'. The 'original method' directly subtracts throw across a shallower horizon from that of the deeper horizon (Chapman & Meneilly, 1991;Petersen et al., 1992), making no assumptions about fault growth style. The 'modified method' subtracts the maximum throw along the horizon from the entire fault surface. This method implicitly assumes that the studied fault grew in accordance with the propagating fault model (Jackson et al., 2017;Rowan et al., 1998). In this study we use the original method to avoid any a priori assumption of the style of fault growth. We then compare our throw backstripping results with depocentre lengths revealed by isochrons; this allows us to assess fault length at different points in the rift history (cf. Jackson et al., 2017;Jackson & Rotevatn, 2013).
The Exmouth Plateau was sediment-starved (i.e. underfilled) during Jurassic extension (Gartrell, 2000;Marshall & Long, 2013). Because of this, accommodation associated with the very earliest stage of extension may have been unfilled, meaning: (1) we cannot determine if any observed fault lengthening occurred by segment linkage or tip propagation (i.e. infilling sediments do not exist to record this); and (2) our measured trace-lengths may thus have been established earlier (i.e. more quickly) than estimated, such that our estimates should be seen as maximum possible values (Jackson et al., 2017;Lathrop et al. 2021).

| Time-structure maps
We first describe the present-day structure of the study area (Figures 7 and 8). We mapped 150 faults over an area of 45 × 39 km that displace the 208.5 Myr base syn-rift horizon (H2, Top Mungaroo); this surface records the total cumulative rift-related strain of the fault network.  (Figures 7c,d). N-S-striking faults are consistently offset across NE-SWstriking faults, indicating the former formed before the latter (Figures 5 and 7a).

| Throw distribution
Maximum throw on the syn-rift fault network varies from 10 to 650 ms TWT (Figure 9a). Throw varies alongstrike on major faults F1-F3, with maximum throw (up to 650 ms TWT) typically occurring near their centres. F I G U R E 7 Time-structure maps (ms TWT) showing the present day structure of the dataset. Faults at each successive interval are mapped using a variance attribute overlay for higher confidence Throw gently decreases to zero at the tips of fault arrays (e.g. F2-a, Figure 9a). Individual segments forming the fault arrays show steep throw gradients towards their tips (e.g. F2-b; see Figure 1 for definitions). Strain is distributed on several splay faults at the ends of fault arrays (Figure 10a for location). Overall, the aggregate throw of major faults produces a broadly coherent, symmetric bell-shaped profile (Figure 9a). In detail, however, throw minima exist along many major faults (e.g. F2-a, F2-b, F3-b), which correlate with along-strike variations (F2, F3; Figure 9a). This observation suggests major faults were previously segmented, at least at this structural level (Figure 9a).
In contrast to the major faults, minor faults typically have a uniform throw value along strike (ranging from 20-80 ms TWT), but exhibit an abrupt decrease in throw towards their lateral tips (Figure 9a). Local throw minima on minor faults also correspond to abrupt changes in the strike, again suggesting a history of fault growth via segment linkage (Figure 7a).

| Rift evolution
In order to determine the growth of the fault network, we now determine the activity of individual fault arrays through the use of isochron maps and throw backstripping.

| Isochron analysis
Within the earliest resolvable time period, SU1 (from 208.5 to 201.3 Ma; i.e. capturing fault activity during the first ca. 7.2 Myrs of rifting), we observe across-fault sediment thickening (ca. 100 ms TWT) in the hanging-wall of all major faults (F1-F3; Figure 10a). These fault-controlled depocentres span the entire present-day fault trace length, indicating that major faults had already reached their near-final length during the first 7.2 Myr of the ca. 85.5 Myr rift history (i.e. only 8% of the total rift duration). Thick (up to 200 ms TWT), mound-shaped features are present in the hanging walls of F1, F3-a and F3-b, which we interpret as the seismic expression of fault scarp-derived talus (Figure 10a) (Bilal et al., 2020;Barrett et al., 2020). The presence of fault scarpderived deposits in the fault hanging walls provide further evidence that the faults were active during this period. These deposits also suggest that fault slip rates outpaced sediment accumulation rates, and that at-surface relief formed which was then eroded.
Thickness variations in SU2 (Figure 10b) show that fault-controlled depocentres persisted adjacent to F1-F3 (from 201.3 to 163.5 Mya; i.e. capturing fault activity during the subsequent 37.5 Myrs of rifting). SU2 depocentres are wider (normal to fault strike) than those in SU1, and show a greater thickness increase up to 3.5 km away from the fault, indicating increased hanging wall F I G U R E 8 Dip direction of the fault network. The rose diagram shows the average strike azimuth and dip direction per fault segment. The rose diagram is split into two halves, whereby the blue coloured faults are east dipping, and the red coloured faults are west dipping flexure (Figure 10b). SU2 is thickest in the middle of the depocentres adjacent to F3-a and F3-b (<170 ms TWT) instead of immediately adjacent to the bounding faults (cf. Figure 10a), which we attribute to the underlying scarpderived sedimentation in SU1 leading to an apparent reduction in the thickness of deposits in SU2 (Figure 10b, Figure 5). Although there is an overall thickness increase across both F3-a and F3-b (ca. 120 ms TWT), hanging wall depocentres are segmented and separated by fault-normal anticlines and indicate that they were soft-linked and separated by a relay zone (Figure 10b). Thickness variations also suggest that segments F2-a and F2-c were active during deposition of SU2, and a lack of thickness changes across segment F2-b suggest that F2 was not hard-linked at this time (Figure 10b).
Within SU3 (from 180 to 143 Mya; i.e. capturing fault activity during the subsequent 37 Myrs of rifting) we observe a thickness increase of 100 ms TWT across all segments of F2. This is in contrast to activity within SU2, where no thickness changes occurred across the faults central segment (F2-b), indicating that by 37 Myrs after the onset of rifting, F2 had hard-linked to form a throughgoing array (Figure 10c). Additionally, thickening across both F3-a and F3-b segments suggest that F3 was now acting as a single, hard-linked fault that was principally growing via the accumulation of displacement at the expense of lengthening (Figure 10c). In summary, F2 and F3 had grown to the maximum lengths after only 33% of the total rift duration, achieving this by segment linkage and relay breaching. This style of growth resulted in the F I G U R E 9 Throw length profiles (segments coloured by Figure 3b) and their map view throw distribution, for (a) the present day fault network, that is the final product of rift history; and (b) the reconstructed, backstripped fault network, where faults are backstripped to 28.5 Myrs (H4; Top Murat) development of abandoned splays in their hanging wall (Figures 7 and 9). SU4 (from 143 to 141 Mya; i.e. capturing fault activity during the subsequent 2 Myrs of rifting) shows a significant thickness increase across F2, indicating this major structure was active at this time. SU5 (from 141 to 138 Mya; i.e. capturing fault activity during the subsequent 3 Myrs of rifting) shows no thickness variations across F1 and F2, suggesting that these structures were inactive by this time (Figure 10e). F3-a also became inactive as evidenced by constant across-fault thickness of SU5. In contrast, the southern part of the fault array remained active (F3-c and F3-d), defining a wide depocentre ca. 4 km across strike. SU6 (from 138 to 123 Mya; i.e. capturing fault activity during the subsequent 15 Myrs of rifting) shows that nearly all the northern faults were largely inactive at this time, as only the southern segments of F2 and F3 are associated with thickness variations. However, these variations appear broad and regional and occur on a length-scale larger than individual faults (Figure 10f). Given that the timing of deposition of SU6 coincides with a period of widespread regional subsidence along the margin (Driscoll & Karner, 1998;Hocking, 1987), we infer that F3-c, F3-d, F2-b, F2-c may be less active than the spatially related thickness changes suggest.

| Fault throw backstripping
Isochrons qualitatively reveal that fault lengths were rapidly established (7.2 Myrs) before significant displacement had accumulated, as fault-controlled depocentres nearly span across the entire length of the present-day fault traces within the earliest resolvable slip increments (Figure 10a,b). This style of growth is consistent with the constant-length model of fault growth, hence the original rather than modified method is considered more appropriate to backstrip the fault network (e.g. Jackson F I G U R E 1 0 Isochron maps showing thickness changes between syn-depositional horizons. The syn-rift unit is divided into six, labelled as SU1-SU6 et al., 2017; see methodology section). Because H3 is not present across the entire fault network due to footwall erosion (see Figure 11a), we backstrip to the next resolvable increment (i.e. H4; 180 Ma) to obtain a view of fault network structure at an intermediate stage of extension (i.e. at 33% of the total rift duration). As a result, we can gain a view of rift-related strain magnitude and patterns after 33% (Figure 9b) and 100% of the rift duration (Figure 9a).
Our results show that by 33% of the rift history, major faults F1, F2 and F3 had accumulated maximum throws of throws of 261, 181 and 365 ms TWT (ca. 261, 183, 372 m) respectively ( Figure 9b); this is approximately half of the present-day throw. In detail, our analysis shows that throw accumulation was relatively uniform (averaging ca. 200 ms TWT) along F1, with an abrupt decrease in throw occurring near a throw minima located 16.3 km along the fault (Figure 9b). This minimum corresponds to a small hanging wall splay and a subtle change in fault strike, thus we infer it defines a relay that was breached before this backstripped interval (i.e. before 33% of the rift duration; Figures 9b and 3b). F1 continues north beyond the seismic dataset, so we cannot tell whether the maximum throw on this structure is located further north; however, extrapolation of the throw gradient at the northern end of F1 suggests this is not the case (Figure 9b). Our results indicate that during the latter part of rifting (33%-100% of the rift duration), F1 accumulated throw near its centre without appreciable lengthening. Our present-day throwlength plot shows a broadly symmetric, bell-shaped profile, with no strong indication of palaeo-segmentation (i.e. displacement minima; Figure 9a). The absence of throw minima indicates that previous locations of relays may become future sites of enhanced throw (e.g. Faure-Walker et al., 2009).
During the initial stage of rifting, F2 exhibited a similar throw distribution as F1, in that throw (average of 130 ms TWT) was fairly evenly distributed along its constituent segments (F2-a, F2-b and F2-d; Figure 9b). The throw minima between F2-a and F2-b suggests that these segments were not geometrically linked at this time. Instead, F2-b developed alongside two small hanging wall splays (pink and brown in Figure 9b). At this time, a small splay (purple F I G U R E 1 1 Our results plotted against global displacement-length data in grey. The growth trajectory represents the final 67% of total rift duration. Data is coloured by the dip map (see Figure 8) whereby blue-coloured faults are east-dipping, and red-coloured faults are westdipping. Data is coloured by the fault map in the inset, where major faults are dark blue (F1), light blue (F2) and green (F3), and remaining minor faults are pink. The overall D-L of fault arrays are indicated as a scatter point, and their constituent-coloured segments are indicated as growth lines (see Figure 9 for mapped fault segments and arrays) in Figure 9b; see Figure 3 for corresponding map) appear to be kinematically interacting with fault segments F2-c and F2-d, given they together define a throw profile that shows a gradual southward increase in throw (Figure 9b). However, the purple splay and F2-d appear to have been abandoned during the latter phase of rifting (33%-100% of rift duration); instead, F2-c gained significant throw (from ca. 80 to 330 ms TWT) and hard-linked with F2 as rifting continued (Figure 9). During the latter phase of rifting, F2-b also accumulated a significant amount of displacement (up to 550 ms TWT) to produce the symmetric, bellshaped throw profile that presently characterises its host array, F2 (Figure 9a).
Backstripping of F3 reveals that during the initial stage of rifting, a series of throw maxima were located towards the northern portions of segments F3-a, F3-b and F3-c. Whereas F3 still displays a generally symmetric throw profile, its constituent segments are defined by a rather irregular, asymmetric throw profile (Figure 9b). We infer that throw minima (e.g. between F3-a and F3-b) record relay zones that were breached before this time-step (Figure 9b). By comparing the throw distribution after ca. 33% of the rift history with that presently observed (Figure 9a), our results show that during the last 67% of rift history, the asymmetric and irregular profile of F3 evolved to have a symmetric, bell-shaped profile, with the throw maxima migrating towards the present centre of the fault.
Backstripping of minor faults show that fault lengths and throw distribution remained relatively similar during the latter 67% of the rift history (Figure 9). A number of minor faults appeared to have died-out before the deposition of H4 (180 Ma; Top Murat Siltstone); however, ductile or so-called continuous deformation, make offset difficult to determine at H4. Our interpretation that the minor faults became inactive <28.5 Myr into the rift event is supported by isochron data, which shows that post-SU2 units do not thicken across them, as well as seismic cross-section evidence that show minor faults tipping-out below H4. Figure 11 shows our fault displacement-length data for all 150 faults mapped in the Thebe fault network, relative to previously collected D-L data (see Figure 2 for references therein). Fault length and displacement measurements for the present-day fault network (Figure 9a), coupled with the length and displacement of the reconstructed fault network at ca. 33% of the total rift duration (Figure 9b), allow us to define how D-L relationships have changed over the last ca. 67% (57 Myr) of rift history. Note that, in this plot, we have converted our throw data in ms TWT to displacement in metres, using an average dip of 60° and borehole-derived time-depth relationships (see Section 3.1). On a log-log scale, our data plot almost three orders of magnitude across the linear trend (n = 1) defined by existing global D-L scaling relationships for normal faults (e.g. Cowie & Scholz, 1992;Marrett & Allmendinger, 1991;Walsh & Watterson, 1988;Watterson, 1986) (Figure 11). D-L data are coloured by dip direction; blue-coloured faults are broadly east-dipping and red-coloured faults are broadly west-dipping (see Figure 8); in the inset the same data is coloured to highlight major faults (F1, F2 and F3) whereby hard-linked geometry are represented as a point and respective constituent segments that make up F1-3 as growth lines. Major, west-dipping (red-coloured) faults are presently somewhat 'over-displaced' (c ~ 0.1) relative to the smaller, antithetic, east-dipping faults (c ~ 0.01) ( Figure 11). Most of the over-displaced faults are short trace-length, yet high-displacement splays abandoned in the hanging walls of faults after relay breaching and segment linkage (Figure 11). Over-displaced faults may be indicative of interpreting faults as segments instead of arrays. For example segment F3-a has a displacement-length ratio of 0.06, although the larger fault array in which it is contained (i.e. F3) is 46-km long and has a throw of 656 ms TWT (ca. 858 m in displacement; Figure 9a). The fault array thus has a lower D-L ratio of 0.02 (Figure 11), demonstrating that D-L (and T-L) ratios are higher for individual fault segments than arrays, although both still plot within the scatter observed in global range D-L data (Figure 11 inset). Similar observations are noted by Peacock & Sanderson (1991), Dawers et al. (1995) and Roberts and Michetti (2004), highlighting how mixing the analysis of fault segments (e.g. Meyer et al., 2002;Roberts et al., 2004) and/or fault arrays (e.g. McLeod et al., 2000), without considering faults kinematics (e.g. relay breaching and tip abandonment), may contribute towards the range of scatter seen in D-L scaling datasets.

| Fault growth trajectory using D-L relationships
Our time-constrained displacement-length relationships reveal a vertical fault growth trajectory, therefore most faults accumulated displacement without lengthening (i.e. the constant fault growth model; Jackson & Rotevatn, 2013;Nicol et al., 2005;Walsh et al., 2003). Quantitative (D-L and fault backstripping) and qualitative (isochron mapping) results both clearly demonstrate that many of the faults established their near-final length relatively early in the ca. 85.5 Myr rift history. Throw backstripping shows that for many of the faults, their lengths were essentially fixed after only ca. 33% of the total rift duration ( Figure 9); this is similar to that documented by Walsh et al. (2002), Jackson et al., (2017) and Rotevatn et al., (2019). In fact, our isochron analysis suggests that fault length establishment may have occurred after as little as 8% of the total rift duration (Figure 10a).
The last ca. 67% (and possibly up to ca. 92%) of the 85.5 Myr rift history was characterised by strain localisation and displacement accumulation on the major faults (F1-F3). Strain localisation onto larger faults that are ultimately active for longer than smaller faults is also seen in the Inner Moray Firth, North Sea (Nicol et al., 1997;Walsh et al., 2003), and the Timor Sea, NW Shelf of Australia (Meyer et al., 2002), and in numerical (e.g. Cowie et al., 1995;Cowie, 1998;Gupta et al., 1998;Naliboff et al., 2020) and physical models (e.g. Ackermann et al., 2001;Mansfield & Cartwright, 2001). Our results demonstrate that the pattern of displacement accumulation and/or strain localisation is itself not straightforward, due to complex fault interactions that result in relay breaching and associated splay abandonment (e.g. Figure 10). Isochrons and throw backstripping reveal a distributed fault network that is initially composed of fault arrays with multiple displacement maxima and minima. Detailed backstripping analyses show that in the early stages of fault development, fault segments exhibited profiles that were broadly symmetric (i.e. throw was greatest at the fault centres and decreases towards their tips). At the end of rifting (captured by present day geometry), throw profiles of fault arrays are also symmetric. However, at an intermediate stage of rifting (33% of total rift duration; Figure 9b), whereas fault segments are symmetric, the overall fault array in which they were contained exhibited an irregular, asymmetric throw profile. During this time, we also observe multiple examples of splay abandonment (e.g. F2-c) related to strain localisation onto more optimally positioned fault segments that lie outside of stress shadows (e.g. F2-b). The alongstrike migration of maximum slip through time may be important when trying to understand multi-fault ruptures that may occur on previously unrecognised faults, producing larger-than-expected earthquakes (e.g. Wesnousky, 1986;Leeder et al., 1991;Goldsworthy and Jackson, 2000;Gupta and Scholz, 2000;Hamling et al., 2017). Our results also suggest that although final lengths may be rapidly established, this may be followed by a period of the systematic along-strike migration of seismogenic slip (Figures 9  and 10).
While major faults experienced displacement accumulation, accounting for ca. 67% of rift history, a number of east-dipping antithetic faults did not slip post-180 Ma and thus remained inactive. These minor, antithetic faults are located between the larger faults and are long (i.e. >20 km), yet have low displacements (<30 m) thus lie significantly 'under-displaced' on D-L profiles (close to c = 0.01 on Figure 11). Although there are no growth strata to directly establish if these structures were lengthening and/or accumulating displacement immediately prior to their death, their long trace lengths, palaeo-segmentation on time-structure maps (Figure 7), thickness variations from isochrons ( Figure 10) and throw minima from throw backstripping (Figure 9), suggest that prior to becoming inactive, these faults were lengthening via segment linkage. We suspect that these long, low-displacement faults may be 'fossilised' in their early stages of development, providing a snapshot of the fault network structure before strain localisation and displacement accrual.
As all the major and minor faults dip westward and eastward, respectively, the structures onto which strain was eventually localised must have been set early in rift history (i.e. as early as the first ca. 8% of the total rift duration). This may be attributed to positive stress feedback processes occurring during early network development, whereby shear stress reduction or 'stress shadow zones' inhibit fault growth and the nucleation of new cracks (e.g. Ackermann and Schlische, 1997;Cowie, 1998;Hu and Evans, 1989), and optimally positioned faults form and slip (i.e. faults with across-strike spacing of 12-14 km). Early establishment of the near-final fault lengths reinforces the suggestion that lateral tips of faults are effectively pinned by stress interactions with neighbouring faults (e.g. Burgmann et al., 1994;Contreras, 2000;Gupta and Scholz, 2000;Willemse et al., 1996). During the initial stages of extension, 'distributed' faulting, where faults grow by predominantly lengthening, may occur until a given rock volume (e.g. a mechanical layer, which could scale from a bed to the crust) is 'strain saturated', as predicted by numerical and physical models (e.g. Ackermann et al., 2001;Cowie & Scholz 1992b;Cowie et al., 1994Cowie et al., , 1995Wu & Pollard, 1995). Thus 'distributed' faulting, as seen in the Central Afar and on the Asal Ghoubbet faults in East Africa (c ~ 0.012 in Manghietti et al., 2015), may characterise relatively immature, active rifts, whereas 'localised' faulting, similar to that seen in fault populations imaged in seismic reflection data (e.g. northern North Sea; McLeod et al., 2000, NW Shelf of Australia;Black et al., 2017), may characterise mature, ancient (i.e. inactive) rifts.
Our observation that minor (E-dipping) faults grow and still appear relatively under-displaced relative to larger faults onto which strain is subsequently localised (with higher displacement of red faults; Figure 11), has only been revealed from studying a complete fault network. Considering the faults imaged in our data span nearly three orders of magnitude, and given that the 150 faults formed in broadly similar rock types in response to the same stress regime, this implies that some of the scatter seen in D-L plots may simply be related to fault maturity and not simply limitations inherent to data collection methods (e.g. Kim and Sanderson, 2005) and/or differences in rock type and causal stress regime (e.g. Cowie & Scholz, 1992;Gillespie et al., 1993;Scholz et al., 1993).
Our findings illustrate the importance of collecting data from entire fault networks developed over large areas if we hope to determine how normal faults grow. Datasets that allow inspection of only a small part of the rift (e.g. small-scale field-based studies), or analysis of faults lacking growth strata, could produce D-L relationships that appear anomalous (i.e. relatively under-or over-displaced) simply due to fault maturity (i.e. the datasets sample only the very young or old faults). Taking the trendline of datasets may produce an exponent value n that may not be representative of geometrical fault scaling relationships. We suspect that this may be why the true value of n is contested (e.g. Clark and Cox, 1995;Cowie & Scholz, 1992;Dawers & Anders, 1995;Gillespie et al., 1993;Marrett & Allmendinger, 1991;Soliva et al., 2005;Torabi et al., 2011;Watterson, 1988;Xu et al., 2006). By studying complete fault networks rather than individual fault arrays, we can better understand how faults grow to accommodate stretching of continental lithosphere.
While our study demonstrates that during the last ca. 67% of rifting faults either grew by displacement accumulation or became inactive, the magnitude of displacement accumulation (i.e. the vertical trajectory on D-L profiles) during this phase does not extend outside of the overall scatter observed in global D-L plots (where c ~ 0.001, Figure 11). This is important, given that recently proposed fault models (constant-length and hybrid) commonly present growth trajectories schematically on D-L plots lacking scales on their axes (e.g. Childs et al., 2017;Rotevatn et al., 2019; Figure 2). Due to a lack of scales, the constant-length and hybrid models do not seemingly take into account the geometric constraints evident within D-L scaling space (i.e. it is unlikely that a fault could have a maximum trace length of 10 km and have only 1 m of displacement, given no such structure has yet been recorded). We suspect that many observations from ancient (inactive) and modern (active) extensional settings that endorse the constant-length model (e.g. Meyer et al., 2002;Tvedt et al., 2016;Jackson et al., 2017;Rotevatn et al., 2019) in fact still lie within the bounds of our proposed geometrical constraints, as our results from offshore Australia show ( Figure 11).

| A conceptual model for faults arrays within an evolving network
Our results show that seismic reflection data can constrain the style and duration of fault array growth during the displacement accumulation stage of crustal-scale, natural fault networks (Figure 11; Meyer et al., 2002;Walsh et al., 2002). In contrast, the preceding stage of fault growth, which is dominated by fault lengthening, is poorly understood (e.g. Rotevatn et al., 2019), principally because the related structures and growth strata are hard to resolve in seismic reflection datasets (Jackson et al., 2017). One interpretation is that immediately prior to the displacement accumulation stage, the fault trajectory on D-L plots inflected from sub-horizontal to sub-vertical outside and, more specifically, below the observational scatter (e.g. Walsh et al., 2002). Scaling data associated with the rupture-slip scaling relationship characterising individual active earthquake slip occur within this region of D-L plots (c ~ 0.0005; Wells and Coppersmith, 2005). An alternative interpretation is that the D-L trajectory may have inflected from sub-horizontal to sub-vertical within the range of the observed D-L scatter. Given that (1) D-L ratios of less than ca. 0.001 are not observed in nature; and (2) physical models (containing faults with displacements of <1 m; e.g. Schlagenhauf et al., 2008) and seismic reflection datasets (imaging faults with displacements of >0.5 km (e.g. Meyer et al., 2002) reveal a switch from dominant fault lengthening to dominant displacement accumulation (i.e. the 'constant-length' model), we present a conceptual model of fault growth which aims to honour both kinematic and geometric constraints imposed by global D-L observations, as well as the recognition that stress feedbacks between arrays play a key process in network evolution (e.g. Ackerman et al., 2001). We stress that Stages 4 and 5 are directly constrained by the results of this study, whereas Stages 1 to 3 occur on shorter temporal and spatial scales than what we can resolve in this study.
Faults initiate as numerous isolated, decimetreto-metre-scale segments that, for example might be fully visible at the scale of individual field exposures (i.e. outcrop-scale faults in Stage 1 in Figure 12; Cowie et al., 2000;Crider & Peacock, 2004;Dawers et al., 1995). Preferential lengthening and linkage of faults along relatively brittle mechanical layers may occur at this stage (Peacock and Sanderson, 1992;Soliva et al., 2006), with faults tipping-out in encasing layers that are unable to maintain a discrete shear fracture, and which instead accommodate strain in a more ductile manner. Fault lengthening mean that these early faults may appear 'under-displaced' when viewed on D-L plots (Stage 2 in Figure 12; Cartwright et al., 1995;Peacock & Sanderson, 1995;Seagall & Pollard 1980). Once the rock volume is sufficiently saturated with fault trace lengths that define the overall pattern of the fault network, strain continues to accumulate on optimally positioned faults (i.e. those that receive positive stress feedback from adjacent structures) gain displacement, whereas those in adjacent stress shadow zones become inactive (Stage 3 in Figure 12; Ackermann et al., 2002;Cowie et al., 1998;Spyropoulos et al., 2002). This phase of strain localisation and displacement accumulation produces a vertical growth trajectory on D-L plots, with the dominant, still-active faults ascending vertically until they reach the upper scaling limit (approximately where c = 0.1). Faults that are not optimally positioned are essentially stranded at various points along this vertical trajectory, the precise position of which is controlled by the exact time at which they die. Stress feedback between the remaining, larger active faults now occur over a larger length-scale, this is kinematically and mechanically the same process that occurred at the smaller, outcrop-scale (Stage 1 and 3 in Figure 12). Stillactive, decametre-to-hectometre scale faults continued to accommodate extensional strain, growing by fault propagation and/or segment linkage (Ackermann et al., 2001;Cowie & Roberts 2001;Kim & Sanderson 2005;McLeod et al., 2000;Soliva et al., 2005). Displacement accumulation on optimally positioned faults mean these faults ultimately form very large, kilometre-scale (in terms of their displacement) structures that span the entire mechanical layer of the brittle crust and that bound crustal-scale (i.e. several tens-of-kilometres-wide) graben and half-graben (i.e. red fault segment on 5; Figure 12; Cowie et al., 1995;Spyropoulos et al., 2002;Walsh et al., 2002) (i.e. the crustal scale faults in Stage 5 in Figure 12). Faults that are located in stress shadow zones of large structures again become inactive, being stranded along vertical growth trajectories and appearing under-displaced in D-L space (i.e. the green segment; Figure 12).
Resolving whether length establishment occurs within, for example a few tens of thousands to a few million years is difficult with seismic reflection data from ancient rifts, given that the limited spatial and temporal resolution of seismic data at depth; it may thus be difficult to learn more about the fault lengthening stage from solely the kinematic analysis of seismic reflection data that image ancient basins. We suggest that future studies should focus on active, shallowly buried (and thus seismically well-imaged) rifts containing fault networks flanked by growth strata (e.g. the Whakatane rift, New Zealand; Taylor et al., 2004), but 3D seismic data rather than grids of 2D data will likely be required. Furthermore, numerical and physical models may provide the spatial and temporal information required to fully investigate the fault growth trajectories unresolved in D-L space.

| CONCLUSIONS
• The timescales over which a fault network evolves have previously not been well constrained, due to a lack of dynamic, kinematic data. • To fill this gap, we studied an extensive fault network (1,200 km 2 ) in the Exmouth Plateau, NW Australia and characterised 150 seismically resolvable faults using displacement-length relationships. • Isochron maps reveal that within the earliest resolvable slip increment of 7. Myrs (out of an overall duration of F I G U R E 1 2 Schematic fault evolution model of how faults may grow within their overall network. The fault trajectory moves through observational D-L scatter, that is they are geometrically constrained. This shows how faults may abide by a linear D-L scaling relationship, but grow via 'constant-length' model trajectories through different length scales, here from outcrop scale (Stages 1-3) to rift scale (3)(4)(5) 85.5 Myrs) the majority of fault lengths have already been established. • Fault displacement backstripping to 28.5 Myrs (33% of total rift duration) shows that the fault network was characterised with distributed faulting. Backstripping reveals along-strike slip migration through time to produce symmetric fault arrays. Throw minima suggest previous segmentation and relay breaching. • Our results allow us to view the last growth trajectory (representing the last 67% of total rift duration) on D-L plots. Faults undergo varying amounts of displacement accrual (i.e. the constant length model). Less (to no) displacement is accrued on faults in stress shadow zones. • Inactive faults located in stress shadow zones contribute towards the large amount of global displacement-length scatter due to differences in fault maturity.