Long‐Term Coupling and Feedback Between Tectonics and Surface Processes During Non‐Volcanic Rifted Margin Formation

Here we present high‐resolution 2‐D coupled tectonic‐surface processes modeling of extensional basin formation. We focus on understanding feedbacks between erosion and deposition and tectonics during rift and passive margin formation. We test the combined effects of crustal rheology and varying surface process efficiency on structural style of rift and passive margin formation. The forward models presented here allow to identify the following four feedback relations between surface processes and tectonic deformation during rifted margin formation. (1) Erosion and deposition promote strain localization and enhance large offset asymmetric normal fault growth. (2) Progressive infill from proximal to more distal half grabens promotes the formation of synthetic sets of basinward dipping normal faults for intermediate crustal strength cases. (3) Sediment loading on top of undeformed crustal rafts in weak crust cases enhances middle and lower crustal flow resulting in sag basin subsidence. (4) Interaction of high sediment supply to the distal margin in very weak crust cases results in detachment‐based rollover sedimentary basins. Our models further show that erosion efficiency and drainage area provide a first‐order control on sediment supply during rifting where rift‐related topography is relatively quickly eroded. Long‐term sustained sediment supply to the rift basins requires elevated onshore drainage basins. We discuss similar variations in structural style observed in natural systems and compare them with the feedbacks identified here.


Introduction
The morphology and structure of rifted passive margins are characterized by a high variability in width, fault geometry, and nature and geometry of the sedimentary succession (Clerc et al., 2017). Studies of active rifts and extensional domains show that their internal structure, history, and dimensions are highly variable (Ruppel, 1995). Contrasting narrow, wide, and core complex modes of extension are predicted to depend on crustal thickness, geotherm, and strain rate (Buck, 1991). Forward numerical models demonstrate the control of rheological stratification, with two end-members types for varying crustal strength: strong crust resulting in narrow rifts with few high angle normal faults (Type 1) and weak crust resulting in wide rifted margins with significant middle/lower crustal flow and core complex style deformation (Type 2; e.g., Huismans & Beaumont, 2011, 2014Svartman Dias et al., 2015). Many aspects of rifts system remain unexplained owing to the complexity of the processes that interact with one another. Weakening processes and structural and rheological inheritance of the lithosphere play an important role in controlling the structural style of deformation (e.g., Brune et al., 2014;Duretz et al., 2016;Huismans & Beaumont, 2002Lavier & Manatschal, 2006;Salazar-Mora et al., 2018;Wu et al., 2015). In this study we focus on the role of surface processes during nonvolcanic rifted passive margin formation.
A range of studies have explored the coupling and feedback between surface processes and fault activity. On short seismic time scales, erosional unloading may contribute to Coulomb stress loading along thrust fault ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2018JB017235
Key Points: • Surface processes enhance plastic strain localization • Surface processes provide important feedback on rifted margin structure • Normal faults distribution, margins width, time of crustal breakup, and stratigraphic records are controlled by these feedbacks Supporting Information: • Supporting Information S1 • Movie S1 • Movie S10 • Movie S11 • Movie S12 • Movie S13 • Movie S14 • Movie S15 • Movie S16 planes leading to increased fault activity . On long-term tectonic time scales, displacement on a single normal fault is enhanced by footwall erosion and hanging wall deposition (Maniatis et al., 2009). Olive et al. (2014) argue that the life span of normal faults is increased in response to surface processes resulting in large fault offsets. Several studies (Buiter et al., 2008;Choi et al., 2013;Olive et al., 2014) indicate that syntectonic erosion and deposition control location and geometry of secondary faults and delay fault migration. Strak et al. (2011) using analog models show that throw rate controls fluvial incision. Similarly, several recent studies have investigated the sensitivity of landscape evolution to normal fault activity (e.g., Demoulin et al., 2018;Petit et al., 2009;Roda-Boluda & Whittaker, 2018).
Several previous studies explore the interaction and feedback of surface processes and extensional tectonics at the scale of the whole rifted margin system. Bialas and Buck (2009) show using forward numerical modeling with simple sediment aggradation that sedimentation may promote narrow rifting. In contrast, Burov & Poliakov (2001 using coupled forward models show that increased surface processes efficiency may enhance distributed deformation (higher crustal thinning and wider margins), while they suggested that erosion of rift flanks in combination with high sedimentation rate in the rift basin induces outward lower crustal flow toward the rift flanks. Andrés-Martínez et al. (2019) show using 2-D coupled forward models with a surface diffusion model for erosion and sedimentation that surface processes enhance strain localization and control margin width and symmetry. High sedimentation rate favors abrupt crustal necking and narrow proximal margins. Efficient mass transport to the distal margin by subaerial processes (low sea level) leads to more symmetric and narrower conjugate margins.
Here we use a high-resolution 2-D thermo-mechanical numerical model (Erdős et al., 2014;Thieulot, 2011) coupled with a mass balancing fluvial erosion, marine deltaic deposition model to investigate the effect of varying rheological conditions, surface process efficiency, and sea level on the structural style of rift and passive margin formation. We first present the numerical modeling approach followed by the model results. We discuss primary feedback relations from the model results and compare these with observations from natural rift systems.

Modeling Approach
We use an arbitrary Lagrangian-Eulerian finite element method for the solution of thermo-mechanical coupled, plane strain, incompressible viscous-plastic creeping flows (Thieulot, 2011;Erdős et al., 2014). The finite element model solves the force balance equations of equilibrium for quasi-static incompressible flows in two dimensions: where μ eff is the effective viscosity, υ is the velocity, ρ the density, and g the gravity acceleration. In equation (1), mass conservation equation with incompressibility constraint (∇υ þ p λ ¼ 0 where p is the pressure and λ a penalty parameter with same dimension as a bulk viscosity) and momentum conservation equation (Stokes equation) are combined (Thieulot, 2011). In addition to the mechanical system, we also solve the thermal evolution in two dimensions using the finite element method: Here T is the temperature, c p is the specific heat capacity, k the thermal conductivity, and H the radiogenic heat production. The last term in the heat balance equation is the temperature correction for adiabatic heating and cooling when material moves vertically at velocity v z , where α T represents the volumetric thermal expansion coefficient.
When the state of stress is below the frictional-plastic yield stress, the flow is viscous and is specified by temperature-dependent nonlinear power law rheologies based on laboratory measurements on "wet" quartz (Gleason & Tullis, 1995) and "wet" olivine (Karato & Wu, 1993). The effective viscosity, μ eff , in the power law rheology is of the general form: where Ė 2 is the second invariant of the deviatoric strain rate tensor 1 2 _ ε ij _ ε ij , A is the preexponential scaling factor, n is the power law exponent, Q is the activation energy, R is the universal gas constant, and V is the activation volume. The factor f is used to scale viscosities calculated from the reference quartz and olivine flow laws, thereby producing strong and weak versions of these materials. Frictional-plastic yielding is modeled with a pressure-dependent Drucker-Prager yield criterion which is equivalent to the Coulomb yield criterion for incompressible deformation in plane strain. Yielding occurs when where J 2 is the second invariant of the deviatoric stress 1 2 τ ij τ ij , lτ ij is the deviatoric stress tensor, ϕ eff is the effective internal angle of friction, psin(ϕ eff )=(p−p f )sin(ϕ), p f is the pore-fluid pressure, and C is the cohesion. With appropriate choices of C and ϕ eff , this yield criterion can approximate frictional sliding in rocks and the effect of pore-fluid pressures. Plastic flow is incompressible. Strain softening is introduced by a linear decrease of ϕ eff (ϵ) from 15°to 2°and C(ϵ) from 20 to 4 MPa with respect to plastic strain (ϵ;  & Beaumont, 2002;Lavier et al., 1999Lavier et al., , 2000; Figure 1). The plastic strain is updated at every time step with the second invariant of the deviatoric strain. The incompressible plastic flow becomes equivalent to a viscous flow if an effective viscosity, μ p eff , for a plastic material is defined such that . Setting the viscosity to μ p eff in regions that are on frictional-plastic yield satisfies the yield condition and allows the velocity field to be determined from the finite element solution of equation (1). The overall nonlinear solution is determined iteratively using μ ¼ μ p eff (for regions of plastic flow) and μ ¼ μ v eff (for regions of viscous flow; Willett, 1992).
The mechanical and thermal evolution are coupled through the temperature dependence of viscosity and density and are solved alternately during each time step. Densities of crust and mantle depend on temperature with ρ(T)=ρ 0 (1−α(T−T 0 )), where α is the coefficient of thermal expansion and ρ 0 is the density at T=T 0 . Parameters are given in Table 1 and in the Supporting Information S1

Model Setup and Boundary Conditions
The models are set up as an idealized rheologically layered lithosphere above a sublithospheric mantle in a 600-km-high and 1,200-km-wide model domain ( Figure 1). The lithosphere consists of a 35-km-thick crust and 85-km mantle lithosphere overlying the sublithospheric mantle. The top 3 km of the crust represents prerift sediments with the same rheology as the upper crust. The Eulerian grid consists of 2,400 and 290 elements in the horizontal and vertical directions, respectively. The distribution of the elements is irregular in the vertical direction, allowing for high resolution in the upper crust with a vertical resolution of Δz=200 m in the first 20 km, Δz=625 m between 20 and 70 km, Δz=1,100 m between 70 and 120 km, and Δz=7,917 m between 120 and 600 km. The resolution in the horizontal direction is 500 m for the entire model domain.
Extensional horizontal velocity boundary conditions (v = ±0.5 cm/year) are applied to the lithosphere, and the corresponding exit flux is balanced by a low velocity inflow in the sublithospheric mantle. The top of the model is a free surface. The sides and base are vertical and horizontal free slip boundaries, respectively ( Figure 1a).
Thermal boundary conditions are specified basal temperature, 1,520°C, and insulated lateral boundaries. The initial temperature field is laterally uniform and increases with depth from the surface (T 0 =0°C) to the base of the crust (initial Moho temperature, Tm=550°C) with a surface heat flux of 55.3 mW/m 2 . Below the Moho, temperature linearly decreases to the base of mantle lithosphere (initially at T = 1,328°C); the temperature of the sublithospheric mantle follows an adiabatic gradient (0.4°C/km). Thermal conductivity linearly increases to 51.46 W·m −1 ·K −1 at 1,350°C (∼125-km depth), corresponding to scaling the thermal conductivity by the Nusselt number of upper mantle convection. The enhanced conductivity maintains a nearly constant heat flux to the lithosphere base and an adiabatic temperature gradient in the sublithospheric mantle (e.g., Pysklywec & Beaumont, 2004).

Rheological Setup
The crust follows a wet quartz rheology (Gleason & Tullis, 1995) with different scaling factors as a way to test the effect of varying crustal strength ( Figure 1). We test four contrasting crustal rheologies, varying scaling factor f uc from 30, 1, 0.1, and 0.02, resulting in thickness of the frictional-plastic upper crust that ranges from 25 km to about 8 km. The viscosity scaling ( Figure 1) represents a simple technique that creates either strong crust or weak viscous crust without recourse to additional flow laws, each with its own uncertainties. The scaling can either be interpreted as a measure of the uncertainty in the flow properties of rocks where flow is dominated by quartz (e.g., wet or dry) or taken to represent variations in thermal state and composition (e.g., Huismans & Beaumont, 2011, 2014. We investigate systems that are characterized by mechanical heterogeneity represented by white noise in the initial strain field (Figures 1 and S1). This approach is designed to represent inheritance of deformation from previous tectonic phases. Specifically, the plastic strain (ϵ) is initialized with white noise that has a Gaussian distribution with a mean value set to 0.3 and with a maximum value 0.8. Inherited weakness is provided by a tapered symmetrical 400-km-wide area in the model center. A small thermal heterogeneity is introduced at the base of the lithosphere in the model center in order to enhance rift localization (see Supporting Information S1 for a detailed description). The orientation of faults and thus the asymmetry of the developing rift is not prescribed and is generated spontaneously.
In order to ensure mass conservation, the average pressure at the bottom part of the model is maintained constant by adjusting the influx of the sublithospheric mantle at the sides. This allows defining an absolute sea level within the model independent from surface displacement. Based on sea level, a water load is implemented in order to fully consider all mass loads on the free surface.

Surface Processes Modeling
Surface processes in the forward model act on the free upper surface ( Figure 2). We implement mass balancing erosion and deposition. A river profile is applied to topography above sea level between each local minimum and maximum similarly to previous studies (e.g., Steer et al., 2011;Willett, 1999Willett, , 2010. Transported sediment can be deposited along a river profile and in lakes above sea level, by deltaic progradation below sea level, or leave the model domain when it reaches the sides. Instantaneous mechanical compaction due to vertical loading of sediments is included with an associated volume and density change (see Supporting Information S1 for detailed description; Albertz & Ings, 2012;Athy, 1930;Chamot-Rooke et al., 1999;Cowie & Karner, 1990;Goteti et al., 2012;Tenzer & Gladkikh, 2014). The change in free-surface topography is a result of the net effect of tectonic uplift or subsidence, erosion, and deposition: where h is the elevation, ė is the erosion rate, _ d is the deposition rate, and _ T the tectonic uplift rate. All river profiles are first computed, and transported sediments are stored at each local minimum or at the "coast line" defined by the sea level. Sediments are then aggraded into the lakes or bypassed down to sea or outside the model domain if they reach sides when lakes are full. Sediments are finally prograded offshore below sea level.
Fluvial erosion is implemented by solving the shear-stress fluvial incision or stream power law (Howard & Kerby, 1983;Howard, 1994;Lague, 2014;Tucker & Whipple, 2002;Whipple & Tucker, 1999;Whipple, Hancock, & Anderson, 2000;Whipple et al., 2000: where ė is the erosion rate, K is a constant called erodibility or erosion efficiency factor, Q w is the stream discharge, and s the topographic slope. The stream discharge Q w represents the total precipitation accumulated on the upstream part of the watershed, Q w =A·P, where P is the precipitation rate and A is upstream drainage area given by A ¼ Li 1:4 À Á 1 0:55 where L i is the length of the drainage area or the river in our case given in kilometers (Hack, 1957;Rigon et al., 1996;Willett, 2010). A sediment transport length function ξ(q) controls sediments deposition along the fluvial profile (Davy & Lague, 2009;Kooi & Beaumont, 1994;Yuan et al., 2019) where q is the stream discharge per unit of river width. The rate of fluvial deposition is controlled by where q s is the sediment river load per unit of river width W, αq corresponds to the sediment transport length, where α (s·m −1 ) >0 allows choosing between detachment-limited (only erosion) and transport-limited (deposition is allowed) end-members. The sediment river load Q s is given at each node by the amount of upstream accumulated sediments and the volume of sediments removed or deposited at this node: where L is the catchment length and the river width W varies as the square root of the water discharge c· ffiffiffiffiffiffi ffi Q w p (Lacey, 1930;Savenije, 2003). There is a river between each local maximum and local minimum above sea level. Each river profile, given by −ė þ _ d river , is computed using a 1-D implicit algorithm (Braun & Willett, 2013). A constant erodibility is used for the entire model. It is assumed that the actual drainage divide in the topography is located between two nodes. The drainage divide is located between the highest node and its lowest neighbor. Consequently, if one side of the drainage divide erodes faster than the other side, the drainage divide can migrate.
In our model, deposition along the river profile only happens when the net change in elevation is positive. Most of the deposition occurs either in local minima above sea level by aggradation in lakes or when sediments are transported beneath sea level by progradation. Each local minimum above sea level is defined as a lake. The elevation of the sill defines the maximum level of sedimentation in the current lake. When a local minimum is full of sediments, then they are bypassed to the next local minimum down to the sea or leave the model domain when they reach sides. Deltaic progradation is based on a characteristic sediment transport length L d , intended to simulate the fall out of a suspension plume (equation (9)).
where q s is the sediment discharge per unit width. To ensure mass balance, q s at the river mouth (river that ends into the sea or into a lake above sea level) corresponds to the integrated net surface change by erosion and deposition along the river profile: All surface processes parameters are summarized in Table 2, and additional detailed explanation of the surface processes model can be found in Supporting Information S1.

Reference Models Without Surface Processes
We first describe the Reference Models M1 to M4 with varying crustal rheologies (R1-R4; with f uc = 30, 1, 0.1, and 0.02) without surface processes and analyze the relief produced for these four cases (Table 3). The end-member without surface process is an important model reference that allows to separate and understand the effect of surface processes on model behavior.

Model M1, Strong Crust, f c =30
Reference Model M1 (Figure 3a) has the strongest crustal rheology with minor decoupling in the viscous lower crust. The tectonic style of deformation is asymmetric and evolves in two phases. Phase 1 Note. Stream power law (see equation (6)), river deposition (see equations (7) and (8)), and hillslope processes and compaction (see Supporting Information S1).
deformation is controlled by one single frictional-plastic shear zone cross cutting the strong crust and upper mantle-lithosphere, while Phase 2 is characterized by ductile necking and symmetric break up. At 5 Ma, most of the deformation is accommodated by a 45°dipping frictional-plastic shear zone that extends into the upper mantle at 55-km depth. Minor secondary localized shear zones develop in the hanging wall of the primary asymmetric normal fault zone. At 10 Ma, continued offset on the large asymmetric shear zone results in rupture of the hanging wall leading to several smaller-scale upper crustal fault blocks on top of the upper mantle-lithosphere bounded by a regional low-angle detachment. Thinning of the mantle-lithosphere and upwelling of the asthenosphere lead to Phase 2 symmetric ductile necking localizing deformation in a narrow region and lithospheric rupture at about 15 Ma. The final margin geometry is highly asymmetric with clear upper and lower plate styles originating from Phase 1. The "lower plate" margin is characterized by a single large offset normal shear zone, several upper crustal allochtonous blocks in contact with the upper mantle-lithosphere, and absent middle to lower crust ( Figure 3). The "upper plate" margin in contrast is characterized by distributed small offset fault zones and a narrow tapered crustal thinning and relief gently dipping toward the distal margin. Model M3 has a weak crustal rheology (Figures 1 and 3). It exhibits a thin about 11-km frictional-plastic upper crust which is fully decoupled from the frictional-plastic upper mantle-lithosphere. This model Note. K2 to 10 −9 , K2-3 to 2.5.10 −9 , and K3 to 5.10 −9 , and sl0 is equivalent to 0 m, sl1 to −500 m, sl2 to −1,500 m.

Reference Models With Surface Processes
A relatively limited range of erodibility (1 order of magnitude in our model; see Tables 2 and 3) results in end-member behavior with either limited or highly efficient (close to complete) erosion. We next present Rift Models M5-M8 with Rheologies R1-R4 including surface processes with intermediate fluvial erosion efficiency (K = K2 = 1×10 −9 m 1−3m s m−1 ) and high sea level fixed at −500 m below the initial reference topography at z=0 m (Table 3).   Model M7 with a weak crust can be compared with the equivalent M3 without surface processes (Figures 1  and 5). Models with weak crust result in very moderate to absent rift flank topography (Figure 4). Given a relatively high sea level, these models have limited potential to produce sediment discharge to the extending margin, which explains the sediment starved character of Model M7. Model M7 exhibits, similar to M3, two characteristic phases of deformation with Phase 1 distributed crustal extension concomitant with narrow mantle-lithospheric necking and Phase 2 progress localization of crustal thinning in the distal margin.

Model Sensitivity to High Erodibility
Models M9 to M12 explore sensitivity to increased erodibility for systems with strong to very weak crustal rheology and high sea level (Figure 6a). For each model, the geometry is shown close to lithosphere breakup (see Movies S1-S16 for model animation). All models show significantly higher sediment input into the basin area throughout their evolution. For the strong crust Case M9, rift flank topography is efficiently eroded with high early sediment export into the basin resulting in a crustal-scale asymmetric half graben.
Most of the strain localizes on a single border fault, and minor conjugate shear zones accommodate hanging wall deformation. Following lithospheric breakup, continued sediment discharge onto the margin creates thick growth fault bordering prograding wedges on the asymmetric conjugate margins. The intermediate strength crust Model M10 shows a similar evolution while evolving in a symmetric rift style. Efficient erosion of the rift flanks on both sides results in high sediment discharge and enhanced offset of the symmetric graben early in its evolution. During and following breakup, continued deposition leads to normal fault bounded growth sequences on new "oceanic" sublithospheric mantle, in this case resulting in more symmetric conjugate margins. Weak crust Model M11 evolves similarly to Model M6 that was characterized by a lower erodibility. A first stage of early distributed crustal extension is followed by progressive necking in the distal margin. Low rift flank topography limits early sediment export to the basin. During its later stage, the large onshore drainage area provides significant higher sediment input into the distal basin area. Late stage prograding sediments bypass early syn-rift deposit and enhance a distal margin crustal-scale shear zone rooting in the midcrust, thus providing significant accommodation space. For very weak crust Model M12, the increased erosion rates result in early sediment input to the distal margin. Deposition in the distal margin interacts with the formation of a large offset low-angle normal fault exhuming midcrust at the base of the sediments. Late in its rift evolution progressive thinning of the distal margin over a wide area results in necking of this sedimentary wedge and formation of a late stage core complex. The distal margin Figure 6. Final architecture for models with surface processes considering (a) high erosion efficiency K3 and a sea level "sl1" at −500 m (Models M9-M12 with Rheologies R1-R4; Table 3) and (b) low sea level "sl2" fixed at −1,500 m and an intermediate erosion efficiency K2 (Models M13-M16 with Rheologies R1-R4; see Table 3). See sections 4.5 and 4.6 for detailed description and Tables 2 and 3 for model parameters.

10.1029/2018JB017235
Journal of Geophysical Research: Solid Earth exhibits also broad platform-type deposition on crustal rafts with very minor upper crustal extension and subsidence resulting from middle and lower crustal removal similar to Model M7.

Model Sensitivity to Low Sea Level
Models M13 to M16 test sensitivity to low sea level at −1,500 m (Figure 6b), with the reference intermediate erodibility and varying crustal rheology (Table 3). Lowering sea level has two main effects.
(1) It increases the amount of onshore upper crust that can potentially be eroded and produce sediment export to the rifted margin.
(2) It reduces at the same time the accommodation space available in the rifted conjugate margin resulting in lower sediment thicknesses in the proximal basin areas and promoting sediment bypass to the distal margin. Increased sediment export can be observed in all four models with low sea level. Interaction between deposition and large offset low-angle normal fault zones in these models leads to large tilted sedimentary basins in the distal margin.

Discussion
We first interpret the first-order factors controlling the structural style of rift and passive margin formation in the models. We then identify and discuss the feedback relations between surface processes and tectonics. Model predictions are compared with observations from natural systems, and finally, limitations of the modeling approach are discussed.

Factors Controlling the Structural Style of Passive Margin Formation
Lithospheric rheological stratification provides the main control on the structural style of rift and passive margin formation. The structural style of passive margins is characterized by the nature and distribution of normal faulting and by its morphology (margin width, onshore topography, and top basement geometry offshore). Models characterized by limited to absent decoupling between the strong frictional-plastic crust and upper mantle layers promote narrow rifts, with most extension accommodated on few large offset listric normal faults, similar to Type 1A margins as described by Beaumont (2011, 2014) and Buck (1991). In contrast, weak crust cases with efficient decoupling between a thin frictional-plastic upper crustal layer and the strong mantle-lithosphere below result in highly depth-dependent extension with crustal extension distributed over a wide area and concomitant narrow rupturing of the lower lithosphere, as in Type 2 margins (e.g., Huismans & Beaumont, 2011, 2014Svartman Dias et al., 2015).
Narrow Type 1 margins as represented by Models 1 and 2 and variations with surface processes are characterized by high angle normal faults and a high crustal strength that allows for high and narrow rift flank topography onshore. In contrast, wide Type 2 margins as in Models 3 and 4, the weak middle and lower crust allows for (1) distributed crustal deformation and fault migration, (2) low or absent rift flank topography, (3) large offset normal faults with hanging wall and footwall rotation and exhumation of midcrust to the surface following the rolling-hinge model (e.g., Lavier et al., 1999), (4) middle and lower crustal flows toward the distal margin, and (5) formation of crustal raft with midlower crust removal and limited upper crustal extension resulting in sag-type subsidence.

Interaction and Feedback Between Surface Processes and Rifted Margin Formation
The forward models presented here allow to identify the following four feedback relations between surface processes and tectonic deformation during rifted margin formation ( These feedbacks have consequences for fault offset, margin width, and time of crustal breakup (Figures 8 and  9). The structural evolution cannot be easily compared between the models as the initial structural inheritance and the nonlinearity of each numerical simulation with different conditions of surface processes introduce variations to timing of initiation and abandonment of individual faults. However, for models with a strong crust, the proximal normal fault is a systematic feature that we use to compare fault offset with changing erodibility. The proximal normal fault in Model 1 without surface processes exhibits 32 km of offset (Figures 8). With increasing surface process, efficiency fault offset increases progressively to 65 km for complete footwall erosion and hanging wall deposition (Figure 8). To first order proximal and distal margins width is controlled by crustal rheology (Figure 9a). Surface processes provide opposite effects on margin width for strong versus weak rheologies. When the crust is strong (R1) or has an intermediate strength (R2), proximal margin width decreases with erodibility owing to efficient strain localization, whereas distal margin width increases with increasing erodibility resulting from efficient hanging wall block rotation along a detachment fault coupled to exhumed mantle. In contrast, when the crust is weak (Rheologies R3 and R4), lower crustal flow toward the rift center due to mass redistribution by surface processes allows for wider proximal and distal margin widths with increasing sediment input. The timing of crustal breakup is relatively insensitive to erodibility for strong and intermediate crustal rheologies (R1 and R2), whereas for weak crustal rheology (R3), the time of crustal breakup increases linearly with increasing erodibility (Figure 9b).
These feedbacks also have consequences for onshore erosion and topography and for offshore sediment supply and stratigraphy. The evolution of the topography onshore controls sediment supply through time. High sediment supply during rifting is enhanced by high onshore paleotopography, large source drainage area, and high erosion rate. In our models, rift-related topography is relatively quickly eroded ( Figure 10). Long-term sustained sediment supply to the rift basins requires elevated onshore drainage basins. Rift flank erosion leads to migration and retreat of the drainage divide and the formation of a rift escarpment. Drainage capture triggers a transient phase of fluvial nick point migration followed by formation of a new drainage divide resulting from regional isostatic uplift. Finally, we distinguish three characteristic stratigraphic architectures: 1. Fault-bounded symmetric or asymmetric basins with growth strata (Figures 7a and 7b), 2. Sag basin deposition on top of crustal rafts with midlower crustal thinning (Figure 7c), and 3. Rollover sedimentary basins with large sections of tilted deposits on top of a flat horizontal detachment fault (Figure 7d).

Comparison With Earlier Studies
The variation of rift width with changing crustal rheology is well established and documented by a range of modeling studies (e.g., Huismans & Beaumont, 2002, 2011, 2014Brune et al.,   (2019), our models exhibit symmetric margins and lack lateral rift migration and large-scale margin asymmetry. The strong localization in these studies results from the combined effects of frictional-plastic strain weakening (included in our models) and viscous strain weakening (not included in our models) that cause the asymmetry (lateral rift migration) as also shown by Huismans and Beaumont (2003).
Several earlier modeling studies explored the interaction and feedback between rifting and surface processes. The feedback of footwall erosion and hanging wall deposition in enhancing fault offset in half grabens described hereby was similarly identified by Olive et al. (2014) and Andrés-Martínez et al. (2019) and suggested by Maniatis et al. (2009). As also shown by these studies, the reduction of the topographic load and flexural force through footwall erosion increase fault life span. The effect of erosion on fault offset is especially effective in the proximal domain of rifted margins where normal faults result in topography above sea level. There appears to be an optimal range of brittle layer thickness (10-20 km) where the feedback of erosion is most pronounced in our models. Lower brittle layer thickness leads to large offset no matter the erosion efficiency, whereas greater brittle layer thickness leads to reduced sensitivity to footwall erosion as also shown by Olive et al. (2014). In the distal margin for strong crust cases with a thick upper crustal brittle layer, normal fault offset is mainly controlled by hanging wall block rotation, and deposition may enhance fault offset.
In contrast to Poliakov (2001, 2003), we do not observe outward midlower crustal flow toward the rift flanks and related stabilization of rift shoulder uplift with high sedimentation rate. Midlower crustal flow in our models occurs generally in the direction of a negative pressure gradient, that is from the elevated rift margins to the subsiding basin center, provided that the crust is sufficiently weak. With weak crustal rheologies, high sedimentation rate enhances midlower crustal flow toward the rift center leading to wide margins and delaying distal margin subsidence.
In contrast with Bialas and Buck (2009), we observe that distal margin width increases with high sediment deposition whatever the crustal rheology. The time of crustal breakup is relatively insensitive to sediment deposition for strong and intermediate crustal rheologies and increases linearly for cases with a weak crust.
Only the proximal domains of cases with strong and intermediate crustal rheology are narrower with efficient surface processes. The effect of surface processes on fault migration and interaction is important, e.g. feedback 2, as demonstrated by the highly different fault evolution with varying surface processes  efficiency. This is consistent with feedback relations inferred by previous modeling studies (Olive et al., 2014;Buiter et al., 2008;Choi et al., 2013).

Natural Cases
We next discuss briefly some natural cases that may provide examples of the inferred four feedback relations between surface processes and tectonics during extensional basins formation.
(2) Many rifts and passive margins as for example the Northern North Sea (Bell et al., 2014;Cowie et al., 2005) and the Gulf of Suez (Gawthorpe et al., 2003) exhibit sets of synthetic basinward normal faults (Figure 11b). We suggest, based on the understanding of the feedback relations developed here, that synthetic basinward faulting may be enhanced by syn-rift deposition.
(3) Sag deposition along the central South Atlantic margin on top of undeformed basement has been explained as resulting from middle and lower crustal thinning (Karner et al., 2003;Moulin et al., 2005;Figure 11c). Middle to lower crustal flow in wide Type 2 margins as identified here may explain the inferred crustal thinning. (4) Tilted sedimentary sections on top of a low-angle detachment as observed in Devonian basins onshore Norway (e.g., Fossen et al., 2016;Seranne & Seguret, 1987) and in the distal part of the wide Gabon margin (Clerc et al., 2017;Figure 11d) are consistent with rollover basins as predicted here in case of high sediment supply in weak Type 2 extensional systems.

Model Limitations
Our modeling approach is inherently 2-D and does not take into account any 3-D variation of tectonic deformation and landscape response. Drainage divide migration (Goren et al., 2014;Braun, 2018) and average topography by river incision (Willett, 2010) could also be modeled with analytical solutions. However, sensitivity tests using a 2-D tectonic model coupled to a 2-D plan form surface process model indicate very similar behavior (Beucher & Huismans, 2016). The model uses a viscoplastic rheology and does not include elasticity. The viscoplastic behavior does, however, exhibit regional and local flexural uplift in response to tectonic and crustal loading, and the model is therefore appropriate for exploring the long-term interactions between surface processes and tectonics. However, the flexural wavelength in our model is rate dependent, while it is rate independent with elasticity. The effect of loading and unloading may therefore generate some differences in term of flexural response and influence fault offset and patterns of fault migration (Olive et al., 2016).
The surface process model does not include short-term variations in precipitation, the effect of storms, and orographic effects, while potentially important are beyond the scope of the present study. Other effects such as sediment blanketing that may affect the thermal structure during rifting and will potentially influence tectonic deformation are not included here and deserve further study.

Conclusion
We use state of the art 2-D forward dynamic models to understand the interaction and feedback between tectonics and surface processes during lithosphere extension and rifted margin formation. The arbitrary Langrangian Eulerian high-resolution 2-D model solves for viscoplastic deformation coupled to a mass balancing fluvial erosion-submarine delta deposition model that allows to bridge a large range of scales from lithosphere-upper mantle scale to the scale of the sedimentary basin. We evaluate competing controls on the structural style of rifted margin formation, associated onshore tectonic morphology, and offshore sedimentary basin architecture and test the sensitivity to crustal strength, fluvial erosion efficiency, and "sea level." Based on our forward models, we conclude the following: 1. The primary control on the structural style of rift and passive margin formation is provided by crustal rheology which is consistent with earlier studies. Type 1 narrow continental margins with strong crust and high coupling with mantle-lithosphere exhibit few large offset listric faults that branch into the mantle. Type 2 wide continental margins with weak crust and high decoupling between upper crust and mantle-lithosphere exhibit distributed deformation and few large offset detachment faults. 2. Footwall unloading by erosion and hanging wall loading by deposition provide a first-order feedback on tectonic deformation by enhancing asymmetric normal fault growth and prolonging fault activity. 3. For intermediate crustal strength, footwall unloading by erosion and hanging wall loading by deposition together with sediment bypass promote the formation of sets of basinward dipping synthetic faults in the proximal margin. 4. Sediment loading on top of undeformed crustal rafts in weak Type 2 margins enhances middle and lower crustal flow resulting in sag basin subsidence. 5. High sediment supply to the distal margin in very weak crust Type 2 systems results in detachment-based rollover sedimentary basins. 6. Erosion efficiency and drainage area provide the first-order control on sediment supply during rifting.