The role of surface processes in basin inversion and breakup unconformity

At divergent plate boundaries, extensional tectonics lead to subsidence, continental rifting and the formation of continental margins. Yet, within this extensional context, transient compressional structures (stress inversion) and phases of uplift (depth inversion) are frequently recorded with no corresponding change in plate motion. Changes in gravitational potential energy during the rifting process have been invoked as a possible source of compressional stresses, but their magnitude, timing and relationship with depth inversions remain unclear. Using high-resolution 2D numerical experiments of the full rifting process, we track the dynamic interplay between the far-field tectonic forces, loading and unloading of the surface via surface processes, and gravitational body forces. Our results show that rift basins tend to localize compressive stresses, they record transient phases of compressional stresses up to 30 MPa and experience a profound depth inversion, 2 km in magnitude, when sediment supply ceases, providing a novel explanation for the breakup unconformity, a well-documented phase of regional uplift typically associated to continental breakup.


INTRODUCTION
Continental divergence leads to normal faults, subsidence, continental rifting and eventually to the formation of conjugate continental margins. In this tectonic setting, it is surprising to frequently observe short transient episodes of basin inversion in the form of i/ surface uplift, ii/ contractional re-activation of normal faults, iii/ long-wavelength, shallow-angle unconformities, and iv/ broad low-amplitude folds (e.g., Withjack et al., 1998;Lundin andDoré, 2002, 2011;Holford et al., 2014). Often, these transient episodes of basin inversion do not correlate with change in plates motion (Withjack et al., 1998;Schlische et al., 2003;Lundin and Doré, 2011), and in the absence of nearby collisional processes, the rift-push force is often invoked as a possible driver of contractional inversion (Dewey, 1989;Bott, 1992;Boldreel and Andersen, 1993;Keleman and Holdbrook, 1995;Withjack et al., 1998;Rey, 2001;Huismans et. al., 2001;Davis and Kusznir, 2002;Pascal and Cloetingh, 2009;Mondy et al., 2018). Here, we use 2D thermo-mechanical numerical experiments of lithospheric scale extension to document the dynamic evolution of continental margins from rifting to drifting. We track through time the stress and strain rate fields, as well as the vertical motion of the pre-rift basement. Results show temporally and spatially variable patterns of extension and compression of similar length scale to that described in nature. In addition, our experiments document a prominent depth inversion of the basement before breakup, from subsidence to uplift, which matches in amplitude and timing the well-documented breakup unconformity often recorded along continental margins (Falvey, 1974;Franke, 2013;Soares et al., 2012).
The magnitude of this inversion depends non-linearly on the basins depth and its timing correlates with the waning of sediment supply due to the increasing distance of the continent.
These results illustrate the importance of surface processes in the dynamics of continental rifting, and may explain some of the complexities of continental breakup (e.g. Péron-Pinvidic et al., 2007).

NUMERICAL EXPERIMENTS PARAMETERS AND SETUP
Using 2D coupled thermo-mechanical experiments, we quantify the evolving stress field during the formation of continental margins. We use Underworld 2.0, an opensource finite-element numerical framework (Moresi et al., 2003(Moresi et al., , 2007Beucher et al., 2019) to solve the Stokes equations to obtain velocity and pressure fields from which strain rate and stress tensors can be derived. Figure 1 summarizes key aspects of our experimental setup. All materials are viscoplastic, with the viscosity dependent on temperature stress, strain, strain rate and melt-fraction.
Plastic behaviour for all rock materials is modelled via a Drucker-Prager yield criterion linearly sensitive to accumulated strain to simulate strain-weakening, which is clipped when the strain reaches 20%. We also impose a stress limiter of 150 MPa and 300 MPa for the crust and mantle respectively, to approximate additional rheological processes such as Peierls creep, and to limit the strength of the lithosphere within observational range (Demouchy, et al., 2013;Zhong and Watts, 2013). The initial geotherm derives from imposing a surface temperature of 20 ºC and a constant temperature of 1353ºC from 140 km to the base of the model. The upper crust and sediments have a radiogenic heat production twice as large as that of the lower crust. This results in an initial temperature at the Moho of 520 ºC. The base of the thermal lithosphere (i.e., isotherm 1350 ºC) initially at a depth of 140 km is let to evolves in a self-consistent manner ( Fig. 1). To capture the decrease in viscosity due to partial melting, we impose a 3-order of magnitude viscosity decrease as melt fraction increases from the solidus to the liquidus. Since our code does not allow melt to escape from its source, we impose a maximum melt-fraction of 30% in the crust and 2% in the mantle. In our model, partial melting also impacts on the geotherm as latent heat energy is consumed in the process. All rheological and thermal parameters are listed in Supplementary Table 1. Extensional deformation is driven by applying a total divergent velocity of 1, 2, 3, or 4 cm/yr, equally partitioned between the left and right vertical walls. Isostasy is achieved by imposing an upward traction at the base of the model matching the initial downward vertical stress. We apply a simple infilling rule to approximate surface processes. When air particles fall below 1km below sea level, fixed through time at elevation 0 m, they become sediment particles, which then record time of 'deposition'. Conversely, when a particle of rock material exceeds a given height, it is "eroded" by switching it to air material. Sedimentation is limited in time by a distance condition on adjacent tracers initially positioned one kilometre apart at the Moho.
Once the distance between two adjacent tracer reaches an imposed threshold (we have tested 30, 60, 80 and 100 km) sedimentation stops. This prevents sediments from accumulating in the ocean far away from the continental source of sediments and allows to investigate how cessation of sediment supply impacts rifting. For comparison we also ran a suite of experiments in which no sediment input was imposed. Finally, we map the Andersonian stress regime throughout the lithosphere by determining which principal stress axes is within 30º of vertical (Zoback, 1992;Mondy et al., 2018).  Figure   SF2 and SF3). Figure 4 shows the evolution of the average depth of the basin's basement (see Suppl. Data Figure SF4 for all other experiments) for a suite of experiments with divergence velocity of 1, 2 and 4 cm/yr, and a maximum distance between adjacent Moho tracers of 30, 80 and 100 km at which sedimentation stops. In experiments with no sediment input, the basement subsides monotonically until breakup occurs. When sediments are supplied, the rate and magnitude of the basement subsidence significantly increases due the weight of the sediment infill (Burov and Poliakov, 2001), and breakup is delayed. In all experiments, we observe that for the same percentage of extension, the magnitude of subsidence increases with the extensional velocity (green squares in Figure 4 show the amount of subsidence at 20% extension). However, because the condition for cessation of sedimentation is meet earlier at faster extensional velocity, the absolute maximum subsidence is smaller for fast extension ( Figure 4). As the conjugated margins move apart beyond the threshold distance for sedimentation (shown as circles on Figure 4), all experiments record a significant inversion from basement subsidence to uplift, proportional to the amount of subsidence and sediment loading. This is particularly visible in experiments with moderate divergence velocities (2 cm/yr) in which the inversion reaches 2 km (Fig. 4). During this inversion, the maximum, median, and average depths of all basins decrease by about 1 to 2 km, and the sediment pile is stretched and thinned by normal faulting and fault block rotation, some basins being completely attenuated during this process. At that stage, the asthenosphere rises to reach its neutral buoyancy level (~3000 m below sea level), continental breakup is reached (triangles in Figure 4), and the maximum basin depth tends to remain steady or slowly decreases.

DISCUSSION
Transient inversions confined to sub-basins have been documented along many continental margins (e.g., Holford et al., 2009;2014;Cloetingh et al., 2008). These transient inversions are associated to either localized contractional reactivation, or epeirogenic uplift, or to the doming of the asthenosphere during rifting which increases a local gravitational potential energy (GPE) in the rift region forcing stress inversion in sub-basins on the adjacent stretched margins (Suppl. Data Figure SF5, Turcotte and Emerman, 1983;Buck, 1991;Huismans et al., 2001;Rey, 2001;Davis and Kusznir, 2002;Mondy et al., 2018). Of particular importance is the well-documented phase of regional uplift and erosion, called breakup unconformities, which typically precedes the onset of seafloor spreading (Falvey, 1974;Braun and Beaumont, 1989;Franke, 2013;Soares et al., 2012). Our experiments reveal a link between the cessation of sedimentation and a regional uplift 1 to 2 km in magnitude. We attribute this inversion to the dynamic imbalance resulting from cessation of sedimentation and surface loading of the thinning lithosphere. This imbalance is relaxed through isostatic uplift which drives the observed depth inversion. Our experiments show that basin inversions are particularly strong for moderate extension rates.
Indeed, during fast extension only shallow basins with a reduced capacity for depth inversion form. In contrast, during slow extension, diffusive cooling becomes important limiting the potential for isostatic uplift and depth inversion.

CONCLUSIONS
We show that in a divergent setting, transient periods of stress and depth inversions are the outcome of competing deep to surface processes altering the dynamic balance at both local and regional scales. Upon lithospheric thinning, the asthenosphere is under the isostatic imperative