Segmentation of Rifts Through Structural Inheritance: Creation of the Davis Strait

Mesozoic‐Cenozoic rifting between Greenland and North America created the Labrador Sea and Baffin Bay, while leaving preserved continental lithosphere in the Davis Strait, which lies between them. Inherited crustal structures from a Palaeoproterozoic collision have been hypothesized to account for the tectonic features of this rift system. However, the role of mantle lithosphere heterogeneities in continental suturing has not been fully explored. Our study uses 3‐D numerical models to analyze the role of crustal and subcrustal heterogeneities in controlling deformation. We implement continental extension in the presence of mantle lithosphere suture zones and deformed crustal structures and present a suite of models analyzing the role of local inheritance related to the region. In particular, we investigate the respective roles of crust and mantle lithospheric scarring during an evolving stress regime in keeping with plate tectonic reconstructions of the Davis Strait. Numerical simulations, for the first time, can reproduce first‐order features that resemble the Labrador Sea, Davis Strait, Baffin Bay continental margins, and ocean basins. The positioning of a mantle lithosphere suture, hypothesized to exist from ancient orogenic activity, produces a more appropriate tectonic evolution of the region than the previously proposed crustal inheritance. Indeed, the obliquity of the continental mantle suture with respect to extension direction is shown here to be important in the preservation of the Davis Strait. Mantle lithosphere heterogeneities are often overlooked as a control of crustal‐scale deformation. Here, we highlight the subcrust as an avenue of exploration in the understanding of rift system evolution.

The Labrador Sea, Davis Strait, and Baffin Bay (Figure 1a) formed due to Mesozoic to Cenozoic divergent motion between Greenland and North America (Abdelmalak et al., 2018;Chalmers & Pulvertaft, 2001;Hosseinpour et al., 2013;Jauer et al., 2019;Peace et al., 2017;Peace, Dempsey, et al., 2018;Wilson et al., 2006). Rifting prior to the opening of the Labrador Sea may have started as early as the Late Triassic to Jurassic, based on ages obtained from dike swarms in southwest Greenland that are interpreted to be related to early rifting (Larsen et al., 2009). Breakup from south to north between Greenland and Canada resulted in oceanic spreading in the Labrador Sea and eventually Baffin Bay (Chian et al., 1995;  (c) Simplified overview of the basements that comprise the NW Atlantic borderlands in a prerifting and breakup configuration modified from Kerr et al. (1997) and St-Onge et al. (2009). (d, e) The NW Atlantic at 60 and 35 Ma, respectively, reconstructed using the model of Matthews et al. (2016) and shown with the calculated extensional directions from Abdelmalak et al. (2012). 1. Rifting south of Davis Strait in the Labrador Sea produced new oceanic crust. 2. Rifting north of Davis Strait in Baffin Bay produced new oceanic crust. 3. A right-stepping segmentation geometry (the Davis Strait) was formed to link Labrador Sea with Baffin Bay. 4. Continental crust is preserved in the Davis Strait during rifting.
The West Greenland-Eastern Canada realm comprises multiple Archean and Proterozoic geological domains, reflecting a complex, multiphase evolution (e.g., Kerr et al., 1997;Grocott & McCaffrey, 2017;St-Onge et al., 2009). The evolution of these domains, including their correlation to a pre-Cretaceous reconstruction is dealt with in detail in van Gool et al. (2002) and St-Onge et al. (2009), as such only the most salient points relevant to this study are reiterated here.
These Archean and Proterozoic domains, and the preexisting structures they contain, likely influenced the Mesozoic-Cenozoic rifting, breakup, and transform system development through the process of structural inheritance (Japsen et al., 2006;Peace et al., 2017;Peace, Dempsey, et al., 2018;Watterson, 1975;Wilson et al., 2006). This previous work has shown that crustal structural inheritance may have controlled the large-scale geometry of breakup and transform systems, the geometry and kinematics of rift-related faulting, and potentially also the location of rifting and breakup-related magmatism. As such it is important to understand the formation of the different basement units that comprise this study area ( Figure 1c). Principally, from north to south the study area herein comprises the following gross tectonic units: the Nagssugtoqidian and Torngat orogens; the North Atlantic Craton and the Nain Province (NP); and the Makkovik Province and the Ketilidian Mobile Belt (Figure 1).
The once continuous Archean North Atlantic Craton is now distributed between Greenland and Labrador (where it is called the NP; St-Onge et al., 2009;Figures 1c-1e). The North Atlantic Craton is bordered to the north and west by segments of Palaeoproterozoic orogenic belts that are tectonically related to the Trans-Hudson Orogen including the Nagssugtoqidian Orogen and Rinkian fold belt on the north side and the Torngat Orogen on the west side of the craton (St-Onge et al., 2009).
The Nagssugtoqidian Orogen ( Figure 1c) is a belt of Palaeoproterozoic deformation and metamorphism in West Greenland considered to have developed simultaneously with the Torngat Orogen in northern Labrador (van Gool et al., 2002). Although the precise spatiotemporal relationship between these orogenic belts is questioned (Scott, 1999), they are interpreted to have formed part of the same Palaeoproterozoic passive margin prior to ocean closure and continental collision with the North Atlantic Craton and NP (van Gool et al., 2002;Grocott & McCaffrey, 2017).
The dynamics of West Greenland rifting and the preservation of the continental Davis Strait is currently a topic of active research, with lithospheric inheritance being discussed as a potential controlling mechanism (Peace et al., 2017;Peace, Dempsey, et al., 2018;Wilson et al., 2006). In this study, we outline a two-phase tectonic history where mantle lithosphere inheritance is generated and then contributes to crustal deformation during subsequent rifting ( Figure 2). We hypothesize that a Palaeoproterozoic collision, which featured the North Atlantic Craton and produced the Nagssugtoqidian Orogen (van Gool et al., 2002), would have left mantle lithosphere scarring during the continental suturing (e.g., Calvert et al., 1995;Vauchez et al., 1997;Holdsworth et al., 2001). Although there is no direct evidence of a mantle structure, deformation during the Nagssugtoqidian Orogen is thought to be on a lithospheric scale, rather than simply crustal scale (Watterson, 1975;Grocott, 1977;van Gool et al., 2002), and as a result we consider that a mantle suture is likely to have been produced (Figure 1; e.g., Biryol et al., 2016;Calvert et al., 1995;Calvert & Ludden, 1999;Hopper & Fischer, 2015;Lie & Husebye, 1994;Mickus & Keller, 1992;Morgan et al., 1994;Steer et al., 1998;Schiffer et al., 2014;Schiffer et al., 2016;Vauchez et al., 1997;Vauchez et al., 1998).  Gool et al., 2002). Line A'-B' in (a) gives an approximate location for the schematic in (b). Annotations: ARF, Arfersiorfik intrusive suite; SIS, Sisimiut charnockite suite; ISB, Itivdleq steep belt; ITZ, Ikertôq thrust zone; NISB, Nordre Isortoq steep belt; NSSZ, Nordre Strømfjord shear zone; SNF, southern Nagssugtoqidian front; p.e.s., present erosion surface. CNO, NNO, and SNO are the central, northern, and southern Nagssugtoqidian Orogen, respectively. (e) We propose this deformation history would leave a mantle scar (highlighted by dashed green lines in c and d). The ITZ is also the proposed location of the suture line in this orogen (van Gool et al., 2002).
Indeed, Watterson (1975) first identified the Palaeoproterozoic Nagssugtoqidian orogenic belt as a lithosphere-scale boundary due to the presence of Cambrian age kimberlites that are crosscut by Mesozoic pseudotachylytes, which was a finding later confirmed by Grocott (1977).
In this study, we model upper crust inheritance and a mantle lithosphere scar that approximates the shape and extent of the suture surrounding the North Atlantic Craton. Below, in a suite of numerical simulations, we analyze the influence of lithosphere inheritance for generating rift tectonics appropriate to the Labrador Sea, Davis Strait, and Baffin Bay.

Methods
The role of three-dimensional (3-D) lithosphere structure in a continental extension tectonic setting similar to that of the Davis Strait is investigated. The models are implemented in a high-resolution 3-D Cartesian box (Figure 3), using the numerical code ASPECT Heister et al., 2017;Kronbichler et al., 2012;Rose et al., 2017), which uses the finite element method to solve the system of equations that describes the motion of a highly viscous fluid. Specifically, we use a nonlinear viscous flow (dislocation creep) and Drucker-Prager plasticity for our model rheology (e.g., Naliboff & Buiter, 2015).

Experimental Setup
The 3-D numerical experiments are conducted within a model domain of 800 km (x axis) by 800 km (y axis) and 600 km vertically (z axis). The computational grid is uniform laterally, but resolution varies vertically with higher resolution prescribed in the top 80 km of the model (from the surface to 80-km depth). Below, the resolution becomes more coarse, with a reduction in resolution between 80 and 180 km, then finally the lowest resolution from 180-km depth to the bottom of the model (supporting information Figure S1). There are 1.7 million active cells in the model, with a horizontal resolution of ∼1 km at the surface.
The 3-D simulations are very computationally expensive, producing 147 million degrees of freedom and needing around 80-GB memory. For most cases, the models used 416 CPUs and took ∼16,000 hr of computational time to generate 12 Myr of deformation on ComputeCanada's Niagara cluster (Loken et al., 2010). Here, three main sections of the North Atlantic Craton suture (south, oblique, and east) are generated into a scar (as shown by red arrows). The scar is applied as a zone of weakness (with a lower angle of internal friction than surrounding material). Abbreviations: UC = upper crust; LC = lower crust; ML = mantle lithosphere.

Governing Equations
In this study, we solve the equations of conservation of momentum, mass, and energy after assuming an incompressible medium with infinite Prandtl number: In the equations above, is the viscosity, . is the strain rate tensor, u is the velocity vector, k is the thermal conductivity, is the density, C p is the thermal heat capacity, H the internal heat production, P the pressure, g gravity, and T the temperature.
The upper crust, lower crust, mantle lithosphere, asthenosphere, and scar are represented by five distinct compositional fields that are advected with the computed flow velocity. For each field c i , this formulation introduces an additional advection equation to the system of equations: Equations (1-4) are solved using the finite element method, where the domain is discretized into hexahedral finite elements and the solution (e.g., velocity, pressure, temperature, and compositional fields) is expanded using Lagrange polynomials as interpolating basis functions (as outlined in Glerum et al., 2018). In this study, we employ second-order polynomials for velocity, temperature, and composition and first-order polynomials for pressure (Q2Q1 elements; Donea & Huerta, 2003). The equations are solved using an iterative Stokes solver (for more details see Kronbichler et al., 2012). The models use a temperature-dependent density in all governing equations, but no pressure dependence, since the model is incompressible. We solve the temperature and composition equation once at the beginning of each time step and then iterate out the Stress exponent 4.0 3.0 3.5 3.5 3.5 Activation energy Note. For angle of internal friction and cohesion, strain weakening occurs over the range 0 to 0.5 (e.g., Brune et al., 2017) and weakens by 50%. UC: upper crust; LC: lower crust; ML: mantle lithosphere; A: asthenosphere; ML scar: mantle lithosphere scar.
solution of the Stokes equation to a solver tolerance of 10 −4 (with a maximum number of nonlinear rheology iterations set to 10).
We use a nonlinear viscous flow (dislocation creep) and Drucker-Prager plasticity for the model rheology and follow a setup similar to previous studies (e.g., Brune et al., 2014;Brune et al., 2017;Huismans & Beaumont, 2011;Naliboff & Buiter, 2015). The viscosity for dislocation creep is defined as follows: where A is the viscosity prefactor, n is the stress exponent (n > 1), . e is the square root of the deviatoric strain rate tensor second invariant, E is activation energy, V is activation volume, and R is the gas exponent (Karato & Wu, 1993;Karato, 2008).
Viscosity is limited through a "yielding" mechanism. Plasticity limits viscous stress through a Drucker-Prager yield criterion, where the yield stress in 3-D is Above, C is the cohesion and is the angle of internal friction. If is 0, the yield stress is fixed and equal to the Von Mises yield criterion. When the viscous stress (2 e ) exceeds the yield stress, the viscosity is rescaled back to the yield surface y = y /(2 . e ) (e.g., Thieulot, 2011). This method of plastic yielding is known as the Viscosity Rescaling Method (Kachanov, 2004;Willett, 1992).
In the models here, strain weakening is implemented for the internal friction angle and cohesion; they are linearly reduced by 50% of their value (from 20 • and 20 MPa; e.g., Bos, 2002) as a function of the finite strain magnitude. This weakening occurs between 0 to 0.5 strain, which is a range used in Brune et al. (2017) rifting study. Other strain ranges for weakening were tested and the findings are presented in supporting information Figures S2-S4.
Compositional material fields (upper crust, lower crust, mantle lithosphere, asthenosphere, and scarring) are each assigned individual values of thermal diffusivity, heat capacity, density, thermal expansivity, and rheological parameters (Table 1). If more than one compositional field is present at a given point (such as in the vicinity of material interfaces), viscosities are averaged with a harmonic scheme (e.g., Glerum et al., 2018). In addition, strain is also tracked as a compositional field.
The rheological setup of these models closely follows that of Naliboff and Buiter (2015). Table 1 outlines the rheological parameters used for the different compositional layers. The upper crust implements a wet quartzite flow law (Rutter & Brodie, 2004), lower crust applies wet anorthite (Rybacki et al., 2006), and the mantle dry olivine (Hirth & Kohlstedt, 2003). All the viscous prefactors described in Table 1 are scaled to plane strain from uniaxial strain experiments (Ranalli, 1995).

10.1029/2019TC005578
An initial reference viscosity of 10 22 Pa·s is applied to each compositional field in the models due to the strain rate dependence of viscosity and the lack of an initial guess for the strain rate for the first time step (Glerum et al., 2018). This initial reference viscosity was changed in the setup of the numerical models and not found to change the outcome of the study. During subsequent time steps, the strain rate of the previous time step is used as an initial guess for the iterative process. The final effective viscosity is capped by a (user-defined) minimum viscosity (set at 10 18 Pa·s) and maximum viscosity (set at 10 26 Pa·s) to avoid extreme excursions and to ensure stability of the numerical scheme. In the models presented here, we apply a viscosity range of 8 orders of magnitude. However, for the majority of models, the viscosity profile stays well within the maximum and minimum cutoffs.

Lithosphere Scarring
In the modeling of a mantle suture, we specify an inherited plane of weakness that has remained over a long period of time (in this case, since the Palaeoproterozoic). There are a number of mechanisms where a mantle lithosphere suture could remain weak over time (e.g., Erdős et al., 2014;Heron et al., 2018;Manatschal et al., 2015;Petersen & Schiffer, 2016), one of which is through grain size reduction of peridotite mylonites at ancient plate boundaries (Bercovici & Ricard, 2014). The mantle lithosphere scar modeled here is 10 km thick, dipping at an angle of 45 • from the horizontal from 32-km depth down to 52 km (Figure 3a), and rheologically weak by having a reduced angle of internal friction compared to the surrounding material (Table 1). Due to the lack of high-resolution geophysical imaging at depth in the region, there is uncertainty in the dip of a mantle structure (or even if there is a heterogeneity present). However, the influence of changing shape and dip angle of generic styles of such weak scars is explored in detail in Heron and Pysklywec (2016) Figure 4 shows the velocity azimuth and magnitude for the Davis Strait using plate reconstruction histories  and the GPlates software (Müller et al., 2018). According to this reconstruction, between 200 and 120 Ma there was no significant extension between Greenland and Eastern Canada. During the early to mid-Cretaceous, extension initiates with an azimuth of between 20 • and 30 • in present-day coordinates (Abdelmalak, 2010). However, it was not until the late Mesozoic/early Cenozoic that there was significant extension in the region. Thus, we identify a "Phase 1" between 75 to 55 Ma as having an average extension velocity of 1 cm/year at an azimuth of approximately 60 • . The azimuth of continental separation rotates anticlockwise in the Cenozoic (Figure 4a), and we identify therefore a "Phase 2" with a higher velocity magnitude at a high angle to the Phase 1 extension direction (Figure 4c). Following the work of , Phase 1 produced the rift shape alongside the spreading in the Labrador Sea and Baffin Bay. In this study, we focus on the initial stages of the rift system and investigate Phase 1 closely. Phase 2 is not thought to have significantly thinned the Davis Strait, although a number of strike-slip crustal features are due to this approximately NE-SW extensional activity .

Extension Rate and Boundary Conditions
To model Phase 1 (Figure 4c), we apply a prescribed boundary velocity on the north and south boundaries, and tangential velocity boundary conditions on the west, east, and base walls of the model, and a free surface on top (Rose et al., 2017). We have modeled the Cartesian 3-D box large enough so that deformation driven from the scarring is not influenced by the tangential boundary conditions (as described below).
The prescribed boundary condition on the north wall is a 0.5-cm/year extension for the lithosphere (120 km) and a return flow of −0.3 cm/year for the bottom 200 km of the box. In between, the velocity tapers from 0.5 to 0 cm/year from 120-to 225-km depth, and from 0 to −0.3 cm/year from 200-to 400-km depth. The reverse is applied to the west wall, with 0.5-cm/year extension for the lithosphere. After extensive testing, we found this velocity profile as a boundary condition to provide stable solutions while maintaining mass balance (meaning no additional volume is added to the box). This prescribed boundary velocity produces an extension rate of 1 cm/year in the lithosphere. This falls within the appropriate velocity magnitude as given in Figure 4.
The free surface allows topography to form and is formulated using an Arbitrary Lagrangian-Eulerian framework for handling motion of the mesh (for more details please refer to Rose et al., 2017). All of the calculations presented here have ∼5,500,000 free surface degrees of freedom.  Matthews et al. (2016). From these values we approximate the two phases of the rift evolution (c).

Thermal Model Setup
An initial temperature field is prescribed (Figure 3a) but is allowed to evolve during the simulation. The initial temperature follows a typical continental geotherm (Chapman, 1986) with no lateral variations. Our initial condition models the late Mesozoic extension of two continental blocks, which collided in the Palaeoproterozoic (Figure 2). Therefore, the closure of the oceanic basin to accrete northern Greenland to the North Atlantic Craton occurred >1 Gyr in the past, and therefore, there are no remaining thermal perturbations from that tectonic event. The temperature equation for calculating the initial geotherm is given as follows: where T o is the temperature at the top of the specific layer, H is the heat production, q is the heat flow through the surface of the specific layer, k is the thermal conductivity, and z is the depth. Table S1 gives the values for the thermal constraints required to generate the geotherm. As described in Naliboff and Buiter (2015), we use a high conductivity in the asthenosphere to maintain the high adiabat in the layer, and to generate a constant heat flux into the lithosphere (Pysklywec & Beaumont, 2004).

Results
Below, we present numerical models of continental extension in the presence of lithosphere inheritance related to West Greenland. Figure 5 shows the impact of crustal and mantle lithosphere inheritance on the evolution of the Davis Strait.

Crustal Inheritance and Mantle Inheritance
To test the applicability of crustal inheritance in generating the first-order tectonics as seen in the Davis Strait, we apply in Model C1 a crustal fault from previous geological studies van Gool et al., 2002;Wilson et al., 2006) that is similar in geometry to the Ikertôq thrust zone (Figures 2 and  5a). After 15 Myr of E-W extension, the surface strain rate in Model C1 indicates that the inheritance does not localize in the region of the scar (Figure 5b). The rifting pattern does not generate the relevant tectonic features, as outlined in the four-point checklist for the Davis Strait.
Model CM1 shows the rifting pattern across the region after 15 Myr of extension in the presence of the crustal scars in Model C1 and a mantle lithosphere inherited structure (Figure 5a). The mantle lithosphere scar represents the Nagssugtoqidian suture separating the North Atlantic craton from the continental material to the north (Figure 3b). In Model CM1, the strain rate replicates the right-stepping segmentation of the rifted conjugate margins (Figure 5b), with the surface tectonics at 15 Myr producing spreading in the north and south of the model (as shown in red in Figure 5c) and preserving continental upper crust material across the oblique section of the suture (as shown in blue in Figure 5c). As a result, Model CM1 meets the four-point checklist for the rift and ocean basin architecture.
In Model M1, the crustal inheritance from CM1 is removed, yet there is little difference in the evolution of the rift system ( Figure 5). In this instance, the tectonics of the region are dominated by the mantle lithosphere structure. Applying only crustal inheritance that mimics the shape of the Nagssugtoqidian suture (e.g., a shallower scar as used in Model M1) is shown in Model C2 (Figure 5a). The tectonic evolution of the region is very different as compared to the similar Model M1. In Model C2, the southern limb of the suture begins to spread and propagates north-south (Figure 5c). Reapplying the mantle lithosphere suture in Model CM2 (Figure 5a) dominates the evolution of the region, as previously shown in Model CM1. Figure 6 shows the rift evolution of the reference Model M1, which enables an interpretation of what is occurring to reproduce the appropriate tectonic patterns of the Davis Strait. The surface strain rate pattern in Figure 6a shows an initial reactivation of the southern limb (A) and the oblique portion of the scar (B) at depth generates localized deformation in the crust. After 4 Myr of extension, there is little activity in the north of the model. However, after 7 Myr, the surface strain rate pattern indicates a localization of deformation to the north of the mantle suture (C, Figure 6a). The eastern limb of the suture does not appear to reactivate.  The thinning and spreading of the upper crust is shown in Figure 6b, with spreading developing first in the south (12 Myr) and then in the north (14 Myr), which is in keeping with geological interpretations of Labrador sea and Baffin Bay Peace et al., 2017). There remains a region between the north and south spreading regions that is preserved, but thinned, continental material, which we describe as being a modeled Davis Strait.

Evolution of the Rift
The mantle lithosphere suture plays a significant role in the development of the southern "Labrador Sea" rift-the southern limb of the structure is perpendicular to the extension direction and as such facilitates the rifting and spreading (Figure 6b, A). However, the oblique portion of the scar transmits strain across it (Figure 6a, B) until it reaches the E-W section of the proposed ancient suture. Instead of leading to deformation on this eastern limb of the scar, the extension diverts to being perpendicular to the extension direction in the north of the model (Figure 6a, C). Figure 7 explores the potential role of obliquity in controlling the rifting pattern. The reference setup for Model M1 has a 45 • angle from the extension axis for the oblique portion of the mantle lithosphere suture. In Figure 7, this angle is changed to be less acute (M70, M65) or more acute (M40, M20) to gauge the range of obliquity at which the reference model can still produce Davis Strait tectonics. For Models M70 and M65, the acute angle cannot maintain the full four-point checklist as a right step segmentation is not generated. However, a narrow region of preserved continental material is still produced (Figure 7).

The Complexity of Obliquity
Decreasing this angle of obliquity in M40 (40 • angle from the extension axis for the oblique portion of the mantle lithosphere suture) maintains the four-point checklist. However, a 20 • obliquity to extension direction mantle suture (Model M20) is not able to propagate strain across it and a north-south rift pattern is produced (Figure 7). Figure 7 further highlights the importance of this oblique portion of the mantle suture in Model M1wide and M1gap. In M1wide, the width of the oblique section is increased (with an angle of 45 • ) and still allows strain to propagate across it. Indeed, the spacing between the north and south spreading regions is increased compared to Model M1 (Figure 7). However, if we remove the oblique portion from M1 altogether (e.g., Model M1gap), we produce north-south spreading. The ability of the suture to transmit strain across the oblique portion is paramount to developing the appropriate rift and ocean basin architecture (Figure 7).

Model Comparison With Gravity Data
For the models presented that satisfy the four-point tectonic checklist (e.g., CM1, M1, CM2, M40, and M1wide), the surface evolution of the models is encouraging ( Figure 6). However, it is important to compare the results of the numerical models with independent estimates of subsurface structure. Figures 8 and  S6 show cross sections of Model M1 after 15 Myr for north and south sections of the rift, as well as across the Davis Strait (lines given in Figure 6c).
The rifting dynamics across the model changes significantly from north to south. In the south of the model, the spreading occurs asymmetrically (Figure 8c), with more pronounced necking to the west of the spreading center. However, in the north (Figure 8a), the rift is more symmetric. The preexisting angled mantle suture promotes this asymmetry in the south, whereas in the north the rift propagates without any inherited features. It appears that in the north the spreading occurs as a result of being perpendicular to the extension direction ( Figure 6). Figure 9a shows a subset of Figure 6b and shows the thinned Davis Strait with variable topography (with a 4-km peak in the west of the model). Figure 9b shows a gravity inversion giving the depth to Moho for the region (as described in the supporting information and Welford & Hall, 2013;Welford et al., 2018). Across the Davis Strait we see some areas of shallow Moho toward Greenland in the east (circled in red). This corresponds to the thinning of the crust across our model Davis Strait (circled in red, Figure 9a). Furthermore, our model could produce decompression melting related to the thin mantle lithosphere as outlined in Figure 9a. It is understood that melting occurred during the Paleogene across the Davis Strait (Larsen et al., 2009). Our angled and thinned mantle lithosphere could be a pathway for such decompression melting patterns, and as a result a potential site of magmatic underplating.

Discussion
Results from our modeling show the impact of mantle lithosphere scarring related to a Palaeoproterozoic orogenic event in the development of the complex Mesozoic-Cenozoic rifting and ocean basin formation between Greenland and Canada ( Figure 5), and highlight the potential role of obliquity in the rift evolution (Figure 7). Although a number of studies have previously modeled oblique rifting in three dimensions (Brune et al., 2014;May et al., 2015;Zwaan et al., 2016;Zwaan & Schreurs, 2017;Brune et al., 2017;Farangitakis et al., 2019), our study shows tectonic features related to West Greenland and offers a new geodynamic explanation for the Phanerozoic rift event in the region. Figure 10 outlines a mechanism for the evolution of the Davis Strait during the Paleocene. First, in the region of the present-day Labrador Sea, there is a reactivation of a mantle lithosphere suture related to the accretion of the North Atlantic Craton. The suture, outlined in green as given by van Gool et al. (2002), in this region is perpendicular to the extension direction generated by the plate motion in the late Paleocene/early Eocene  Funck et al., 2012;N2, Gerlings et al., 2009;F12, Funck et al., 2012;and S13, Suckro et al., 2013), used to assess the reliability of the gravity inversion results. . The model rifting first to the south of the Davis Strait follows the geological history of the region, with Labrador Sea spreading occurring before Baffin Bay spreading in the north (Abdelmalak et al., 2018;Peace et al., 2017).

Mechanism
The middle panel of Figure 6 shows the angled portion of the suture that connects the NP and how the north of the North Atlantic Craton plays an important role in the evolution of the Davis Strait (gray region). The obliquity of the suture to the late Paleocene extension direction does not permit the Davis Strait to achieve breakup in the same way as in the Labrador Sea. Although our model Davis Strait undergoes extensive thinning (transtension; Figure 8), the oblique suture delays, and ultimately prohibits spreading. Stress is transmitted across the oblique suture, creating the elevated region that becomes the Davis Strait.
This stress transfer follows the mantle suture until it becomes parallel to the extension direction on Greenland ( Figure 10). Despite the presence of a weak region of mantle lithosphere in the east of the model, the rift propagates north perpendicular to the extension direction into the Baffin Bay region (Figure 6). In our models there is no structural inheritance to the north of our Davis Strait; because of this, the Labrador Sea and Baffin Bay spreading patterns are very different (Figure 9a).

The Influence of the Mantle Lithosphere in Tectonic Processes
It has been noted by many workers that the Mesozoic rifting and Cenozoic margin and basin formation in West Greeland crosscuts basement orogenic belts and cratons (e.g., Buiter & Torsvik, 2014;Larsen & Rex, 1992;Tappe et al., 2007). In particular, the North Atlantic Craton was split by the Labrador Sea, meaning a thin sliver is now located to the west in the Torngat region on Labrador (see Figure 1). This crosscutting relationship prompts the question as to why the Labrador Sea ocean basin opened where it did and not further west at the surface in the Torngat Palaeoproterozoic belts. A possible answer provided by this study is that an east dipping Palaeoproterozoic suture ( Figure 2) would have been located at mantle depths many tens of kilometers inboard of the surface trace of the western margin of the North Atlantic Craton. If this mantle feature localized extensional strain in the crust directly above, as demonstrated by our models, then it would be entirely feasible that a strip of the North Atlantic Craton could end up on the Labrador side of the newly formed ocean basin ( Figure 1). Indeed, this provides a mechanism by which structural inheritance by unseen mantle structures influences upper crustal deformation patterns and creates crustal slivers. In this case promoting cross-cutting narrow margins by necking of the overlying crust (e.g., Wenker & Beaumont, 2018), where the extension direction is perpendicular to the preexisting mantle scar.
The Palaeoproterozoic Nagssugtoqidian orogenic belt to the north of the North Atlantic Craton was first identified as a persistent (>2.5 Gyr) tectonic lineament by Watterson (1975), who regarded the boundary as a lithosphere-scale structure due to the presence of Cambrian age kimberlites that are crosscut by Mesozoic age pseudotachylytes (Grocott, 1977). Subsequent investigation of brittle deformation in exposures of the Nagssugtoqidian Orogen adjacent to the Davis Strait by Wilson et al. (2006) revealed a two-phase model for fault development that is compatible with the development of the Mesozoic to Cenozoic continental margin offshore (Chalmers et al., 1993;Oakey & Chalmers, 2012). Wilson et al. (2006) found that the Phase 1 generally N-S trending normal faults were kinematically compatible with the early NW-SE opening of the Labrador Sea-Davis Strait-Baffin Bay seaway in the Cretaceous to Paleocene. Phase 2 faults are strike-slip and thrust structures that are spatially confined to ductile shear zones within the Nagssugtoqidian (such as the Norder Isortoq shear zone; Figure 2) and explained by partitioning of the wrench deformation that formed the Eocene Ungava transform system (via preexisting structures; Wilson et al., 2006).
The main Palaeoproterozoic shear zones identified as part of the Nagssugtoqidian Orogen continue offshore and control the primary depocenters and later transpressional deformation in the Davis Strait region Peace et al., 2017). Early Cretaceous synrift fault patterns show generally margin parallel NNW-SSE trends in Labrador Sea and Baffin Bay; however, in the Davis Strait region the faults show a broad, diffuse pattern with the main faults rotated clockwise relative to the overall margin trend (Alsulami et al., 2015;Chalmers et al., 1993;Oakey & Chalmers, 2012;Peace et al., 2017). This pattern is compatible with the Davis Strait forming as a zone of transtensional deformation, under local ENE-WSW extension in a right-stepping transfer zone from Labrador Sea into Baffin Bay. In this scenario the Davis Strait was a primary structure formed prior to the change in spreading direction in the Eocene. The coincidence of the Davis Strait transtensional zone with the offshore continuation of the Nagssugtoqidian orogenic shear zones led Wilson et al. (2006) to suggest that deformation was "strongly influenced by basement fabrics such that this region experienced complex 3-D strain." We now suggest this crustal inheritance, which is clearly expressed in the fault patterns, the depositional history of the basins and the overall crustal thickness , was in effect a passive response to the oblique deformation controlled by the mantle scar beneath. The overall right step of the margin, which set up the oblique extensional deformation zone which becomes the Davis Strait was a first-order response to the locus of stretching deformation seeking to follow the proposed mantle scar where it was oblique to the stretching direction. Further east, where the proposed mantle scar related to the ancient suture becomes perpendicular to the overall stretching direction is the point (south end of Baffin Bay) where the locus of deformation resumed its NNW-SSE direction.
Inheritance in rift evolution is commonly attributed to the influence of preexisting crustal structures such as Proterozoic mobile belts in the East African rift and the Red Sea (Corti et al., 2007;Daly et al., 1989;Jarrige et al., 1986). There have been fewer studies on the effects of a pre-existing mantle heterogeneity (e.g., Schiffer et al., 2018). Autin et al. (2013) tested models to explain the evolution of the Gulf of Aden that included an oblique mantle weakness, concluding that observed offset spreading centers were in this case more likely to have been produced by an inherited Mesozoic basin. We suggest that on the west Greenland margin, in the absence of an obvious inherited basin, a mantle scar induced a similar offset spreading center and that offset zone became the Davis Strait.
This study presents a new, deep origin of the inheritance that may drive deformation in a region where only crustal processes have previously been suggested Peace, Dempsey, et al., 2018). It should be noted that it is indeed unexpected that applying a North Atlantic Craton mantle suture (Figure 3b) in the presence of an extension field that is relevant in velocity and orientation to the Paleogene (Figure 4) would produce appropriate rift dynamics for the Davis Strait system ( Figure 6). However, the study here complements a growing body of work that highlights the potential of the mantle lithosphere to play an important role in tectonic processes (Pysklywec & Beaumont, 2004;Babuška & Plomerová, 2013;Jourdon et al., 2017;Salazar-Mora et al., 2018;Phillips et al., 2018;Balázs et al., 2018;Heron et al., 2019).
The tectonic mechanism presented here may also be applicable to other oblique and transform margins overlying significant crustal and lithospheric terrane boundaries (e.g., Equatorial Atlantic of Africa and South America; North West Shelf Australia; and Southern and East coast of Africa).

3-D Modeling
We present 3-D numerical models that are 800 km × 800 km and have a crustal resolution of 1 km. There are distinct advantages to using such models over 2-D simulations, as discussed in Le Pourhiet et al. (2018). However, there are drawbacks related to these higher dimension models. For instance, due to computational expense, we are unable to model dynamically the full evolution of the region. That is, simulate the continental collision that produced the Nagssugtoqidian Orogen and subsequent hypothesized mantle lithosphere sutures in the Palaeoproterozoic, then organically generate the Mesozoic-Cenozoic rifting as a result of far-field plate motion (e.g., Naliboff & Buiter, 2015;Salazar-Mora et al., 2018).
By manually implementing such a mantle scar as an initial condition we negate a lot of the geological history of the region. However, our hypothesis that ancient tectonic activity could produce weak lithospheric structures that remain dormant over long timescales before reactivation is well established (e.g., Vauchez et al., 1997;Holdsworth et al., 2001). Indeed, this study is important as it applies well-established theories regarding mantle lithosphere inheritance (e.g., Bercovici & Ricard, 2014) to a regional geological feature. Here, a mantle lithosphere structure can generate appropriate deformation related to the Davis Strait and follows a number of previous studies highlighting the importance of the mantle lithosphere in tectonic processes (e.g., Balázs et al., 2018;Babuška & Plomerová, 2013;Pysklywec & Beaumont, 2004;Heron et al., 2019;Heron et al., 2015;Hopper & Fischer, 2015).
The plate motion, which we apply here as a boundary condition, is an important part of the history of the region. Figure 4 shows the relative velocities and orientation of the plate motion over the course of the rift . In our modeling, we have fixed the extension velocity and orientation for 15 Myr in order to approximate Phase 1 of the rift history ( Figure 4). Our modeled Davis Strait region is susceptible to rifting and indeed thins throughout the simulation, which is in keeping with geophysical interpretation of the region (Figure 9b; Funck et al., 2007;Suckro et al., 2013). If we allow our reference case Model M1 to deform for longer than 15 Myr, the modeled Davis Strait thins further before joining up to the north and south spreading zones after 19 Myr ( Figure S5).
Due to numerical complexity and computational expense, it is difficult to apply a time-dependent extension velocity covering the whole rift sequence (Figure 4). However, the extension velocity and orientation used here fall within the estimation for Phase 1. As outlined in  and shown in Figure 1b), the four-point checklist for the rift evolution of the region has already been satisfied at the end of Phase 1 (60 Ma, 15-20 Myr after extension is initiated). The rotation of the extension axis to approximately north-south in Phase 2 ( Figure 4) has an impact on the fault orientation and kinematics but not on the overall geometry of breakup . As a result, modeling only Phase 1 (e.g., 1 cm/year at 15 Myr) is appropriate for our study.

Parameter Analysis
In testing the robustness of our study, we explored the parameter space surrounding these 3-D numerical models of extension of continental lithosphere finding that the choice of rheological parameters is important to the development of appropriate Davis Strait tectonics. Schiffer et al. (2016) interpret mantle lithosphere scarring on the continental margin of East Greenland to be of higher density than the surrounding mantle material, with Petersen and Schiffer (2016) providing modeling on the topic. In our study, through changing our mantle lithosphere scar from an area of weakness to being stronger than the surrounding material, we were unable to produce any focusing of strain that would allow a Davis Strait-type geometry rift to develop. However, a number of studies have discussed the weakening impact of tectonic processes on the lithosphere to facilitate continental rifting (Dunbar & Sawyer, 1988. The subduction of crustal material into the mantle through ancient processes could increase volatiles to the lower lithosphere, weakening the seismically imaged scarred material (Petersen & Schiffer, 2016;Pollack, 1986).
We also studied the strain range over which material is weakened (e.g., Figure S3). In Models M2-M5 we used Model M1 setup and changed the strain range for weakening to different values used in recent studies (e.g., Huismans & Beaumont, 2011;Naliboff & Buiter, 2015;Brune et al., 2013;Salazar-Mora et al., 2018). We found some differences between the results with regards to the evolution of the rift (e.g., Figures S3 and  S4); however, they all satisfied the required four-point checklist for Davis Strait tectonics ( Figure S3). Allken et al. (2012) showed the amount of strain weakening is an important factor in controlling the mode of rift interaction in brittle-ductile coupled crustal systems. Although the parameters used in the main manuscript are in keeping with the rest of the community (e.g., Brune et al., 2017), the work presented here (alongside the previous work of Allken et al., 2012) highlights the difficulty of modeling strain weakening due to the unconstrained nature of the values for different rheologies.

Conclusions
For the first time, numerical simulations show that rifting of lithosphere with a preexisting mantle structure can reproduce first-order features that resemble the Labrador Sea, Davis Strait, Baffin Bay continental margins, and ocean basins ( Figure 6). The results offer a new mechanism for rifting in the region, focusing on the role of ancient mantle lithosphere suturing rather than or in addition to crustal inheritance ( Figure 5). The obliquity of the suture to the extension direction is important for the tectonic evolution of the region and generates a segmented rift pattern (Figure 7). This study supplements a growing body of work that is posing questions on the fundamentals of inheritance and shows that we should be looking deeper than the Moho for controls on the tectonic style of lithosphere-scale deformation.