Morphotectonic Evolution of Passive Margins Undergoing Active Surface Processes: Large‐Scale Experiments Using Numerical Models

Extension of the continental lithosphere can lead to the formation of rifted margins with contrasting tectonic and geomorphologic characteristics. Many of these characteristics depend on the manner extension spatially distributed. Here we investigate the feedback between tectonics and the transfer of material at the surface resulting from erosion, transport, and sedimentation and discuss how they influence the rifting process. We use large-scale (1, 200 × 600 km), high-resolution (1 km) numerical experiments coupling a 2-D upper-mantle-scale thermo-mechanical model with a planform 2-D surface processes model. We test the sensitivity of the coupled models to varying crust-lithospheric rheology and erosional efficiency. We confirm that the development and long-term support of topography is dependent on the strength of the coupling between the crust and the mantle lithosphere. Strong coupling promotes high topography as the integrated strength of the lithosphere is sufficient to support the additional stress. Weak coupling results in the stress being relaxed via viscous flow in the middle/lower crust and leads to more subdued topography. Erosion and transport of sediment modulates this behavior but has only minor effect on the overall structure of the rift. High erosion efficiency counters the development of high topography and creates complex landscape morphologies while low erosion efficiency allows for longer standing high topography and results in more simple landscape morphologies. The transfer of mass between the continent and the basin alters the stress field at the onshore-offshore transition and facilitates the development of faults, increasing their offsets and keeping them active over a longer period.


Introduction
The extension and subsequent breakup of the continental lithosphere is responsible for the formation of a wide range of rifted passive margins, with contrasting tectonic and geomorphological characteristics, ranging from the high-elevation cases (south-eastern South Africa, Red Sea rift) characterized by sharp sub-aerial topography and prominent escarpments, to margins with low elevations (southern Australia, eastern Africa) with more subdued topography (e.g., Gallagher & Brown, 1997;Gilchrist & Summerfield, 1990) (Figure 1). A variety of mechanisms have been proposed to explain the existence of high-elevation topography onshore, including depth-dependent extension (e.g., Huismans & Beaumont, 2011, lateral heat flow (Steckler, 1981), induced asthenospheric convection (e.g., Huismans et al., 2001), plastic necking and flexural response to unloading (e.g., Braun & Beaumont, 1989), or magmatic underplating (e.g., McKenzie, 1984). These mechanisms can be invoked to explain a wide range of topographic wavelengths depending on the strength of the middle and lower crust which control the coupling between upper and lower parts of the lithosphere.
From the onset of rifting, the topography created is subject to erosion, with varying degree of efficiency depending on the regional climatic context and the nature of the eroded material. Evidence of sediment, several kilometers thick, infilling passive margin sedimentary basins (e.g., Whittaker et al., 2013) highlight a strong generic link between onshore and offshore domains ( Figure 1) and point to the importance of onshore topography as the source of material transferred into the basins. It also raises important questions on how the redistribution of mass by surface processes may affect the evolution of rifting.  (>500 m). Elevation data are from the ETOPO1 data set (https://www.ngdc.noaa.gov) while offshore sediment thicknesses are from Whittaker et al. (2013). Four examples of passive margins morphologies with associated cross sections: (b) Area off Morocco with section modified after Osmundsen and Redfield (2011). The breakup is dated around 175 to 190 Ma. Morphology shows relatively subdued topography with low elevations and a smooth transition to the oceanic domain. (c) Area off Namibia with section modified after Osmundsen and Redfield (2011). Although the age of breakup is relatively old (approximately 127 Ma), the morphology shows a sharp escarpment with high elevation and a sharp transition to the oceanic domain. (d) Area from the Red Sea rift in Arabia showing first-order morphology with a sharp escarpment and high elevations. (e) Area of the Ethiopian margin across the Afar depression. The margin is at a syn-rift stage and presents high elevations.

Geochemistry, Geophysics, Geosystems
10.1029/2019GC008884 coupled response to tectonics, climate, and surface processes (Whittaker, 2012). In some cases, redistribution of mass by surface processes are sufficient to alter the internal dynamics of orogens (Avouac & Burov, 1996;Erdős et al., 2014;Stolar et al., 2006;Steer et al., 2014). Recent numerical modeling works have investigated the large-scale interactions between tectonics and surface processes during rifting and have highlighted the importance of the crustal strength in modulating topography and the subsequent mass redistribution (Andrés-Martínez et al., 2019;Huismans, & Beaumont,2005;Theunissen & Huismans, 2019). Rivers have long been recognized as the main agents of long-range transport of sediment between continents and oceans (e.g., Whipple, 2009). They are also active agents of erosion with river incision primarily relating to the shear stress exerted by the flowing water and erodibility of the river bed (e.g., Howard, 1994). The lower limit for river erosion is ultimately controlled by the sea level (although local base levels can also be established away from the ocean where rivers flow into large rivers or lakes). Supply limited rivers may lead to rapid erosion and thus rapid rock uplift in local areas (e.g., Zeitler et al., 2001) while at a larger scale the transfer of mass modifies the balance between gravitational and tectonic forces. In rift systems, the sediment deposits offshore extend over hundreds of kilometers and commonly reach thickness of several kilometers (Figure 1) which can translate into additional vertical stresses of tens of mega-pascals, comparable to tectonic forces. Moreover, erosion and sedimentation can increase the throw rate (Maniatis et al., 2009) and the maximum throw of faults (Olive et al., 2014) and can favor narrow over distributed deformation within the crust (Bialas & Buck, 2009;Buiter et al., 2008). Climate also plays an important role as it controls the amount of water available and modulates the discharge in the rivers. It also strongly influences chemical alteration of rocks, development of vegetation, soil evolution, and so forth. Topography resulting from tectonic uplift can alter the wind and atmospheric circulation patterns that in turn can affect the pattern of precipitation. Change in precipitation patterns may in turn cause heterogeneous distribution of erosion with the potential for changing the stress distribution and thus the structural evolution of a system (Beaumont et al., 1992;Stolar et al., 2006;Willett, 1999;Whipple, 2009).
Common approaches for understanding topography and basin evolution include mass balance studies between source areas and sedimentary basins (e.g., Anell et al., 2010;Guillocheau et al., 2012), geomorphological studies of the source areas (e.g., Hall & Bishop, 2002), cosmogenic nuclides (e.g., Cockburn et al., 2000), and thermochronological studies (e.g., Cowie et al., 2000;Sacek et al., 2012;Wildman et al., 2016). These techniques can record different processes (surface uplift, rock uplift, and rock exhumation) over a range of different timescales, and resolving the apparent conflicts between methods can be challenging. Understanding the coupled tectonic-morphological response of passive margin development requires controlled experiments only possible through modeling of the entire onshore-offshore system. In the following, we investigate the feedback between tectonics and the transfer of material at the surface in the context of passive margin formation using 2-D thermo-mechanical models coupled to a state-of-the-art 2-D surface processes model. We explore how mass redistribution during continental rifting and passive margin formation affects the structural style and conversely, how different modes of rifting and their associated characteristic topographies lead to a different surface process response.

Tectonics
We use the geodynamics code Sopale (Fullsack, 1995) which solves for 2-D thermo-mechanically coupled viscous-plastic plane-strain deformation of an incompressible continuum. The equations solved are where is the stress tensor, P is pressure, is density, g is acceleration owing to gravity, x i are the spatial coordinates, T is temperature, C p is heat capacity, k is thermal conductivity, t is time, v is the velocity vector, and H is the radiogenic heat production per mass unit. The first two equations solve for the scalar pressure field (P) and the velocity vector field (v) over the spatial coordinates (x i ) with being the density of the Radiogenic heat production A R 0.9 × 10 −6 W∕m 3 materials, g the gravity which allows for buoyancy forces, and the stress tensor. The last equation is the heat equation: T is temperature, C p is heat capacity, k is thermal conductivity, t is time, and H is the radiogenic heat production. The last term accounts for adiabatic heating or cooling when material is advected vertically at velocity v 2 . The thermo-mechanical coupling is carried out using temperature-dependent density and viscosity. The relation between density and temperature is calculated as follows: is the coefficient of thermal expansion and is set to 2.0 × 10 −5 • C −1 . Table 1 reports the densities used for the materials at T 0 = 0 • C.
The model is composed of materials (crust, mantle lithosphere, and sublithospheric mantle). When the state of stress is below the frictional-plastic yield criterion, the flow is viscous and follows a strain-rate, pressure and temperature-dependent nonlinear rheology with the effective viscosity eff specified as where A is the pre-exponential scaling factor, n is the stress exponent,̇I ∕ 2 is the second invariant of the deviatoric strain rate tensor , Q is activation energy, V * is activation volume, and R is the universal gas constant. A, n, and Q derive from laboratory experiments and vary according to mineralogy and water content (Wet Quartz, Olivine) (Karato & Wu, 1993). The parameter f is a scaling factor chosen to represent viscous weakening or strengthening with respect to the reference flow law. This scaling is a simple way to address the large uncertainty in the effective viscosity of earth materials (e.g., Huismans & Beaumont, 2003).
Brittle (frictional-plastic) failure is pressure-dependent according to a Drucker-Prager yield criterion which in 2-D is equivalent to the Coulomb yield criterion: where is the second invariant of the deviatoric stress, / ij is the deviatoric stress tensor, eff is the internal angle of friction given as P sin( e ) = (P − P ) sin( ) for pore-fluid pressure P f , and C is the cohesion. A linear decrease of the effective internal angle of friction from 15 • to 2 • for accumulated strain values between 0.5 and 1.5 is used to account for strain softening in frictional-plastic shear zones (Huismans & Beaumont, 2003).

Surface Processes 2.2.1. Fluvial Erosion and Deposition
The top free surface of our 2-D thermo-mechanical model is subjected to the planform Landscape Evolution Model FastScape (Braun & Willett, 2013) which includes fluvial erosion, fluvial deposition, and hillslope processes. As for most landscape evolution models, the erosion rateĖ is assumed to be dependent on the local slope S and the upstream drainage area A. The model solves the classical stream power law equation in two dimensions:Ė where K f is an erosion efficiency factor which depends on lithology, climate, and local hydrology. The values of the parameters m and n are still largely debated and vary depending on the climatic and tectonic context (see discussion in Croissant & Braun, 2014). Estimates of the ratio m∕n are generally determined from the gradient of slope versus area plots and the assumption of a steady state profile and spatially uniform uplift. m∕n has been shown to vary between 0.35 and 0.6 for regions where these assumptions are correct (e.g., Whipple, 2009). It is reasonable to take a value for n = 1 as the basal shear stress exerted on the river bed varies linearly with channel slope. We use an intermediate value for m = 0.4 and n = 1 (see Table 1).
The second term of the equation incorporates the effects of sediment on incision (Davy & Lague, 2009). It allows for deposition in the channels at a rate which is taken to be proportional to the sediment discharge Q s but inversely to water discharge Q w .
Hillslope processes are modeled using linear diffusion by assuming that transport is linearly proportional to slope: where H∕ t represents the change in elevation H over time t and K d is the coefficient of diffusive transport. Contrary to nature, there is no realization of actual channels in the numerical models. All points of the model free surface are subjected to fluvial erosion and hillslope processes. Numerical fluvial erosion is effectively the dominant process, and hillslope processes could be discarded. However, a small amount of diffusion is beneficial to numerical stability.

Offshore Sedimentation
FastScape has been modified to include offshore sedimentation. The sediment fluxes at the river mouths are converted to sediment volumes and transferred into the basin depending on distance and amount of available accommodation space. For the offshore transport of sediment, we used the same approach as in Sacek et al. (2012) which is similar to hillslope diffusion, but with a depth-dependent coefficient, K m , to account for the decrease in energy available for transporting sediment with increasing water depth: where H∕ t represents the change in elevation H over time t and where K * m is a constant, H sea is sea level, and H c is a characteristic depth.

Coupling FastScape With Sopale
We have developed a module that allows coupling between the free surface of the thermo-mechanical model (S) and the 2-D FastScape surface (F). The coupled models represent a T-model (or 2.5-D model). For every tectonic time step dt S , the S surface (effectively a 1-D line) is tectonically advected (both vertically and horizontally) by the velocity field V which derives from the mechanical solution ( Figure 2). The new positions of the points define a new surface that can be approximated using a cubic spline. The spline is then used to derive the uplift at the node locations of the higher-resolution Fastscape surface. This interpolation step is necessary as FastScape uses a static computational grid and does not allow for direct horizontal tectonic advection. As Sopale returns a 2-D velocity field (V (x,z) ), the surface is only tectonically advected in the x and z directions, and no tectonic advection occurs in the y direction (V (y) = 0.). FastScape calculates erosional and depositional fluxes and updates the 2-D surface elevations. The new surface elevations are averaged along a 1-D line which is used to update the tectonic model free surface (S) and the material field (essentially adding sediment where needed) before a new tectonic solve is initiated. The evolution of the free surface is thus determined from both the erosional/depositional fluxes and the tectonic advection. The numerical stability of the landscape evolution model typically relies on shorter time steps than the tectonic model (≈100 to 500 yr). Therefore, the Sopale tectonic time step dt S is divided into sub-time steps dt F . The Fastscape surface is conserved in memory at every time steps to ensure the conservation of the geomorphic features.

Model Design
Our coupling approach is similar to the one recently used by Theunissen and Huismans (2019). The main difference resides in the use of a 2-D plan landscape evolution model which explicitly uses the stream-power law and allow for predicting river networks and landscape morphologies. A 1-D fluvial erosion model (e.g., Andrés-Martínez et al., 2019;Theunissen & Huismans, 2019;Willett, 2010), although sufficient to estimate net amount of erosion and sedimentation, does not allow for tracking the position of the rivers heads nor the location of escarpments and cannot be used to assess rates of rift flank retreat. Furthermore, 2-D plan models are more suitable for studying transient evolution of the topography as they do not assume that the landscape has reached a steady state. On the other hand, tectonics models are only sensitive to the mean topography (e.g., Willett, 2010) and in so, 1-D fluvial erosion models remain a valid option for studying the effect of coupling with erosion/deposition on a large scale. It is however worth noting that as 3-D thermo-mechanical models become more ubiquitous, thanks to the increasing availability of fast and reliable solvers, the community faces an increasing need for robust and reliable workflow coupling state-of-the-art surface processes models and 3-D thermo-mechanical models; we believe this work is part of such effort.

Thermo-mechanical Model Setup
The thermo-mechanical model ( Figure 2) extends down to 600 km depth over a 1, 200 km wide region. The lithosphere includes a 35 km thick crust and a 90 km thick mantle lithosphere. The crust has a wet quartz rheology (WQz) (Gleason & Tullis, 1995) and a reference density of 2, 800 kg∕m 3 while the mantle lithosphere and the sublithospheric mantle have a wet olivine rheology (WOl) (Karato & Wu, 1993) and a reference density of 3, 300 kg∕m 3 . A scaling factor f (equation (5)) is used to define a moderately weak (f = 1) and a strong (f = 30) wet quartz-based crustal rheology. For a nominal strain rate and an initial geotherm, changing f from f = 1 to f = 30 increases the viscous strength of the crust and the depth of the brittle/viscous transition with a resulting thicker frictional layer and a thinner viscous layer. As the model evolves this translates into more or less efficient coupling between the crust and the mantle lithosphere. The wet olivine rheology is scaled to differentiate the sublithospheric mantle (f = 1) from the mantle lithosphere (f = 5) to account for a mechanically stronger mantle lithosphere (Burov & Watts, 2006) (Figure 2 and Table 1). A rectangular (12 km wide by 10 km high) strain-weakened seed (WOl rheology, = 2 o ) is placed at the center of the model, in the frictional-plastic uppermost part of the mantle lithosphere, in order to initiate localization of deformation. The model does not account for elastic behavior, and there is no short-term elastic flexure of the lithosphere. However, the viscous-plastic component of the rheologies and the layering of the lithosphere allows for a non-local, viscous-plastic flexural isostatic response to the changing loads.
The Eulerian grid consists of 600 elements in the horizontal and 200 elements in the vertical dimension. The horizontal resolution is 2 km for the entire domain. The distribution of the elements is irregular in the vertical direction, and the resolution varies with depth, with 100 elements in the first 25 km (250 m resolution), 50 elements between 25 and 125 km depth (2 km resolution), and 50 elements between 125 km and 600 km (5 km resolution). Each cell of the Eulerian grid is initialized with nine Lagrangian (material) particles that allow for the tracking of the materials through space and time.
Lithospheric divergence is achieved at the lateral boundaries by symmetric outward-directed horizontal velocities, which sum to give a full divergence rate of V = 1.5 cm/yr. The rift forms through thinning and necking of the lithospheric layers and results in vertical advection of the underlying asthenospheric mantle. There is no "active" forcing from the asthenosphere. A material inflow is prescribed in the sublithospheric mantle with a velocity scaled to maintain the models' isostatic balance. The vertical sides are free to move in the vertical direction. Finally, a free slip condition is applied at the base of the models.
The temperature field is initialized to a laterally uniform steady state with temperatures at the Moho and the base of the lithosphere initially at 550 and 1330 • C, respectively. No heat flux is allowed across the vertical side boundaries. Internal heat production is uniformly distributed in the crust and absent in the mantle. Note that we do not account for the potentially lower diffusivity of newly deposited sediment but use the same diffusivity used for the crust . All thermal parameters are listed in Table 1.

Landscape Evolution Model
The landscape evolution model surface dimension is 1, 200 km × 200 km with a uniform resolution of 500 m (2,400 × 400 cells) (Table 1). In all the following experiments, the side boundary conditions of the landscape evolution model are set to periodic for the north/top and south/bottom sides (water and sediments leaving the domain re-enter the model from the opposite side), and the west/left and east/right sides are set to be at base level. The sea level is initialized to a value H 0 = −300 m and remains constant through the duration of all experiments. The initial surface elevation (relative to sea level) is thus at 300 m which is consistent with the mean elevation of continents. Sediments deposited under water are diffused using the zero-depth diffusion coefficient K * m = 100 m 2 ∕yr and a characteristic depth of 1, 000 m. We vary the erosion efficiency by changing the fluvial erosion constant K within 1 order of magnitude (10 −6 to, 10 −5 m 1−3m s m−1 ) between models.

Results
We present a series of model snapshots at t ≈ 4, 7, and 10 Myr for both studied cases. Figures 3 and 4 show the evolution of the topography and the structural evolution of the rift while Figures 5 and 6 illustrate the distribution and amount of erosion/sedimentation together with the stress field (second invariant of the stress tensor) within the lithosphere. We chose to use the stress field as it provides an overall view of the strength of the lithosphere (a similar figure showing the strain rate is available as additional material).
We define a set of metrics to characterize the morphological and structural evolution of the models (Figure 7). We focus on the evolution of the onshore domain, and therefore, only elevations above sea level are used in the calculations. (Z mean ) and (Z max ) record the mean and maximum elevations with respect to sea level. The mean slope (S mean ) corresponds to the mean of the absolute value of the topographic slope over the subaerial part of the model domain. The topographic relief (R max ) is measured across the model as the difference between the minimum and maximum elevations along-strike the rift (i.e., between rivers and hilltops). Finally, the mean rock uplift (U) and mean erosion rates (E) illustrate the relative control and dynamics of tectonics and surface processes through time.

Model 1: Strong Crust Rheology
Model 1 (Figures 3 and 5) is similar to the "Type I" margins of Beaumont (2011, 2014) in which the crust has a wet quartz viscous rheology (Gleason & Tullis, 1995) and scaling factor f = 30 is efficiently coupled to the mantle lithosphere ( Figure 2).
Deformation starts with distributed pure shear before localizing in the central part of the model. The deformation occurs along symmetric conjugate shear zones that develops in the brittle part of the strong crust while underneath the ductile mantle lithosphere is rapidly necking (Figure 3). A small asymmetry can be observed as strain-weakening results in one of the shear zone being activated before the other (see , for discussion). At 7.0 Myr the deformation has thinned the lithosphere to approximately 20 km. Lithospheric breakup occurs at 9 Myr.
The resulting continental margins show thinning from 35 to 0 km across a narrow region which is less than 100 km wide while necking of the lithospheric mantle occurs over a region that is twice as wide. The shear zones accommodate the rotation of tilted blocks in the basin. The topography at the surface shows elevated rift flanks with elevation over 5, 000 m ( Figure 3). A sharp transition from onshore to offshore occurs along the main shear zones.
Model 1A surface evolution (Figure 3a) shows development of asymmetric rift flanks which define a clear water divide: Steep rivers drain small watershed areas delimited by short wavelength interfluves into the basin while, on the opposite sides, rivers show more gradual slopes and drain larger areas toward hinterland domains. The source of sediment arriving into the basin is limited to the small watersheds flowing into the basin where the maximum of erosion is recorded (Figure 5a). Slopes on the basin side of the flanks control the escarpments retreat, and the diminution in sediment production is only partially compensated by an increasing drainage area. However, as erosion is low, the morphology remains controlled by tectonics, and the escarpments conserve their asymmetry through time. The shorelines remain close to the base of the escarpments and define a sharp transition between sub-aerial and subaqueous topography. Transport of sediment from sources to sinks occurs on a short distance as most sediment is deposited in the proximal part of the basin (Figure 5a). Model 1B (Figure 3b) shows the effect of increased erosion on the morphology. The asymmetry of the rift flanks is less pronounced, and the slopes on both sides of the rift flanks reduce as fast as they are eroded (c, f) Evolution of the distance between the top of the escarpment/rift shoulder and the shore line. We only consider the left part of the rift and extract the location of the maximum elevation for elevation greater than 300 m (we consider that there is no rift shoulder otherwise). Position of the shoreline is estimated by extracting the location where the mean topography intersect with sea level.
( Figure 8). The headwaters of rivers flowing into the basin interweave with those belonging to the rivers flowing toward hinterland areas. The water divide wanders in a broad region and moves as the headwaters of rivers flowing to the basins penetrate further inland. The resulting morphology is more complex than in the "Low Erosion" case: Watersheds are more elongated and wider with interfluves showing greater wavelengths. The area affected by erosion is much wider (approximately 100 km at 7.0 Myr, Figure 5b). Distance between escarpments and shorelines increases as sediment progrades ( Figure 8) and large flat areas make a smooth transition from sub-aerial to subaqueous areas. Previously submerged areas are uplifted above sea level and are then eroded. There are thus large areas previously dominated by sedimentation that are reworked by erosion (see Figure 5 at 7.0 Myr). Sediment arriving into the basin fills up the accommodation space by prograding from the rivers' mouths toward the distal domain. Lateral distribution of material is important, and sedimentation tends to fill up the accommodation space before reaching more distal blocks in the basin.
The large amount of sediment deposited on the offshore part of the margin is accompanied by an increase in offset (≈+50%) along the shear zones (Figure 3b, at 10.0 Myr). The deformation tends to localize more, and the thinning of the lithosphere is delayed in the distal part of the rift (Figure 3b, at 7.0 Myr). The effect of the mass transfer from the onshore to the offshore is seen on the stress field: The high stress values, at the base of crust, below the rift shoulders and the main shear zones, are shifted toward the center of the basin. As the blocks rotate more than for Model 1A, basement highs appear in the basin. They are subjected to erosion and become a secondary source of sediment (Figure 3b, at 10.0 Myr).

Morphology
Evolution through time of the distance between the top of the rift flank and the center of the model (Figures 8a and 8b) shows horizontal advection of Model 1A rift flanks at a rate of 0.8 cm∕yr which corresponds to the half-extension rate setup for the model. The positions of the rift shoulders remain stable through time. The rift flanks in Model 1B migrate at a much faster rate (1.4 cm∕yr), illustrating the effect of increased erosion on the escarpment retreat rate (0.7 cm∕yr) (Figures 8a and 8b). Similarly, the positions of the shoreline relative to the top of the escarpments remain constant during the evolution of Model 1A (Figure 8c), while in the case of Model 1B those distances increase significantly at an approximate rate of 1.2 cm∕yr as a consequence of the basin being filled with sediment (Figure 8c).

10.1029/2019GC008884
Model 1A mean (Z mean ) and maximum (Z max ) elevations (Figures 7a and 7b) increase significantly around 2 Myr in response to a significant increase of the uplift rate (U) to a maximum value of ≈ 0.8 mm∕yr (Figure 7e). The low erosion efficiency and high uplift rates allow for rapid increase of the mean slope (S mean ) to a value of ≈ 6 • at 5 Myr (Figure 7c) while the maximum relief (R max ) quickly reaches ≈ 4 km (Figure 7d). The mean slope (S mean ) only shows a relatively small increase after 4 Myr while R max rapidly decreases to 1 km and then remains constant (Figures 7c and 7d). Erosion (E) is essentially equal to 0 after 5 Myr. The uplift rate (U) shows on average negative values (Figures 7e and 7f) consistent with both the mean (Z mean ) and maximum elevations (Z max ) that continue to decrease (Figures 7a and 7b).
Model 1B mean (Z mean ) and maximum (Z max ) elevations (Figures 7a and 7b) show the same response to high uplift rates (Figure 7e) as for Model 1A. However, the high erosion efficiency limits the mean elevation (Z mean ) to ≈ 1.3 km, 0.3 to 0.4 km lower than for Model 1A. The maximum elevation (Z max ) also reaches a slightly (0.3 to 0.4 km) lower value. The combination of high erosion rate (E) and high uplift (U) between 2 and 4 Myr results in an increase of the mean slope (S mean ) which reaches a maximum value of ≈ 11 • at 5 Myr before rapidly decreasing. The maximum rate of uplift (U) occurs at the same time as for Model 1A (2 Myr). The maximum relief (R max ) shows a similar evolution, with values higher than Model 1A and a maximum of ≈ 4.5 km. The maximum relief (R max ) decreases much slower than for Model 1A after ≈ 4 Myr, as a consequence of erosion being more efficient. Slopes, relief, and uplift rates result in higher erosion rates (E ≫ 0.1 mm∕yr), active over a longer period (2 − 8 Myr).

Model 2: Weak Crust Rheology
Model 2 (Figures 3 and 5) has a wet quartz viscous rheology, 30 times weaker than Model 1 (WQz f = 1). The thermal properties of the model and the initial temperature field are identical to Model 1. The model is characterized by widespread viscous horizontal decoupling between the crust and the mantle lithosphere ( Figure 2). The geometry of the margin contrasts with that of Model 1 in that the crust extends over a wider area into the continent. The upper part of the crust and the mantle lithosphere have a plastic rheology. The deformation localizes along two frictional-plastic shear zones as in Model 1. The shear zones rapidly flatten into the viscous part of the crust forming low angle detachments, where they thin and neck the lower crust resulting in uplift of the footwall areas and slow subsidence of the central keystone area. With time deformation migrates to the center of the model resulting in secondary frictional shears rupturing the keystone block concomitant with necking of mantle lithosphere. The crust ruptures shortly after the viscous mantle lithosphere layer. The final geometry of the rift shows wider passive margins and lower topography in comparison with Model 1. The basin depth is significantly less than in the strong rheology case.

Effect of Erosion and Deposition
Model 2A (Figures 4 and 6a) is the "Low Erosion" case. The rift flanks are less elevated than for Model 1, and their morphology does not show significant asymmetry between the sea and hinterland side. Early in the model evolution the slopes are sufficient to define a clear water divide with rivers flowing toward the basin or hinterland domains. The water divide then becomes harder to localize as topography decreases. There is no real morphological difference between rivers flowing to the basin or to the hinterland. The rift flanks are rapidly eroded to sea level. The amount of sediment that reaches the basin is limited by the low relief onshore and is strongly controlled by the sea level at the start of the experiment.
The "High Erosion" case (Model 2B, Figure 4b) does not show any significant difference with Model 2A in terms of tectonics, and the stress field within the lithosphere is identical. As erosion is more efficient, the landscape is quickly eroded away, and the surface rapidly flattens and reaches sea level. The amount of sediment is again limited by the low elevated topography undergoing erosion. The total amount is similar to that of Model 2A, but sediment reaches the offshore domain earlier in the rift evolution, before mechanical and thermal subsidence had significantly deepened the basin. As a result the accommodation space is less and sediment progrades further into the distal domain (Figure 6b). Rift flanks are much less elevated and the position of the drainage divide and the shoreline become poorly defined after approximately 5 Myr (Figures 4b and 8) as the topography approaches sea level. The high erosion rate does not significantly affect the position of the highest elevation before 4 Myr.

Morphology
Model 2A mean elevation (Z mean ) reaches a maximum value of 0.5 km at ≈ 3 Myr and then slowly decreases before stabilizing around ≈ 0.3 km at 6 Myr ( Figure 7a). The maximum elevation (Z max ) shows the same evolution with a maximum value of ≈ 0.8 km at ≈ 2 Myr. The mean slope (S mean ) presents typical values below 1 • over the entire duration of the model. A maximum relief (R max ) value of ≈ 0.5 km is recorded at In comparison to Model 2A, Model 2B has higher maximum elevations (Z max up to ≈ 1 km) (Figure 7b). It also shows steeper slopes (S mean ), up to 0.3 • at ≈ 3 Myr, which translate into higher reliefs (R max ) with amplitudes approximately twice those of Model 2A (up to ≈ 1 km, Figures 7c and 7d). After ≈ 4 Myr, the uplift rate is insufficient to counter the ≈ 0.1 mm∕yr erosion rate which results in decreasing mean elevation (Z mean ) and relief (R max ).

Discussion
The rheology of the crust determines the strength of the coupling with the strong upper mantle and dictates the amplitude and steepness of the topography and thus the potential for erosion and sediment transport. Erodibility and amount of precipitation determine the timing and efficiency of the mass transfer between onshore and offshore domains. From our results, it is clear that tectonics controls the overall geometry of the rift and that erosion only weakly modulates its evolution.

Origin of the Topography
The formation of rift flanks is a consequence of regional isostatic compensation of mass that is redistributed during the necking of the lithosphere as previously proposed by Braun and Beaumont (1989), Fletcher and Hallet (1983), and Zuber and Parmentier (1986). Rift flanks uplift occurs via flexure of lithosphere in response to mechanical or thermal unloading of the lithosphere during extension. The amplitude of the resulting topography depends on the flexural strength of the lithosphere which is affected by extension. It is then expected that a strong lithosphere would have a strong flexural response and high rift flank topography, while a weaker lithosphere would lead to reduced flank topography.
As the flexural strength of the lithosphere reflects its integrated brittle, ductile, and elastic strengths (e.g., Watts & Burov, 2003), the ability to sustain topographic heights and steep slopes depends on the strength of the coupling between the brittle crust and the strong upper mantle (e.g., Burov & Watts, 2006). In the case of Model 1, efficient coupling between the strong upper crust and upper mantle lithosphere results in efficient stress transfer and significant uplift of the rift flanks. As the integrated strength of the lithosphere is high, the flexural response of the lithosphere leads to the formation of prominent flank topography with a clear asymmetry between steep slopes dipping toward the basin and gentler slope dipping toward hinterland domains. Rift flank uplift occurs while the two major boundary faults are slipping which suggests that the flexure of the lithosphere is a direct consequence of the mechanical unloading of the lithosphere along the faults. In the Model 2 case, the weak middle to lower crust relaxes the stresses, reducing the flexural strength of the lithosphere and in so limits uplifts of the rift flanks. The role played by the thermal weakening of the lithosphere due to upwelling of hot mantle lithosphere is difficult to quantify but is strongly dependent on the rate of extension.
The strength envelope of the lithosphere is the main control on the amplitude of the rift flank uplift. A strong crust coupled to a strong upper mantle lithosphere is expected to promote high topography while any decoupling between the upper crust and upper mantle lithosphere owing to the presence of a weak layer will relax the stresses and will result in more subdued topography. This highlights the importance of the regional thermal state and preexisting crustal composition at the time of rifting.

Topography and Erosion
Rift flank topography controls the drainage pattern on a regional scale and exerts a first-order control on the sediment fill of the basin. A water divide rapidly separates rivers flowing toward the basin from streams draining the hinterland domain. Erosion efficiency, controlled by climate and function of lithology, modulates the response to tectonic deformation. If erosion is not efficient enough to balance the uplift, the water divide remains stable through time, the sizes of the catchments remain constant, and erosion is limited. A natural example is found in the lake Baikal rift system where thermochonological analysis shows that only approximately 100 to 200 m of erosional denudation has affected the flanks adjacent to the central rift segment since the onset of rifting (van der Beek, 1997). The water divide is clearly established at the edge of the rift flanks, controlling the shape of the drainage area. On the other hand, high erosion efficiency allows rivers to rapidly cut through the uplifting rift flanks, and the morphology is then determined by surface processes. Topography is the direct expression of the competing effect of erosion and rock uplift. The models presented here show that for most of the time, topography is not at steady state (U∕E ≠ 1) and is either increasing, due to tectonic (U), or lowering in response to erosion (E). Models with low erosion efficiency (Models 1A and 2A) converge to steady state as shown by evolution of elevation, slope, and relief which reach stable values after 5-6 Myr. In the cases with high erosion efficiency (Models 2A and 2B), elevation, slope, and relief converge toward zero and will eventually achieve a steady state as topography reaches sea level. It is, however, important to point out that in this case steady state implies close to zero rock uplift and zero erosion, which for passive margins corresponds to the period of post-rift tectonic quiescence. 5.2.1. Sediment Flux and Onshore to Offshore Transition Sediment production is essentially limited by the ability of the rivers to erode the landscape created by tectonic uplift. The amount of sediment that reaches the basin also depends on the size of the drainage area that connects to the basin itself. The size of the catchments flowing toward the sea increases through time as the escarpments defined by the rift flanks migrate inland. The basin is thus continuously fed by erosion of the escarpment and subsequent transport by the rivers. The strength of the lithosphere controls both the development of the rift flanks and also affects the subsidence and the amount of accommodation space in the basins. High rift flanks may provide large amounts of material through erosion, but they also create natural topographic barriers for large river systems draining the continent. However, the link between high topography and high sediment thicknesses is not always evident in nature: Large rivers such as the Amazon, the Mississippi, or the Ganges have brought enormous amounts of sediment to the margins from distances greater than 1, 000 km (e.g., Hori & Saito, 2007).

Surface Processes and Tectonic Feedback
The tectonic evolution and the general geometry of the system is conserved as we go from no surface processes to increased erosion and sedimentation. Tectonics is the main factor controlling the style of deformation. Transfer of mass from the onshore to the basin changes the lithostatic pressure gradient within the crust and results in local isostatic adjustments. The main feedback of surface processes on tectonic deformation is that loading by sediment results in increased offsets and rotations of the rigid hanging walls along the main fault/shear zones. The sediment load locally increases the vertical stress which facilitates yielding of the crust underneath. Accommodation space increases as the hanging wall rotates and subsides. The accommodation space is progressively filled by sediment. Overall, the steepness of the fault/shear zone increases with increasing erosion and sediment loading owing to the combined effect of footwall uplift and hanging wall subsidence and rotation. This is in agreement with Olive et al. (2014) who showed that some of the large offsets observed on rift faults have been facilitated by the unloading and loading of the normal fault footwalls and hanging walls through erosion and sedimentation. How sediment volume is accommodated depends again on the rheology of the lithosphere. A strong lithosphere tends to localize deformation and is more likely to accommodate sedimentation within some relatively local depressions. This is where the additional loads have the potential to affect tectonics by keeping faults active, increasing their offsets and delaying the formation of more distal faults (Olive et al., 2014). In the case of a weak lithosphere, extension is much more distributed, the depth of the basin tends to be lower, and there is less accommodation space available. If the energy available for transport is sufficient, we expect sediment to be distributed over longer distance. In that case the load is too thinly spread to have a real impact on the structure of the margin.

Comparison With Natural Cases
Our predicted morphologies are consistent with observations of real passive margins as shown by Petit et al. (2007) for the Dhofar margin where the flexural rigidity and climatic-driven erosion efficiencies are likely responsible for the contrasting morphologies flanking the Gulf of Aden. In a context where the monsoon has established a strong climatic contrast, areas with high precipitations (high erosion efficiency) show morphologies with low relief, flat low-lying sediment plains, no significant retreat of the escarpment, and low flexural response (Petit et al., 2007). In contrast, semi-arid areas with fewer precipitations (low erosion efficiency) have morphologies revealing a high component of flexural rebound, significant escarpment retreat (25-30 km), and bedrock exposure at the foot of the escarpment (Petit et al., 2007). They concluded that although the lithospheric elastic thickness is responsible for the large-scale flexural response and large-scale retreat of the escarpment, erosion is the main factor controlling the landscape morphology and that transport of eroded material to the foot of the escarpment can potentially affect its evolution (Petit et al., 2007).

Model Limitations
A systematic comparison of the passive margins morphologies worldwide is beyond the scope of this study, and the limited set of models presented here is not meant to capture the morphological variability of passive margins. Such work would require constraining numerical models with careful observations of the morphologies onshore, estimates of spatial and temporal distribution of erosion, and broader availability of the sediment records offshore. This is arguably challenging as the present topography is often the only criteria available to modelers. However, the perspectives of such a work could offer interesting insights into the evolution of passive margins under different climatic/erosional contexts.
The models presented have tested the sensitivity of landscape morphology related to rifting dynamics by varying the strength of the lithosphere and the erosion efficiency K. The limitations of our models reside both in the approach used to model the surface processes and in the tectonic model driving the boundary conditions. The model configuration used in this study does not capture the complexity of natural systems and must be seen as deliberately simple experiments used to investigate the interaction between surface processes and tectonics. Here we discuss, in a qualitative manner, the factors that affect the setup of our models.
The T-Model approach used in our study is one of the main model limitations. It does not allow investigation of the effects of tectonics on sediment routing in 3-D as expected from the interaction of multiple faults and the subsequent formation of relay ramps (e.g., Athmer & Luthi, 2011). Local interaction between erosion, faulting, and sedimentation is essentially a 3-D problem which requires a high-resolution surface process model coupled to a 3-D thermo-mechanical model. The resolution of the tectonic model is limited and affected by our choice to include the sublithospheric mantle in the model domain required to resolve isostasy. The main drawback is that faults are resolved at a resolution of 2 km.
The study of the interactions between tectonics and surface processes is essentially a timescale problem. Major earthquakes aside, the structural evolution of tectonic settings occurs over timescale of hundreds of thousands to millions of years while erosion typically act on much shorter timescales (hundreds to thousands of years depending on climatic forcing and erodibility). The temporal and spatial variability of erosion and transport efficiency is likely to result in tectonics and surface processes being out of sync most of the time. Rapid deposition of large quantity of sediment offshore might be sufficient to alter fault dynamics but depend on the distance between sources and sinks (likely to be short in the case of high elevated rift flanks but potentially very large for continental scale rivers), the duration of transport (which depends on the capacity of the rivers to carry sediment loads), and the rate of tectonic subsidence of the basin (function of the tectonics and the thermal field). We do not enforce mass balance in our experiments, and material is allowed to leave the domain from the left and right sides. This probably conforms to what happens in natural systems where eroded material can potentially be transported away from the rift. Sediments eroded from the flanks facing toward the rift end-up being deposited in the main basin. The effect of erosion and sedimentation on the tectonic model are averaged over the width of the model which effectively assumes a maximum distance of transport of 200 km along the rift axis. That width has been chosen to limit the size of the computational domain while capturing what we believe to be a reasonable section of a rift. Our Model 1B represents an end-member case where erosion and deposition are relatively localized and where distance and duration of transport between eroded areas and basins are short. In that case, sediment accumulation does show an effect on faults dynamics but mostly on the proximal faults and it is unlikely that sedimentation would affect the evolution of more distal faults as the rift evolves. What about long-term effect? It has been suggested that the erosion of the rift shoulder and the subsequent loading of the nearby basins would alter/create horizontal stress gradients in the crust which would be relaxed through viscous flow in the lower part of the crust (Gartrell, 2000). How important that effect could be remains unclear and would require more targeted work using numerical model with better resolved lower crust. The problem of the balancing of lateral loads across the lithosphere also require more efficient ways to monitor the isostatic response of the systems which remains challenging for large geodynamics models.
Our models show symmetric margins similarly to Theunissen and Huismans (2019) but in contrast to models from Brune et al. (2014) and Andrés-Martínez et al. (2019). Huismans and Beaumont (2003) showed that the large-scale asymmetry observed in numerical models is likely the result of feedback between the dominant rheology (in a layered lithosphere context) and the parameterization of the strain softening. Softening of the dominant rheology via a combination of viscous and frictional-plastic strain softening promotes rift asymmetry but is highly dependent on the extension rate. The extension rate chosen in our experiments and the lack of viscous strain softening is probably responsible for the lack of asymmetry in our models. This is clearly a limitation of our work as many real passive margins around the world exhibit large-scale crustal asymmetry (e.g., Brune et al., 2014).
The extension rate chosen for our experiments is probably in the higher range of rates observed in nature (Brune et al., 2016). The tectonic response to erosion and the potential feedback due to sedimentation would be different for a slower rift as the ratio of tectonic uplift to erosion/sedimentation would be different. The amount of accommodation space available for the sediment would also likely be different which would affect the way sediment is distributed within the rift. It should therefore be noted that the landscape predicted by the models is not only the result of a specific choice of K value but also extension rate.
Our thermo-mechanical models use a visco-plastic representation of the lithosphere. Elasticity was considered to be unimportant for studying the phenomenon explored in this paper. However, it is worth mentioning that the importance of elasticity remains a matter of debate in the community. Olive et al. (2016) showed that contrary to visco-elasto-plastic model visco-plastic models such as the one used in this study are strongly dependent on the product of extension rate and viscosity. This has consequence on faulting development and may affect the resulting topography. In our case the product is high and will tend to favor rigid block offsets.
The flexural response to loading and unloading through erosion and deposition is through the visco-plastic behavior of the strong upper crust and upper mantle lithosphere. It has been recognized that large sediment fans such as the Amazon delta may flex the lithosphere by a few kilometers over distance of several 100 km (Watts et al., 2009). However, this largely depends on the strength of the underlying lithosphere at the time of deposition. Using a visco-elasto-plastic rheology could potentially change the wavelength of the zone affected by regional isostatic subsidence and could also intensify the uplift and thus the erosion onshore. Burov and Poliakov (2003) used a visco-elastic-plastic rheology in their study of the erosional forcing of basin dynamics and inferred a widening of the basin with increased erosion. We do not observe this behavior which could be related to large-scale flexural response of the lithosphere. Additionally, it has been shown that including elasticity may help develop more realistic fault geometries (e.g., Lavier et al., 2000), with dip angles satisfying the Coulomb failure criterion. In the model presented here, the mesh resolution and the geometry of the weak zone used to initiate plastic deformation are the main factors controlling the initial angle of the shear zones. The observed angle is usually close to 50 • and corresponds to the Roscoe angle ( = 45 • ± Φ 2 , with the angle of friction, Φ = 15 • , see Table 1) and is in the range of accepted mechanical values (Kaus, 2010). Higher tectonic model resolution will affect the distribution and localization of strain and may result in a different structural expression of faulting (Kaus, 2010).
An important assumption made in the surface process model FastScape is that river erosion can be modeled using an effective water discharge. The variability of water discharge in both time and space is a common feature of natural systems, and studies have shown that discharge distributions may strongly influence the dynamics of the stream-power law (e.g., Lague et al., 2005;Tucker, 2004). Furthermore, the linear dependency of incision rate to slope (n = 1) used as a computational convenience does not reflect highly variable climatic settings where nonlinearities affect interactions between erosion and tectonics (Lague et al., 2005).
Sediments accumulated in the basin have the mechanical and thermal properties of the upper crust. The models do not explicitly account for an ocean layer, which would affect the normal stress on the free surface in the submarine domain. We do not model the effects associated with porosity, water content, or compaction, which may affect the tectonic response but are beyond the scope of our study. Furthermore, although accumulation of sediment limits the cooling of the central part of our model, there is no thermal blanketing effect that may result from lower thermal-conductivity sediments. This is an important factor that may affect the absolute amount and duration of subsidence in the basin (Cacace & Scheck-Wenderoth, 2016).

Conclusions
We present a series of forward numerical models coupling a large-scale two-dimensional thermomechanical model to a plane-view 2-D surface processes model. The model is used to investigate the importance of surface processes in the dynamics of rift systems. The model results show that the structural style and dynamics of rifting are mostly controlled by the rheology of the lithosphere. Efficiency of erosion controls the landscape morphology and modulates the tectonic response in the crust.
• The development and long-term support of topography is dependent on the strength of the coupling between the crust and the mantle lithosphere. Strong coupling promotes high topography as the integrated strength of the lithosphere is sufficient to support the additional load. Weak coupling results in the stress being relaxed via viscous flow in the lower part of the crust and leads to more subdued topography. • The development of the rift is controlled by the tectonics. The structure of the rift reflects the rheology of the upper crust while erosion and sedimentation only weakly modulate its evolution: -High erosion efficiency translates into fast landscape response to tectonic uplift: It limits the development of high topography and creates complex landscape morphologies. -Low erosion efficiency allows for longer standing high topography and results in morphologies showing sharp relief with relatively high peaks and deeply incised streams. -The nature of the onshore-offshore transition depends on both the rheology of the crust and the erosion efficiency: Spatial variation in the accommodation space depends on whether the deformation is localized or distributed and whether the basin is over or underfed with sediment. -The transfer of mass between the continent and the basin enhances fault development, increases offsets, and keeps them active over a longer period.
These model results make an important contribution toward a better understanding of the complex interactions and feedback between surface processes and tectonics in rift settings. The main challenge is to transpose those experiments into 3-D at high spatial resolution in order to simulate the range of mechanical and thermal interactions.