Exploring the Interactions Between Rift Propagation and Inherited Crustal Fabrics Through Experimental Modeling

Continental rifting is a geodynamic process that involves the breakup of the crust and may eventually evolve to seafloor spreading. Although it is often assumed to be a product of orthogonal divergence, continental rifting may result from oblique extension, and in several cases, it is related to the rotation of plates or crustal blocks about a vertical axis. This implies the occurrence of rifts with straight but not parallel margins and rift axis‐parallel gradients in extension velocity and amount of strain. The effects of rift propagation through the continental crust has only recently started to be addressed, and even less investigated is the interaction between rift propagation and inherited crustal fabrics. We have studied this issue by carrying out a series of analog laboratory experiments. Modeling results suggest that inherited linear discontinuities in the model brittle upper crust that are oriented above a threshold angle (≥45°) from the orthogonal to the rift axis have the ability to interact with rift‐related structures. Depending on their orientation, such inherited discontinuities are reactivated as either transfer zones or rift‐bounding faults. We compare our models with four natural prototypes in which rift propagation is likely to have occurred (the Trans Mexican Volcanic belt, the Gofa Province, the Gulf of Suez, and the Kenya Rift). For each case study, we show how the kinematics of propagating rifts is comparable to our model results, and we provide insights into how rift‐related deformation may interact with inherited crustal fabrics.


Introduction
Rifting is one of the most important geodynamic processes that has shaped the Earth through crustal and lithospheric thinning, which may ultimately result in the breakup of continental plates. Whereas the process is commonly assumed to be cylindrical or symmetric in many rift settings, some examples show a rift axis-parallel gradient of deformation associated with rift propagation, which implies variations in extension rates along rift margins (Figure 1a; e.g., . These examples include the Red Sea, the Pacific Galapagos Rise, the Woodlark Basin in the SW Pacific Ocean, SE of Papua New Guinea, the South China Basin, and part of the Main Ethiopian Rift . Recent studies (e.g., Molnar et al., 2019; have shown how strain features in these settings can be explained by the rotation of crustal blocks around an Euler pole, as a result of the so-called "rotational tectonics" (e.g., . Despite some pioneering attempts (e.g., Benes & Scott, 1996;Mart & Dauteuil, 2000;Souriot & Brun, 1992), the issue of rift propagation generated by rotation of crustal blocks has been addressed only recently in analog modeling. In particular,  showed that rotational rifting obtained by an expanding foam stretching an overlaying package of viscous and brittle material can properly reproduce a first-order triangular or "V-shaped" rift propagation pattern, which shows similarity with the above-mentioned natural cases worldwide. The seminal studies by Molnar et al. (2017Molnar et al. ( , 2018Molnar et al. ( , 2019 reproduced lithospheric scale rift propagation induced by rotational tectonics, investigating how this process can influence and drive the formation of micro-continents (Molnar et al., 2018) and how propagating rifts may interact with weak zones in the lithospheric mantle and lower crust (Molnar et al., 2017(Molnar et al., , 2019(Molnar et al., , 2020. In fact, many rift systems (both with orthogonal extension and rotational propagation style) have been shown to interfere with inherited structures (e.g., Autin et al., 2013;Brune et al., 2017;Corti, 2012;Corti et al., 2007;Molnar et al., 2017Molnar et al., , 2018Molnar et al., , 2019Pongwapee et al., 2019;Rotevatn et al., 2018;Zwaan & Schreurs, 2017). The relationships between the orientation of rift extension direction and that of lithospheric-and crustal-scale inherited structures are of primary interest since inherited fabrics might affect rift genesis and evolution in multiple ways and at various scales (e.g., Corti et al., 2018). (b) Schematic cartoon of the modeling apparatus. Two moving plates are rotated using a stepper motor. Deformation was distributed using a basal rubber sheet (RS; green rectangle). Inherited structures (one or two sets; S1 and S2) were reproduced within the brittle upper crust and placed at defined angles with respect to an arbitrary "model north" (N m ) corresponding to the direction of σ 3 , which is approximately taken as orthogonal to the direction of rift propagation (Table 2). (c) Top-view cartoon showing the rotational opening imposed on the model. (d) Sketch illustrating the adopted monitoring strategy: a fixed camera acquired high-resolution top-view photos, while a moving camera acquired high-resolution photos at predetermined positions (R1 0 to R1 7 and R2 0 to R2 7 ) for elaborating photogrammetric digital elevation model (DEM) reconstructions. The correct scaling of DEMs was achieved through the use of locally georeferenced markers (M 1 to M 4 ).
In this study, we implement the approach of previous analog modeling studies by presenting new laboratory experiments investigating the mutual effect between the propagation of crustal extension and inherited fabrics. In particular, we focus on upper crustal brittle fabrics, whose influence on the deformation pattern resulting from rift propagation was not considered in previous experimental work. We apply our modeling results to the Trans-Mexican Volcanic Belt (TMVB), Mexico, as well as to other propagating rifts worldwide, namely, the Gofa Province in Ethiopia, the Gulf of Suez, part of the Red Sea-Gulf of Aden system, and the Kenya Rift (KR).

Modeling Strategy
Our experiments were designed by adopting a parametric approach, which consisted of the systematic variation of the orientation of brittle discontinuities (e.g., mimicking faults and fractures association). This modeling was carried out to evaluate if and how rift propagation may interact with crustal-scale inherited fabrics. The structure of the models was analyzed quantitatively by means of photogrammetric digital elevation model (DEM) construction, digital particle image velocimetry (PIV), and semiautomatic fault pattern quantification. This procedure provided the basis for comparing model results with the selected natural prototypes.

Model Setup, Materials, and Scaling
The conceptual model of rifting (Figure 1a) was achieved through a modeling apparatus consisting of two moving plates (MPs) fixed at their external tips and rotating around two pivot points. This model setup allowed the progressive opening and consequent rift propagation within the overlying analog brittle/ductile continental crust, whereby the minimum horizontal principal stress trends roughly orthogonal to rift propagation (Figures 1b and 1c). The movement was obtained by dragging the plates with a stepper motor. A linear opening velocity of 10 mm/hr was applied at the left-hand side of each MP (see the red dot in Figures 1b  and 1c). Due to the rotational movement of the MP, the velocity gradually decreased toward the pole of rotation along the propagating rift axis. The boundaries of the MP generated a velocity discontinuity (VD), whose effect was dissipated by employing a basal elastic rubber sheet (RS; green rectangle in Figure 1b), the stretching of which allowed the distribution of deformation and avoided the formation of faults with unrealistic lengths and vertical throw. The RS is a natural latex elastic band 150 mm wide and 0.35 mm thick that was attached to the MP for 10 mm on each side, using high-resistance double-sided tape and then secured with American tape. The area of interest of each model corresponded to the center of the RS, where boundary effects were minimized. The models reproduced a two-layer, brittle-ductile continental crust (e.g., Bonini et al., 1997). We simulated the upper crust (UC) using a 1-cm-thick layer of a quartz and K-feldspar sand mixture (80:20 proportion % in weight; see Montanari et al., 2017 andDel Ventisette et al., 2019) ( Figure 1c). The lower crust (LC) was simulated using a 1-cm-thick layer of polydimethylsiloxane (PDMS) mixed with corundum powder in the proportion of 1:1 (% in weight) to obtain the correct density for the LC. All the details of material properties are reported in Table 1. The total thickness of the models (LC + UC) was 2 cm, over a total modeling area of about 70 × 35 cm and a more restricted area of interest of approximately 65 × 15 cm.
Models were scaled down to natural conditions in order to achieve correct length, dynamic, rheological, and kinematic similarity (Hubbert, 1937;Ramberg, 1981;Weijermars & Schmeling, 1986). In particular, the geometrical scaling ratio was~6.7 × 10 −7 , such that 1 cm in the models simulated 15 km in nature. Scaling of stress, viscosity, and velocity resulted in a natural extension rate of~1 cm/yr. Details of the scaling ratios are reported in Table 1, and the detailed scaling procedure is reported in Appendix A.

Description of the Experimental Series
We present the results of 10 experiments ( Table 2; raw data of these models are freely available to download from Maestrelli et al., 2020 at the link http://doi.org/10.5281/zenodo.3724666) out of 30 model runs for testing the relationships between rift propagation and inherited fabrics, as well as the reproducibility of results. To this aim, we introduced in the undeformed model one or two sets of discontinuities (referred to as S1 and S2) with specific azimuth angles (α and β, respectively) with respect to the orthogonal to the direction of rift propagation (Figure 1b), which defines the arbitrary "north" in the experiments (N m ) and that

10.1029/2020TC006211
Tectonics approximately corresponds to the initial "instantaneous stretching direction" (ISD). The discontinuities are represented by subvertical artificial dilatation zones obtained by "cutting" the sand pack with a knife secured to adjustable guides. Reorientation of particles and variation of their original compaction cause local weakening of the sand mixture (e.g., Bellahsen & Daniel, 2005), which has the ability to focus deformation when the discontinuities are favorably oriented under the local stress field induced by model deformation. Discontinuities were introduced in the right half of the model with an average spacing of 4 cm (slightly varying depending on their orientation). In each model, we varied the orientation of discontinuities as reported in Table 2. Model RP-1 had no discontinuities and was used as a "reference model." Models RP-2 and RP-3 tested only one set of discontinuities in order to avoid possible mutual effects resulting from the presence of two sets. Models RP-4 to RP-9 tested the effect of rift propagation when two distinct sets of discontinuities trending at various angles are present in the experiments.

Monitoring and Analysis of Deformation
The evolution of model surfaces was constantly monitored through the automatic acquisition of highresolution top-view photos with 120-s time steps. The finite topography was quantified using photogrammetric techniques (e.g., Donnadieu et al., 2003). Using structure from motion software (Agisoft Photoscan®), photographs taken from multiple perspectives of the model surface were used to generate a 3-D point cloud Note. S1 and S2 represent the two sets of discontinuities introduced artificially in the layer simulating the upper crust to simulate inherited crustal fabrics. Azimuth was measured in a clockwise direction to the arbitrary "north of the model" (N m ), that is, the orthogonal to rift propagation direction corresponding to the instantaneous stretching direction (ISD). ) v* = ε*l* = 14,000 Note. The asterisk (*) denotes the ratio between model and nature for a given parameter. Characteristics of the granular material represented by the Qz-K-feldspar sand mixture are derived from Montanari et al. (2017) and Del Ventisette et al. (2019). The range for the internal friction coefficients is calculated considering the "peak friction" and "stable friction" values of the granular mixture. and from that a DEM ( Figure 1d). The use of markers placed at fixed and locally geo-referenced positions in the model system ensured that all calculated DEMs were on the same scale and in the same coordinate system (more details on this procedure are reported in Maestrelli et al., 2020). Quantitative analysis of models was performed using the free software FracPaQ (Healy et al., 2017) developed for MATLAB™, which allows the quantification of parameters that describe structural patterns in the models (e.g., fault trace length and orientation, estimated structure density [E D ] and intensity [E I ], and dilation tendency of structures). Digital PIV analyses were performed on selected model top-view photos to evaluate displacement vectors and quantify deformation using the free software PIVLab for MATLAB™ developed by Thielicke and Stamhuis (2014).

Model RP-1: Reference Model
Model RP-1 was run without introducing discontinuities (Table 1) in order to obtain a reference model allowing observation of the basic system of rift development. From t = 0 to t = 80 min, the opening of the MP caused incipient deformation within and outside of the area of interest, which is indicated by the darker area in Figure 2. After t = 80 up to t = 120 min (Figure 2a), small graben systems started to develop and link together, propagating toward the pole of rotation. Fault throw decreased from the proximal, left-hand (western) side of the model (i.e., the area opposite to the pole of rotation), which underwent more finite and incremental deformation toward the eastern distal pole of rotation, where only shallow grabens were formed. Overall, the cumulative throw of graben systems led to a structural configuration similar to a

10.1029/2020TC006211
Tectonics "wide rift" setting (e.g., Buck, 1991). This deformation can be referred to as a "deformation wedge" propagating toward the pole of rotation and increasing the number, length, and vertical throw of the faults (Figures 2a-2f), which trend on average east-west (see Figure 1b). These observations are supported by slope analysis performed on the DEM obtained for the final deformation stage (Figure 3b).
Gradients in the degree of deformation were detected and highlighted by carrying out a slope analysis (Figure 3b) of the model DEM ( Figure 3c) using ArcGIS® software, where the slope attribute highlights fault scarps. Steeper fault scarps localize in the leftmost sector of the experiment, where the amount of extension is greatest and where fault density is highest. More to the right (i.e., closer to the pole of rotation), vertical fault throw is smaller (Figure 3c), similar to the slope of the fault scarps. This is particularly evident by comparing topographic profiles obtained from the model DEM (Figure 3d). Profile a-a′ shows a gradual decrease of elevation from the distal area (i.e., close to the pole of rotation) toward the left-hand side of the model,

Models RP-2 to RP-9: Inherited Fabrics 3.2.1. Models RP-2 to RP-4: Inherited Discontinuities Trending N30°and N125°M
odels RP-2 and RP-3 were designed to test the role of a single set of discontinuities, trending at angles of α = 30°(Set S1) and β = 125°(Set S2), respectively. Model RP-4 was then designed to test the two sets in combination (see Table 2).
In Model RP-2 (Set S1 with α = 30°), the general development did not differ significantly from the evolution of Model RP-1, showing incipient normal faulting in the proximal area of the model, which progressively propagated toward the rotation pole (right-hand side of the model) to form a "V-shaped" subsiding wedge in the center of the model. The wedge tip migrated progressively with the opening of the metal plates. Model deformation showed a clear gradient in terms of fault throw and fault number (both larger for faults on the proximal part of the model). The average fault trend (red line in Figure 5a) is N m 90°. No evidence of reactivation of Set S1 discontinuities was observed ( Figure 5a).

10.1029/2020TC006211
Tectonics Model RP-3 (Table 2 and Figure 5b) tested Set S2 and showed a similar evolution to Models RP-1 and RP-2. Larger faults, with the largest vertical throw, reached a maximum length of~9 cm (Figure 5b), as in Model RP-2. Model RP-3 showed clear evidence for reactivation of the S2 discontinuities, which occurred in a scattered and discontinuous manner during deformation, in the form of incipient small graben systems trending N m 125°( Figure 5b).
Model RP-4 tested the presence of Sets S1 and S2 together and showed an overall evolution similar to Models RP-2 and RP-3 ( Figure 5c). The average trend of structures (N m 95°) marks the prominent role of rift-related faults. As in Model RP-3, the S2 weak zones were reactivated (red lines in Figure 5c, trending N m 125°), whereas the S1 discontinuities (with α = 30°) were not reactivated, similarly to Model RP-2.
Notably, both in Models RP-3 and RP-4, inherited structures to the left experience stronger reactivation than those to the right (closer to the distal area) and a more marked reactivation in the topmost area (in planview) than those in the bottom portion of the models. Model RP-5 ( Figure 6a) tested two sets of discontinuities with α = 15°(Set S1) and β = 165°(Set S2). In this model, the "V-shaped" propagation of the graben systems is extremely well developed. None of the S1 and S2 sets was reactivated during model deformation. Even though no evident reactivation was visible (i.e., slip along inherited structures acquiring new throw), they may have subtle effects on the propagation of rift-related structures since the latter are slightly segmented and shifted laterally from their straight propagation trajectory in correspondence with the inherited structures.
Model RP-6 tested the effect of discontinuities trending at lower angles (Set S1 with α = 45°and S2 with β = 135°). This model was repeated three times since it did not provide unequivocal indications on fault reactivation (we show here two out of all the performed attempts; Figures 6b and 6c). In fact, two models out of three showed no reactivation of discontinuity Sets S1 and S2 (Figure 6c shows one of the two models with no reactivation), while one model showed a slight reactivation of both sets (Figure 6b), with no continuous faults and small amount of subsidence between grabens generated by reactivated fault couples. Nonetheless, as in Model RP-5, the discontinuities may not be reactivating but they seem to be exerting some control over the rift architecture.
Model RP-7 ( Figure 6d) tested discontinuity sets trending at α = 60°(Set S1) and α = 120°(Set S2). In the lefthand side proximal area of the model, the cumulative deformation was accommodated by the development of N m 90°-trending (on average) graben systems. Both S1 and S2 weak zones were largely reactivated by rift-related deformation, causing the development of a "chessboard" pattern of reactivated structures.
(a-f) Top-view photos of the final stage of Models RP-5 to RP-9 with different orientation of discontinuities with respect to the propagating rift (symbols as in Figure 6). Panels (b) and (c) show two models testing the effect of inherited discontinuities trending at 45°(Set S1) and 135°(Set S2). In Model RP-6a (b), discontinuities were only slightly reactivated, while in Model RP-6b (c), no reactivation of the two sets was observed. Rose diagrams legend as in Figure 5.

10.1029/2020TC006211
Tectonics Similarly, Model RP-9 showed the continuity of deformation proceeding from the rift-related structures in the area of incipient deformation (left-hand side of the model in Figure 6f) to the area of reactivation of S1 discontinuities (with α = 90°). A "V-shaped" depression did not develop, and strain was accommodated along the reactivated weak zones. In the area between two reactivated discontinuities (marked by red lines in Figure 6f), we observed parallel rift-related structures. These latter show less variability in their orientation than those observed in the left-hand part of the model (where no preexisting discontinuities were cut), together with a longer continuity. Besides, no evidence of fault reactivation was observed for S2 discontinuities, which trend at a high angle to the direction of rift propagation (β = 180°).
Overall, the behavior of the reactivated discontinuities (i.e., Models RP-6a, RP-7, RP-8, and RP-9) was similar. In particular, a pair of conjugate faults developed in weak zones. Therefore, when reactivated, a single sand cut invariably localized a couple of faults delimiting a graben rooted at the cuts' deep end.

Rift Propagation Versus Inherited Fabrics
Distribution and chronology of fault formation reflect the rotational nature of the model setup, implying a gradient of deformation propagating toward the pole of rotation. Fault displacement is also consistently higher for structures developed in the proximal regions of the models (where cumulative deformation is larger) and smaller in distal areas adjacent to the rotation pole, which experience lower strain. Velocity vectors appear orthogonal to the main graben faults (see PIV analysis in Appendix B), suggesting a dominant dip-slip component, while the reactivated structures can have a strike-slip component that is a function of their trend.
Overall, our reference Model RP-1 is in agreement with the modeling results shown by , who reproduced the effect of rotational tectonics with a different yet comparable model setup. A main difference occurs due to the presence of the RS, which we employed to distribute deformation, while Zwaan et al.'s (2020) models aimed to localize extensional faults at the center of the model using a viscous "seed." Compared to their setup, the presence of inherited brittle discontinuities in our experimental series represents a further difference. Nonetheless, the role of inherited fabrics was also addressed by Molnar et al. (2017Molnar et al. ( , 2018Molnar et al. ( , 2019Molnar et al. ( , 2020, who investigated rotational extension at a lithospheric-scale, highlighting the key role of crustal and lithospheric weak zones in controlling the formation of microcontinental blocks and testing their influence on surface deformation patterns. Despite a different aim and model setup, the overall surface deformation compares well with our models, as both show a "V-shaped" system of extensional faults with a clear decrease in the amount of extensional strain toward the rotation pole. Furthermore, Molnar et al. (2018) document how weak-zones are crucial to favor (or not) the detachment of microcontinents, highlighting that the orientation of inherited fabrics at the lithospheric scale is a key factor controlling this process.
In our models, many~N m 90°-trending rift-related faults (i.e., the black lines in Figure 6) were influenced by the reactivation of discontinuities. Discontinuity Set S2 (with α = 125°) has a noteworthy effect: These structures interacted with the rift-related faults during their propagation, deflecting the trajectory of the latter (Figure 7a). Reactivated structures have often acted as transfer zones, exerting a hard linkage, forcing the propagation of rift-related faults to shift laterally (yellow arrows in Figure 7a) and connecting the graben boundary faults. In addition, discontinuities trending at a low-angle to the direction of rift propagation were reactivated as highly dilatant structures, localizing all of the extensional strain, as observed in Model RP-8 ( Figure 6e).
Finally, there is a clear relationship between the trend of discontinuities and their reactivation: Only favorably oriented discontinuities can be reactivated by extension (Figure 7b). Specifically, in the experiments reported here, only discontinuities ranging between α = 45°and β = 135°are favorably oriented for reactivation. A further corollary is that the discontinuity angles tested in Model RP-6 (i.e., Set S1 trending α = 45°a nd Set S2 trending β = 135°; Figures 5b and 5c) likely represent "limit angles" (LA) for reactivation, since out of the three experiments performed with this setup, two did not experience reactivation, and the other one (i.e., Model RP-6a, Figure 6b) showed only little evidence of reactivation. This behavior is confirmed by the observation that no reactivation was observed for discontinuities oriented with an angle α < 30°1 0.1029/2020TC006211 Tectonics and β > 135° (Figure 7b). Nonetheless, as shown by Models RP-5 and RP-6b, high-angle inherited structures that are not reactivated (i.e., that are not developing visible slip and generating throw) can also exert some control on the architecture of rift-related structures. In particular, inherited structures can force rift-related faults to shift laterally from their propagation trajectory, inducing some "soft linkage" (Figure 7c), conversely to what is observed when inherited structures are fully reactivated (e.g., Figure 7a). This contributes to an increase in rift segmentation and shows that high-angle inherited structures (α < 45°a nd β > 135°) may act as accommodation zones for rift propagation (Figures 7b and 7c).
Further consideration can be given to the gradient of reactivation observed on inherited structures (e.g., Models RP-3, RP-4, RP-6a, and RP-7; Figures 5b, 5c, 6b, and 6d). Both a proximal-distal gradient and an along-strike gradient of reactivation can be observed (e.g., Model RP-3; Figure 8a). Both are controlled by the distance from the pole of rotation and by the trend of the inherited structures, respectively. These gradients are particularly evident if we extend the observation outside the area of interest of the models (Figure 8a). The proximal-distal gradient directly depends on the distance of a certain structure from the pole of rotation: the greater the distance, the greater the degree of reactivation, in agreement with the gradient of deformation imposed by the rotational nature of the model setup (Figure 8b). The along-strike gradient is controlled by the trend of the inherited structures, since the lower the angle, the greater the degree of reactivation (Figure 8b). A decrease in the discontinuity trend angle implies that the left-hand side of a given inherited structure is closer to the proximal area, while its right portion is farther. This implies an

10.1029/2020TC006211
Tectonics increase in the along-strike gradient (i.e., the reactivation is stronger in its proximal portion than in its distal one; Figures 8a and 8b). For example, the along-strike gradient is zero for inherited structures that are orthogonal to the direction of rift propagation, while it coincides with the proximal-distal gradient for structures trending N m 90°. These considerations do not apply to inherited structures trending <45°and >135°(i.e., above and below the LA; Figures 8b and 9c) since they were not reactivated. In these cases, inherited structures acted only as accommodation zones. Phillips et al. (2019) showed how inherited structures trending 45°to 90°to the direction of extension affected the localization (and the orientation) of rift-related structures in the North Sea, while structures trending at high-angles do not influence the formation of new faults, although they can contribute to basin segmentation. This consideration is in line with our observations, indicating that high-angle discontinuities can interact with rift-related structures to generate accommodation zones. Furthermore, the observation that inherited structures strongly influence rift-related fault propagation is in agreement with other natural examples. For example, Phillips et al. (2016) and Schiffer et al. (2019) report that inherited basement structures had a major influence on the evolution of rifting in offshore Norway and in the North Atlantic, respectively. This influence on rift development has been observed in analog experiments. Corti et al. (2007) demonstrated that inherited structures are able to induce oblique-slip kinematics on newly forming rift-related structures, explaining the Z-shaped structural pattern observed in some natural rifts. Finally, a similar influence was demonstrated in the analog experiments of Bellahsen and Daniel (2005) and Bellahsen et al. (2006), which showed structural patterns similar to those observed in our models (c.f. Figure

Comparison With Natural Cases 4.2.1. Trans-Mexican Volcanic Belt
The motivation for the present study was to investigate both the role of tectonic inheritance and active tectonics in specific regions of the TMVB, which is a large-scale, NW-SE trending, >1,000 km long volcano-tectonic feature extending through central Mexico. Ferrari (2004) considered that continental rifting propagated southeastward as a result of progressive tearing of the Cocos slab. As a result of slab tearing and sinking, material from the asthenosphere was forced to migrate upward, inducing mafic magmatism, which was emplaced progressively from NW to SE in agreement with the direction of continental rift propagation (Figure 9b). In the TMVB, several inherited structures have been interpreted to interact with Miocene-Quaternary faults and to have localized eruptive centers (Garduño-Monroy et al., 2009;Gómez-Vasconcelos et al., 2020). A variety of observations (e.g., focal mechanisms, borehole data, and geological indicators) indicate a range of local extension directions throughout the TMVB (Heidbach et al., 2016). In central Mexico, Quaternary structures trending E-W and NNE-SSW suggests an average N-S Quaternary extension direction (e.g., Garcıa-Palomo et al., 2000; Figure 9c). The large-scale crustal extension observed in our experiments has similar deformation gradients to those observed in the TMVB (Figure 9a). Ferrari et al. (2012) showed how the crustal thickness doubles from WNW toward ESE, reaching a maximum value below Mexico City (Figure 9b). A similar increase was also observed in our experiments, where the analog crust thickens toward the pole of rotation (e.g., see section b-b′ along Model RP-1; Figure 3d). Similarities were also observed when comparing extensional tectonic structures, which trend roughly parallel to the direction of rift propagation. Nonetheless, equivalents to~N-S-trending structural systems (e.g., the Colima Rift and the Penajmillo Graben; Figure 9b) were not observed in our experiments. These systems might be related to larger-scale N-S trending lithospheric inherited features (e.g., Ferrari et al., 2012) that were not simulated in our setup. Furthermore, fault density is higher in the area where the crustal extension initiated, while less deformation is observed toward the SE, where no NW-SE trending structures are observed at the surface. This is in agreement with the pattern of deformation reproduced by our models (compare with fault density and intensity analysis in Figure 4b). Nevertheless, erupted surface deposits from recent volcanic centers obscured or covered evidence of tectonic structures. Volcanic centers show an older-to-younger trend toward the SE, suggesting that continental rifting propagated southeastward (Ferrari, 2004) (Figure 9b). This kinematic model is in good agreement with our analog experiments. In the central area of the TMVB, the Morelia-Acambay Fault System (MAFS) shows important similarities with our Model RP-4 (cf. Figure 9c and Figure 9d). Here, a system of roughly E-W-trending Miocene-Quaternary grabens that are aligned with historical seismicity (Figure 9c) are dissected by older (~30 Ma) NW-SE trending faults, which show evidence of right-lateral reactivation (Garduño-Monroy et al., 2009).
In order to compare our models with nature, we have rotated Model RP-4 to fit the geographic north (N; see Figure 9). The trend of faults in the MAFS is comparable to the S2 discontinuities in Model RP-4 (anisotropies trending β = 125°corresponding to N140°in nature) and Models RP-6a,b (anisotropies trending β = 135°i n the model corresponding to N150°in nature).
Models RP-6a,b showed that discontinuities with β = 135°are occasionally reactivated during extension, while discontinuities trending β = 125°are always reactivated. This is in good agreement with the setting of the MAFS, which has been interpreted as a reactivated fault system (Garduño-Monroy et al., 2009). Model results suggest that fault reactivation can be related to the propagation of rift-related extensional faults. Furthermore, Model RP-4 (as well as Model RP-3) showed that inherited faults are often reactivated with a strike-slip component, which may offset the E-W trending rift segments (see Figure 7a).
During progressive southeastward propagation of crustal extension, the TMVB developed several large-scale volcanoes and calderas, some of which are the target of geothermal exploration.
Inherited crustal fabrics have also been reported in the easternmost sector of the TMVB, where they are represented by~NE-SW and NW-SE-trending regional faults (Campos-Enriquez & Garduño-Monroy, 1987). The reactivation of inherited faults by continental rifting is an important issue, as these features have the ability to focus the ascent of magmas and related hydrothermal fluids. These structures may have played an important role in the localization of volcanic and caldera complex systems (e.g., the Acoculco caldera complex and the Los Humeros volcanic complex; Avellán et al., 2020; Figure 9b)  By rotating Model RP-4 15°counterclockwise, we obtain α = 30°and β = 135°for discontinuities S1 and S2, respectively. This model showed that discontinuity Set S2 was reactivated by crustal extension induced by the rift propagation, implying that the discontinuities trending around N140°in this area might have been reactivated. This result accords well with the observation that the NW-SE trending faults in the Acoculco caldera system show the highest alteration produced by the circulation of hydrothermal fluids (Liotta & WP4 Working Group, 2019). A further observation based on our modeling results is that the interaction between rift-related faults and inherited fabrics may provide an efficient process for producing intersecting active structures, which provide favorable conditions for localizing volcanic centers. This is confirmed by fault dilation tendency (T D ; Ferrill et al., 1999) analysis (see Appendix C) on inherited faults performed on our models, showing that these structures are prone to be reactivated by the imposed stress field.

Gofa Province
A further case of comparison for our models is the Gofa Province (Moore & Davidson, 1978), north of the Chew Bahir basin and west of the southern segment of the Main Ethiopian Rift (Figure 10). Here, a system of normal faults trending roughly SW-NE has been hypothesized to result from the rift propagation that resulted from local block rotation induced by the Somalian plate motion (Philippon et al., 2014). Furthermore, the role of NW-SE-trending inherited structures in the basement has been invoked to explain the zigzag deformation pattern acquired by the basins (Corti, 2009;Philippon et al., 2014). The rift-related structures (white lines in Figure 10) are offset by inherited faults (yellow lines), documenting that the southern portion of the rifted area is wider than the northern sector, as often occurs in propagating rift settings. The inherited faults trend on average 45-50°to the regional extension direction, with an ISD trending N105°(black arrows in Figure 10; Philippon et al., 2014), which is similar to present-day GPS vectors that indicate an extension direction trending N100°(green arrows in Figure 10; Kogan et al., 2012). Seismicity in the area is localized along both rift-related and inherited structures, and the few available focal mechanisms show dip-slip to transtensional solutions. Assuming an orthogonal ISD to the direction of rift propagation (and coincident with N m ), the trend of reactivated structures is comparable with the trends tested in Models RP-6 and RP-7. As discussed above (Figure 7), Model RP-6 suggests that 45°is a limiting angle for reactivation, and Model RP-7 indicates that discontinuities trending 60°to the ISD are always reactivated. This observation therefore supports the idea that inherited structures trending 45-50°to the ISD could have been reactivated during rift propagation, influencing the development of the structural pattern of rift-related faults in the Gofa Province. In particular, Philippon et al. (2014) suggest that the NW-SE-trending structures could have acted as transfer zones, a conclusion confirmed both by seismicity and by our modeling. The presence of transfer zones is documented not only in Models RP-6a and RP-7 (Figures 6 and 10) but also in many other models experiencing reactivation (see Figure 7).

Gulf of Suez
A third natural example is the Gulf of Suez, Egypt (Figures 11a and 11b), which represents the~350 km long northern continuation of the propagating, and partly oceanized, Red Sea. The evolution of the Red Sea-Gulf of the Aden system is driven by triple-junction dynamics related to the African and Arabian plates. Bosworth et al. (2005) hypothesized a three-stage evolutionary model, the last of which implies rotation around an Eulerian pole that provided the basis of analog modeling experiments addressing the Red Sea with a rotational tectonic setup (e.g., Molnar et al., 2017Molnar et al., , 2018Molnar et al., , 2020. The Gulf of Suez represents an interesting case for comparison since inherited structures in the Precambrian basement have been mapped and hypothesized to control the pattern of rift-related faults during northwestward propagation of rifting. Inherited structures were identified by McClay and Khalil (1998) as~N-trending faults, fractures, and shear zones in the Esh El Mallha area and in the Hadaid block (EEM and HB in Figure 11a), while Younes and McClay (2002) documented the reactivation of the Ekma and Durba faults as transfer zones (Ekma-Durba Transfer Fault System, EDTFS in Figure 11a). The reactivation of these structures (for the EDTFS also supported by the localization of historical seismicity; Figure 11a) under an N60°-directed extension (Younes & McClay, 2002) compares well with the setting of Model RP-3 (Figure 11c), which tested S2 discontinuities with β = 125°. The latter represents the average angle between the ISD and the inherited faults in the three areas of interest. Model RP-3 showed that such inherited structures fall in the experimental reactivation field (Figure 7b) and thus could have been favorably oriented for reactivation. These structures often act as transfer zones (cf. areas marked by yellow dashed circles in Figure 11c), which displace the pattern of rift-related faults and show a good correlation between our models and these natural examples. A deviation from this favorable comparison is represented by the absence of the steep rift margins observed in the Gulf of Suez 10.1029/2020TC006211 Tectonics system, which is not reproduced in our models. This is likely due to the experimental setup architecture, which was designed to distribute deformation and avoid localization of deformation at the VD. In this regard, a stronger influence of the VD may have improved this similarity.

Kenya Rift
The Keny Rift is one of the major rifts segments of the East African Rift System (EARS) (e.g., Morley, 1999a), and it has propagated from the Turkana depression toward the Tanzania Divergence to the south (Figures 12a and 12b; Nyblade & Brazier, 2002). Focal mechanisms in the northern sector show dip-slip to transtensional solutions aligned with the main rift-related structures. A similar pattern occurs in the southern sector, but with a higher complexity due to the interaction with the Tanzania Divergence. The southward propagation was likely induced by the clockwise rotation of the Somalian Block (e.g., Bonini et al., 2005;Fernandes et al., 2004;Nocquet et al., 2006) (Figure 12c). The KR also shows some curvature as a result of interaction with the Tanzanian Craton, whose rotation is also partly responsible for the northward propagation of the Tanganyika Rift. The areas surrounding the KR have been shown to be largely shaped by the presence of inherited basement structures, with a wide range of orientations. Many of these structures are parallel to the major NW-SE-trending Aswa Shear Zone (Le Turdu et al., 1999;Morley, 1999b;Muirhead & Kattenhorn, 2018). Inherited faults associated with this major lineament interact with the KR. For example, the Kordjya Fault System (α ≈ 65°) in the Magadi Basin is interpreted to result from the interaction between rift-related faults (white lines in Figure 12c) and inherited faults (yellow lines in Figure 12c) (Muirhead & Kattenhorn, 2018). Model RP-7 shows similarities with the KR (Figure 12d), as rift-related faults reactivated inherited fault segments with α ≈ 60°. The interaction between rift-related and inherited faults formed structural patterns that are similar to the Kordjya Fault system in the Magadi Basin (cf. Figure 12c and Figure 12e). Furthermore, seismicity localized on some of the inherited structures

Concluding Remarks
We performed rotational extensional experiments in order to investigate the interaction between propagating rifts and preexisting brittle crustal fabrics. The orientation of preexisting discontinuities was varied relative to the extension direction. Our models provide general and specific insights into the rotational rift propagation process and the possible structural evolution of some specific natural systems. We summarize our main findings below: 1. Our models reproduced rift propagation in terms of density of structures, number and depth of grabens, and the overall crustal thinning toward the pole of rotation. Models are comparable in general terms with existing processes described in the literature. 2. The results of modeling provide new insights into the reactivation of inherited structures, suggesting that the reactivation field is delimited by an angle ≥45°to the extension direction. 3. Depending on their orientation, reactivated inherited faults often acted as transfer faults or oblique-slip faults. High-angle inherited faults, although not reactivated, may also act as accommodation zones. These features have the ability to offset rift propagation segments (i.e., exerting respectively hard and soft linkage effects), generating fault patterns similar to some natural examples. 4. Gradients of reactivation were defined as a function of the distance from the pole of rotation and the trend of the inherited structure. Figure 11. (a) The Gulf of Suez (Egypt), which is the (b) northern part of the Red Sea-Gulf of Aden rift system. Rift propagation is related to the final evolutionary stage of the system, implying the counterclockwise rotation (b) of the Arabian Plate and the opening of the two rifts (Bosworth et al., 2005). Structures in panel (a) are mapped after Moustafa and El-Raey (1993) and Younes and McClay (2002), while orange dots show historical earthquakes (time interval 1900-2020). Seismicity and focal mechanisms are plotted for earthquakes with M ≥ 3 from the USGS Earthquake Catalog available at https://www.usgs.gov/natural-hazards/earthquake-hazards/earthquakes. (c) Model RP-3 showing the reactivation of inherited discontinuities as transfer zones (trending at 125°to the direction of σ 3 ) and displacing the pattern of rift-related faults. Model RP-3 has been rotated to fit the natural setting.

10.1029/2020TC006211
Tectonics 5. Large-scale similarities are found between our models and the TMVB, as well as the local-scale setting of the MAFS and the Los Humeros and Acoculco areas. Here, inherited discontinuities were reactivated as transfer faults and may interact with rift-related faults. 6. Our models also show similarities to other continental rifts, particularly the Gofa Province, the Gulf of Suez, and KR. In all these cases, inherited structures have strongly interacted with propagating rift segments, often acting as transfer zones.

Appendix A: Model Scaling and Analog Material Characterization
Models were scaled to achieve a correct geometric, dynamic, rheological, and kinematic similarity (Hubbert, 1937;Ramberg, 1981) (Table 1). Specifically, the length and thickness (geometrical similarity) of models were scaled using a scaling ratio (defined as l* = l model /l nature and h* = h model /h nature ) of 6.67 × 10 −6 , such that 1 cm in the model corresponds to 15 km in nature. The modeling procedure was developed in normal gravity so that the gravity scaling ratio (g*) between the model and nature was equal to 1. The density of the mixture of granular material (quartz and K-feldspar) and silicone and corundum (1,440 kg m −3 ) implied a density scaling ratio (ρ*) of 0.53, considering a natural rock density of~2,700 kg m −3 leading to a stress scaling ratio (σ* = ρ*g* l*) of 3.53 × 10 −7 . The coefficient of internal friction (μ) requires to  Morley (1999a), Morley and Ngenoh (1999), and Vétel et al. (2005) and structures in panel (c) from Muirhead and Kattenhorn (2018). The trend of regional extension direction (σ 3 ) is from Muirhead and Kattenhorn (2018). Orange dots in panel (a) show historical earthquakes (time interval 1900-2020). Seismicity and focal mechanisms for earthquakes with M ≥ 3 are plotted from the USGS Earthquake Catalog available athttps://www.usgs.gov/ natural-hazards/earthquake-hazards/earthquakes. (d) Model RP-7 shows inherited faults (α = 60°) with a trend similar to those observed in the Magadi Basin (α ≈ 65°). Model RP-7, particularly the area displayed in panel (e), shows a fault reactivation pattern similar to the one observed in the Magadi basin and the Kordjya Fault System. Model RP-7 has been rotated to fit the natural setting.

10.1029/2020TC006211
Tectonics be similar both in the analog material and the natural rock, implying a μ*~1. In order to obtain dynamic scaling of models, cohesion ratio (c*) needs to equal the stress ratio (σ* = c*), a condition that is fulfilled as c* and σ* share the same order of magnitude (10 −7 ) ( Table 1).
Scaling of the ductile layer simulating the LC implies the consideration of both viscosity and the strain rate. The LC was simulated with a mixture of PDMS and corundum sand. Corundum was added to PDMS in a quantity necessary to reach the correct density value of 1,440 kg m −3 . In order to obtain a proper density of the PDMS-corundum mixture comparable to the density of the model UC (1,440 kg m −3 ), we empirically estimated a proportion of PDMS-corundum of~1:1 in weight. The viscosity of this material was measured using a coni-cylinder viscometer apparatus.
The viscosity of the model ductile lower crustal layer is relatively high, leading to a rather low brittle/ductile strength ratio. This will influence to some extent both the spacing of faults, which is expected to result in a larger number of faults and a wider region of deformation, and surface topography relief, which will be somewhat less amplified.
Through progressive steps, we calculated the shear-stress versus strain rate curve for the PDMS-corundum material ( Figure A1). We obtained the following equation using the regression line (R 2 = 0.99): Y ¼ 2 × 10 −6 X 1:24 where Y ¼ _ ε ( _ ε being the strain rate) and X = τ (τ being the shear stress): _ ε ¼ 2 × 10 −6 τ 1:24 In our case, n = 1.24, and thus we have assumed the PDMS + corundum mixture as a near-Newtonian material. Viscosity is therefore extrapolated at model strain rate ( _ ε mod ) according to the following equations: η ¼ τ=_ ε mod _ ε mod ¼ 2 × 10 −6 τ 1:24 _ ε mod was calculated from our models as the maximum engineering strain rate considering the linear deformation and the velocity at the proximal, left-hand side of our model setup (see Figure 2).
The dynamic similarity of the viscous layer can be effectively defined by the dimensionless ratio R m of gravitational to viscous forces (Ramberg, 1981), termed "Ramberg number" (R m ) by Weijermars and Schmeling (1986): where ρ d and z d are the density and thickness of the ductile layer, g is the gravitational acceleration, η is the viscosity, and _ ε is the strain rate given by the ratio between the velocity v and the thickness of the ductile layer. In nature, R m can be estimated to be about 1.6 (considering ρ n = 2,900 kg m 3 , g = 9.81 m s −2 , h = 15 km, η n = 10 22 Pa s, and v n = 12.5 mm yr −1 ), while in the Model R m can be estimated at 1.5 (considering ρ m = 1,440 kg m 3 , g = 9.81 m s −2 , h = 0.01 m, η m = 1.69 × 10 5 Pa s, and v m = 20 mm hr −1 ; Table 1).

Tectonics
The dynamic similarity of the brittle layer can be estimated using the ratio S m of gravitational to frictional forces in model and nature (Corti et al., 2003), which is essentially equivalent to the Smoluchowsky number (S m ) defined by Ramberg (1981): where ρ b and z b are the density and thickness of the brittle layer, g is the gravitational acceleration, c is the cohesion, and μ is the internal friction coefficient. In nature, S m can be estimated to be about 1.6 (considering ρ n = 2,700 kg m 3 , g = 9.81 m s −2 , μ = 0.6, h = 15 km, and c o ≈ 10 7 Pa), while in the model, S m can be estimated at~1.5 (considering ρ n = 1,440 kg m 3 , g = 9.81 m s −2 , μ = 0.6, h = 0.01 m, and c o ≈ 6 Pa). Since both nature and models share very similar R m and S m numbers, the requirement of dynamic similarity is fulfilled.
The Reynolds number, R, is another dimensionless number (ratio of inertial to viscous forces) relevant for the modeling, which is given as (Ramberg, 1981, notations as above) However, R e numbers cannot be equivalent in model and nature, even though a dynamic similarity with natural conditions can be obtained for analog models deformed under a natural gravitational field whether R e < <1 and inertial forces can be considered as negligible, which leads to the "quasi-static" approximation (Hubbert, 1937;Ramberg, 1981). More specifically, the lack of a strict equivalence of Reynolds numbers in model and nature can be disregarded when inertias are negligible, that is when (Ramberg, 1981) "acceleration, in terms of rate of change of velocity, but certainly not in terms of force per unit mass in a body-force