Development of 3‐D Rift Heterogeneity Through Fault Network Evolution

Observations of rift and rifted margin architecture suggest that significant spatial and temporal structural heterogeneity develops during the multiphase evolution of continental rifting. Inheritance is often invoked to explain this heterogeneity, such as preexisting anisotropies in rock composition, rheology, and deformation. Here, we use high‐resolution 3‐D thermal‐mechanical numerical models of continental extension to demonstrate that rift‐parallel heterogeneity may develop solely through fault network evolution during the transition from distributed to localized deformation. In our models, the initial phase of distributed normal faulting is seeded through randomized initial strength perturbations in an otherwise laterally homogeneous lithosphere extending at a constant rate. Continued extension localizes deformation onto lithosphere‐scale faults, which are laterally offset by tens of km and discontinuous along‐strike. These results demonstrate that rift‐ and margin‐parallel heterogeneity of large‐scale fault patterns may in‐part be a natural byproduct of fault network coalescence.


Introduction
The majority of continental rifts form in previously deformed lithosphere. They preferably exploit preexisting zones of weakness such as mobile belts and avoid stronger regions such as cratons (Dunbar & Sawyer, 1988). The East African Rift and the North Sea Rift, for instance, are classical examples where rifts are guided by inherited weaknesses related to past collisional structures (Chorowicz, 2005;Phillips et al., 2016). These lithosphere-scale strength heterogeneities are capable of controlling rift kinematics on regional scale (Glerum et al., 2020); however, to which degree crustal fault networks are controlled by inherited weaknesses is still unclear.
Fault network characteristics have been elucidated through seismic studies of rifted margins. Recent observations inferred multiple characteristic phases of extension that partition rifted margins into distinct rift-normal structural entities (e.g., Cowie et al., 2005;Faleide et al., 2010;Franke, 2013;Mohriak & Leroy, 2013;Reston, 2009). As summarized in Péron-Pinvidic et al. (2013), these distinct phases of extension can be tied to three widely observed rifted margin domains: (1) The proximal domain corresponds to areas of low crustal thinning with deformation patterns characterized by Andersonian-type normal faults, including high angle segments that generate series of tilted blocks and half-grabens basins with moderate accommodation space.
(2) The necking domain marks a transition to more significant crustal thinning and an increase in accommodation space, with deformation typically focused along large-scale detachment faults. (3) The distal domain records the final stages of rifting prior to breakup and may include zones of hyperextended crust, exhumed mantle lithosphere, and magmatic products. While this conceptual model reconciles the common observations across many extensional systems, both rifts and rifted margins exhibit large variations in width, structural pattern, sedimentation, and magmatism. 2-D numerical investigations suggest that many of these variations reflect the initial rheological structure and the extension velocity (e.g., Brune et al., 2014;Huismans & Beaumont, 2011;Jammes & Lavier, 2016;Lavier & Manatschal, 2006;Naliboff & Buiter, 2015), but they fail to capture the rift-parallel component of deformation.
Indeed, many rifts and rifted margins exhibit significant along-strike variability in terms of structures, sedimentary architecture, and magmatic products (e.g., Autin et al., 2013;Corti et al., 2018;Osmundsen et al., 2016;Péron-Pinvidic et al., 2015, Tsikalas et al., 2008. These variabilities often lead to the definition of boundaries between distinct margin segments (for a detailed discussion, see Osmundsen & Péron-Pinvidic, 2018). However, their origin and geological significance remain obscure. The formation of such along-strike heterogeneity likely reflects multiple competing factors, including variations in the initial lithospheric structure (e.g., inheritance), reactivation of strength heterogeneities (e.g., shear zones), and spatial and temporal changes in the magnitude and obliquity of plate driving forces (e.g., Heine et al., 2013;Manatschal et al., 2014).
Within this context, here we use high-resolution 3-D simulations of continental rifting to study fault network evolution and coalescence. In particular, we conduct model runs with different rates of extension, magnitudes of strain softening, and geothermal gradients to assess patterns of along-strike fault interaction across a wide range of conditions. The modeling results demonstrate that distributed fault networks in the early phase of rifting give rise to rift-parallel structural heterogeneity, including lateral offsets of lithospheric-scale fault complexes as observed in many rifts and rifted margins (e.g., Osmundsen & Péron-Pinvidic, 2018). Furthermore, we demonstrate that the predicted fault network characteristics reproduce observational data across the range of tested physical parameters.

Numerical Methods
We model the 3-D thermal-mechanical evolution of continental extension using the open-source and CIG-supported finite element code ASPECT version 2.0.1-pre (Bangerth et al., 2018;Heister et al., 2017). Additional information on how to reproduce our experiments is contained within the supporting information (Text S2).
Velocity and pressure are solved for using the incompressible Boussinesq approximation, where the Stokes equations are defined as Above, u is velocity, μ is viscosity, _ ε is the deviatoric strain rate, Pis pressure, ρ is density, and g is gravitational acceleration.
Thermal evolution is modeled through the advection-diffusion equation: where C p is the heat capacity, T is temperature, t is time, K is thermal conductivity, and H is the rate of internal heating. Following the Boussinesq approximation, the density varies linearly as a function of the reference density (ρ 0 ), thermal expansivity (α), reference temperature (T 0 ), and temperature: Compositional fields are used to track and advect distinct lithologic domains (e.g., rock types) and other time-dependent quantities such as plastic strain. Each field requires solving an additional advection equation.
The constitutive behavior combines nonlinear viscous flow with brittle failure (see Glerum et al., 2018 and Text S1 for further details), with viscous flow following a dislocation creep formulation: Above, σ ′ eff is the second invariant of the effective stress, A is the viscous prefactor, n is the stress exponent, _ ε eff is the second invariant of the deviatoric strain rate (effective strain rate), Q is the activation energy, P is pressure, V is the activation volume, T is temperature, and R is the gas constant.
Brittle (plastic) behavior follows a Drucker-Prager yield criterion formulation, where the yield stress in 3-D is a function of the cohesion (C), angle of internal friction (ϕ), and pressure (P): To help localize deformation and account for geologic observations of strain localization, we reduce the cohesion and friction angle linearly as a function of accumulated plastic strain, which is derived from the second invariant of strain rate in regions undergoing plastic (brittle) deformation.
The procedure for calculating the viscosity at every point follows the viscosity rescaling method, which first compares the predicted effective stresses from viscous flow and plastic failure. If the viscous stress exceeds the plastic yield stress, the viscosity is reduced so the effective stress exactly matches the plastic yield stress (for further details, see Glerum et al., 2018).
Throughout the model domain we use quadratic elements (Q2) to solve the advection-diffusion equation for temperature and composition, while the Stokes equation is solved on elements that are quadratic for velocity and continuous linear for pressure (Q2Q1). To reduce the effect of diffusion of compositional field interfaces through time, averaging of material parameters between compositional fields follows a "maximum composition" approach (Glerum et al., 2018).
Nonlinearity introduced by the constitutive model is resolved using standard Picard iterations for the velocity and pressure with temperature and composition values from the current time step. We use a maximum time step of 20,000 (for an applied extension velocity = 10 mm/yr) or 40,000 (velocity = 5 mm/yr) years to limit numerical instabilities during advection and to improve the nonlinear convergence behavior.

Model Design
The numerical experiment design ( Figure 1) follows our goal of elucidating the relationship between distributed normal faulting and rift-parallel heterogeneity. The 3-D model domain is 500 × 500 × 100 km, respectively, in the horizontal (X,Y) and vertical (Z) directions ( Figure 1a).
The numerical resolution is 1.25 and 2.5 km, respectively, above and below 50 km depth. This transition reflects our focus on brittle failure, which always occurs above 50 km depth in our simulations. Significantly, we do not use ASPECT's adaptive mesh refinement capabilities to track brittle shear bands, as their width and the associated accumulation of finite deformation is inherently mesh dependent. We also note that additional models presented in the supporting information have a coarser resolution (2.5-5 km), but nevertheless align with our conclusions.
Deformation is driven by constant outward horizontal velocities prescribed on the model sides, which have no velocity constraints in the vertical directions. Inflow at the model base balances lateral outflow, while the upper boundary is a free surface, and the front and back walls are free-slip. The free surface advects according to a vertical velocity projection, which stabilizes the Stokes solution (Rose et al., 2017). The outward velocity (full extension rate) is 5 or 10 mm/yr.
We follow the method of Chapman (1986) and our previous work (Naliboff et al., 2017) to produce an initial temperature profile characteristic of the continental lithosphere. The surface heat flow is specified as 55 or 60 mW/m 2 , which we use to systematically calculate the temperature with depth using the thermodynamic properties of each layer within the lithosphere ( Figure 1a, Table 1). The resulting Moho and basal temperature, respectively, are 600°C and 1340°C (for a surface heat flow of 55 mW/m 2 ) or 700°C and 1540°C (for The lithosphere contains three distinct compositional (lithologic) layers ( Figure 1a, Table 1), representing the upper crust, lower crust, and mantle lithosphere. Respectively, the flow laws for these layers are wet quartzite (Gleason & Tullis, 1995), wet anorthite (Rybacki et al., 2006), and dry olivine (Hirth & Kholstedt, 2003). The initial friction angle and cohesion for all layers are, respectively, 30°and 20 MPa, which linearly weaken by a factor of 2 or 4 between plastic strains of 0.5 and 1.5 (Naliboff et al., 2017;Naliboff & Buiter, 2015). Combining these properties with the initial geotherm produces a strength profile containing two ductile layers within the crust (Figure 1b), which prevent lithospheric-scale detachment faults from forming near the onset of rifting.
We localize deformation in the center of the model domain by prescribing an initial zone of plastic strain (e.g., "Damaged" region in Figure 1a) with values randomized between the brittle strain weakening range (0.5-1.5). This randomized damage leads to an initial phase of distributed normal faulting (Figure 2), consistent with the findings of Naliboff et al. (2017). However, we note that a higher degree of initial plastic strain is required to localize deformation in the model center relative to Naliboff et al. (2017), which is likely due to the coarser resolution (1.25 km vs. 250 m) and our application of the randomization at quadrature points versus a single value over the entire element as in Naliboff et al. (2017). Without sufficient initial plastic strain for the given strain weakening magnitudes, simulations often exhibited extended pure shear flattening of the entire lithosphere for extended periods and in some cases propagation of deformation toward the model edges. However, even at the comparatively coarse resolution of 1.25 km grid-spacing, complex spatial-temporal finite deformation patterns develop along the margin axis as deformation progresses.
Notably, Duclaux et al. (2020) also successfully used a zone of randomized plastic strain to initiate distributed faulting in high-resolution 3-D simulations of oblique continental extension, albeit with a higher rate (strain bounds of 0.1 and 0.5) and magnitude (7.5X for friction, 5X for cohesion) of strain weakening. We  Table 1) for the upper crust, lower crust, and mantle lithosphere. The model length (along strike) is either 20, 100, or 500 km. A 200 km wide zone of randomized initial plastic strain (i.e., damage, lightly colored regions) weakens the lithosphere and localizes deformation in this region. Extension is driven by outflow on the side walls, which is balanced by inflow at the base and deformation at the top free surface. (b, c) Initial strength and temperature profile for a surface heat flow of 55 or 60 mW/m 2 and a strain rate of 1 × 10 −14 1/s.   (Gleason & Tullis, 1995); Lower crust: wet anorthite (Rybacki et al., 2006); Mantle: dry olivine (Hirth & Kholstedt, 2003). * The friction angle and cohesion decrease linearly by a factor of 2 or 4 between plastic strain values of 0.5 and 1.5.
further note that prescribing offset large-scale weaknesses can also produce significant margin heterogeneity (e.g., Le Pourhiet et al., 2017), but it does not capture the formation of heterogeneity due to initial distributed fault network evolution and effectively represents the lithosphere's response to longer wavelength inheritance. In contrast, our approach mimics the effect of inherited smaller-scale fault and fracture networks (km to tens of km) that reactivate at the start of rifting and coalesce into larger faults segments.

Model Limitations and Interpretation
Numerical results must be interpreted within the context of model limitations and design choices. Given that the largest simulations presented here require~500-1,000 cores to run over multiple days, we conducted a large number of 2-D and coarse 3-D simulations to ensure production runs would exhibit reasonable convergence behavior. In some cases, coarse-to moderate resolution 3-D tests unpredictably crashed when solving the Stokes linear system if the asthenosphere was included. Similarly, we found the compositional field tracking accumulated plastic finite deformation also became difficult to interpret at later rifting stages, in part due to diffusion. In a follow-up work we have employed and developed new methods to address these issues and also taken advantage of a recently developed Newton Solver (Fraters et al., 2019) to solve the nonlinear system of equations more efficiently and to a greater accuracy. However, here we purposefully stop our simulations prior to breakup as this is the point up to which we are confident in the numerical accuracy and applicability to modern observations. While the results are applicable to observations from both rifts and rifted margins, additional processes during the exhumation and breakup stage of rifting (e.g., magmatism and hydrotherm alteration) may introduce further along-strike heterogeneity.

Rift Localization Through Time
In order to provide a reference case for examining the effects of multiple physical parameters on rift-parallel heterogeneity, we first consider the temporal evolution of a model with an extension velocity of 10 mm/yr, an initial surface heat flow of 55 mW/m 2 , and a 2X brittle strain weakening (Figure 2). Consistent with both geologic observations and previous numerical studies, widely distributed normal faults accommodate the initial phases of deformation, which localize through time into a narrower rift zone containing larger faults.
Both active (Figure 2a) and finite (Figure 2b) deformation patterns after 8 Myr of extension clearly illustrate the initial phase of rifting, where the majority of active faults spans tens of km in length and are roughly evenly distributed across a 100-150 km wide rift zone ( Figure 2a). Significantly, along-strike heterogeneity already exists at this stage of rifting, with faults laterally offset (<10 to tens of km) along the length of the domain and localized rotation of fault orientations to accommodate overlapping segments. Brittle finite deformation (Figure 2b) largely mimics the active deformation patterns, although the width of individual fault segments appears wider due to fault migration and the deformation along transfer zones appears to connect distinct active structures along strike.
Continued extension to 12 Myr (Figures 2c and 2d) produces narrowing of the rift zone (~50 km wide) and along-rift linking of discrete fault segments, with multiple active faults spanning 100-200 km in length (Figure 2c). Nearly all fault segments exhibit some degree of curvature along their length, which often coincides with either fault deactivation, reactivation, or segmentation of specific structures. While active faults appear as discrete localized structures, accumulated deformation ( Figure 2d) reveals a complex pattern of relatively wide (up to~10 km) finite shear zones, which contain significant overprinting of the initial rifting phase. Despite this overprinting, along-strike offsets between discrete shear zones can still be distinguished. The processes of fault deactivation, activation, segmentation, rotation, and rift localization through which these structures develop over time is clearly illustrated when viewed in animation sequences (Movies S1a-S1c).
Our results provide critical constraints on possible model configurations. We find that the development of along-strike segmentation requires a model extent of at least 100 km in the rift-parallel direction, unless additional initial structural heterogeneity is implemented. In a model extending only 20 km in the y direction (Movies S2a-S2c), the deformation patterns effectively mimic those of a 2-D simulation (i.e., in-plane deformation), while extending the model length to 100 km (Movies S3a-S3c) enables some degree of along-strike segmentation. The exact longitudinal length scale required to produce along-strike segmentation likely depends on a number of factors, including brittle layer thickness, fault strength, and structural inheritance.

Kinematic and Rheological Controls on Rift Structure
To assess the sensitivity of observed rift heterogeneity development to our model parameters, we conducted a series of simulations that examine a slower rate of extension, larger magnitude of brittle strain weakening, and a higher geothermal gradient (Figure 3). Reducing the extension velocity from 10 mm/yr (Figures 2 and  3a, Movies S1a-S1c) to 5 mm/yr (Figure 3b, Movies S4a-S4c) has only minimal impact on the first-order rift structure after equivalent amounts of total extension, with the exception of higher effective fault viscosities at lower velocities and less prevalent transfer structures between rift distinct fault segments.
In contrast, maintaining a velocity of 5 mm/yr while increasing the magnitude of brittle strain weakening from 2 (Figure 3b, Movies S4a-S4c) to 4 (Figure 3c, Movies S5a-S5c) produces a distinct change of the rift regime, with two separate rift basins after 20 Myr of extension. These rift basins initiate as crustal faults at the edges of the initial plastic strain zone, which then evolve to two full rift basins separated by a stable, 100 km wide horst block. While the faults at the outer edges of each rift basin run nearly continuously along the model edges, the inner-rift border faults exhibit significant along-strike segmentation that is defined by sharp changes in fault orientation. The majority of deformation is accommodated along the basin-bounding faults. Nevertheless, second-order intrarift faults with a spacing of tens of km exhibit similar patterns of along-strike heterogeneity.
Maintaining a brittle strain weakening factor of 4 and increasing the initial surface heat flow to 60 mW/m 2 produces another significant change in rift style, with deformation distributed across a 300 km wide domain (Figure 3d, Movies S6a-S6c). The restriction of faulting to the upper crust and the relatively flat Moho reflect a predominantly ductile deformation with pronounced crust-mantle decoupling (Buck, 1991). Higher deformation rates (lower viscosities) occur in the central region with randomized initial strength perturbations. Faults in this region exhibit similar qualitative patterns of along-strike segmentation, linkage, and rotation as previous models, albeit with larger fault spacing (up to 25 km between major faults). Faults in the outer segments of the rift deform at a lower rate and also show greater continuity and less segmentation along their length, with select faults comprising nearly the entire model length.
Additional testing at lower resolutions (Movies S7-S12) for a wider range of extension velocities (5-20 mm/yr) reveal qualitatively similar results. However, certain results also highlight the mesh-dependency of plasticity, as the multiple rift basin mode of deformation transitions to a continuous rift basin with significant changes in orientation along strike (Movies S8a-S8c). This transition reflects as (c) but with an initial surface heat flow of 60 mW/m 2 . All models are depicted after the same amount of total extension. As in Figure 2, the color scale is selected to highlight regions of active deformation, although here the viscosity contrast spans 2 orders of magnitude and thus reveals smaller faults and transfer structures masked by the strain-rate cutoff in Figures 2a and 2c. The spatial extent, scale bars, and density contour for the Moho are the same as those described in Figure 2. See supporting information for model animations.

Geophysical Research Letters
NALIBOFF ET AL.
that narrower faults at higher-resolutions enhance strain localization, which allows deformation to initially localize at the damage zone boundaries and eventually produce two rift basins (Figure 3c, Movies S5a-S5c). Recently proposed methods (e.g., Duretz et al., 2020) provide a solution for producing mesh-independent and iteratively stable plastic behavior, which can be incorporated into future simulations.
Furthermore, we again note the importance of resolution and strain-softening within the context of a recently published study (Duclaux et al., 2020), which used a similar approach to look at normal fault development under different extension obliquities. Duclaux et al. (2020) demonstrate that fault length systematically decreases with increasing obliquity, but at low obliquities (16°) and equivalent extension magnitudes obtain fault lengths (>500 km in some cases) significantly greater than those obtained here (tens to <300 km). We suspect this is both a function of the lower strain softening magnitude used here relative to Duclaux et al. (2020), as well as the fact that oblique extension tends to form relatively ordered en echelon fault patterns. The mechanical thickness of the crust and additional numerical decisions likely also affect the average fault length. Similarly, numerous model parameters and design choices likely also affect the average lateral offsets between fault segments. Here, lateral offsets are typically on the order of 10-50 km, while models (e.g., Balázs et al., 2018;Le Pourhiet et al., 2017) with initial localized inhomogeneities (i.e., larger-scale inheritance) produce offsets greater than 100 km in some instances.

Geologic Interpretation
Our numerical simulations show that along-strike structural variability arises without any segmentation prescribed through explicit inherited structure or changes in rift propagation, obliquity, or linkage. These results highlight that a first-order margin-scale segmentation may be naturally generated in all extensional settings without a need for specific inherited lithospheric anisotropies.
The evolution of our numerical fault networks follows the classical interpretations derived from the Gulf of Suez and the North Sea (Gawthorpe & Leeder, 2008, Cowie et al., 2005, where initially small-scale faults grow and interact before they successively link to form a mature fault network with large-scale normal faults. The observed lengths of these mature normal faults reach 100-150 km in East Africa (Accardo et al., 2018;Corti et al., 2018;Hodge et al., 2018), the Baikal Rift (Petit & Déverchère, 2006), and the North Sea (Bell et al., 2014), which correlates well with the maximum lengths reached in our models. A comparison of fault displacement versus length data between our models and natural systems shows close similarity in terms of values and variability ( Figure 4). However, no clear correlations or trends appear between the varied parameters and the displacement-length relationship. This may reflect that the initial plastic strain slightly overprints and skews the displacement-length relationship, distinct trends manifest in a time-dependent analysis, or additional yet to be determined factors.
Inherited structures are proven to influence rift evolution in various settings. In portions of the North Sea, for instance, orogenic-collapse related shear zones partly controlled the early rift phases (Lenhart et al., 2019;Phillips et al., 2016). However, the influence of this inheritance decreases drastically during later rift phases, and the final rift geometries are often interpreted to be more related to the thermal and structural effects of the prior phases of rifting deformation, than to any initial inheritance. Indeed, observations from the Mid-Norwegian margin (Péron-Pinvidic & Osmundsen, 2018) suggest that along-strike margin segmentation reflects changes in the extent of different margin domains (e.g., necking, distal), which are in turn defined by distinct bounding faults. In other words, along-strike changes in the structural role of the fault network (e.g., necking vs. exhumation) define distinct margin segments. While our models only reach the necking stage of rifting, the initiation of domain segmentation is evident from the rift-parallel fault heterogeneity and arises without requiring complex preexisting structures.
Finally, our models show that the overall 3-D geometries are determined by complex fault evolution, which includes multiple phases of activity (fault-initiation, fault-activation, fault-linkage, fault-crosscutting, faultdeactivation, and possibly fault-reactivation; as shown in 2-D by Naliboff et al. (2017). The 3-D approach highlights that the geometries are furthermore dependent on the length scale of the faults and their linkage patterns, which produces segmentation on length scales ranging from tens to hundreds of km (i.e., segment lengths) and average lateral offsets on the order of~10 to 50 km (Figures 2 and 3).

Conclusions
The 3-D thermal-mechanical numerical models presented in this study reveal that distributed normal fault networks naturally evolve to produce rift-parallel structural heterogeneity over length scales of hundreds of km. In the case of moderate fault weakening, lithospheric-scale faults are offset by tens of km along the margin length. Stronger fault weakening may lead to the formation of multiple rift basins, while higher geothermal gradients promote formation of wide and distributed fault networks. These findings suggest that along-strike variations in rift and rifted margin structure may naturally arise from the complex evolution of fault networks, without any large-scale prerift anisotropies, with further heterogeneity arising from variable lithospheric structure and plate driving forces.