Influence of Zones of Pre‐Existing Crustal Weakness on Strain Localization and Partitioning During Rifting: Insights From Analog Modeling Using High‐Resolution 3D Digital Image Correlation

Pre‐existing crustal structures are known to influence rifting, but the factors controlling their influence remain poorly understood. We present results of digital image correlation that allows for the surface strain analysis of a series of analog rifting experiments designed to test the influence of the size, orientation, depth, and geometry of pre‐existing crustal weak zones on strain localization and partitioning. We apply distributed basal extension to crustal‐scale models consisting of a silicone weak zone embedded in a quartz sand layer. We vary the size and orientation (α‐angle) of the weak zone with respect to the extension direction, reduce the thickness of the sand layer to simulate a shallow weak zone, and vary the geometry of the weak zone. Our results show that at higher α‐angle (≥60°) both small‐scale and large‐scale weak zones localize strain into graben‐bounding (oblique‐) normal faults. At lower α‐angle (≤45°), small‐scale weak zones do not localize strain effectively, unless they are shallow. In most models, we observe diffuse, second‐order strike‐slip intra‐graben structures, which are conjugate and antithetic under orthogonal and oblique extension, respectively. Generally, the observed spectrum of rift faulting styles (from discrete fault planes to diffuse fault zones, from normal to oblique and strike‐slip) highlights the sensitivity of rift architecture to the orientation, size, depth, and geometry of pre‐existing weak zones. Our generic models are comparable to observations from many natural rift systems like the North Sea and East Africa, and thus have implications for understanding the role of structural inheritance in rift basins.

Most previous studies have focused on first-order attributes like the nature of pre-existing weak structures such as discrete versus distributed structures (e.g., Bellahsen & Daniel, 2005;Bonini et al., 2015;Deng et al., 2017;Tong et al., 2014), and the orientation (obliquity) of the inherited structure with regards to the extension direction (e.g., Agostini et al., 2009;Autin et al., 2013;McClay & White, 1995;Michon & Sokoutis, 2005;Molnar et al., 2019;Withjack & Jamison, 1986). With the progress made in seismic imaging of deeply buried structures located in crystalline basement, higher-order attributes such as the size, depth, and shape of pre-existing structures are becoming evident (e.g., Kolawole et al., 2021;Osagiede et al., 2020;Phillips et al., 2016;Wrona et al., 2020); such attributes have received little attention in modeling studies, thus little is known about how they influence the evolving rifts, especially in the context of pre-existing structures like crustal shear zones (here referred to as weak zones). Shear zones represent weak zones in nature because they are often characterized by lower shear strength (due to phyllosilicate mineralization and the preferential orientation of shear fabrics) compared to their bounding undeformed host rocks (e.g., Ramsay, 1980;White et al., 1980).
In this study, we investigate the controls that pre-existing crustal weak zones have on the strain localization process and development of rift-related structures during extension. We achieve this through a series of extensional analog experiments that test how the (a) size, (b) depth, and (c) geometry of pre-existing crustal weak zones (specifically shear zones) affect their propensity to influence younger rift faults. We deploy state-of-the-art stereoscopic (3D) digital image correlation (DIC) technique that allows us to quantitatively assess the evolution of the model surface deformation and structural pattern at high resolution. Our results are of generic significance and have implications for understanding how pre-existing weak zones selectively influence younger rift faults in natural rift systems.

Experimental Setup
The experimental apparatus can be considered a pure shear extension experiment ( Figure 2). A basal neoprene foam block with dimensions of 50 × 50 × 12 cm is first compressed by 5 cm (i.e., down to an initial width of 45 cm) between a fixed and mobile wall, with open side walls, before placing the rock analog materials on top of it. Subsequently, we start the model run by slowly moving the mobile wall at a constant rate of 0.005 mm/s, and releasing the compressed foam by a total of 5 cm, translating theoretically into c. 11% of Representative plan-view sketch of our weak zone geometry oriented at α = 90° to the extension direction represented by the black arrow, and the focus area used for the time-series digital image correlation (DIC) processing. (c) Representative cross-sectional view, indicating our model layering and distributed deformation at the base of the models. (d) Experimental deformation rig-setups, recording-setups, and processing-setups. (e) Hypothetical rheological layering and associated strength profile for the homogeneous (pristine)-crust and weakened-crust, respectively, in our experimental series A (not drawn to scale). (f) The 3D geometry and dimensions of the weak zone used in the different models to simulate a range of natural weak zone geometries (note the map view insert of the weak zone in Model E). extension during each model run. Although this small magnitude (11%) of extension is a limiting constraint of the setup, it however allows us to investigate the early stage of rifting in detail. Due to some localization of extension at the model boundaries, the effective extension in the model domain is a bit less (c. 9 ± 1%) and distributed homogeneously as verified by benchmarks ( Figure A1). The use of a basal foam in analog model studies has the advantage of simulating distributed extension with a constant basal velocity gradient ( Figure A1), compared to the use of rigid base plates that typically serve to strongly localize basal extension (e.g., Schlagenhauf et al., 2008;Zwaan et al., 2019).
A side effect of using elastic materials to impose basal kinematic boundary conditions is the association of transverse contraction along with the longitudinal extension. This "Poisson effect" may lead to a switch in simulated tectonic setting from crustal extension to strike-slip if the Poisson's ratio is close to 50%, as often happens with a basal rubber sheet, for example (Bahroudi et al., 2003;Zwaan et al., 2019). However, in our model setup, lateral contraction of the basal foam block is c. 16% of its longitudinal extension (i.e., 0.8 cm contraction associated with 5 cm extension, Figure A1) and therefore suitable for simulating crustal extension with only a minor contribution of transverse contraction.
In this study, our focus is on the surface deformation of the brittle crust; we do not explicitly model a ductile lower crust, which has been sufficiently addressed by previous analog modeling studies (e.g., Allemand & Brun, 1991;Bellahsen et al., 2003;Brun, 1999;Zwaan et al., 2019). Our reference model setup therefore simulates an old and cold stable (e.g., Brun, 1999;Zwaan et al., 2019) or a highly coupled (e.g., Dyksterhuis et al., 2007) crustal/lithospheric setting. That is, our model setup does not capture the vertical rheological stratification that may characterize the crust in relatively young and/or hot crustal settings, but rather focuses on the lateral mechanical anisotropy induced by relatively local pre-existing crustal weak zones (e.g., shear zones; see Vauchez et al., 1998). We also acknowledge the presence of lithospheric-scale and/ or mantle weaknesses that may have complex effects on rift structures (see e.g., Molnar et al., 2019;Zwaan et al., 2021a); however, the focus of the present work is on crustal weaknesses.

Rock Analog Materials
As a rock analog material for the brittle upper crust, we used a mix of natural-colored and a few percent (4.7 vol%) of black-colored, dry, quartz sand (type G12, Rosenau et al., 2018). This mixture provides an appropriate visual contrast that allows for the digital correlation of recorded images (e.g., White et al., 2001). The grain size is 100-400 μm, with an average of 240 μm. The bulk density of the sand, which is sieved onto the basal foam block during model preparation is about 1700 kg/m³. The sand exhibits a frictional-Coulomb  Rosenau et al. (2018). b Rudolf et al. (2016). plastic behavior, with static and dynamic friction coefficients of 0.69 and 0.55, respectively, and a cohesion in the order of a few tens of Pascal (see Table 1; Rosenau et al., 2018). Quartz sand like this has been widely used to represent brittle upper crustal rocks in numerous analog models (Adam et al., 2005;Del Ventisette et al., 2019;Klinkmüller et al., 2016;Lohrmann et al., 2003;Panien et al., 2006;Ritter et al., 2016;Schellart & Strak, 2016;Schreurs et al., 2006Schreurs et al., , 2016. We use viscoelastic objects made of Polydimethylsiloxane (PDMS) or silicone oil to represent pre-existing weak zones, simulating kilometer-scale shear zones within the brittle crust. While several methods permit simulation of weak zones into sand-based analog models (e.g., silicon seed; Le Calvez & Vendeville, 2002;Zwaan et al., 2019, pre-cut;Bellahsen & Daniel, 2005;Tong et al., 2014;Zwaan et al., 2021a), a major advantage of using silicone oil is that it can be easily molded to produce a range of sizes and shapes that reflects a diverse range of natural weak zone geometries. The silicone oil has a density (ρ) of c. 960 kg/m 3 , and a zero-shear viscosity (ɳ) of c. 2.24 × 10 4 Pa s (Table 1; Rudolf et al., 2016). As a Maxwell viscoelastic fluid, it has a linear viscous rheology (Newtonian) at relatively low strain rates as realized in our experiments, changing to a nonlinear viscous rheology at much higher strain rates than achieved in this study (Rudolf et al., 2016). We choose this material for the weak zone because it can respond to the extension applied at the model base while maintaining its shape and height over the experimental time scale.

Scaling
To achieve analog models that are applicable to nature, an adequate geometric (length), kinematic (time), and dynamic (stress) scaling between model and nature must be established (Hubbert, 1937;Mulugeta, 1988;Ramberg, 1981; see also review by Schellart and Strak [2016]). For proper dynamic scaling and similarity in the brittle regime, the following equation, relating material strength (here cohesion C) and gravitational stresses (or overburden pressure), should be satisfied: where C*, ρ*, g*, l*, and σ* are the model versus nature ratios (called "scaling factors") for cohesion, density, gravity, length, and stress, respectively. Since the experiment is carried out under normal gravitational field in the laboratory, g* is 1, and Equation 1 reduces to: With setting up the model at laboratory scale we impose a geometric scaling factor (length ratio l*) in the order of 10 −5 to 10 −6 , that is, 1 cm in the model represents c. 1-10 km in nature. The density ratio ρ* is c. 0.7 (assuming an average density value of 2,400 kg/m 3 for sedimentary and crystalline rocks in the upper crust). Therefore, dynamic similarity (from Equation 2) is achieved with a material satisfying a cohesion ratio (C*) that equals a stress ratio (σ*) of 7 × 10 −6 to −7 . Accordingly, our quartz sand with an average cohesion of a few tens of Pascal will be able to simulate brittle crustal rocks with cohesion in the range of approximately 10-100 MPa. This range is well within the range of cohesion values for brittle crustal rocks (e.g., Klinkmüller et al., 2016;Schellart, 2000;Schellart & Strak, 2016, and references therein).
We note that even if we consider that our models are brittle-dominated and display time-independent behavior, the rheology of the weak zone material needs consideration in the context of its viscous behavior, as the expected strain rate defines how weak the crust is in the weak zone region compared to the normal crust. The viscous strength of any flowing material is the viscosity times the strain rate. For a strain rate of 10 −5 /s as in our experiments (see Appendix B), it follows that the weak zone material strength is about 0.1 Pa, which is 100-10,000 times lower than the strength of the quartz sand layer that increases linearly with depth from cohesion values near the surface to about 700-800 Pa at the base of the 4 cm thick models. Therefore, we consider the contribution of the weak zone itself to the integrated strength of the model crust to be negligible. The integrated model crustal strength (i.e., the area beneath the strength profile in Figure 2e) in the area underlain by the weak zone is proportional to the squared thickness of the layer above the weak zone (see Zwaan et al., 2020, for a geometric derivation of this scaling). Accordingly, for 1-cmthick and 2 cm-thick weak zones in an overall 4 cm thick layer model crust, the integrated strength in the weak zone area is reduced to ca. 56% and 25%, respectively, of the integrated strength of the normal (pristine or sand-only) model crust.

Model Series Design
We present a reference model (R), and 10 main models grouped into five main series (A-E; Table 2). The reference model consists of a 4-cm-thick sand layer (without a pre-extension weak zone) directly coupled to the extending foam. Series A and B focused on testing the extent to which the size (small vs. large, respectively) and the orientation of pre-existing weak zones control the pattern of strain localization and partitioning in the overlying cover. In Series C, we use the thickness of the overburden brittle layer as a proxy to examine the influence of shallow versus deep burial of pre-existing weak zones on the extent of strain localization in the overlying cover. Series D and E test the influence of the overall 3D geometry of pre-existing weak zones on deformation patterns during crustal extension. Details of the (a) orientation angle α, of the weak zone with respect to the applied extension direction, measured clockwise (positive angle) or anti-clockwise (negative angle) ( Figure 2b) (b) dimension and geometry of the weak zone, and (c) integrated strength in the area underlain by the weak zone, with respect to the normal model crust, for each model are provided in Table 2. In this study, models where the weak zone is orthogonal to the extension direction (i.e., α = 90°) are referred to as orthogonal extension models, whereas those with oblique weak zone versus extension direction angles (i.e., α = 60° and 45°) are referred to as oblique extension models. Note that in Series E model, the weak zone was curvilinear implying that about half of the weak zone is orthogonal to the extension direction (α = 90°) and half is oblique to the extension direction (α = 60°).

Surface Deformation Monitoring
We used a stereoscopic pair of two 12-bit, 29-megapixel monochrome CCD (charge-coupled device) cameras (LaVision Imager XLite 29M) mounted c. 1 m above the model surface with oblique viewing angles ( Figure 2d). Images were recorded at an image frequency of 0.05 Hz, corresponding to one image every 20 s, which, given the extension velocity (0.005 mm/s) of our models, represents one image per 0.1 mm of wall displacement.
To process the recorded images, we deployed a DIC technique, also known as particle imaging velocimetry (PIV), that allows us to monitor the surface deformation at high spatial and temporal resolution (e.g., Adam   White et al., 2001). We used commercial LaVision Davis 10 software employing a least-squares-method (LSM) to calculate the three-component deformation field between successive images. DIC processing was done at increments comparing every 10th image, representing 1 mm of sidewall displacement. We find that this interval is not only computationally convenient, but also provides adequate information on the evolution of strain within our models at a high spatial and sufficient temporal resolution. To avoid boundary effects, and further optimize computation time, DIC processing on the raw images was restricted to the central part of the model (area of interest) by means of a rectangular mask (see Figures 2b and 3).
In the DIC results, we present in this study are based on time-series of the 3D cumulative surface displacement field, where incremental displacements are summed up in a Lagrangian reference frame. From the 3D surface displacement field, the in-plane (2D horizontal) cumulative normal strain (ε n ) and shear strain (ε s ) were derived. The ability to decompose strain into normal and shear components highlight the key strength of the DIC analysis deployed in this study, in comparison to the numerous previous analog models where deformation analysis has been mostly based on qualitative (rather than quantitative) visual inspection of model surface deformation at a much lower spatial and temporal resolution. The computed strain distribution allows us to quantitatively assess the geometry and pattern of strain distribution of extension-related structures in our models at unprecedented resolution.

Model Results
The results are presented first in Sections 3.1-3.4 based on the different experimental series. Within each section, we first present the final stage (c. 11% extension) of the model evolution by means of the surface fault pattern and the DIC-derived surface normal strain (ε n ). We also present the DIC-derived surface shear strain (ε s ) of the reference model, and Series A and B models, highlighting the oblique nature of first-order and secondary strike-slip structures. Finally, we provide cumulative model surface subsidence profiles illustrating the temporal evolution of rift topography. In Section 3.5, we briefly highlight the temporal aspects of structural evolution, and in Section 3.6, we quantify the obliquity of normal faulting in selected models. All data underlying this study and additional results of the DIC processing (    = graben axis, t = time in seconds. Note that an apparent "northern" direction is imposed for the purpose of description.

Reference Model R: Brittle-Only (Homogeneous, Isotropic Layer, No Weak Zone)
The final stage (c. 11% extension) of the reference model surface is characterized by the occurrence of faults only along the upper and lower edges of the model due to boundary effect resulting from the discontinuity between the sand layer and the model sidewalls ( Figure 3a). There is no localization of deformation in the central portion of the model. Instead, distributed, diffuse normal strain bands (c. 10% strain), orthogonal to the extension direction develop (Figure 3b) with no associated shear strain (Figure 3c). We interpret these diffuse strain bands as products of the distributed basal extension boundary condition (i.e., brittle layer directly on basal foam) in our model setup, and are equivalent to the distributed faulting reported in previous analog models by both Schlagenhauf et al. (2008) and Zwaan et al. (2019). Diffuse strain is plausibly a consequence of the larger grain size of the sand used in our setup (our mean grain size of 240 μm vs. their c. 120 μm) and the higher amount of extension applied in their models (our c. 11% vs. their ≥13%). Apart from the diffuse strain bands, there is no significant shear strain localization in the model (Figure 3c). The result of the reference model is important in this study because it allows us to differentiate the inherent model observations that are due to (a) the absence of a pre-existing weak zone and (b) edge and/or basal boundary conditions associated with the overall model setup. The diffuse strain bands are not considered or interpreted further in this study.

Experimental Series A and B: Small-Sized Versus Large-Sized Weak Zones With Variable Orientation
The model surface deformation after the final stage (c. 11% extension) of Models A1 (small weak zone, α = 90°), A2 (small weak zone, α = 60°), B1 (large weak zone, α = 90°), B2 (large weak zone, α = 60°), and B3 (large weak zone, α = 45°) are characterized by surface displacement along normal faults that bound model rift-related graben (herein referred to as graben-bounding faults) and which form directly above and parallel to the underlying weak zone (Figures 4a, 4f, 5a, 5f, and 5k). The graben-bounding faults are well constrained as localized, high normal strain zones (≥20%, Figures 4b, 4g, 5b, 5g, and 5l). The final stage of Model A3 (small weak zone, α = 45°) differed from other models in that a rift-related graben structure did not develop (Figures 4k and 4l). This suggests that the weak zone in this model does not weaken the brittle layer sufficiently to allow for strain localization and surface displacement on graben-bounding faults under moderate obliquity.
In the orthogonal extension models (A1 and B1), the graben-bounding faults are through-going with an overall orientation that is perpendicular to the extension direction, except around the fault tips where the fault geometry slightly deflects outward (Figures 4a, 4b, 4d, 5a, 5b, and 5d). The fault tip deflection is most likely due to the absence of the weak zone on the outer edges of the model, thereby making the sand layer thicker and the graben slightly wider at the edge. In Model B1, intra-graben faults that are parallel but antithetic to the main graben-bounding faults also develop, resulting in the formation of a central horst separating double marginal grabens (Figures 5a-5d). In the oblique extension models A2, B2, and B3, the resulting graben-bounding faults are generally oblique to the extension direction, except at the lateral tips of the faults where their geometry deflects, becoming near-orthogonal to the extension direction (Figures 4f-4i, 5f, 5g, 5i, 5k, 5l, and 5n). The normal strain component shows that in some of the models, specifically models B1 and B2, there are zones of approximately zero-strain (little or no deformation) in the immediate footwall and/or hangingwall of the graben-bounding faults (dark blue zones in Figures 5b and 5g). These zones are here referred to as strain shadows (Figures 5d and 5i) (e.g., Ackermann & Schlische, 1997;Cowie, 1998;Deng et al., 2017;Gupta & Scholz, 2000;Soliva et al., 2006).
The shear strain component of the deformation provides insight into the role of oblique-normal-slip and strike-slip strain components during the model run. In the orthogonal extension models (models A1 and B1), the pattern of the shear strain distribution along the graben-bounding faults is largely chaotic in the central part of the faults, and more systematic only at the fault tips (Figures 4c and 5c). The chaotic shear strain pattern reflects pure-dip slip in the central part of the graben-bounding faults. Some systematic oblique-slip is indicated at the fault tips which is consistent with lateral contraction imposed on the model by the basal foam. Within the graben (especially Model A1), there are poorly developed, diffuse shear strain zones (c. ±0.3%). These were not obvious in the normal strain pattern indicating their pure strikeslip nature. We interpret these as sets of conjugate strike-slip faults formed in response to the small lateral contraction imposed by the basal foam in an overall pure shear regime (see appendix; Figure C1). These poorly developed conjugate shears within the graben are typically not observed in natural rifts, and thus are considered here as basal boundary effects, and are not considered further in this paper.
In oblique extension models A2, B2 and B3, the overall shear strain distribution is generally similar and consistent with strain partitioning in transtensional kinematics (see appendix; Figure C1). The graben-bounding faults are consequently characterized by a minor component of dextral shear motion (Figures 4h, 5h, and 5m). Within the graben, en échelon, sigmoidal zones of sinistral (i.e., antithetic) shear developed (especially obvious in models A2 and B3). These intra-graben structures did not accommodate vertical surface displacement as they are not apparent in the normal strain component of the model deformation, indicating that they are purely strike-slip structures. We interpret them as antithetic Riedel (R') shears related to the simple shear component of transtension. These R' shears are even seen in Model A3 (Figure 4m), although no prominent graben developed (see Figures 4k and 4l). The overall en échelon arrangement of the R' shears is parallel to the underlying weak zone, whereas the individual R' shears are near-orthogonal to the extension direction.
The width of the rift-graben is approximately the same (c. 4 cm) for models A1 and B1, irrespective of the size of the underlying weak zone (compare Figures 4a and 5a), whereas it is slightly narrower with decreasing α-angle (compare Figures 5a-5k). The width of the rift-graben is unaffected by the height of the weak zone (which determines the thickness of the brittle layer above the weak zone) because the graben width is not only dependent on the thickness of the brittle layer above the weak zone and the dip of the graben-bounding faults, but also on the width of the weak zone (Corti, 2004) (see appendix; Figure D1). Time-series cross-sectional profiles of the models show that the subsidence is symmetric within the rift-graben (Figures 4e, 4j, 5e, 5j, and 5o), except for Model A3 where the profiles are mostly flat because no prominent graben developed in the model (Figure 4o). Whereas the model surface outside the model rift-graben generally subsides at rather constant rates (equidistant horizontal profile segments) consistent with distributed thinning, the rift-graben floors show indications of variable subsidence rates over the model run.
For example, in the rift-graben floor of Model A1, the subsidence rate shows a gradual increase from the onset of the model run until the ca. 8000 s subsidence profile line, then decreasing from 8000 up to 10,000 s line (Figure 4e). The absolute and maximum magnitudes of subsidence vary between models, generally increasing with increasing size of the underlying weak zone (compare Figures 4e and 5e) and decreasing with decreasing α-angle (compare Figures 5e, 5j, and 5o). The overall subsidence is positively correlated to the maximum vertical displacement accommodated by the graben-bounding faults. Interestingly, the subsidence profile also show zones of relative footwall uplift with magnitudes that positively correlate with amount of vertical displacement of the associated graben-bounding faults (Figures 4e, 5e, 5j).  = time in seconds. Note that an apparent "northern" direction is imposed for the purpose of description.

Experimental Series C: Thin Brittle Layer
We here present an additional model (Model C; small weak zone, α = 45°) with a brittle layer thickness of 1.5 cm. That is, the only difference between the setup of models C and A3 is the thickness of the brittle layer.
In the final stage of extension, Model C is characterized by surface displacement and strain localization along graben-bounding faults that form directly above and parallel to the underlying weak zone (Figures 6a-6c). This observation provides insight on the effect of the depth of the weak zone, since Model A3 with a thicker overburden brittle layer did not localize deformation (compare Figures 6a and 4k). The graben structure in Model C is c. 1 cm wide and narrower than those formed in Series A. Furthermore, the normal strain distribution along the graben-bounding faults is characterized by several local maxima (Figures 6b and 6c). The location of these strain maxima coincides with the termination zone of the diffuse strain bands on the obliquely oriented graben-bounding faults, and thus could be considered a boundary effect. Time-series cross-sectional profiles of the model show that the subsidence is symmetric within the rift-graben (Figure 6d).  Figure 7 shows the results of the final stage (c. 11% extension) of three models (models D1, D2, and E) with different weak zone geometries, simulating a range of cross-sectional and plan-view geometries that pre-existing weak zones may exhibit in nature. Overall, the deformation patterns of these models are more complex compared to the models in Series A-C, in which we used geometrically simple, cuboid-shaped weak zones.

Experimental Series D and E: Variable Weak Zone Geometry
In Model D1 with an upward-pointing triangular prism-shaped weak zone, the normal strain distribution across the graben structure is quite asymmetrical (Figures 7a and 7b). The "southern" boundary of the graben is characterized by a through-going graben-bounding fault, with a zig-zag geometry consisting of relatively longer segments that are parallel to the underlying weak zone, and shorter jogs that are subparallel to the extension direction ( Figure 7c). The zig-zag geometry reflects growth by the (at-surface) linkage of initially isolated weak zone-parallel fault segments. However, the strain distribution pattern along the "northern" boundary of the graben is more complex. Unlike the southern boundary, the northern boundary is characterized by a wider zone of faulting, with a dominant right-stepping en échelon faulting style (Figures 7a-7c). The en échelon fault set consists of individual fault segments which are oriented sub-perpendicular to the extension direction but is overall aligned parallel to the underlying weak zone. Time-series cross-sectional profiles show asymmetric subsidence within the rift-graben (Figure 7d).
In Model D2 with a broader, deformed semi-cylindrical-shaped weak zone, normal strain distribution on the borders of the graben structure is near-symmetrical (Figures 7e and 7f). This differs from the asymmetric strain distribution observed in Model D1. In contrast to the discrete graben-bounding faults observed in Series A-C models, the graben in Model D2 is bounded by relatively wide (c. 2.5 cm) deformation (fault) zones on both margins (Figures 7e and 7f). The wide fault zones are characterized by variable fault styles, including; subparallel fault sets, anastomosing faults (sensu Peacock et al., 2016), and right-stepping en échelon faults (within the "northern" graben-bounding fault zone, similar to Model D1) (Figure 7g). Within the graben, discrete intra-graben fault develops parallel to the underlying weak zone. The intra-graben fault accumulates higher amount of deformation (c. ≥20% strain) compared to the graben-bounding fault zone (Figure 7f). The time-series cross-sectional profiles of the graben subsidence are near-symmetrical, unlike that of Model D1 (Figure 7h). Also, the symmetry of the subsidence profile is characterized by a "dual slope" and differs from the "single slope" symmetry of Series A-C models (compare Figure 7h with profiles in Figures 4-6).
Model E contains a weak zone that has a deformed semi-circular cross-section, but curvilinear planform geometry. The normal strain distribution is similar to Model D2, with wide (c. 2.5 cm) fault zones bounding the rift-graben structure (Figures 7i and 7j). The wide fault zones also exhibit variable fault styles (Figure 7k). Like Model D2, the intra-graben fault is parallel to the underlying weak zone, and accumulates higher amount of deformation compared to the graben-bounding fault zone ( Figure 7j). However, the development of the intra-graben fault is only limited to the eastern portion of the graben structure that is oblique to the extension direction, and does not extend to, or develop in, the western portion that is perpendicular to the extension direction. The time-series cross-sectional profiles of the graben subsidence are symmetrical, and somewhat similar to that of Model D2 (Figure 7l).

Temporal Evolution of Model Structures
To illustrate the temporal evolution of the model structures, we briefly present the early (2.2% extension), intermediate (6.7% extension), and final (11.1% extension) stages of normal strain evolution of some Series A and B models (Figure 8). The timing of faulting and amount of strain accommodated by the faults varies in these models. For example, whereas the graben-bounding faults in Model A1 (small weak zone, α = 90°) accommodated only c. 5% strain at the early stage of extension, its equivalent in Model B1 (large weak zone, α = 90°) already accommodated up to 15% strain by the same stage (compare Figures 8a and 8c). Similarly, the graben-bounding faults in Model B2 (large weak zone, α = 60°) accommodated more strain than in Model A2 (small weak zone, α = 60°) at the same stages (compare Figures 8b and 8d). Furthermore, comparing models A1 and A2 show that whereas graben-bounding faults initiated early (during the early stage of the model run; 2.2% extension) as through-going faults in Model A1, they initiated later as mostly isolated, en échelon faults in Model A2 (compare Figures 8a and 8b). The en échelon faults are left-stepping and aligned parallel to the underlying pre-existing weak zone (Figure 8b).
There is a two-phase sequential evolution of the rift-related faults in Model B1. The first phase is characterized by the initiation and accumulation of normal strain (vertical displacement) on the graben-bounding faults during the early stage of the model deformation, and the second phase is characterized by the initiation and accumulation of normal strain on antithetic intra-graben faults during the intermediate to late stages (Figure 8c). A similar, but less-well developed two-phase development of the faults is also observed in Model B2 (Figure 8d).

Quantifying the Normal Versus Shear Strain on Model Structures
To quantify the obliqueness of normal faults we analyze the relative amount of normal versus shear components of strain accommodated by the graben-bounding faults and discrete intra-graben faults (when formed). The finite normal and shear strain are plotted along the profiles that are parallel to the model extension direction (Figure 9). The plots highlight the overall style of the graben-bounding faults, consisting mainly of either narrow (single) graben-bounding faults (models A1, A2, B1, B2, B3, and C), wide graben-bounding fault zones (models D2 and E), or both (Model D1) (Figure 9). In Model A3, both the normal and shear strain lines are flat because strain did not localize in the model domain (Figure 9i).
In the orthogonal extension models (models A1, B1, and the orthogonal half of Model E), the strain accommodated by the graben-bounding and intra-graben faults is 100% normal strain and zero shear strain consistent with pure normal faulting (Figures 9a-9c). In the oblique extension models (models A2, B2, B3, C, D1, D2, and the oblique half of Model E), the strain accommodated by the graben-bounding and intra-graben faults compose of both normal and shear components with a variable relative ratio between the models (Figures 9d-9h, 9j-9k). In Model A2, the finite normal strain on the graben-bounding fault is c.  Figures 9d and 9h). Consequently, the normal and shear strain on the wide graben-bounding fault zone and intra-graben fault in models E and D2 are similar.
In summary, in the orthogonal extension models, shear strain does not contribute to the total strain accommodated by the graben-bounding and intra-graben faults, except at the fault tips as observed in the strain maps described earlier (see Figures 4c and 5c). These faults are orthogonal to the extension direction. In the oblique extension models, when α = ±60° (low obliquity), shear strain account for only c. 20% (one-fifth) of the total strain accommodated by the graben-bounding faults and 40% (two-fifth) of the total strain accommodated by the intra-graben fault. When α = 45° (moderate obliquity), the amount of shear strain on the graben-bounding faults increases to c. 25% (one-fourth) of the total strain.

Comparison With Previous Analog and Numerical Modeling Studies
Observations from our models are broadly consistent with previous studies. Here, we briefly compare our key observations with observations from previous modeling studies and highlight the similarities and differences, where applicable.
First, the varying influence of a pre-existing crustal weak zone because of its orientation or depth compares well with previous models. For example, the ability of a weak zone to localize deformation is strongly 10.1029/2021TC006970 18 of 30 influenced by the orientation of the crustal weak zone with respect to the model-extension direction, that is, it localizes more when it is orthogonal (i.e., optimally oriented) and less with increasing obliquity (e.g., analytical solution by Ranalli and Yin [1990]; and analog models by Zwaan and Schreurs [2017], Molnar et al. [2019], Zwaan et al. [2021aZwaan et al. [ , 2021b). Under oblique extension, the graben-bounding faults tend to develop as characteristic en échelon faults which is comparable to observations in several oblique rift modeling studies (Agostini et al., 2009;Clifton et al., 2000;McClay & White, 1995;Tron & Brun, 1991;Van Wijk, 2005;Withjack & Jamison, 1986;Zwaan et al., 2021a). Also, our observation that a shallow weak zone (thin crustal cover) localizes more strain than a deeper weak zone agrees with previous analog model results by Sokoutis et al. (2007). The observation that the graben width narrows due to decreasing α-angle (increasing obliquity) of the weak zone agrees with previous analog modeling studies (e.g., Clifton et al., 2000;Tron & Brun, 1991;Zwaan et al., 2016). This is plausibly due to the increasing steepness/dip of the graben-bounding faults as they accommodate oblique-slip (e.g., Tron & Brun, 1991;Zwaan et al., 2016).
Second, the two-phase temporal evolution of rift-related faults (early graben-bounding faults vs. late intra-graben structures) in our models compares well with the early boundary faults and later internal faults observed in the orthogonal and low to moderate oblique rift analog models of Agostini et al. (2009). This reflects the progressive migration of strain from the borders of the rift-graben toward the center of the graben (Agostini et al., 2009). However, in our models, the intra-graben (internal) faults are less well developed and thus generally poorly expressed in the model surface (e.g., Figures 5a and 5f) compared to those generated by Agostini et al. (2009) and Philippon et al. (2015). This difference is likely due to a combination of factors, including (a) the difference in the model setup and boundary conditions (e.g., brittle-ductile in their models vs. brittle-only in our model), (b) the intrinsic strain-localization properties of the brittle material (coarse quartz sand in our model vs. fine-grained K-feldspar powder in their models), and (c) the amount of applied extension (c. 11% in our model vs. approximately double in their models). The later development of intra-graben faults in our orthogonal model B1 also led to the formation of a double marginal graben separated by a central horst (Figures 5a-5d). Similar deformation patterns where double marginal grabens develop have been previously described in numerous analog studies (e.g., Allemand & Brun, 1991;Corti, 2004Corti, , 2012Keep & McClay, 1997;Schreurs et al., 2006;Zwaan et al., 2019).
Third, the partitioning of deformation (in terms of both fault orientation and sense of fault slip) between the graben-bounding faults and intra-graben structures in our low (α = 60°) and moderately (α = 45°) oblique extension is consistent with several previous models. In terms of fault orientation, the graben-bounding faults are oriented obliquely to the extension direction, whereas the intra-graben structures (consisting of en échelon R' shears) are approximately orthogonal to the extension direction. Such variability in the orientation of the graben-bounding-structures versus internal-structures have been observed in many previous low to moderate oblique-rifts analog models (e.g., Agostini et al., 2009;Autin et al., 2010;M. Bonini et al., 1997;Corti, 2008;McClay & White, 1995;Philippon et al., 2015;Sani et al., 2019;Tron & Brun, 1991;Withjack & Jamison, 1986;Zwaan et al., 2016).
In terms of the sense of slip, an advantage of the high-resolution DIC technique we deploy is that it allows even low strains (not captured by visual inspection) to be quantitatively decomposed into the normal-and shear-components based on the 3D displacement field. Thus, both the normal and shear sense of slip on the graben-bounding-structure and internal-structure are constrained in our oblique extension models (Figures 4, 5, and 9). Our data show that the graben-bounding faults are oblique to the extension direction, and accommodate a major (80%-75%) normal-component and minor (20%-25%) shear-component of strain, indicating they are oblique-normal slip faults, rather than pure dip-slip. Similar oblique-normal sense of slip on faults striking oblique to the extension direction has been interpreted qualitatively in previous analog models (e.g., Agostini et al., 2009;Corti, 2012;Corti et al., 2007;Ghosh et al., 2020;Molnar et al., 2019;Tron & Brun, 1991;Withjack & Jamison, 1986). The extension-orthogonal, en échelon internal structures in our oblique extension models are pure strike-slip, differing from the dip-slip (vertical displacement) interpreted in earlier analog (e.g., Agostini et al., 2009;Corti, 2008;Philippon et al., 2015) and numerical models (e.g., Brune, 2014;Duclaux et al., 2020). This difference can again possibly be explained by the difference in the amount of applied extension, which is comparably small (about half) in our study. For example, Agostini et al. (2009) showed that at low to moderate obliquity, the internal-rift faults only started accommodating vertical displacement and thus detectable surface expression after c. 12.5% of extension, which is higher than the total extension in our models. Similar oriented en-échelon internal structures but interpreted as normal faults are detected in numerical simulations at about 4% (Duclaux et al., 2020) to 20% (Brune, 2014) extension. Duclaux et al. (2020) describe the rotation of these structures in a synthetic sense with respect to the overall transtensional shear sense which would imply these normal faults have an antithetic strikeslip component as seen in our analog models. The simulations of Brune (2014) suggest a set of strike-slip faults, where one is sub-parallel to the normal faults and one more oblique. We here infer that the internal structures possibly initiate as pure strike-slip (R' shear) structures (largely invisible to traditional monitoring techniques), which will later rotate and evolve to mainly dip-slip faults as extension increases (c. >10% extension) and strain migrates to the center of the rift-graben. It is however noteworthy that material properties controlling strain localization (i.e., strain weakening parameters) may significantly differ between nature, analog, and numerical models and thus the evolution may only qualitatively follow the path suggested here. In any case, the application of the DIC technique in our analog models thus provides new insights into the early stage development of these internal rift structures that are not resolvable in previous analog modeling studies using traditional visual inspection techniques (e.g., Agostini et al., 2009;Philippon et al., 2015).
Finally, the along-strike slip variation that characterizes the extension-orthogonal graben-bounding faults in our orthogonal extension models, i.e., pure dip-slip at the center of the faults and opposing sense of obliqueslip at the fault tips resemble observations from previous analog-( Corti et al., 2013;Philippon et al., 2015) and numerical- (Maniatis & Hampel, 2008) models. However, in the previous models the opposing sense of oblique-slip kinematics at the fault tips converges toward the center of the individual faults (see Philippon et al., 2015), whereas in our model, it diverges away from the fault center (Figures 4 and 5d). While convergent slip along single faults can be understood in terms of stretching with volume conservation (i.e., Poisson effect), divergent slip as observed in our models is interpreted as a boundary effect. More specifically, this likely relates to the inward drag of the deforming footwall block because of its direct frictional coupling to the laterally shortening basal foam while the hangingwall block, underlain by silicone, is decoupled from the foam and behaves like a more rigid block, resisting lateral shortening (see also, Zwaan et al., 2019).

Controls on Strain Localization Above Pre-Existing Weak Zones
Our models provide insights into the ways in which pre-existing weak zones may influence the structural style and drive strain localization during the early stage of continental extension. During extension, the presence of an underlying pre-existing weak zone results in the development and localization of strain on graben-bounding faults and fault zones that form in the brittle cover, in all but one of the simulated scenarios (Model A3) (Figures 4-7).
In experimental Series A and B, in which the size and orientation of the weak zone were varied, we observe that strain preferentially localized on graben-bounding faults in the brittle cover when the weak zone is oriented at ≥60° to the extension direction, irrespective of the size of the weak zone (Figures 4a-4j and 5a-j). When the weak zone is oriented 45° to the extension direction, the presence of the large weak zone (Model B3) resulted in strain localization along graben-bounding faults, whereas the small weak zone (Model A3) did not (Figures 4k and 5k). This suggests that the control exerted by the size of pre-existing weak zones becomes increasingly more important with increasing obliquity of the weak zone with respect to the regional maximum extension direction. As the angle between the orientations of pre-existing weak zone and the extension direction reduces (α ≤ 45°), smaller-scale weak zones are less likely to locally perturb the regional stress field and localize strain, whereas, larger-scale weak zones may still localize strain at a much lower angle of orientation with respect to the extension direction. Ranalli and Yin (1990) presented detailed analytical solutions that demonstrate that the critical differential stress (σ 1 -σ 3 ) required for strain localization and eventual faulting on a plane parallel to a pre-existing strength anisotropy (i.e., weak zone) is dependent on: (a) the material parameters and layering, which controls the yield stress curve/strength profile, and the integrated strength in the vicinity of the anisotropy, compared to that of the homogeneous, isotropic column of the simulated crust (as quantified in Table 2, see also Figures 2e) and (b) the orientation of the strength anisotropy, that is, there is a maximum orientation with respect to the principal extension direction (as a function of the integrated strength), beyond which strain localization is no longer possible along the anisotropy. Comparing the results of Models A3 and C provides insight on the possible role of the depth of pre-existing weak zones on normal fault evolution, since the size (small) and orientation (α = 45°) of the weak zone are the same, but Model C has a thinner brittle layer. In contrast to the lack of graben-bounding faults in the model with the thicker brittle layer (Model A3; Figures 4k-4o), strain localization on graben-bounding faults occurs above the weak zone in the thinner brittle layer (Model C; Figures 6a-6d). This suggests that a shallow pre-existing weak zone is more likely to localize strain in the brittle cover compared to a deep weak zone. Again, this is in strong agreement with Ranalli and Yin (1990), where they showed that critical differential stress increases with depth, implying that a lower critical differential stress is required for strain to localize in the vicinity of shallow weak zones compared to deeper weak zones.
Our model observations correlate well with observations in many natural rifts such as the northern North Sea Rift (e.g., Osagiede et al., 2020;Phillips et al., 2016Phillips et al., , 2019, Taranaki Basin, New Zealand (Collanega et al., 2019), and East African Rift (e.g., Daly et al., 1989;Kolawole et al., 2021;Morley, 1995), where some underlying pre-existing shear zones influenced the location, segmentation, and geometry of subsequent rift-related structures, whereas others had limited or no influence. For example, to the south of the western branch of the East African Rift System, the Rukwa-Malawi Rift segments preferentially developed along the large-scale NW-SE-trending inter-cratonic Ubendian mobile belt consisting of amalgamated Precambrian shear zones, whereas the N-S-trending Kenya Rift cross-cut the similar, but relatively smaller-scale NW-SE-trending Aswa Shear Zone (see Figure 1b) (e.g., Daly et al., 1989;Morley, 1995;Theunissen et al., 1996). More specifically, our work suggests that this selective influence of pre-existing weak zones on strain localization during extension is controlled by a complex interplay between the orientation, size, and depth of occurrence of the weak zones ( Figure 10). This conclusion supports the previous suggestion by Phillips et al. (2016) where only Devonian shear zones with vertical thickness of 1-2 km (and not thin c. 100 m ones) were preferentially exploited and influenced younger rift faults in the northern North Sea. Similarly, Osagiede et al. (2020) report that the Utsira Shear Zone (with a vertical thickness >3 km) locally perturbed the regional stress field during the Middle Jurassic-Early Cretaceous rift phase in the northern North Sea, thus influencing the geometry and growth of the cover rift faults, whereas the relatively thinner Heimdal Shear Zone (<1 km of vertical thickness) had no influence on cover faults (see Figure 1d).

Influence on the Development and Timing of Rift Structures
In the models in which the underlying weak zone induces strain localization, the resulting graben is broadly parallel to the weak zone. However, the graben-bounding faults exhibit a range of styles from discrete through-going graben-bounding faults that mimic the orientation of the weak zone (especially in orthogonal extension models), to individual en échelon fault segments that may be oblique to the weak zone (in some oblique extension models) (compare Figures 8a and 8b). Overall, the correlation between the location, segmentation, and orientation of graben-bounding faults and underlying pre-existing weak zones indicate that the growth and geometry of rift normal faults can be strongly influenced by inherited crustal weak zones (e.g., Collanega et al., 2019;Osagiede et al., 2020).
Our quantitative model results provide new insights into how pre-existing weak zones control the timing of fault development in rifts. We find differences in the timing of faulting and the amount of strain accommodated by graben-bounding faults, as a function of either the orientation, size, or the depth of the weak zone. Graben-bounding faults nucleate earlier and accommodate more normal strain in orthogonal extension models compared to oblique extension models, where a complementary part of the strain is accommodated by obliqueslip and intra-graben strike-slip (cf. models A1 and A2; Figures 8 and 9). For models with the same weak zone orientation, the graben-bounding faults nucleate earlier for the thicker and/or shallower weak zone and accumulate more strain causing more pronounced strain shadows (cf. Models A2 and B2; Figures 8b, 8d, and 5g).

Influence on the Style of Bounding Fault System
We observe that the 3D geometry of the underlying weak zone dictates whether strain localizes on either (a) narrow, discrete faults, or (b) wide, diffuse fault zones, which ultimately bound rift-related graben. In many natural rift systems, for example, the North Sea, pre-existing weak zones (e.g., ductile shear zones) exhibit a range of map-view geometries, from largely linear (e.g., Hardangerfjord Shear Zone) to curvilinear (e.g., Utsira Shear Zone) (Figure 1a). In this study, we did not attempt to directly simulate any specific natural 3D weak zone geometry due to the complexity of these geometries in nature. However, we used four simplified, generic weak zone geometries in our models (see Table 2 and Figure 2f). All models in Experimental Series A, B, and C, where rift graben developed, were characterized by near-symmetrical strain localization on very narrow zones at the margins of the graben, resulting in a single or segmented graben-bounding faults (Figures 4-6 and 9). In all these models, the same weak zone geometry that simulates a block-like weak zone was used, and thus provides evidence of a striking relationship between the geometry of the weak zone and the distribution of strain.
Conversely, in Experimental Series D and E where we used weak zone geometries that simulated arched weak zones, strain at the graben margins was more diffuse, resulting in a more structurally complex rift architecture characterized by the development of fault zones (and not a single-fault) on both margins (Figures 7 and 9). These fault zones are characterized by a range of fault styles, including, near-parallel-, en échelon-, and anastomosing fault sets. Furthermore, the observed differences in the strain distribution in experimental Series D and E appear to also be related to whether the arched geometry is "tight" (Model D1), or "gentle" (models D2 and E) (Figure 7). These results suggest that rift architecture is at least partly controlled by the 3D geometry of the pre-existing weakness.
The Rukwa segment of the East African Rift System provides a natural example of where pre-existing structures impose a first-order control on the across-strike variation in the style and magnitude of rift-related strain (Kolawole et al., 2021). Kolawole et al. (2021) observe that the wide (>12 km) Chisi Shear Zone localized the development of a distributed fault zone containing a cluster of en échelon-faults, that bounds the basin to the southwest, whereas the discrete Lupa Fault bounds the basin to the northeast. These two styles of strain distribution, with a wide fault zone to the southwest versus a narrow discrete fault to the northeast, is particularly similar to the style of the bounding faults in our Model D1.

Influence on Strain Partitioning Under Oblique Extension
We have also observed significant strain partitioning in the oblique extension models, suggesting that obliquenormal-slip and strike-slip are important mechanisms in the development of oblique rifts (e.g., Figures 4h  and 5m). A well-documented example is the NW-SE-trending Gulf of California, Mexico. Withjack and Jamison (1986) compared the pattern of faulting and strain partitioning in their analog and analytical models to the normal to oblique-normal, and strike-slip faulting pattern observed in the Gulf of California. Based on this comparison, they inferred that the pattern of faulting and strain partitioning observed in the Gulf of California best fits with an oblique rift system where the rift trend is c. 30° to the relative displacement direction. The Gulf of California is therefore a highly oblique rift system where fault slip data show that E-W extension is accommodated by both oblique-normal slip (majorly dip-slip with minor dextral component) dominating the rift margin faults, and strike-slip dominating the intra-Gulf domain (see Bonini et al., 2019).
Our observations of oblique-normal slip on graben-bounding faults in our oblique extension models imply that in extensional settings, normal faults that are oblique to the regional extension direction most likely accommodate deformation by both dip-slip and strike-slip displacements. Although the dip-slip component of such faults may be greater, they are, however, not pure normal faults, but rather oblique-normal faults. This is the case in the graben-bounding faults of the NE-trending rift graben of the Parnaíba Basin, Brazil, which developed as a result of the brittle exploitation of the underlying NE-SW-trending Transbrasiliano Shear Zone under N-S regional extension (Figure 1c) (see de Castro et al., 2016). Similar minor strike-slip component is reported in several graben-bounding fault systems that are oblique to the regional extension direction in other natural rift settings, including the Northern Main Ethiopian Rift (e.g., Boccaletti et al., 1998;Corti, 2009), the Gulf of Aden (Dauteuil et al., 2001), the Gulf of California (e.g., Bonini et al., 2019;Withjack & Jamison, 1986), and the Mohns Ridge, Norwegian Sea (Dauteuil & Brun, 1996). We note that, in outcrop-based fault analysis, the horizontal offset of markers and kinematic indicators such as slickenlines allows for the discrimination of oblique-normal faults from pure normal faults. However, in seismic reflection-based analysis, it is more challenging to constrain the strike-slip component of displacement on rift faults due to the absence of obvious kinematic indicators, and as such only the dip-slip (throw) component is quantified. Therefore, care should be taken when using results of such analysis to scale fault length-displacement relationship, as they may likely skew such relationship more toward under-displaced faults.
In contrast to settings where faults that are oblique to the regional extension direction accommodate deformation by both dip-slip and strike-slip displacements, other natural examples have been observed where faults striking oblique to the regional extension are pure dip-slip (normal) faults, for example, the Rukwa Rift segment of the East African Rift (Morley, 2010), the Main Ethiopian Rift (Corti et al., 2013), and the Southern Malawi Rift (Williams et al., 2019). To explain the latter, Morley (2010) proposed that the presence of pre-existing weak zones in the brittle crust can result in the local re-orientation of the regional stress field, such that the maximum horizontal compressive stress lies sub-parallel to the weak zone. This explanation has been supported by the interpretation of analog models by Corti et al. (2013) and Philippon et al. (2015) partly based on re-analysis of earlier analog models. This contrasting influence of pre-existing weak zones raises an important question of which conditions favor strain partitioning versus stress re-orientation in extensional settings. Although answering this question is beyond the scope of this current study, we note that Philippon and Corti (2016) suggested that a moderate obliquity threshold of 45° could mark the transition from stress re-orientation to strain partitioning in a divergent setting, while our model results suggests that strain partitioning may already occur at lower obliquity. The difference could lie in the different experimental methods used (centrifuge vs. normal gravity modeling) including the respective boundary conditions and/or material properties. Detailed field as well as high-resolution analog modeling and numerical modeling studies will be key to addressing this question.

Conclusions
In this study, we have investigated how the variability in the size, orientation, depth, and overall geometry of pre-existing weak zones influence the structural style and pattern of strain localization during rifting. We achieved this through a series of extensional analog models consisting of weak zone made from silicone oil that is placed inside a layer of quartz sand, mimicking a brittle crust with inherited weak zone. We deployed DIC technique to monitor the progressive surface deformation, allowing the cumulative horizontal displacement to be constrained at a high resolution. Our results have implications for improving the understanding of the role of inherited structures, specifically crustal shear zones on rifts and rift-related faulting.
Our key findings may be summarized as follows: 1. The presence of pre-existing weak zones in the brittle crust reduces the integrated strength of the crust that may facilitate strain localization in the vicinity of the weakened crust, resulting in the formation of rift graben. The graben is bounded by new faults whose geometries are influenced by the orientation of the pre-existing weak zones when the weak zone is oriented at ≥ 60° to the extension direction. The scale of the weak zone becomes an important factor when it is oriented at a much lower angle (≤45°) to the extension direction, in which case only large-scale weak zones may effectively weaken the crust and localize strain. Additionally, shallow weak zones are more likely to influence the pattern of deformation in the cover during subsequent rifting compared to deep weak zones. These observations underscore how different properties of a pre-existing weak zone may interplay to control its influence during rifting. 2. The timing of faulting and the amount of strain accommodated by graben-bounding faults are influenced by the orientation, size, or the depth of pre-existing weak zone. Graben-bounding faults nucleate earlier and accommodate more normal strain when the weak zone is: (i) oriented orthogonal to the extension direction compared to oblique orientation, (ii) large-scale compared to small-scale, and (iii) shallow compared to deeply buried. 3. Strain is mainly accommodated by pure dip-slip on normal faults, when the weak zone is orthogonal to the extension direction (i.e., orthogonal rift). Whereas, when the weak zone is oblique to the extension direction (i.e., oblique rift), strain is partitioned (increasing partitioning with increasing obliquity) into: (i) graben-bounding faults characterized by oblique-normal slip motion (i.e., major dip-slip, with minor component of strike-slip), and (ii) intra-graben domain dominated by strike-slip structures during the early stage of rifting. 4. Under oblique extension, intra-graben faults that are orthogonal to the extension direction, initiates as strike-slip, antithetic Riedel shears. 5. Our model results provide insights that may help to interpret observations in many natural rift systems such as the North Sea Rift and the East African Rift. n m (where E n is the matrix index in x and E m in y direction) is approximated by taking the average between the neighboring vectors    , 1 y E V n m and    , 1 y E V n m in the respective direction and multiplying it with the vector spacing E d : Figure A1. Benchmarks verifying the basal boundary condition of distributed, non-plane strain pure shear: The upper row shows final displacements after 50 mm of externally applied wall motion both parallel (y, left) and normal (x, right) to the extension direction for the benchmark model X (basal foam only) and reference model R (sand layer 40 mm thick). Displacement profiles are linear and parallel suggesting first-order quantitative transfer of homogeneous deformation from the basal foam to the sand layer. Middle row shows the corresponding longitudinal strains Eyy-y and Exx-y. Note that strain localization is unavoidable at the boundaries but the central area of interest is homogeneously extended longitudinally by c. 8%-10% and shortened transversally by c. 1.5% resulting in a "Poisson's ratio" of c. 16%. Appendix C Figure C1. Geometrical relationship of structures associated with pure shear and simple shear kinematics, and how they correlate with the major structures observed in our orthogonal and oblique models based on longitudinal and shear strain analysis (modified from Sylvester, 1988).