Bridging spatiotemporal scales of normal fault growth using 1 numerical models of continental extension 2

11 Continental extension is accommodated by the development of kilometre-scale normal 12 faults, which grow by accumulating metre-scale earthquake slip over millions of years. 13 Reconstructing the entire lifespan of a fault remains challenging due to a lack of observational 14 data of appropriate scale, particularly over intermediate timescales (104-106 yrs). Using 3D 15 numerical simulations of continental extension and novel automated image processing, we 16 examine key factors controlling the growth of very large faults over their entire lifetime. 17 Modelled faults quantitatively show key geometric and kinematic similarities with natural fault 18 populations, with early faults (i.e., those formed within ca. 100 kyrs of extension) exhibiting 19 scaling ratios consistent with those characterising individual earthquake ruptures on active faults. 20 Our models also show that while finite lengths are rapidly established (< 100 kyrs), active 21 deformation is highly transient, migrating both alongand across-strike. Competing stress 22 interactions determine the overall distribution of active strain, which oscillates locally between 23 localised and continuous slip, to distributed and segmented slip. These findings demonstrate that 24 our understanding of fault growth and the related occurrence of earthquakes is more complex 25 than that currently inferred from observing finite displacement patterns on now-inactive 26 structures, which only provide a spatialand time-averaged picture of fault kinematics and 27 related geohazard. 28

www.geosociety.org/pubs/ft20XX.htm, or on request from editing@geosociety.org. Figure 1. The top row shows the strain rate invariant (s -1 ) in the upper crust (5 km depth), documenting active deformation patterns within the last resolvable time increment (10, 5 and 2.5 Myrs for models ran at 2.5, 5, and 10 mm yr -1 respectively). The bottom row shows their corresponding fault length property extracted from the active deformation field. The spatial extent covers the high-resolution (625 m) portion of the model domain. Figure 2. Geometric statistics of time-dependent fault properties. A) Fault D-L evolution for the modelled fault network that experienced 5 mm yr -1 extension. D-L geometries of observable faults are plotted in grey, where different shades correlate to different datasets (references therein). B) The number of active faults throughout model evolution all models. C) The average active length of a given population throughout model evolution. Note that all three models output 50 timesteps, however the age duration for models deformed at 2.5 mm yr -1 , 5 mm yr -1 and 10 mm yr -1 are 10 Myrs, 5 Myrs and 2.5 Myrs respectively.

Governing equations
We use the open-source, mantle convection and lithospheric dynamics code ASPECT Where # is the velocity, , is the viscosity, -̇ is the second deviator of the strain rate tensor, 0 is pressure, 1 is density, and 2 is gravitational acceleration.
Temperature evolves through a combination of advection, heat conduction, shear heating, and adiabatic heating: Where 3 ! is the heat capacity, 6 is temperature, 7 is time, : is thermal diffusivity, and < is the rate of internal heating. Respectively, the terms on the right-hand side correspond to internal head production, shear heating, and adiabatic heating.

Rheological Formulation
Rheological behaviour combines nonlinear viscous flow with brittle failure (see Glerum et al., 2018). Viscous flow follows dislocation creep, formulated as: Above, B ## $ is the second invariant of the deviatoric stress, C is the viscous prefactor, F is the stress exponent, -̇# # is the second invariant of the deviatoric strain rate (effective strain rate), G is the activation energy, H is pressure, I is the activation volume, 6 is temperature, and J is the gas constant Brittle plastic deformation follows a Drucker Prager yield criterion, which accounts for softening of the angle of internal friction (K) and cohesion (3) as a function of accumulated plastic strain: The initial friction angle and cohesion are 30 and 20 MPa respectively, and linearly weaken by a factor of 2 as a function of finite plastic strain. We localise deformation in the center of the model by prescribing an initial plastic strain field. Rather than a single weak seed (e.g., Lavier et al., 2000;Thieulot, 2011;Huismans et al., 2007), or randomised distribution (e.g., Naliboff et al., 2017;Naliboff et al., 2020;Duclaux et al., 2020), the initial plastic strain field is partitioned into 5 km coarse blocks which are randomly assigned 0.5 or 1.5. This method results in statistically random but pervasive damage using an adjustable wavelength (i.e., the block size) which allows for the adjustment of spatial distribution without changing the overall initial damage/strain. The data file (composition.txt) containing the initial distribution of plastic strain and additional compositional field data, as well as the python script to generate this data (composition.py) are located in the supplementary data set. Viscosity prefactor (A*) 8.57 x 10 -28 Pa -n m -p s -1 7.13 x 10 -18 Pa -n m -p s -1 6.52 x 10 -16 Pa -n m -p s -1 The extension rate (i.e., the prescribed outward velocity) is 2.5, 5 and 10 mm/yr (Fig. 1).

Compiling and Running ASPECT
The model in this study were run with ASPECT version 2. The input image can be a depth slice of any field documenting strain (e.g., the finite plastic strain, strain rate invariant, or components of the velocity vector). Our results use horizontal, vertical and diagonal) of pixels (E) (Dillencourt et al., 1992). Noise is filtered out by removing labelled pixels which are smaller than 20. In the case that branching remains, a post-processing script breaks up any remaining branches by locating euclidean distance transformation anomalies which arise as a result of bifurcation, and use the locality as a mask to split labels into two discrete fault segments.