Modes of Propagation of Continental Breakup and Associated Oblique Rift Structures

V‐shaped oceanic propagators are widespread around the world. Their geometry combined with magnetic anomalies associated with their opening shows at the first order that ridge propagation in the third dimension occurs by pulses. In this study we use 3D thermomechanical numerical models to show how oblique kinematic boundary conditions control both the intracontinental rift development and the oceanic ridge propagation. To do so, we apply a shortening velocity boundary condition in the direction perpendicular to the extension for “strong” and “weak” crustal rheologies. Numerical model results highlight the finding that three ridge propagation modes can occur. For low out‐of‐plane velocities (12% to 15% of the extension rate), the ridge propagation is fast (>1.5 cm year−1) and straight. Higher shortening velocities (15% to 17%) lead to a ridge propagation by pulses alternating between fast propagation (~1.5 cm year−1) and stalling phases. Finally, for higher velocities (17% to 20%) a ridge jump propagation mode occurs, localizing a new spreading center between 100 and 200 km far from the initial ridge. We also show that ridge propagation phases are associated with dip‐slip‐dominated deformation, while stalling phases are dominated by strike‐slip deformation. These deformation regimes are marked by structure reorientation, while kinematic boundary conditions remain constant. We discuss these results in terms of plate tectonic reconstructions and regional geological studies.


Introduction
Magnetic anomalies record the past direction and rate of seafloor spreading (Vine & Matthews, 1963). As such, they are our best constrain for global plate tectonic reconstructions (Müller et al., 2008(Müller et al., , 2016Seton et al., 2012). In perfectly cylindrical condition, the orientation of the first magnetic anomaly is parallel to the rifted margin of continents. However, in some areas the first magnetic anomalies intercept the rifted margins of continents at an acute angle instead of being parallel to it. These remnants are precious as they allow us to quantify the rate at which continental breakup propagates along the strike of the mid-oceanic ridges (Courtillot, 1982). Indeed, using only basic kinematic considerations, it is clear that if the propagation occurs at the same velocity as the spreading rate, the angle between the margin and the magnetic anomaly should be roughly 45°. Faster propagation leads to the parallelization of the magnetic anomaly with the margin of the oceanic basin matching with the 2D cylindrical model. Slower propagation produces magnetic anomalies that normalize to the rifted margin, until the limit case of a transform margin (no ridge propagation at all).
In a map, these phases of sluggish ridge propagation produce V-shaped margins. V-shaped propagators are found in many places on Earth (Figure 1), and seafloor age map compilation by Müller et al. (2008) is often sufficient to unravel the first-order features of the timing of continental breakup at the regional or more global scale. Examples displayed in Figure 1 show propagators separated by very long segments where magnetic anomalies are parallel to the margin. This indicates that pulsation in the rate of continental breakup propagation first introduced by Courtillot (1982) in the Gulf of Aden is a major feature of divergent plate boundaries. Despite the well-documented geological record of these changes of rate, and their implications for plate tectonics and kinematic reconstructions or the thermal regime of oblique margins, the dynamics of the propagation at long/geodynamic time scale remain poorly studied or understood.
Kinematic interpretations of V-shaped propagators are based on models established in the late 1970s and early 1980s, which are still in use today. The first model, called the rotational extension (Hey et al., 1980;Molnar et al., 2017;Mondy et al., 2018;Zwaan et al., 2020), assumes that seafloor spreading occurs near  (Ryan et al., 2009) with oceanic crust ages of (a) the Gulf of Aden (Fournier et al., 2010;Leroy et al., 2012), (c) the South China Sea (Le Pourhiet et al., 2018;Li et al., 2014), (e) the Red Sea (Almalki et al., 2016;Bosworth et al., 2020), and (g) the North Atlantic Ocean (Klitgord & Schouten, 1986). Main transform fault zones are reported with thick black lines. GPS velocities are from ArRajehi et al. (2010) and Reilinger and McClusky (2011). Graphs with red curves represent the ridge propagation in (b) the Gulf of Aden, (d) the South China Sea, (f) the Red Sea, and (h) the North Atlantic Ocean interpreted from the oceanic seafloor ages distribution along strike. Topographic and bathymetric maps made with GeoMapApp (www.geomapapp.org)/CC BY.
the Euler pole and is therefore described by very large variations in the opening rate in the direction of the ridge. In this case, the first magnetic anomaly forms parallel to the rifted margin, and then the linearly increasing rates of opening with distance from the Euler pole produce a series concentric V's over time Figure 2c. This model infers that very rapid changes in a tectonic regime must occur at the tip of the V-shaped ocean domain. The second kinematic model, published by Vink (1982), considers that the divergence remains constant in the direction of the rift axis, but it is variably partitioned between the oceanic accretion and the continental rifting ( Figure 2b). This model applies only when magnetic anomalies indicate that the Euler Pole is far enough so that the rate of divergence is constant at the regional scale. This is generally the case for a distance greater than 3,000 km or 30°from the Euler pole. The Vink model allows the formation of V-shaped propagators when the propagation rate is similar to the oceanic accretion rate. For faster propagation, magnetic anomalies are parallel to the margin ( Figure 2a). The Vink model is the basis for Courtillot's (1982) proposal (Figure 2d), who interpreted the change in angle between magnetic anomalies and the continental margins along the direction in the Gulf of Aden as the geological record of pulsating continental breakup propagation. It should be noted here that these two kinematic models are not contradictory but rather complement each other as two end-member cases. However, the prediction of Vink's and Courtillot's models in terms of magnetic anomalies and associated tectonic structures at the front of the propagator are different and therefore permit the description of different natural examples of V-shaped propagators. The Woodlark Basin only opens 500 km from the Euler Pole. It displays a strong segmentation, variations from extensional to strike-slip structure (Holm et al., 2016;Mondy et al., 2018), and only the youngest ridge ages at the tip that conform with characteristics of a rotational extension opening basin. In contrast, the South China Sea that opened at about 4,500 km from its pole of rotation (Briais et al., 1993;Le Pourhiet et al., 2018;Mazur et al., 2012) with very little variations in the spreading rate along the strike would better illustrate the Vink kinematic model. Dynamically speaking, on the one hand, the rotational opening model has been reproduced by a number of dynamic simulations based on either analogue materials (Molnar et al., 2017;Zwaan et al., 2020) or numerical methods (Mondy et al., 2018); on the other hand, almost no properly scaled dynamic simulations, including numerical (Allken et al., 2012;Beaussier et al., 2019;Le Pourhiet et al., 2017;Liao & Gerya, 2015) or analogue models (Ding & Li, 2016;Zwaan et al., 2016), reproduces the kinematics proposed by Vink (1982) at the lithospheric scale. With these boundary conditions and a realistic viscosity in the asthenospheric mantle, lithospheric necking should propagate approximately 10 times faster than plates diverge (Allken et al., 2012;Le Pourhiet et al., 2018), resulting invariably in the fast propagation kinematic displayed in Figure 2a. Therefore, to produce a kinematic solution, which resembles Vink's (1982) model ( Figure 2b)  (Vink, 1982), (b) slow ridge propagation (Vink, 1982), and (c) scissor opening (Hey et al., 1980). Dynamic models of the formation of oceanic V-shaped propagators and continental rifts: (d) strong ribbon in the continental crust slowing down ridge propagation (Courtillot, 1982), (e) offset propagators (Le Pourhiet et al., 2017), and (f) shortening perpendicular to extension direction (Le Pourhiet et al., 2018). and a number of geological observations (Figure 1), some additional ingredients must therefore be added to the dynamic simulations to stop or slow down the propagation in order to form a V-shaped oceanic propagator.
The obvious candidate to slow down the propagation of continental breakup is the presence of a strong rheological heterogeneity, the so-called strong ribbon (Figure 2d). This hypothesis was first proposed by Morgan and Parmentier (1985) based on the theory of crack propagation. It was supported by small-strain numerical simulations of the propagation of continental rifting published by Van Wijk and Blackman (2005). The results of this first study are consistent with the results presented by Allken et al. (2012): the propagation rate of a neck (continental rift) in a frictional layer (upper crust) lying above a viscous layer (lower crust) decreases with the viscosity of the lower layer. Many similar analogue or numerical experiments (Allken et al., 2012;Brune et al., 2017;Corti, 2012) have proved that increasing the viscosity of the lower crust resulted in a slower lengthening of grabens in the upper crust. However, at the lithosphere scale, there is nothing like a high-viscosity asthenosphere to delay or slow the propagation of continental breakup. At that scale, the strength of the lower crust has been shown to enhance fast lithospheric necking and tends to increase the rate of propagation of continental breakup (Benes & Scott, 1996;Le Pourhiet et al., 2018).
Oblique continental rifting controlled by the presence of preexisting viscous weak zones oblique to the extension direction creates en échelon rift basins, which are separated by the transfer zones that accommodate strain differences between segments (Ammann et al., 2017;Brune, 2014;Duclaux et al., 2020;Mart & Dauteuil, 2000;Mcclay & White, 1995;Molnar et al., 2017). These isolated basins form small segments of oceanic accretion before complete lithosphere continental plate breakup is achieved. As they grow in width, they also grow in length, propagating from the basin center toward both the extremities. Yet their propagation is limited by the compressional stress buildup that arises from the growth of a nearby segment. Le Calvez and Vendeville (2002) proposed a model setup to study in detail the dynamic of linkage between en échelon segments. This setup consists of two offset weak zones that grow laterally in map view and link or not depending on spacing, frictional layer thickness, and viscosity of the underlying layer. Allken et al. (2012) presented the result for a similar set of numerical experiments and showed the importance of the rate of softening in the frictional layer to produce a transform fault. More recently, Le Pourhiet et al. (2017) have extended this model to study continental breakup at the lithosphere scale (Figure 2e). At that scale, the linkage between segments can be delayed by tens of millions of years. Meanwhile, two asymmetric V-shaped propagators form, as shown in Figure 2e. Full linkage occurs only after oblique necking of the mantle lithosphere is reached, showing similarities with the initial conditions of the Ammann et al. (2017) study in which an oblique weak zone is imposed. These asymmetric propagators better fit the observations in the Gulf of Aden or in the equatorial Atlantic (Le Pourhiet et al., 2017) than the original propagator models do (Figures 2a and 2b), but they do not apply to all geological examples of V-shaped propagators.
To produce oblique extension, pure shear boundary conditions in map view similar to the Alboran Sea setup in Le  (i.e., a stretching velocity in one horizontal direction and a shortening velocity in the other horizontal direction) are a good alternative to the introduction of oblique heterogeneities in cylindrical boundary conditions. Indeed, strain localization causes the formation of structures that can rotate internally within the pure shear far-field strain and that favor the emergence of an instable four-plates system that stabilizes into stable three-plates systems (Gerya & Burov, 2018). It has been demonstrated that pure shear conditions produce V-shaped propagators both during the reorganization of oceanic plates (Gerya & Burov, 2018) and during continental breakup (Le Pourhiet et al., 2018). This results together with the analogue experiments of Ding and Li (2016), who introduced a gradient in gravitational potential energy facing the propagation direction, suggests that applying a small amount of shortening or compression out of the plane of the main driving extension slows down the continental rifting propagation rate. One should note that in 3D the transition from a transtensional to transpressional regime in the pure shear map setting is related to change in the surface of the model and that as long as the surface remains in the dilatational regime, no shortening-related structures form (Le . GPS velocity vectors plotted relative to the fixed Arabian plate in Figure 1 clearly show frontal displacement toward an actively propagating oceanic ridge existing in active propagating continental rifts of the Gulf of Aden and Red Sea. One might suspect the influence of the Afar plume in the first case and the vicinity of the Mediterranean subduction zone in the second case. In both cases, the frontal shortening is not associated with shortening tectonic structures, which could "fossilize" its existence. Here we therefore intend to focus our attention on this still largely unexplored direction. We adopt the same boundary conditions as in Le Pourhiet et al. (2018) (Figures 2f and 3) in order to analyze the sensitivity of the continental breakup propagation rate to both the rate of ridge-parallel shortening and lower crustal rheology. We first extract from the simulations synthetic propagation rates and expected magnetic anomalies in opening basins that correspond to data that are widely available and relate them to large-scale plate deformations. We then focus on the tectonic structures that accommodate the frontal shortening in the models in order to discuss what are key tectonic observations that might help to decipher past kinematics and lithospheric dynamics from seismic lines and thermal data sets.

Modeling Approach
In order to model the long-term deformation of the lithosphere, we use pTatin3D (May et al., 2014(May et al., , 2015, a highly scalable, massively parallel implementation of the finite element method that employs an arbitrary Lagrangian-Eulerian discretization together with the material point method to solve the conservation of momentum for an incompressible fluid: where η is the nonlinear effective viscosity, _ ε the strain rate tensor, P the pressure, ρ the density, g the gravity acceleration vector, and v the velocity vector of the fluid. This Stokes flow is coupled with the time-dependent heat conservation equation: where T is the temperature, v the velocity vector of the fluid, κ the thermal diffusivity, H the heat source (here only radioactive heat sources are considered), and Cp is the heat capacity. The material density is considered to vary with temperature such as where ρ 0 is the initial material density and α the thermal expansion coefficient.
The Stokes problem is solved using a Q 2 -P 1 discretization, while the heat equation is discretized on Q 1 elements. pTatin3D uses a free-surface boundary condition that can dynamically evolve due to deformation and surface where h is the altitude and k the diffusion coefficient. The equation is solved using a first-order explicit time integration.

Rheological Model
To simulate the long-time behavior of the lithosphere, we use temperature-and pressure-dependent nonlinear rheologies. Brittle parts of the lithosphere follow the Drucker-Prager pseudo-plastic yield criterion adapted to continuum mechanics: where C is the cohesion, ϕ the friction coefficient, P the pressure and _ ε II the second invariant of the strain rate tensor. To simulate the friction drop in deforming faults, we use a simple linear decrease of the friction angle from 30°to 5°with accumulated plastic strain from 0 to 1. However, laboratory experiments show that under high confining pressures (>1 GPa), brittle behavior changes to plastic behavior (e.g., Kameyama et al., 1999;Precigout et al., 2007). We account for that change by limiting the Drucker-Prager yield stress to a maximum deviatoric stress of 400 MPa according to the calibration performed by Watremez et al. (2013).
Viscous parts of the lithosphere are simulated using an Arrhenius flow law for dislocation creep: where A, n, and Q are material-dependent parameters (see Table 1); R is the gas constant; and V is the activation volume. Quartz (Ranalli & Murphy, 1987) Anorthite (Rybacki & Dresen, 2000) Olivine (Hirth & Kohlstedt, 2003) A MPa −n s

Initial Conditions
Our experiments represent a 1,200 (in the x direction) × 600 (in the z direction) × 250 (in the y direction) km 3 domain discretized with 512 × 256 × 128 elements, respectively ( Figure 3a). In the y direction, two thirds of the 128 elements are concentrated in the first third of the 250 km, giving approximately 8 km 3 (2 × 2 × 2 km 3 ) elements in the first 80 km. Four rheological layers subdivide the model domain to represent the upper and lower crust (each of 20 km thick), the lithosphere mantle (80 km thick), and the asthenosphere (170 km thick). The Moho is located at 40 km depth, and the lithosphere-asthenosphere boundary is at 120 km depth ( Figure 3a). The lithosphere mantle and asthenosphere are both simulated using a dry olivine flow law (Hirth & Kohlstedt, 2003) (Table 1).
Depending on the model, the crust rheology may vary. The rheology of the lower crust has been proven to be of first-order importance for the strain localization and propagation in 3D (Allken et al., 2012;Brune et al., 2017;Corti, 2012;Le Pourhiet et al., 2014. Therefore, in this study we use two different crustal lithologies denoted as weak and strong (Figure 3(c)). The weak crust corresponds to an upper and lower crust modeled with a quartz rheology (Ranalli & Murphy, 1987), while the strong crust corresponds to an upper crust modeled with a quartz rheology and the lower crust with an anorthite rheology (Rybacki & Dresen, 2000). We modeled those two rheological layering with the same boundary and initial conditions in order to assess the role of the lower crust on the evolution of the rifting and continental breakup.
The initial temperature field is computed with a steady state analytical solution (Turcotte & Schubert, 2002) for a T y¼0 of 0°C, an incoming mantle heat flux q m of 20 mW m −2 , a radiogenic heat production H of 1.2 × 10 −6 W m −3 , a characteristic radiogenic layer y p of 40 km, and a conductivity of 3.3 W m −1 K −1 . This analytical solution gives 610°C at 40 km (depth of the Moho; Figure 3c) and 1300°C at 120 km (depth of the lithosphere-asthenosphere boundary; Figure 3c). Then from 120 to 250 km, we compute a linear increase of temperature for an adiabatic gradient of 0.5°C/km ( Figure 3c). Although this second part of the geotherm is not steady state, the cooling by diffusivity is very slow (less than 2°C Myr −1 for the maximum cooling rate), and it allows keeping reasonable conductivity values in the asthenosphere (3.3 W m −1 K −1 ).
The initial radiogenic heat production is set as an exponential decay of heat production with depth according to Turcotte and Schubert (2002), as follows: for a surface production H 0 of 1.2 × 10 −6 W m −3 .
In order to initiate localization on one side of the modeling domain, we define a weak zone in which a random amount of initial plastic strain ranging from 0 to 1 is applied to the particles located in a cube of coordinates x ¼ 550 to x ¼ 650, y ¼ 0 to y ¼ 100, and z ¼ 0 to z ¼ 100 instead of the 0.0 to 0.03 background value.
Since we use plastic strain softening, the initial amount of plastic strain will decrease the initial friction angle in the concerned area as where ϕ is the friction angle; ϕ 0 the initial maximum friction angle (30°); ϕ ∞ the minimal friction angle reached after softening (5°); ε max and ε min are the bounds of maximum and minimum plastic strain amount, respectively, between which softening is applied; and ε p the actual plastic strain computed from Equation 6.

Boundary Conditions
The boundary conditions to solve the heat conservation equation are T ¼ 0°C at the top boundary (y ¼ 0 km), T ¼ 1440°C at the bottom boundary (y ¼ 250 km), and null heat fluxes on all vertical boundaries of the domain. To solve the conservation of momentum, we impose a free surface on the top boundary, free-slip boundary conditions on the face of normal z at coordinates x, y, z ¼ 0, and Dirichlet boundary conditions on the normal component of velocity on all other vertical sides. An extension velocity v x ¼ 5 mm year −1 is imposed on both faces of normal x and a shortening velocity v z is imposed on the face of normal z at coordinates x, y, z ¼ 600 (Figures 3a and 3b). This shortening velocity may vary between models (from 1.2 to 2.0 mm year −1 ). It introduces a part of constriction in the model domain that intrinsically modifies the ratio between normal stretching and thinning as compared to a cylindrical model. In this range of shortening rates, propagation rates are indeed of the same order of magnitude (less than 1.5 cm year −1 ) of the plate velocity imposed in the model (1 cm year −1 ), and it is for this range of velocities that we expect the highest dependence of deformation pattern to the rheology of the lithosphere (Le Pourhiet et al., 2018).

Postprocessing
In order to compare the numerical model results with natural examples and known geological and geophysical data, we computed an equivalent of seafloor magnetic anomalies. Magnetic anomalies measured on the seafloor are used to model seafloor spreading rates and oceanic crust ages in oceans (Müller et al., 2008). In our models, the formation of oceanic crust by partial melting of the mantle and magma crystallization is not implemented. Therefore, we compute an age of mantle exhumation below the 800°C isotherm, which can be roughly approximated as the Curie temperature of a ferromagnetic material (Rushbrooke & Wood, 1958). This age corresponds to the time in the model at which the mantle cooled below 800°C. In the models, we use 0 Myr as the beginning, and time goes forward from this initial value, while geological time is counted from present (0 Ma) to past. Therefore, the age of the mantle exhumed during the latest stages of the models is greater than the age of the mantle exhumed in the earlier stages. Those age values that need to be interpreted as the greater are also the younger.
Figures 8, 10, and 12 and supporting information Figures S11 and S12 display the deformation regime computed from the stress tensor orientation, also known as the regime stress ratio (Brune, 2014;Brune & Autin, 2013, and reference therein). The color intensity is related to the strain rate second invariant value.
Figures 4, 5, and 7-9 show the plastic and/or viscous strain deformation. Those values represent respectively the cumulated finite strain computed with Equations 6 and 7. Therefore, the plastic strain displays the deformations related to the frictional (or brittle) parts of the lithosphere, while the viscous strain shows the deformations related to the dislocation creep flow laws in the "ductile" parts of the lithosphere.
Supporting information Figure S1 shows the β factor of the crust computed as where h t¼0 c is the initial thickness of the crust at the beginning of the simulation (h c ¼ 40 km, t ¼ 0 Myr) and h t¼n c is the thickness of the crust at a given time n. Thus, the β factor value represents the thinning ratio of the crust at a given time. Figure 6a displays the position of the tip and the ocean propagation with respect to time. The position of the tip is measured as the higher z coordinates of the mantle-continent transition along the strike of the oceanic ridge. Therefore, it represents a particular area at the frontier between a plate boundary (the oceanic ridge) and an intracontinental deformation zone (the rift). As discussed in the numerical results section 4, this area may evolve with time in a triple junction or remain as a plate boundary to intracontinental deformation transition. The ridge propagation (Figure 6b) is the time derivative of this position, computed as where Vrp is the ridge propagation velocity, z t ¼ n and z t ¼ n − 1 are the z coordinates at time t n and t n -1 ,

10.1029/2020JB019906
Journal of Geophysical Research: Solid Earth respectively. Figure 6c shows the tectonic forces evolution through time. Fx and Fz represent the forces normal to the faces of normal x and z, respectively. To compute forces, we first compute the deviatoric stress tensor as where η is the viscosity and _ ε the strain rate tensor. We then compute the traction along the faces of normal x and z, respectively, as and where x ! and z ! are the unit vectors in the x and z directions, respectively. We then compute the tectonic forces per linear meter as . Map view of nine models with a strong and weak lower crust and varying v z velocities displaying the plastic strain and the mantle exhumation age. The plastic strain represents all the deformations that occurred under the regime described by Equation 6. The mantle exhumation age is set when the mantle cools below the isotherm 800°C.
where Lx and Lz are the length of the domain along the x and z axes, respectively.

Dynamic of Continental Breakup Propagation
Numerical models are aimed at investigating the role of both the crust rheology and the intensity of rift-parallel shortening applied in our simulations. We show here that if we impose a sufficient velocity v z (perpendicular to the extension direction), the propagation rate of the continental breakup is not For the same amount of shortening velocity applied in the direction perpendicular to the extension, the main difference between strong and weak lower crust models resides in the propagation mode. Furthermore, for large shortening velocities, weak and strong lower crust models tend to show the same behavior, forming oblique rifts (Figures 4 and 5). In order to highlight the principal results of this study, we first focus on the ridge propagation mode as a function of rift-parallel inflow (obliquity).

Rheology of the Crust and Sensitivity to the Shortening Velocity
Our numerical modeling results show that the rift-parallel shortening velocity and the crust rheology strongly affect the propagation of the rift and its structure (Figures 4-6). With a strong lower crust, the minimum shortening velocity necessary to form a strike-slip conjugate and isolate a new microplate is 17% of the extension velocity (v z ¼ 1.7 mm year −1 ; Figures 4e and 5e, model 3). With a weak lower crust, the rift propagation starts to be sensitive to the shortening velocity when the latter reaches 12% of the total extension velocity (Figure 4b and 5b, model 1). The influence of this third boundary condition becomes very important for 15% of the extension velocity (Figure 4d and 5d, model 2). For both lower crust rheologies (strong and

10.1029/2020JB019906
Journal of Geophysical Research: Solid Earth weak), the shortening velocity represents a threshold, across which pure extension gives place to a transtensional and strike-slip deformation regime. A strong lower crust favors the formation of three tectonic plates joining at a ridge-ridge-fault triple junction because it allows for a more localized strain than do weaker crust rheologies (Figures 4e, 4g, 4h 5e, 5g, and 5h). Margins are narrow, and few shear zones accommodate the crustal thinning (as shown by the β factor varying from 1 to 5 in less than 50 km; supporting information Figure S1). In contrast, a weak lower crust rheology favors a more distributed deformation (Figures 4b, 4d, 4f, 5b, 5d, and 5f). Wide rifts are formed, and many shear zones develop to accommodate the crustal thinning (β factor of~3 along 300 km from the exhumed mantle; supporting information Figure S1). This distribution prevents the formation of three distinct plates, and the effects of the out-of-plane boundary condition are visible in the rift evolution and propagation. Finally, for both weak and strong lower crusts, a too large amount of shortening (20% for the strong case [ Figures 4h and 5h] and 17% for the weak case [ Figures 4f and 5f]) results in the formation of an oblique continental rift followed by oblique spreading.

Ridge Propagation
The numerical models highlight three types of rifts propagation modes (Figure 4 and Figure 6a). The first mode is a very fast one (from 3 to 1.5 cm year −1 ; Figure 6b) that results in an oceanic opening almost simultaneously along the rift axis. Models with a strong lower crust and a shortening velocity of less than 17% of the total extension velocity (1.7 mm year −1 ) display this evolution (Figures 4a, 4c, and 6).

Journal of Geophysical Research: Solid Earth
The second propagation mode is characterized by a succession of fast propagation and stalling phases ( Figure 6). This mode is principally visible in the weak lower crust models for shortening velocities less than 1.7 mm year −1 (Figures 4b and 4d). Due to the weak lower crust, a long phase of continental rifting occurs (~40 Myr) before ridge propagation takes place (Figure 6d). The propagation velocity starts at about 1.5 cm year −1 before slowing down in~10 Myr to less than 5 mm year −1 (Figure 6b). This propagation velocity drop can be related to crustal thickness variations and competition between extension velocity and shortening velocity. Indeed, the initial continental rift phase efficiently thins the crust along 200 to 300 km (in the z direction) where the velocity field is principally perpendicular to the rift. Therefore, once crustal breakup occurs, the propagation along this already thinned crust is fast. However, once the rift tip reaches a less thinned crust, the propagation stops or slows down to allow the continental crust domain in front of the tip to thin enough before resuming the propagation (supporting information Figures S4, S7, and S9). During the whole simulation time (~120 Myr), models were able to produce two or three stalling phases, each one lasting less than the previous one. This feature can be explained by the decreasing crustal thickness with time in front of the rift propagator.
Finally, a third propagation mode reveals a "jump" of the spreading center from the tip of the stalling rift to a new basin. Models with strong and weak lower crust both present this propagation mode for shortening velocities larger than 1.7 mm year −1 (Figures 4e, 4g, 4h, and 6a). As in the two previous propagation modes, the rift starts with a fast propagation (~1.5 cm year −1 ; Figure 6b) before stalling. But instead of propagating again from where it has stopped, a new spreading center forms 200 to 300 km farther (Figures 4e and 6a). The main characteristic of this propagation mode resides in the deformation regime between the onset of the stalling phase and the formation of a new spreading center. As soon as the along-strike oceanic propagation stops,  Figure 3. The superimposed colors represent the stress-inferred deformation regime (Brune, 2014;Brune & Autin, 2013). The intensity of the color scale is related to the intensity of the strain rate tensor second invariant.

Journal of Geophysical Research: Solid Earth
the intracontinental deformation at the front of the ridge switches from a dip-slip-dominated regime to a strike-slip-dominated one. During this deformation regime, the rifting continues but localizes on an offset branch and leads to the formation of two offset oceanic segments and an ultraoblique margin (Figures 4e,  4g and 4h). Figure 6c shows the tectonic forces Fx (Equation 14a) and Fz (Equation 14b) evolution along boundaries of normal x and z, respectively. Along the faces of normal x, the Fx force corresponds to the extensional forces. Strong lower crust models display initial extensional forces (from 15 to 17 TN m −1 ) greater than weak lower crust ones (from 10 to 12 TN m −1 ) at the beginning of the experiments and through time. The Fx force displays a small dependency to the rift-parallel shortening magnitude at the beginning of the simulations. However, since the rift-parallel shortening influences the strain localization, the force evolution differs between models. At the first order, the Fx force decreases with the lithosphere thinning through time for all models. The force decrease is directly related to the thinning and propagation rates of the models: The faster the lithosphere thins, the faster the force decreases. Then, once the continental breakup occurred along the whole length of the modeled domain, the force reaches or tends to approach a plateau, representing the strength of the lithosphere.

Tectonic Forces Evolution
Along the face of normal z, a first-order observation is that the Fz force is not dependent on the rift-parallel shortening velocity but on the strength of the lithosphere. Indeed, the models with a strong lower crust display an initial state between 4 and 5 TN m −1 , while the models with a weak lower crust show an initial state at~3 TN m −1 . During the deformation of the lithosphere, the Fz force shows few variations related to strain localization events in the models except for models with a strong lower crust and a rift-parallel shortening of 1.2 to 1.5 mm a −1 . These two models display a fast force increase associated with the fast ridge propagation and a plateau between 8 and 9 TN m −1 when the ridge had propagated all along the domain. In absolute values, the Fz force is 2 to 5 times lower than the extensional force Fx.

Characteristic Structures Associated With Propagation of Breakup
In order to decipher the geological record of the different propagation modes and at the same time to gain insight on the transition from straight mode of propagation to partitioned mode, we have chosen to detail the evolution of models with rates of shortening 1.2 mm a −1 (model 1), 1.5 mm a −1 (model 2), and 1.7 mm a −1 (model 3).

Weak Lower Crust, 12% of Shortening Velocity
The first model (model 1) is characterized by a weak lower crust and a shortening velocity v z ¼ 1.2 mm year −1 (12% of the total extension velocity). Continental rifting lasts 40 Myr before continental breakup occurs at the southern boundary of the model (Figures 7a and 8a). During this initial phase, normal faulting is perpendicular to the extension direction. The weak lower crust results in strain distribution over a wide continental rift (Figures 7a and 7b, cross sections 2 and 3, and Figures 8a and 8b). At 50 Myr, the continental breakup reaches z ¼ 200 km, and a localized strain band replaces the wide strain distribution in the exhuming mantle between z ¼ 0 and z ¼ 200 km (Figures 7b, cross section 3, and Figure 8b), while in the rest of the model continental rifting is still active (Figure 7b, sections 1 and 2). At 60 Myr, the domain of exhumed mantle reaches z ¼ 300 km, and the oceanic rift enters a second phase characterized by short episodes of propagation, and stalling takes place. Propagation first slows down for 10 Myr, while the stretching of the continental crust increases in a narrow domain and finally localizes as crustal breakup occurs, allowing oceanic rift propagation to continue (Figures 7c and 7d). During this phase, a minor strike-slip crustal deformation initiates at the edge of the propagator tip (coordinates: 300 < z < 400; 700 < x < 800; Figure 7d, cross section 2, and Figure 8d) that accommodates the differential movement between the stalling ridge and ongoing crustal extension away from the ridge. Then, from 80 to 100 Myr, the ridge propagation resumes, and as the mantle is exhumed, the active deformation in the crust progressively vanishes (Figures 7e-7g and 8e-8g). The age map of mantle exhumation displays a typical pattern where the oldest ages are decreasing in the direction of propagation (Figure 7g). At the end of the model, the angle between the oceanic rift and the continental margins varies along strike. This angle variation developed during the slowing down of the propagation phase. The direction of normal shear zones changes in time and along the rift axis. The first faults formed during the continental rift phase have been rotated during the first stalling phase (between 50 and 60 Myr; Figures 7c and 7d). Then, as the rift propagates, the shear zones formed in the crust continue to rotate and align with the rifted margins.

Weak Lower Crust, 15% of Shortening Velocity
The second model (model 2) also presents a weak lower crust but with a velocity v z ¼ 1.5 mm year −1 (15% of the total extension velocity). Here again continental rifting takes place during approximately 40 Myr before reaching the crustal breakup and mantle exhumation (Figures 8a, S5a, and S6a). Like in the previous model, a wide rift develops, but deformation appears to be more localized. At 50 Myr, the ridge is propagating northward (up to z ¼ 200), and normal faulting is reduced (Figures 9b and 10b). At this stage, the propagation starts to slow down, and the ridge propagation enters a stalling phase until 80 Myr. During these 30 Myr, the exhumed mantle domain widens, forming stripes of exhumation ages parallel to the direction of the extension but perpendicular to the continental margins located at the edge of the propagator (Figures 9b-9d). This stalling phase results in drastic modification of the rift architecture. The rifted margins located around the tip, and normal faults rotate until they form an angle of 90°with the active ridge (Figures 9c-9e and 10c-10e). As a result, the rifted margins located at the propagator tip are striking parallel to the extension direction. This rotation is accompanied with the development of a left-lateral strike-slip fault zone in the crust (Figures 9c-9e and 10c-10e). Transcurrent deformation first develops at around 50 Myr as ridge propagation stops ( Figure 10b) and evolves as a highly localized crustal-scale strike-slip shear zone (Figure 10e). The rotated normal faults of the continental rift reach an obliquity of~45°with respect to the stretching direction and ultimately form releasing bends aligned in the strike-slip fault direction. During this phase, the continental crust is thinned under a transtensional strain regime, exhuming the lower crust in elongated dome structures parallel to the stretching direction as described on the field (e.g., Jolivet et al., 2004) and reproduced previous modeling studies (e.g., Le Pourhiet et al., 2012, 2014. Around 80 Myr, the propagation resumes and the exhumed mantle domain opens obliquely to the stretching direction along the former strike-slip zone (Figures 9e and 10e). At 90 Myr, the ridge propagation slows down once again, and a second left-lateral strike-slip structure develops (Figures 9f and 10f). The most remarkable result is that the rift orientation change associated with the development of strike-slip deformation occurs spontaneously without involving changes in plate kinematics (Figures 9 and 10).

Strong Lower Crust, 17% of Shortening Velocity
The third model (model 3) involves a strong lower crust and a velocity v z ¼ 1.7 mm year −1 (17% of the total extension velocity). In this model, the continental rifting phase stands for less than 30 Myr before crustal breakup occurs (Figure 6d). Crustal deformation is also more localized than in the previous models (Figures 11 and 12). At 30 Myr, the ridge reaches z ¼ 200 km (Figures 11a and 12a), and its propagation slows down. This propagation rate change coincides with the formation of two new deformation zones ( Figure 12a): a right-lateral strike-slip structure to the left and a left-lateral transtensional shear zone to the right. At 40 Myr, the ridge propagation enters in a stalling phase and stops (Figure 11b). The two branches localize now all the intracontinental deformation at the front of the ridge and a ridge-rift-fault triple junction that rapidly evolves in ridge-ridge-fault forms at the tip of the propagator (Figure 12b). The strike-slip motion in both branches intensifies and becomes predominant (Figures 12b and 12c) until a new exhumed mantle domain appears at 60 Myr (800 < x < 900; 400 < z < 500; Figure 11d, cross section 1). Between the initial propagator and this new basin, a thin band of continental crust remains before the two exhumed mantle domains join at 70 Myr (Figure 11e, cross section 2). As soon as the new exhumed mantle domain forms, the deformation regime on the right branch changes once again from a strike-slip regime to a transtensional extension ( Figure 12e). As in the second model (weak lower crust, v z ¼ 1.5 mm year −1 ), the strike-slip deformation in the right branch accommodates the rotation of the continental margin previously located at the edge of the initial propagator. The rifted margin not only rotates but also thins during the transtensional deformation phase, before being torn apart when both exhumed mantle domains join (Figure 11e, cross section 2). The finite rotation results in a rifted margin oriented parallel to the extension direction, while it was formed perpendicularly to it. This process may be responsible for the initiation of a transform margin. Here again, the lower crust exhumes preferentially in areas where strike-slip deformation is dominant. In the left branch, the right-lateral strike-slip deformation is active from 30 Myr to the end of the simulation (Figures 11 and 12). This

Journal of Geophysical Research: Solid Earth
localized vertical shear zone progressively offsets the former rifted margin structures (Figure 11, cross section 3) and forms a vertical margin where the transition from continent to exhumed mantle is sharp (Figures 11e and 11f, cross section 4). Compared with that in the previous models, the lithosphere-scale strain partitioning is a main characteristic of model 3 and results in the formation of a triple junction and three distinct tectonic plates. As in model 2, the change in deformation regime and rift orientation is provoked by the stalling phase, while kinematic boundary conditions are kept constant.

Geodynamic Interpretation of Rift-Parallel Inflow
Lithospheric modeling aims at understanding the dynamics that lead to the emergence of new localized plate boundaries at the regional scale. These plate boundaries are of three types (fault, ridge, and trench), which can be studied separately. But plate tectonics usually deal with triple junctions and the velocity field surrounding the triple junction is not cylindrical in a map; it is rather a combination of pure and simple shear components. The geodynamic pertinence of noncylindrical boundary conditions makes therefore no doubts if one wants to study the emergence of a triple junction.
In terms of modeling, Le  have shown that transtensional tectonic structures (A-type domes) emerging from pure shear boundary conditions are similar to those produced by imposing simple shear boundary conditions. Yet pure shear boundary conditions have the advantage over simple shear boundary conditions because they allow the formation of large strike-slip structures without having to impose a step function in the boundary conditions of the models, that is, without prescribing the location of the strike-slip faults or their occurrence at all. They are therefore useful to studying the emergence or not of a localized strike-slip plate boundary rather than an oblique rift.
Finally, as for all models run with kinematic boundary conditions, one needs to discuss the forces needed to deform the models. Our results and especially that of the force calculation are unexpected, which needs to be discussed.
For 3D deformation patterns to emerge within cylindrical boundary conditions, one must include a form of inheritance in the initial conditions. Models having a large area with random noise (Duclaux et al., 2020;Naliboff et al., 2020) demonstrate that the offset of oblique structures in cylindrical conditions is not very  Figure 3. The superimposed colors represent the stress-inferred deformation regime (Brune, 2014;Brune & Autin, 2013). The intensity of the color scale is related to the intensity of the strain rate tensor second invariant.
In this study, we used the Le Pourhiet et al. (2018) kinematic boundary conditions to impose obliquity. As for pure shear boundary conditions, we showed that accounting for sufficient rift-parallel shortening leads to the formation of a triple junction between three tectonic plates in models with a strong lower crust. The coupled rheology of the crust seems to be a first-order constraint to produce a large-scale strike-slip fault in extensional contexts as also shown in analogue models (Zwaan et al., 2019).
In terms of tectonic forces, the models show that the forces are more dependent on the rheology of the lithosphere than on the rift-parallel shortening magnitude. In a weak lower crust lithosphere, the force associated with the rift-parallel shortening (Fz) is in the same order of magnitude as the force generated by the gravitational potential energy (Bird et al., 2006). For a strong lower crust lithosphere, the gravitational potential energy may not be sufficient to explain such rift-parallel shortening, and an active tectonic process like a subduction or a collision may be required in the vicinity of the rift to allow rift-parallel shortening. However, the force generated by this rift-parallel shortening is 10 times lower than the extensional force and likely represents the resistance of the strike-slip fault developing in strong lower crust models, which is several hundreds of kilometers long. Such a force may be lowered by a more efficient strain softening in the brittle crust.

Ridge Migration and Comparison With Previous Modeling Work
Models presented in this study show that a small amount (12% to 20%) of rift-parallel shortening leads to alternation of ridge propagation phases and ridge stalling phases. These alternating phases are each associated with strain regime changes. During ridge propagation, the deformation in the continental crust located at the front of the propagator is dip-slip dominated, while during stalling phases the deformation in the continental crust is strike-slip dominated. According to our model results, the ridge propagation occurs under three different mechanisms (Figure 6 and 14). For low rift-parallel shortening velocities (12% to 15% of the extension rate; Figures 4a-4c), the ridge propagation is fast (between 1.5 and 3.5 cm year −1 ; Figure 6b) and straight. Therefore, in this case, the ridge propagation velocity exceeds the total extension rate (1 cm year −1 ). This propagation mode has been documented in crustal-scale models involving a low-viscosity layer beneath a plastic layer (Allken et al., 2012), where the lower viscosities lead to faster propagation. The same behavior is observed in lithospheric-scale models with cylindrical boundary conditions (Le Pourhiet et al., 2018), where free-slip boundary conditions prevent deformation in the rift-parallel direction (which is equivalent to a boundary opposing no resistance to propagation) but with the mantle acting as the viscous layer controlling the propagation rate.
Higher rift-parallel shortening velocities (15% to 17%; Figure 4d) lead to ridge propagation by pulses alternating between fast propagation (between 0.5 and 1.5 cm year −1 ; Figure 6b) and stalling phases. This propagation mode obviously requires that the ridge propagation stall at some point. Until now, this behavior has been modeled only during soft linkage of en échelon propagators (Ammann et al., 2017;Le Pourhiet et al., 2017;Liao & Gerya, 2015) and a rift-parallel inflow boundary condition model (Le Pourhiet et al., 2018, and this study). The common feature between these models is that during the stalling phase the velocity field at the tip of the propagator changes and leads to a strike-slip-dominated deformation. This mode of propagation was also proposed by Courtillot (1982) in the Gulf of Aden to explain the obliquity between seafloor magnetic anomalies and the passive margin. In Courtillot's (1982) model, a hard node in the lithosphere structure would explain why the ridge propagation stopped. All these thermomechanical and conceptual models show that to enter in a stalling ridge phase, a noncylindrical parameter is required. It can be rheological (Courtillot, 1982) or kinematic (Ammann et al., 2017;Le Pourhiet et al., 2017Liao & Gerya, 2015).
Finally, for higher velocities (17% to 20%; Figures 4e-4g) a ridge jump propagation mode occurs. This ridge propagation mode also requires a ridge stalling phase, but rather than resuming from the ridge tip, a new spreading center localizes between 100 and 200 km away from the ridge (Figures 6a and 11d). This propagation mode is particularly well observed in models evolving toward three distinct tectonic plates (Figures 4e,  4g, and 4h), therefore for highly oblique kinematic conditions. A propagation by jump was also obtained in analogue models with rotational extension coupled with an oblique weak seed (Molnar et al., 2017) where only highly oblique weak zones lead to jumps of spreading centers. Effectively, a weak zone oblique to the ridge propagation direction has the same effect as a strong zone or an inflow in the rift-parallel direction. In the first case, it will guide the deformation, while in the second case, it will force the deformation to propagate elsewhere.
These different propagation modes are highly dependent on the lithosphere's ability to deform in the rift-parallel direction, whether it is due to rheological heterogeneities (e.g., Courtillot, 1982;Molnar et al., 2017) or to velocity field variations (e.g., Ammann et al., 2017;Le Pourhiet et al., 2017.

Comparison With Natural Cases
Although the four oceanic basins presented in Figure 1 are not at the same time and space scales, they all share similar features between each other and with our models. All of them are well-documented examples of ridge propagation characterized by pulses of alternating fast propagation along straight margins and stalling phases, which often corresponds to the location of oblique and transform margins. Our numerical simulations show that during propagation phases, the deformation occurs orthogonal to the rift axis, whereas oblique extension and strike-slip deformation characterize stalling phases. Initiation of strike-slip deformation results in the reorientation of preexisting rift normal faults. As stretching continues, an oblique continental rift forms, which may result in the formation of a transform margin or oblique continental breakup propagation (Figure 14). A first-order observation is that as in the natural systems described in this study, none of the models involving rift-parallel shortening produce shortening structures (as thrusts or folds). However, the noncylindrical component of the velocity field is well expressed in the GPS velocities of active oblique rift systems (Gulf of Aden and northern Red Sea basins).
Among the four oceanic basins in Figure 1, the South China Sea is probably the most compelling example of secondary oblique continental breakup propagation following a stalling phase. Figure 1b outlines three V-shaped propagators that indicate that continental breakup propagated from east to west by pulses of rapid propagation, separated by two periods of stalling. The longest period occurred from 32 to 23 Ma, where magnetic anomalies in the propagator region form at a 60°angle with the continental margin indicating that the propagation rate is much lower than the spreading rate. Le Pourhiet et al. (2018) proposed that rift-parallel inflow could explain this first phase of stalling using a model similar to the one with a weak crust and 1.2 mm year −1 rift-parallel inflow. They then posit that this long stalling period ends when both the spreading and propagating direction of the mid-oceanic ridge changed by 15° (Sibuet et al., 2016) in the South China Sea (Figure 1b). The origin of this rift-parallel inflow has been proposed as the expression of gravitational potential energy gradient (e.g., Ding & Li, 2016;Le Pourhiet et al., 2018) between the developing basin and the higher Indochina to the west. Bird et al. (2006) showed that a topographic load can propagate vertical forces of 1 to 2 TN m −1 over hundreds of kilometers. Here with a more thorough exploration of the effect of shortening on the propagation of continental breakup, we show that increasing the shortening rate to 1.5 mm year −1 reproduces the 15°rotation of the spreading direction reported by Sibuet et al. (2016) and the propagation of a younger oceanic propagator self-consistently. Moreover, associated with the strike-slip deformation, lower crust domes exhume in the continental margin as in the South China Sea. This finding questions whether magnetic anomalies do record a change in far-field conditions or simply the localization of the deformation on an oblique continental rift branch. This result may have significant implications on plate reconstruction and geodynamic interpretations. Indeed, changes of structural directions are common during continental rifting. These phases coincide with a distribution of the deformation in the crust in front of the oceanic propagator as in the Atlantic realm (Enachescu, 2006;Mohn et al., 2015;Nirrengarten et al., 2018;Sibuet et al., 2007), in the Gulf of Aden (Fournier et al., 2010;Leroy et al., 2012), or in the northern Red Sea (Almalki et al., 2016;Bartov et al., 1980;Bosworth, 2015Bosworth, , 2020, where major transform fault zones are segmenting the oceanic basins at locations where stalling phases occurred. Previous studies proposed that these stalling phases can be related to several processes acting as an opposing force to the rift propagation in the undeformed continental lithosphere like the presence of high-strength continental domains or inherited structures in the lithosphere (Courtillot, 1982;Morgan & Parmentier, 1985;Nirrengarten et al., 2018). As a result, they are generally interpreted as indicating varying kinematic conditions or inheritance. This hypothesis was already questioned by a numerical model of rift linkage at an oblique continental margin, which showed that with increasing strain localization the tectonic regime changes from initially dip-slip dominated to strike-slip dominated (Le Pourhiet et al., 2017). Here we show that strain partitioning and change in deformation regime can occur at the tip of oceanic propagators under constant boundary kinematic conditions (Figures 9 and 11). This implies that a relationship seems to exist between major transform margins and the rift stalling phase during ocean opening. In the North Atlantic Ocean, we observe that the Newfoundland-Azores-Gibraltar fault zone and the Charlie-Gibbs fault are the two major transform faults that separate segments of different spreading ages (Nirrengarten et al., 2018, and references therein). In the same way, the Gulf of Aden (Figure 1a) opened almost synchronously between the Alula-Fartak fault zone and the Shukra el Sheik fault zone but was stalling along the Shukra el Sheik fault zone for 9 to 12 Ma before propagation restarted (Fournier et al., 2010;Leroy et al., 2012). However, in our models, we were not able to produce oceanic transform faults segmenting the ridge. Thermomechanical models tend to show that magmatism and serpentinization are essential to produce oceanic transform faults (e.g., Ammann et al., 2017;Gerya, 2012).
The South China Sea and the Gulf of Aden are both good examples to illustrate the rifting of a weak continental lithosphere (Figures 4b, 4d, 7, and 9) due to the geodynamic context in which those two oceanic basins have formed. The Gulf of Aden is located near an important mantle thermal anomaly (the Afar hot spot; e.g., Fournier et al., 2010;Leroy et al., 2012), while the South China Sea opening followed a long-term subduction associated with magmatism (Pubellier et al., 2003;Wang et al., 2013). The Gulf of Aden presents the advantage of being an active plate boundary monitored with GPS stations. The GPS velocities with respect to Arabia fixed (ArRajehi et al., 2010) display the eastward motion (i.e., parallel to the ridge) of the continent at the tip of the oceanic propagator and an oblique northeastward motion for Africa fixed. The oceanic basin and continental margins display stalling phases and propagation phases, and the rifting phase is characterized as a wide oblique rift. Figure 13c shows the relative velocity field in model 1 (1.2 mm year −1 rift-parallel inflow), which displays the oblique velocity at the tip of the propagator. In the Gulf of Aden, the force that triggers the rift-parallel inflow as displayed by GPS velocities could come from the differential gravitational potential energy between the high topography above the Afar hot spot and the opening Gulf of Aden.
However, the strain localization and the evolution of the deformation in a stronger lithosphere (Figures 4e,  4g, 4h, and 11) results in the development of a localized strike-slip fault and the formation of a microplate. These models display similarities with the northern Red Sea and its two branches, the Gulf of Aqaba to the east and the Gulf of Suez to the west (Almalki et al., 2016;Bartov et al., 1980;Bosworth, 2015Bosworth, , 2020. In the northern Red Sea, the Sinai block behaves as a third tectonic microplate separated from Arabia by the Levant transform fault and from Africa by an active oblique rift (with respect to the Red Sea). Although a continental rift is not a complete plate boundary, it is the initiation. The GPS velocities with respect to Arabia display the southward motion of the Sinai block along the Levant Fault. Therefore, regarding the absence of shortening structures in the northern Red Sea, the GPS velocities showing the rift-parallel displacement, and the presence of the transform continental fault of the Levant, we propose that a rotational extension model (Molnar et al., 2017;Smith, 1993;Zwaan et al., 2020) may not be the best candidate to explain the dynamics of the northern Red Sea. However, the formation of the triple junction between a continental transform fault, an active oblique rift, and a ridge is well reproduced in our strong lower crust models when the rift-parallel inflow is greater than 16% of extension velocity (Figure 14b, model 3). The geodynamic origin of this regional rift-parallel inflow arises from the contrast between the fast motion of the Arabian plate driven by subduction in Zagros and Makran and the slow northward motion of Africa in collision with the European plate.
Moreover, the lithosphere-scale strain partitioning between a highly localized strike-slip fault and a very oblique rift leads to the development of two types of transform margins with extremely different kinematic, dynamics, and thermal history.
The classical model for transform margin formation involves the development of transform faults during the intracontinental rifting phase (Basile, 2015;Francheteau & Le Pichon, 1972;Mascle & Blarez, 1987). The evolution from rifting to oceanic spreading should preserve the continental structures and develop oceanic transform faults along the previous continental ones. This model has major implications on the development of the margin structure and temperature. It implies that oceanic ridge migrates along the transform margin as a triple junction at half the spreading rate (for large offset transform margin, it may take several million years), heating the continental margin and its basins (Rüpke et al., 2010). It also implies that the strike-slip faults responsible for the steep structure of the transform margins (Mercier de Lépinay et al., 2016) are active at the beginning of the rifting phase and that the transform margin was born with characteristics observable today on the fossil transform margins. However, models of this study show different characteristics and evolution for transform margins.
On the one hand, models with high obliquity display the typical kinematic context of a transform margin (Basile, 2015;Francheteau & Le Pichon, 1972). In this condition, one of the oblique continental rift 10.1029/2020JB019906 Journal of Geophysical Research: Solid Earth branches ( Figure 11, cross sections 1 and 2, and Figure 12) reaches continental breakup, and the system evolves into two en échelon ridges separated by a zone of soft linkage that evolves into an oblique rift with strike-slip deformation similar to those in previous studies (Ammann et al., 2017;Laetitia Le Pourhiet et al., 2017;Liao & Gerya, 2015). However, structural properties of these modeled margins correspond to hyperthinned continental crust but not steep margins along a localized transform fault. On the other hand, simulations with high obliquity and a strong lower crust produce steep margins along a pure strike-slip fault (e.g., Figure 11f, cross section 4) that display the structural characteristics of a transform margin as a slope break, a marginal ridge close to the slope break, and a strike-slip fault separating the continental crust from the exhumed mantle (or oceanic crust) (Mercier de Lépinay et al., 2016). However, these strike-slip margins form in the atypical kinematic context for a transform margin. From an initially hyperthinned margin developed during the initial extension phase (until 30 Myr), a steep margin develops (Figures 11e and 11f, cross section 4) as the deformation increases along the strike-slip fault. This strike-slip shear zone puts in contact a 50 km thick continental lithosphere with the asthenosphere. As a result, the formation of these steep margins starts 10 Myr after the continental breakup, meaning that the margin geometry is acquired after the intracontinental rifting phase. These results highlight the finding that oblique and transform margins can acquire their final structure several million years after continental breakup due to the strike-slip deformation.

Conclusion
Three-dimensional thermomechanical numerical models integrating a rift-parallel inflow show the following: 1. The 3D boundary conditions can be interpreted as an absolute velocity field that can be converted to a relative velocity field comparable with natural cases of triple junctions involving continental transform faults. 2. At least three propagation modes can exist. A fast propagation resulting in a straight oceanic basin (two plates), an alternation of propagation and stalling phases (extension opposing velocity of 12% to 15% for the weak lower crust) and propagation with a "jump" of the spreading center (opposing extension velocity over 17% for the strong and weak lower crusts). 3. Continental breakup propagation slows down when rift-parallel inflow is applied (Le Pourhiet et al., 2018), and we show that during these stalling phases the continental margins and tectonic structures rotate and the deformation regime changes from dip-slip dominated to strike-slip dominated, while boundary kinematic conditions are kept constant. 4. There might be a relationship between stalling phases, transform margin formation, and the subsequent location of large oceanic transform fault zones due to margin rotation. 5. Two modes of transform margin can develop: steep margin associated with strike slip faults and hyperextended ones associated with an offset spreading segment similar to that in Le Pourhiet et al. (2017) 6. All oblique margins might not emerge from a single geodynamic process, and it could be worthwhile in the future to revisit the available data using the soft linkage oblique margin and the strike-slip oblique margin models as a classification.

Data Availability Statement
pTatin3D is an open-source code, which is available publicly at https://bitbucket.org/ptatin/ptatin3d/src/ master/. The version of the code used in this study, an example options file to reproduce the results, and the models results are available publicly at http://doi.org/10.5281/zenodo.3967184. Maps from Figure 1 were made with GeomapApp and post-processing with Paraview.