# Toward the calibration of 2D thermomechanical simulations of magma poor passive continental margins: method, validation and case example

Paul Perron, Laetitia Le Pourhiet, Anthony Jourdon, Tristan Cornu, Claude Gout

<figure><img src="image_1.png" /></figure>

## To cite this version:

Paul Perron, Laetitia Le Pourhiet, Anthony Jourdon, Tristan Cornu, Claude Gout. Toward the calibration of 2D thermomechanical simulations of magma poor passive continental margins: method, validation and case example. Comptes Rendus. Géoscience, 2024, 356 (S2), pp.1 - 22. 10.5802/cr-geos.258. hal-04587010

HAL Id: hal-04587010

https://hal.science/hal-04587010

Submitted on 24 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Toward the calibration of 2D thermomechanical simulations of magma poor passive continental margins: method, validation and case example

Paul Perron, Laetitia Le Pourhiet, Anthony Jourdon, Tristan Cornu, Claude Gout

<figure><img src="image_1.png" /></figure>

## To cite this version:

Paul Perron, Laetitia Le Pourhiet, Anthony Jourdon, Tristan Cornu, Claude Gout. Toward the calibration of 2D thermomechanical simulations of magma poor passive continental margins: method, validation and case example. Comptes Rendus. Géoscience, 2024, 356 (S2), pp.1 - 22. 10.5802/cr-geos.258. hal-04587010

HAL Id: hal-04587010

https://hal.science/hal-04587010v1

Submitted on 24 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Paul Perron, Laetitia Le Pourhiet, Anthony Jourdon, Tristan Cornu and Claude Gout

**Toward the calibration of 2D thermomechanical simulations of magma poor passive continental margins: method, validation and case example**

Published online: 22 May 2024

**Part of Special Issue:** Geodynamics of Continents and Oceans – A tribute to Jean Aubouin

**Guest editors:** Olivier Fabbri (Université de Franche-Comté, UMR CNRS 6249, Besançon), Michel Faure (Université d'Orléans-BRGM, UMR CNRS 7325, Institut des Sciences de la Terre, Orléans), Jacky Ferrière (Université de Lille, faculté des Sciences), Laurent Jolivet (Sorbonne Université, ISTeP, UMR 7193, Paris) and Sylvie Leroy (Sorbonne Université, CNRS-INSU, ISTeP, Paris)

https://doi.org/10.5802/crgeos.258

<figure><img src="image_3.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Paul Perron, Laetitia Le Pourhiet, Anthony Jourdon, Tristan Cornu and Claude Gout

**Toward the calibration of 2D thermomechanical simulations of magma poor passive continental margins: method, validation and case example**

Published online: 22 May 2024

**Part of Special Issue:** Geodynamics of Continents and Oceans – A tribute to Jean Aubouin

**Guest editors:** Olivier Fabbri (Université de Franche-Comté, UMR CNRS 6249, Besançon), Michel Faure (Université d'Orléans-BRGM, UMR CNRS 7325, Institut des Sciences de la Terre, Orléans), Jacky Ferrière (Université de Lille, faculté des Sciences), Laurent Jolivet (Sorbonne Université, ISTeP, UMR 7193, Paris) and Sylvie Leroy (Sorbonne Université, CNRS-INSU, ISTeP, Paris)

https://doi.org/10.5802/crgeos.258

<figure><img src="image_3.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/
Research article

Geodynamics of Continents and Oceans – A tribute to Jean Aubouin

# Toward the calibration of 2D thermomechanical simulations of magma poor passive continental margins: method, validation and case example

Paul Perron<sup>a, *</sup>, Laetitia Le Pourhiet<sup>a, *</sup>, Anthony Jourdon<sup>b, *</sup>, Tristan Cornu<sup>b</sup> and Claude Gout<sup>b</sup>

<sup>a</sup> Sorbonne Université, CNRS-INSU, Institut des Sciences de la Terre Paris, ISTE P UMR 7193, F-75005 Paris, France

<sup>b</sup> TotalEnergies, CSTJF, Avenue Larribau, 64000, Pau, France

URL: https://www.idref.fr/086801902

Current address: Ludwig-Maximilian University, Department of Earth and Environmental Sciences, Munich, Germany (A. Jourdon)

E-mail: laetitia.le_pourhiet@ sorbonne-universite.fr (L. Le Pourhiet)

**Abstract.** In recent decades, geodynamic modelers have aimed to comprehend key factors governing continental rifting, such as the extension rate, lithospheric thickness, Moho mechanical coupling, and mantle convection's thermal influence. While prior models offered insights into rifting processes, they lacked the calibration to specific Earth regions. Introducing heterogeneities into the model does in some cases help to calibrate the simulation results to a geological data from a specific region. Acknowledging structural inheritance as a form of kinematic forcing in the models, and recognizing the challenge of anticipating and identifying all inherited geological structures present before rifting, a new modeling approach was devised. This method integrates a new kinematic module into the pTatin2D code, allowing for calibrating numerical simulations with regional geological and geophysical dataset over time while solving for mechanical balance using Stokes flow to ensure that crustal deformation remains consistent with mantle dynamics. By calibrating against a 2D cross-section extracted from the final state of a 3D model, we show that the approach predicts thermal history and deformation paths beyond calibration points. In particular, the thermo-mechanical feedback can help mitigate some uncertainties in the deformation path. Applied to Iberia–Newfoundland margins, the method demonstrates effectiveness in real-case scenarios, aligning with previous reconstructions by incorporating faults and lower crustal flow.

**Keywords.** Thermomechanical modelling, Rifted margins, Kinematic reconstruction.

Manuscript received 3 August 2023, revised 3 April 2024, accepted 9 April 2024.

## 1. State of the art

Thermo-mechanically coupled simulations of the lithosphere and asthenosphere tectonic evolution

are of both academic and industrial interest, as they offer a physics-based understanding of the tectono-thermal evolution and sedimentary infilling of basins. However, reconstructing the polyphase stretching and thinning modes during the history of rifting, leading to continental breakup and seafloor spreading, remains challenging due to

<footnote>* Corresponding author.</footnote>
Research article

Geodynamics of Continents and Oceans – A tribute to Jean Aubouin

# Toward the calibration of 2D thermomechanical simulations of magma poor passive continental margins: method, validation and case example

Paul Perron<sup>a, *</sup>, Laetitia Le Pourhiet<sup>a, *</sup>, Anthony Jourdon<sup>b, *</sup>, Tristan Cornu<sup>b</sup> and Claude Gout<sup>b</sup>

<sup>a</sup> Sorbonne Université, CNRS-INSU, Institut des Sciences de la Terre Paris, ISTE P UMR 7193, F-75005 Paris, France

<sup>b</sup> TotalEnergies, CSTJF, Avenue Larribau, 64000, Pau, France

URL: https://www.idref.fr/086801902

Current address: Ludwig-Maximilian University, Department of Earth and Environmental Sciences, Munich, Germany (A. Jourdon)

E-mail: laetitia.le_pourhiet@ sorbonne-universite.fr (L. Le Pourhiet)

**Abstract.** In recent decades, geodynamic modelers have aimed to comprehend key factors governing continental rifting, such as the extension rate, lithospheric thickness, Moho mechanical coupling, and mantle convection's thermal influence. While prior models offered insights into rifting processes, they lacked the calibration to specific Earth regions. Introducing heterogeneities into the model does in some cases help to calibrate the simulation results to a geological data from a specific region. Acknowledging structural inheritance as a form of kinematic forcing in the models, and recognizing the challenge of anticipating and identifying all inherited geological structures present before rifting, a new modeling approach was devised. This method integrates a new kinematic module into the pTatin2D code, allowing for calibrating numerical simulations with regional geological and geophysical dataset over time while solving for mechanical balance using Stokes flow to ensure that crustal deformation remains consistent with mantle dynamics. By calibrating against a 2D cross-section extracted from the final state of a 3D model, we show that the approach predicts thermal history and deformation paths beyond calibration points. In particular, the thermo-mechanical feedback can help mitigate some uncertainties in the deformation path. Applied to Iberia–Newfoundland margins, the method demonstrates effectiveness in real-case scenarios, aligning with previous reconstructions by incorporating faults and lower crustal flow.

**Keywords.** Thermomechanical modelling, Rifted margins, Kinematic reconstruction.

Manuscript received 3 August 2023, revised 3 April 2024, accepted 9 April 2024.

## 1. State of the art

Thermo-mechanically coupled simulations of the lithosphere and asthenosphere tectonic evolution

are of both academic and industrial interest, as they offer a physics-based understanding of the tectono-thermal evolution and sedimentary infilling of basins. However, reconstructing the polyphase stretching and thinning modes during the history of rifting, leading to continental breakup and seafloor spreading, remains challenging due to

<footnote>* Corresponding author.</footnote>
several geological uncertainties, such as inheritance (rheology, tectonic structures, magmatism...), stress direction, and strain rate. These last two pieces of information are constrained by global plate tectonic reconstructions [e.g., Seton et al., 2012]. Yet, the lack of magnetic anomalies before seafloor spreading [Vine and Matthews, 1963] requires indirect quantification of the deformation using approximate restoration methods [e.g., Crosby et al., 2008, Sutra et al., 2013].

Passive margins result from the process of thinning and breakup of the continental lithosphere leading to the formation of new oceanic lithosphere [McKenzie, 1978]. On Earth, conjugated passive margins present various morphologies [e.g., Sapin et al., 2021]: They can be symmetrical as in the case of Pelotas and Walvis margins [Blaich et al., 2011] or asymmetrical and hyperextended in the case of Iberia–Newfoundland margins [Sutra et al., 2013].

Many dynamic numerical models in the literature attempt to explain the variability of structural architecture of passive margins, either in two dimensions [e.g., Buck, 1991, Brune et al., 2014, Gueydan and Précigout, 2014, Huismans and Beaumont, 2014, Dias et al., 2015, Brune et al., 2017, Tetreault and Buiter, 2018] or in three dimensions [e.g., Brune and Autin, 2013, Heine and Brune, 2014, Le Pourhiet et al., 2017, 2018, Gouiza and Naliboff, 2021, Jourdon et al., 2020, 2021, Neuharth et al., 2021]. These models present “general concepts” so as to comprehend the evolution of rifted margins and none aims to be directly applicable to real examples of passive margins. Indeed, passive continental margins evolve self-consistently with thermo-rheological laws of continental crust and mantle boundary and initial conditions. Radically different geometries of basins can be obtained simply by varying some of the input parameters, notably the rate of extension, the initial lithospheric thickness, mechanical coupling or decoupling at the Moho, and the thermal contribution of mantle convection. Moreover, a lot of complexity can emerge from running simulations in 3D, due in particular to (1) the obliquity between the strike of the mantle neck and that of the crustal structures [e.g., Brune and Autin, 2013, Heine and Brune, 2014, Liao and Gerya, 2015, Le Pourhiet et al., 2017, Ammann et al., 2018, Jourdon et al., 2021, Neuharth et al., 2021], or (2) the rate of propagation of the mantle neck [e.g., Le Pourhiet et al., 2018, Jourdon et al., 2020].

Although the influence of inheritance on the architecture and tectonic evolution of passive margins is well accepted in the community [Taylor et al., 1999, Manatschal et al., 2015, Brune et al., 2017], highlighting it through forward thermomechanical numerical modelling remains difficult, particularly during the initiation and evolution of rifting and lithospheric breakup. In most cases, the approach consists in trying different hypotheses, i.e. testing several simulations with initial imposed inherited “anomalies” such as faults, plutons and/or cratons in order to reproduce the present-day configuration of passive margins consistent with the geological observations [Chenin and Beaumont, 2013, Manatschal et al., 2015, Duretz et al., 2016, Balázs et al., 2018, Wenker and Beaumont, 2018, Chenin et al., 2019, Petri et al., 2019, Gouiza and Naliboff, 2021]. This modelling workflow remains fastidious and does not necessarily lead to the reconstruction of a particular passive margin architecture which can be quantitatively calibrated with data.

Alternatively, the kinematic modelling approach as favoured in Jeanniot et al. [2016] permits the model predictions of a particular continental margin being quantitatively calibrated against geological and geophysical observations. Utilizing a linear Stokes solver, this fast approach combines the instant flow field with a consistently fixed area of pure shear deformation (representing rift location, width, and extension rate) and a constant buoyancy flux at each phase of deformation (10 or 12 phases are needed). Material and heat are then transported using the calculated flow field considering small time steps that respect CFL condition.

However, despite the adaptability of Jeanniot et al. [2016] parameterization allowing calibration for complex passive margins, these simulations lack essential rheological feedback. The absence of complete thermomechanical coupling leads to mantle dynamics that are inconsistent with the ongoing crustal thinning annihilating natural transition from passive to active rifting during necking. Moreover, pure shear deformation disregards faulting and flexural deformation of the lithosphere during the early stages of continental rifting. This leaves a lot of free independent parameters, which could be constrained by introducing a physical model. This possibility of overfitting impairs the predictive power of the preferred model outside of the calibration points.
several geological uncertainties, such as inheritance (rheology, tectonic structures, magmatism...), stress direction, and strain rate. These last two pieces of information are constrained by global plate tectonic reconstructions [e.g., Seton et al., 2012]. Yet, the lack of magnetic anomalies before seafloor spreading [Vine and Matthews, 1963] requires indirect quantification of the deformation using approximate restoration methods [e.g., Crosby et al., 2008, Sutra et al., 2013].

Passive margins result from the process of thinning and breakup of the continental lithosphere leading to the formation of new oceanic lithosphere [McKenzie, 1978]. On Earth, conjugated passive margins present various morphologies [e.g., Sapin et al., 2021]: They can be symmetrical as in the case of Pelotas and Walvis margins [Blaich et al., 2011] or asymmetrical and hyperextended in the case of Iberia–Newfoundland margins [Sutra et al., 2013].

Many dynamic numerical models in the literature attempt to explain the variability of structural architecture of passive margins, either in two dimensions [e.g., Buck, 1991, Brune et al., 2014, Gueydan and Précigout, 2014, Huismans and Beaumont, 2014, Dias et al., 2015, Brune et al., 2017, Tetreault and Buiter, 2018] or in three dimensions [e.g., Brune and Autin, 2013, Heine and Brune, 2014, Le Pourhiet et al., 2017, 2018, Gouiza and Naliboff, 2021, Jourdon et al., 2020, 2021, Neuharth et al., 2021]. These models present “general concepts” so as to comprehend the evolution of rifted margins and none aims to be directly applicable to real examples of passive margins. Indeed, passive continental margins evolve self-consistently with thermo-rheological laws of continental crust and mantle boundary and initial conditions. Radically different geometries of basins can be obtained simply by varying some of the input parameters, notably the rate of extension, the initial lithospheric thickness, mechanical coupling or decoupling at the Moho, and the thermal contribution of mantle convection. Moreover, a lot of complexity can emerge from running simulations in 3D, due in particular to (1) the obliquity between the strike of the mantle neck and that of the crustal structures [e.g., Brune and Autin, 2013, Heine and Brune, 2014, Liao and Gerya, 2015, Le Pourhiet et al., 2017, Ammann et al., 2018, Jourdon et al., 2021, Neuharth et al., 2021], or (2) the rate of propagation of the mantle neck [e.g., Le Pourhiet et al., 2018, Jourdon et al., 2020].

Although the influence of inheritance on the architecture and tectonic evolution of passive margins is well accepted in the community [Taylor et al., 1999, Manatschal et al., 2015, Brune et al., 2017], highlighting it through forward thermomechanical numerical modelling remains difficult, particularly during the initiation and evolution of rifting and lithospheric breakup. In most cases, the approach consists in trying different hypotheses, i.e. testing several simulations with initial imposed inherited “anomalies” such as faults, plutons and/or cratons in order to reproduce the present-day configuration of passive margins consistent with the geological observations [Chenin and Beaumont, 2013, Manatschal et al., 2015, Duretz et al., 2016, Balázs et al., 2018, Wenker and Beaumont, 2018, Chenin et al., 2019, Petri et al., 2019, Gouiza and Naliboff, 2021]. This modelling workflow remains fastidious and does not necessarily lead to the reconstruction of a particular passive margin architecture which can be quantitatively calibrated with data.

Alternatively, the kinematic modelling approach as favoured in Jeanniot et al. [2016] permits the model predictions of a particular continental margin being quantitatively calibrated against geological and geophysical observations. Utilizing a linear Stokes solver, this fast approach combines the instant flow field with a consistently fixed area of pure shear deformation (representing rift location, width, and extension rate) and a constant buoyancy flux at each phase of deformation (10 or 12 phases are needed). Material and heat are then transported using the calculated flow field considering small time steps that respect CFL condition.

However, despite the adaptability of Jeanniot et al. [2016] parameterization allowing calibration for complex passive margins, these simulations lack essential rheological feedback. The absence of complete thermomechanical coupling leads to mantle dynamics that are inconsistent with the ongoing crustal thinning annihilating natural transition from passive to active rifting during necking. Moreover, pure shear deformation disregards faulting and flexural deformation of the lithosphere during the early stages of continental rifting. This leaves a lot of free independent parameters, which could be constrained by introducing a physical model. This possibility of overfitting impairs the predictive power of the preferred model outside of the calibration points.
## 2. Why calibrating a thermomechanical simulation?

This study aims at developing an alternative approach allying kinematics calibration and physically consistent thermo-mechanics evolution. Combining these two approaches increases the predictive power of the calibrated models in terms of tectono-thermal history of sedimentary basins of passive continental margins. More precisely, the objective is to get thermal feed-back from the simulations of lithosphere-asthenosphere boundaries (LAB) through time by the calibration of the crustal scale features (basins, faults and the Moho depth and geometry) known from geophysical and geophysical data (e.g., 2D seismic profiles). The overall goal is therefore to be able to produce complex real passive margins formed through multiphase continental rifting and break-up that can be calibrated to subsurface data, in particular seismic data at a minimal numerical cost.

The most important hypothesis behind our work is that, while we admit that structural inheritance plays an important role in introducing complexity in the system, we do not want to include structural inheritance explicitly in the models because the initial geometry is even less known than the final one.

In other words, we posit that mechanical heterogeneities, present in the crust at the onset of rifting, or 3D effects, like propagating ridge or strike-slip system, can be replaced by adequate kinematic forcing. By kinematic forcing, we mean that we are imposing “boundary conditions”, i.e. Dirichlet conditions, within the model domain to turn on and off upper crustal extension when and where the data justifies it. Compared to a regular self-consistent thermo-mechanical approach, we therefore intend to limit the complexity embedded in the initial conditions (see Section 3.1) and commit the “academic scientific crime” to enforce the complexity using extremely rigid boundary conditions within the model domain (see Section 3.2). Compared to a regular kinematic approach, the small domains that can deform freely in the crust do not follow pure-shear flattening but may be faulted at smaller length scale according to the thermo-rheological profile of the lithosphere. Most importantly, the deformation of the LAB during rifting, which cannot really be constrained by seismic data on passive margin, is completely self-consistent with the better constrained deformation imposed

in the crust. This allows for better constraining the mantle contribution to the tectono-thermal history of basins of continental lithosphere during rifting. To perform the study, we use the thermo-mechanical numerical code pTatin2D, which solves for non-linear variable viscosity and density Stokes equation coupled with conservation of heat including advection, diffusion and production, and with simple diffusion-based surface evolution model. Appendix A provides details on the partial differential equations involved.

In order to test the relevance of coupling self-consistent thermo-mechanical modelling approach with a kinematic modelling approach, a new kinematic module has been implemented into the two-dimensional thermomechanical code pTatin2D, allowing to force the location and range of extensional deformation over time in the crust. First of all, sensitivity tests are performed to examine the influence of each parameter that is involved in the kinematic calibration of the new module in pTatin2D. We then test the predictivity of this 2D approach by benchmarking it with a cross-section extracted from a complex 3D model. Finally, we present the application of the new kinematic module to the well-constrained geological case study of Iberia–Newfoundland conjugate rifted margins.

## 3. The kinematic module

### 3.1. *Model setup: rheological parameters and initial conditions*

Our experiments represent a 1200 km (in the x direction) × 250 km (in the z direction) domain discretized with 90 × 40 elements, respectively (Figure 1a). This low-resolution together with mesh refinement towards the surface permits fast computing while capturing large-scale mantle dynamics and crustal-scale tectonics. The model domain is kept to a minimal complexity. It is only subdivided in two rheological layers: the crust simulated by a diorite flow law [Carter and Tsenn, 1987] and the mantle simulated using olivine flow law [Goetze and Evans, 1979]. The initial Moho depth is calibrated by the user. The initial Lithosphere Asthenosphere Boundary (LAB), that is the 1300 °C isotherm, is initially lying flat at 120 km depth (Figure 1a). The full rheological parameters are compiled in Table 1 and its caption.
**2. Why calibrating a thermomechanical simulation?**

This study aims at developing an alternative approach allying kinematics calibration and physically consistent thermo-mechanics evolution. Combining these two approaches increases the predictive power of the calibrated models in terms of tectono-thermal history of sedimentary basins of passive continental margins. More precisely, the objective is to get thermal feed-back from the simulations of lithosphere-asthenosphere boundaries (LAB) through time by the calibration of the crustal scale features (basins, faults and the Moho depth and geometry) known from geological and geophysical data (e.g., 2D seismic profiles). The overall goal is therefore to be able to produce complex real passive margins formed through multiphase continental rifting and break-up that can be calibrated to subsurface data, in particular seismic data at a minimal numerical cost.

The most important hypothesis behind our work is that, while we admit that structural inheritance plays an important role in introducing complexity in the system, we do not want to include structural inheritance explicitly in the models because the initial geometry is even less known than the final one.

In other words, we posit that mechanical heterogeneities, present in the crust at the onset of rifting, or 3D effects, like propagating ridge or strike-slip system, can be replaced by adequate kinematic forcing. By kinematic forcing, we mean that we are imposing “boundary conditions”, i.e. Dirichlet conditions, within the model domain to turn on and off upper crustal extension when and where the data justifies it. Compared to a regular self-consistent thermo-mechanical approach, we therefore intend to limit the complexity embedded in the initial conditions (see Section 3.1) and commit the “academic scientific crime” to enforce the complexity using extremely rigid boundary conditions within the model domain (see Section 3.2). Compared to a regular kinematic approach, the small domains that can deform freely in the crust do not follow pure-shear flattening but may be faulted at smaller length scale according to the thermo-rheological profile of the lithosphere. Most importantly, the deformation of the LAB during rifting, which cannot really be constrained by seismic data on passive margin, is completely self-consistent with the better constrained deformation imposed

in the crust. This allows for better constraining the mantle contribution to the tectono-thermal history of basins of continental lithosphere during rifting. To perform the study, we use the thermo-mechanical numerical code pTatin2D, which solves for non-linear variable viscosity and density Stokes equation coupled with conservation of heat including advection, diffusion and production, and with simple diffusion-based surface evolution model. Appendix A provides details on the partial differential equations involved.

In order to test the relevance of coupling self-consistent thermo-mechanical modelling approach with a kinematic modelling approach, a new kinematic module has been implemented into the two-dimensional thermomechanical code pTatin2D, allowing to force the location and range of extensional deformation over time in the crust. First of all, sensitivity tests are performed to examine the influence of each parameter that is involved in the kinematic calibration of the new module in pTatin2D. We then test the predictivity of this 2D approach by benchmarking it with a cross-section extracted from a complex 3D model. Finally, we present the application of the new kinematic module to the well-constrained geological case study of Iberia–Newfoundland conjugate rifted margins.

**3. The kinematic module**

**3.1. Model setup: rheological parameters and initial conditions**

Our experiments represent a 1200 km (in the x direction) × 250 km (in the z direction) domain discretized with 90 × 40 elements, respectively (Figure 1a). This low-resolution together with mesh refinement towards the surface permits fast computing while capturing large-scale mantle dynamics and crustal-scale tectonics. The model domain is kept to a minimal complexity. It is only subdivided in two rheological layers: the crust simulated by a diorite flow law [Carter and Tsenn, 1987] and the mantle simulated using olivine flow law [Goetze and Evans, 1979]. The initial Moho depth is calibrated by the user. The initial Lithosphere Asthenosphere Boundary (LAB), that is the 1300 °C isotherm, is initially lying flat at 120 km depth (Figure 1a). The full rheological parameters are compiled in Table 1 and its caption.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Initial model setup for calibration. (a) Geometry of the model, with dimensions and boundary conditions for all simulations. (b) Initial temperature state of the lithosphere for all the simulations presented. (c) The kinematic module presenting the different inputs parameters **W**, W₀, X₀, vₓR, vₓL and **zvel** (see part 3.2 for full description). (d) Yield-strength envelopes (Mohr-Coulomb representation) of the simulated materials calculated for a strain rate of 10⁻¹⁴ s⁻¹ and a stress limited at 400 MPa [value from Watremez et al., 2013]. Δσ represents the second invariant of deviatoric stress.</figcaption></figure>

<table><caption>**Table 1.** Rheological parameters used for the Dislocation Creep Arrhenius' Laws of each simulated rheologies</caption><thead><tr><th></th><th>Lithology</th><th>Sediments (Quartz)</th><th>Crust (Diorite)</th><th>Mantle/Asthenosphere (Olivine)</th><th>Units</th></tr></thead><tbody><tr><td>ρ<sub>0</sub></td><td>Density</td><td>2000</td><td>2850</td><td>3300</td><td>kg/m<sup>3</sup></td></tr><tr><td colspan="6">Creep parameters</td></tr><tr><td>n</td><td>Stress exponent</td><td>3</td><td>3.05</td><td>3</td><td></td></tr><tr><td>A</td><td>Pre-exp. factor</td><td>6.80 × 10<sup>-6</sup></td><td>6.3 × 10<sup>-2</sup></td><td>7.0 × 10<sup>3</sup></td><td>MPa<sup>-n</sup>·s<sup>-1</sup></td></tr><tr><td>Q</td><td>Activation energy</td><td>156</td><td>275</td><td>510</td><td>kJ·mol<sup>-1</sup></td></tr><tr><td colspan="6">Brittle parameters Intact “0” and fully damaged “∞”</td></tr><tr><td>φ<sub>0</sub></td><td>Friction</td><td>20</td><td>30</td><td>30</td><td>°</td></tr><tr><td>φ<sub>∞</sub></td><td>Friction</td><td>20</td><td>30</td><td>10</td><td>°</td></tr><tr><td>Co<sub>0</sub></td><td>Cohesion</td><td>2</td><td>20</td><td>20</td><td>MPa</td></tr><tr><td>Co<sub>∞</sub></td><td>Cohesion</td><td>1</td><td>10</td><td>20</td><td>MPa</td></tr><tr><td>Ref.</td><td></td><td>R & M</td><td>C & T</td><td>G & E</td><td></td></tr></tbody></table>

R & M stands for Ranalli and Murphy [1987], C & T for Carter and Tsenn [1987] and G & E for Goetze and Evans [1979]. All the parameters are described in Appendix A and the parameters that are not listed in Table 1 have fixed values for all lithologies: α = 3 × 10⁻⁵ K⁻¹ (thermal expansion), β = 1 × 10⁻¹¹ Pa⁻¹ (adiabatic compressibility), κ = 5 × 10⁻⁷ m²·s⁻¹ (erosional diffusive coefficient), κ = 1 × 10⁻⁶ m²·s⁻¹ (heat diffusivity), εmin = 0, εmax = 1, Cp = 1000 J·kg⁻¹·K⁻¹ (heat capacity), H = 3.33 × 10⁻¹⁰ W·kg⁻¹ (heat production).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Initial model setup for calibration. (a) Geometry of the model, with dimensions and boundary conditions for all simulations. (b) Initial temperature state of the lithosphere for all the simulations presented. (c) The kinematic module presenting the different inputs parameters **W**, $W_0$, $X_0$, $v_{xR}$, $v_{xL}$ and **zvel** (see part 3.2 for full description). (d) Yield-strength envelopes (Mohr-Coulomb representation) of the simulated materials calculated for a strain rate of $10^{-14}$ s$^{-1}$ and a stress limited at 400 MPa [value from Watremez et al., 2013]. $\Delta\sigma$ represents the second invariant of deviatoric stress.</figcaption></figure>

<table><caption>**Table 1.** Rheological parameters used for the Dislocation Creep Arrhenius' Laws of each simulated rheologies</caption><thead><tr><th></th><th>Lithology</th><th>Sediments (Quartz)</th><th>Crust (Diorite)</th><th>Mantle/Asthenosphere (Olivine)</th><th>Units</th></tr></thead><tbody><tr><td>ρ<sub>0</sub></td><td>Density</td><td>2000</td><td>2850</td><td>3300</td><td>kg/m<sup>3</sup></td></tr><tr><td colspan="6">Creep parameters</td></tr><tr><td>n</td><td>Stress exponent</td><td>3</td><td>3.05</td><td>3</td><td></td></tr><tr><td>A</td><td>Pre-exp. factor</td><td>6.80 × 10<sup>-6</sup></td><td>6.3 × 10<sup>-2</sup></td><td>7.0 × 10<sup>3</sup></td><td>MPa<sup>-n</sup>·s<sup>-1</sup></td></tr><tr><td>Q</td><td>Activation energy</td><td>156</td><td>275</td><td>510</td><td>kJ·mol<sup>-1</sup></td></tr><tr><td colspan="6">Brittle parameters Intact “0” and fully damaged “∞”</td></tr><tr><td>φ<sub>0</sub></td><td>Friction</td><td>20</td><td>30</td><td>30</td><td>°</td></tr><tr><td>φ<sub>∞</sub></td><td>Friction</td><td>20</td><td>30</td><td>10</td><td>°</td></tr><tr><td>Co<sub>0</sub></td><td>Cohesion</td><td>2</td><td>20</td><td>20</td><td>MPa</td></tr><tr><td>Co<sub>∞</sub></td><td>Cohesion</td><td>1</td><td>10</td><td>20</td><td>MPa</td></tr><tr><td>Ref.</td><td></td><td>R & M</td><td>C & T</td><td>G & E</td><td></td></tr></tbody></table>

R & M stands for Ranalli and Murphy [1987], C & T for Carter and Tsenn [1987] and G & E for Goetze and Evans [1979]. All the parameters are described in Appendix A and the parameters that are not listed in Table 1 have fixed values for all lithologies: $\alpha = 3 \times 10^{-5}$ K$^{-1}$ (thermal expansion), $\beta = 1 \times 10^{-11}$ Pa$^{-1}$ (adiabatic compressibility), $\kappa = 5 \times 10^{-7}$ m$^2$·s$^{-1}$ (erosional diffusive coefficient), $\kappa = 1 \times 10^{-6}$ m$^2$·s$^{-1}$ (heat diffusivity), $\epsilon_{min} = 0$, $\epsilon_{max} = 1$, $Cp = 1000$ J·kg$^{-1}$·K$^{-1}$ (heat capacity), $H = 3.33 \times 10^{-10}$ W·kg$^{-1}$·s$^{-1}$ (heat production).
The model upper boundary behaves as a free-surface subjected to erosion and sedimentation, which allows the development of sedimentary basins. These surface processes are modelled using Culling [1965]'s law with a constant diffusivity ($\kappa_e$) of $10^{-6}$ m²·s⁻¹ accounting for local erosion/deposition processes. We do not introduce a source term, which can be used to account for local carbonates production or out-of-plane sediment infill [e.g. Jourdon et al., 2018b, Perron et al., 2021]. The basal boundary is prescribed an inflow to compensate for the outflow imposed by the extensional velocity varying along the lateral boundaries. The temporal variation and implementation of these boundary conditions are detailed in the next section.

The thermal boundary conditions are fixed at the top and base of the model to 0 °C and 1400 °C respectively, and a null heat flux is assumed at the model lateral boundaries (Figure 1b). A radiogenic heat production is only imposed in the crust with a value of $3.33 \times 10^{-10}$ W·kg⁻¹, we neglect it in the mantle.

## 3.2. Boundary conditions

The kinematic module allows the user to impose regions where deformation can occur or, on the contrary, regions acting as rigid blocks. These choices are arbitrarily based on previously collected geological evidence. This approach is particularly useful in the case of abandoned rifts due to multiphase rifting, out-of-plane source of localisation like a propagating rift segments [Le Pourhiet et al., 2017, Neuharth et al., 2021] or strike slip system juxtaposing two margins segments with different histories [Jourdon et al., 2020, 2021]. For this type of problem, one could argue that 3D simulations would be more appropriate, but data on crustal scale, which encompass the whole margin are usually 2D seismic lines rather than 3D seismic blocks. Moreover, running 3D simulations would defeat the relatively low computational cost of the new approach we want to elaborate.

This module, newly implemented in pTatin2D, is largely inspired by the work of Jeanniot et al. [2016] which imposed different delimited areas deformed by pure shear at lithospheric scale during each given lapse of time of the simulation. Here, instead of imposing pure shear, we let these parts of the crustal domain free to deform according to their rheological parameters and the current thermal field. In this

zone of mechanically consistent deformation, our simulation can produce strain localisation at a smaller scale than the one imposed by the user. Therefore, the kinematic module imposes a rigid behaviour down to a depth **zvel** outside of a domain of width $W(t)$.

The kinematic parameter **zvel**, corresponds to the fixed thickness of the lithosphere where the variable velocities $v_{xL}$ and $v_{xR}$ are imposed (Figure 1c). This boundary condition is similar to that proposed by Gorczyk et al. [2007]. We investigate the impact of **zvel** on the evolution of the rifted margin architecture, by varying its thickness from 0 to 40 km (Figure 2). We find that driving the model by only fixing the velocity on the surface nodes ($zvel = 0$) is not satisfactory because the very shallow crust and sediments have not a sufficient frictional strength to transmit the deformation to the whole crust. It results in the formation of large artefactual shear strain rate at the top of the simulation and very little drag is transmitted to the crust (Figure 2). When **zvel** includes at least 5 km of crust, kinematic forcing works fine and further increase of **zvel** show negligible impact on the architecture of the rifted margins.

Each phase of deformation $P^i$ is defined by a start time $t_i$ and an end time $t_e$, and a set of user-defined kinematic parameters. The kinematic parameters define the geometry of the rigid domains at the beginning of the phase $P^i$ and the velocity applied at the limit during that phase (Figure 1c). In details the user must provide:

* **W** is a user input and corresponds to the initial width of the delimited zone of the lithosphere, where the deformation is not imposed by user down to depth **zvel**. $W_0$ used in the computation represents the half of **W**;

* $X_0$ represents the localisation of the centre of the domain **W**.

* $v_{xL}$ is the velocity imposed on all the points of the domain located above **zvel** which abscisa is

$$x \leq X_0 - W_0 + v_{xL}(t - t_i). \quad (1a)$$

* $v_{xR}$ is the velocity imposed on all the points of the domain located above **zvel** which abscisa is

$$x \geq X_0 + W_0 + v_{xR}(t - t_i). \quad (1b)$$
The model upper boundary behaves as a free-surface subjected to erosion and sedimentation, which allows the development of sedimentary basins. These surface processes are modelled using Culling [1965]'s law with a constant diffusivity ($\kappa_e$) of $10^{-6}$ m$^2$.s$^{-1}$ accounting for local erosion/deposition processes. We do not introduce a source term, which can be used to account for local carbonates production or out-of-plane sediment infill [e.g. Jourdon et al., 2018b, Perron et al., 2021]. The basal boundary is prescribed an inflow to compensate for the outflow imposed by the extensional velocity varying along the lateral boundaries. The temporal variation and implementation of these boundary conditions are detailed in the next section.

The thermal boundary conditions are fixed at the top and base of the model to 0 °C and 1400 °C respectively, and a null heat flux is assumed at the model lateral boundaries (Figure 1b). A radiogenic heat production is only imposed in the crust with a value of $3.33 \times 10^{-10}$ W.kg$^{-1}$, we neglect it in the mantle.

## 3.2. Boundary conditions

The kinematic module allows the user to impose regions where deformation can occur or, on the contrary, regions acting as rigid blocks. These choices are arbitrarily based on previously collected geological evidence. This approach is particularly useful in the case of abandoned rifts due to multiphase rifting, out-of-plane source of localisation like a propagating rift segments [Le Pourhiet et al., 2017, Neuharth et al., 2021] or strike slip system juxtaposing two margins segments with different histories [Jourdon et al., 2020, 2021]. For this type of problem, one could argue that 3D simulations would be more appropriate, but data on crustal scale, which encompass the whole margin are usually 2D seismic lines rather than 3D seismic blocks. Moreover, running 3D simulations would defeat the relatively low computational cost of the new approach we want to elaborate.

This module, newly implemented in pTatin2D, is largely inspired by the work of Jeanniot et al. [2016] which imposed different delimited areas deformed by pure shear at lithospheric scale during each given lapse of time of the simulation. Here, instead of imposing pure shear, we let these parts of the crustal domain free to deform according to their rheological parameters and the current thermal field. In this

zone of mechanically consistent deformation, our simulation can produce strain localisation at a smaller scale than the one imposed by the user. Therefore, the kinematic module imposes a rigid behaviour down to a depth **zvel** outside of a domain of width $W(t)$.

The kinematic parameter **zvel**, corresponds to the fixed thickness of the lithosphere where the variable velocities $v_{xL}$ and $v_{xR}$ are imposed (Figure 1c). This boundary condition is similar to that proposed by Gorczyk et al. [2007]. We investigate the impact of **zvel** on the evolution of the rifted margin architecture, by varying its thickness from 0 to 40 km (Figure 2). We find that driving the model by only fixing the velocity on the surface nodes ($zvel = 0$) is not satisfactory because the very shallow crust and sediments have not a sufficient frictional strength to transmit the deformation to the whole crust. It results in the formation of large artefactual shear strain rate at the top of the simulation and very little drag is transmitted to the crust (Figure 2). When **zvel** includes at least 5 km of crust, kinematic forcing works fine and further increase of **zvel** show negligible impact on the architecture of the rifted margins.

Each phase of deformation $P^i$ is defined by a start time $t_i$ and an end time $t_e$, and a set of user-defined kinematic parameters. The kinematic parameters define the geometry of the rigid domains at the beginning of the phase $P^i$ and the velocity applied at the limit during that phase (Figure 1c). In details the user must provide:

* **W** is a user input and corresponds to the initial width of the delimited zone of the lithosphere, where the deformation is not imposed by user down to depth **zvel**. $W_0$ used in the computation represents the half of **W**;

* $X_0$ represents the localisation of the centre of the domain **W**.

* $v_{xL}$ is the velocity imposed on all the points of the domain located above **zvel** which abscisa is

$$x \leq X_0 - W_0 + v_{xL}(t - t_i). \quad (1a)$$

* $v_{xR}$ is the velocity imposed on all the points of the domain located above **zvel** which abscisa is

$$x \geq X_0 + W_0 + v_{xR}(t - t_i). \quad (1b)$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Influence of initial thickness of zvel kinematic parameter on the architecture of rifted margins. Sediments and plastic strain are shown on markers with their respective colormaps (stratigraphy and Plastic strain). The second invariant of strain rate is represented using contour in its respective colormap. Note a constant velocity $vx = 0.8$ cm/yr is applied and that the initial width of W is 200 km for all simulations.</figcaption></figure>

The deforming zone widens with time as

$$W(t) = L_x - (2W_0 + (v_{xL} + v_{xR})(t - t_i)), \quad (2)$$

with $L_x$ the full length of the simulation domain, to ensure that the structures formed consistently with mechanical equilibrium during one phase do not exit the computational domain. While the code can handle asymmetric stretching velocities, i.e. $v_{xL} \neq v_{xR}$, all the simulations presented herein use $v_{xL} = v_{xR}$ so that we introduce $v_x = (v_{xL} + v_{xR})$.

### 3.3. Influence of the strength of the crust and stretching rate using our parametrisation

For the reconstruction, it is interesting to get the initial basin spacing as a feature of the simulation using a large W, when possible because this will constrain quite well the rheological parameter. The spacing between normal faults is known to depend, at first order, on the thickness of the resistant layer of the lithosphere [e.g., Spadini and Podladchikov, 1996], that is, the depth to the brittle ductile transition. Using only one lithological layer for the whole crust, this depth is controlled by the viscous strength of the crust and the strain rate. In our approach, we use an ad hoc non-dimensional parameter $A_{scale}$ which

modifies the Arrhenius creep law to control the viscous strength of the crust. As the power law parameter, A, is in MPa$^{-n}$.s$^{-1}$ in Table 1 (n being the creep law exponent), a $A_{scale}$ value of $10^6$ produces the laboratory creep law, a lower value produces a lower viscosity. The mean stretching rate is controlled by the value of $vx/W$.

In most cases, using a large W (200 km, Figure 3) permits to calibrate $A_{scale}$ with the natural graben spacing expected in the early rifting phase if there is not too much influence of inheritance. As expected from all the literature cited in the introduction, weak crust leads to metamorphic core complexes, detachment systems or wide rift systems. Low extension rates lead to diffuse extension, while strong crust together with intermediate to high extension rates tend to lead to narrow rift systems. Figure 3 is intended as a catalogue for the user to pick the parameters to perform the initial reconstruction step. Supplementary Figure S1 provides a more comprehensive phase-map of possible initial rift geometry as a function of $A_{scale}$ and W than what could possibly be displayed in Figure 3. Figure 4 displays the outcome of an intermediate case with $A_{scale} = 10^6$ with W = 200 km and $v_{xR} + v_{xL} = 0.8$ cm/yr. Comparing with $A_{scale} = 10^7$ in Figure 3, one can see that intermediate strength crust (Figure 4) leads to the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Influence of initial thickness of zvel kinematic parameter on the architecture of rifted margins. Sediments and plastic strain are shown on markers with their respective colormaps (stratigraphy and Plastic strain). The second invariant of strain rate is represented using contour in its respective colormap. Note a constant velocity $vx = 0.8$ cm/yr is applied and that the initial width of W is 200 km for all simulations.</figcaption></figure>

The deforming zone widens with time as

$$W(t) = L_x - (2W_0 + (v_{xL} + v_{xR})(t - t_i)), \quad (2)$$

with $L_x$ the full length of the simulation domain, to ensure that the structures formed consistently with mechanical equilibrium during one phase do not exit the computational domain. While the code can handle asymmetric stretching velocities, i.e. $v_{xL} \neq v_{xR}$, all the simulations presented herein use $v_{xL} = v_{xR}$ so that we introduce $v_x = (v_{xL} + v_{xR})$.

### 3.3. Influence of the strength of the crust and stretching rate using our parametrisation

For the reconstruction, it is interesting to get the initial basin spacing as a feature of the simulation using a large W, when possible because this will constrain quite well the rheological parameter. The spacing between normal faults is known to depend, at first order, on the thickness of the resistant layer of the lithosphere [e.g., Spadini and Podladchikov, 1996], that is, the depth to the brittle ductile transition. Using only one lithological layer for the whole crust, this depth is controlled by the viscous strength of the crust and the strain rate. In our approach, we use an ad hoc non-dimensional parameter $A_{scale}$ which

modifies the Arrhenius creep law to control the viscous strength of the crust. As the power law parameter, A, is in MPa$^{-n}$.s$^{-1}$ in Table 1 (n being the creep law exponent), a $A_{scale}$ value of $10^6$ produces the laboratory creep law, a lower value produces a lower viscosity. The mean stretching rate is controlled by the value of $vx/W$.

In most cases, using a large W (200 km, Figure 3) permits to calibrate $A_{scale}$ with the natural graben spacing expected in the early rifting phase if there is not too much influence of inheritance. As expected from all the literature cited in the introduction, weak crust leads to metamorphic core complexes, detachment systems or wide rift systems. Low extension rates lead to diffuse extension, while strong crust together with intermediate to high extension rates tend to lead to narrow rift systems. Figure 3 is intended as a catalogue for the user to pick the parameters to perform the initial reconstruction step. Supplementary Figure S1 provides a more comprehensive phase-map of possible initial rift geometry as a function of $A_{scale}$ and W than what could possibly be displayed in Figure 3. Figure 4 displays the outcome of an intermediate case with $A_{scale} = 10^6$ with W = 200 km and $v_{xR} + v_{xL} = 0.8$ cm/yr. Comparing with $A_{scale} = 10^7$ in Figure 3, one can see that intermediate strength crust (Figure 4) leads to the
<figure><img src="image_6.png" /></figure>

**Figure 3.** Influence of the parametrised strength of the crust on the architecture of continental rift and symmetry of the rifted margin. Sediments and plastic strain are shown on markers with their respective color scales. The second invariant of strain rate is plotted using contour lines and its color scale. We use a wide (W = 200 km) thermo-mechanically controlled zone, so that the width of the margin and the type of rifting depends on the rheology of the crust [wide rift for a weak crust, narrow rift for a strong crust, e.g. Buck, 1991], which is parametrised with a single parameter *A*<sub>scale</sub> in our approach, and on the stretching rate W/vx.

formation of an initial graben which is then broken
in two parts while very strong crust, with its very low
level of initial necking, does not lead to the formation
of any such graben.

As mentioned earlier, it is preferable when a re-
construction can start with a large **W**. This permits
to better constrain the rheology of the crust (A<sub>scale</sub>)
using the initial basin spacing. Yet, as stated in
the introduction, heterogeneities in the crust may
sometimes lead to more localised modes of exten-
sion [Le Pourhiet et al., 2004, Huet et al., 2011a,b] or
to more diffuse modes if many heterogeneities are
present [e.g., Petri et al., 2019]. This can be dealt with
by using a small initial **W** (initially localised deforma-
tion, see next section) or a weaker on average crustal
strength (small A<sub>scale</sub>) in the calibration.

3.4. *Impact of the width* W

As mentioned before, the kinematic parameter **W** corresponds to the initial width of the zone of the lithosphere with mechanically consistent deformation during one phase (Figure 1c). We tested the impact of **W** on the architecture of rifted margins for values of 40, 120 and 200 km (Figure 4).

After 2.5 Myr, corresponding to an amount of ex-
tension of 20 km (i.e., 2.5 Myr × 0.8 cm/yr), the mod-
els with a small **W** exhibit a high level of necking
which delimit a central horst rather than a central
graben in the model with a large **W** of 200 km. Af-
ter 15 Myr and an amount of extension of 120 km, the
rifted margins with initially small **W**, that is, with a
central horst, have turned into asymmetric margins
<figure><img src="image_6.png" /></figure>

**Figure 3.** Influence of the parametrised strength of the crust on the architecture of continental rift and symmetry of the rifted margin. Sediments and plastic strain are shown on markers with their respective color scales. The second invariant of strain rate is plotted using contour lines and its color scale. We use a wide (W = 200 km) thermo-mechanically controlled zone, so that the width of the margin and the type of rifting depends on the rheology of the crust [wide rift for a weak crust, narrow rift for a strong crust, e.g. Buck, 1991], which is parametrised with a single parameter *A*<sub>scale</sub> in our approach, and on the stretching rate W/vx.

formation of an initial graben which is then broken
in two parts while very strong crust, with its very low
level of initial necking, does not lead to the formation
of any such graben.

As mentioned earlier, it is preferable when a re-
construction can start with a large **W**. This permits
to better constrain the rheology of the crust (A<sub>scale</sub>)
using the initial basin spacing. Yet, as stated in
the introduction, heterogeneities in the crust may
sometimes lead to more localised modes of exten-
sion [Le Pourhiet et al., 2004, Huet et al., 2011a,b] or
to more diffuse modes if many heterogeneities are
present [e.g., Petri et al., 2019]. This can be dealt with
by using a small initial **W** (initially localised deforma-
tion, see next section) or a weaker on average crustal
strength (small A<sub>scale</sub>) in the calibration.

3.4. *Impact of the width* W

As mentioned before, the kinematic parameter **W** corresponds to the initial width of the zone of the lithosphere with mechanically consistent deformation during one phase (Figure 1c). We tested the impact of **W** on the architecture of rifted margins for values of 40, 120 and 200 km (Figure 4).

After 2.5 Myr, corresponding to an amount of ex-
tension of 20 km (i.e., 2.5 Myr × 0.8 cm/yr), the mod-
els with a small **W** exhibit a high level of necking
which delimit a central horst rather than a central
graben in the model with a large **W** of 200 km. Af-
ter 15 Myr and an amount of extension of 120 km, the
rifted margins with initially small **W**, that is, with a
central horst, have turned into asymmetric margins
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Influence of W kinematic parameter (initial width of the thermo-mechanically controlled zone) on the architecture of continental rift and symmetry of the rifted margin for an intermediate strength ($A_{scale} = 10^6$), a constant velocity $v_{xL} = v_{xR} = 0.4$ cm/yr and **zvel** equal to 40 km. Smaller **W** may favour higher level of necking that results in final asymmetry.</figcaption></figure>

while the simulation with a central graben has re-
mained symmetric. The origin of this behaviour, ex-
cept for the initial forcing, is purely dynamic. Models
with initial central horsts develop two mantle necks
that compete until one horst evolves into an active
upwelling and inhibits the activity on the other horst.
This leads to the formation of an asymmetric mar-
gin. For W equal to 40 km, the left margin is more
extended (about 100 km) than its right counterpart
(about 60 km) and features two sub-basins instead
of one to the right. For W equal to 120 km, it is the
right margin which extends over 160 km with two
sub-basins. The left margin is smaller with 60 km in
length. It is difficult to control which side will break
in a completely dynamic simulation. The direction
of asymmetry indeed mostly relies on small numeri-
cal errors amplified by the non-linearities of the rhe-
ology. With our new kinematic module, we can stop
the simulation and pick-up the neck that was pre-
ferred in nature. Note that for values of W greater
than 200 km with this rheological configuration, the
rifted margins are symmetrical and 110 km long be-
cause the initial set up featured a low level of necking
that resulted in a single mantle neck at the centre of
the simulation.

The user and reader must keep in mind that de-
creasing **W** at constant $v_x$ lead to an increase of the

strain rate and that there are also some trade-offs
with the strength of the crust. Nonetheless, the rule
of thumb with the shallow/high and deep/low ini-
tial level of necking is robust at intermediate strength
and the impact of **W** is less important for both very
strong and very weak crusts as outlined in Figure 5.
A more comprehensive mapping of the parameter
space can be found in Supplementary Figure S2.

3.5. Workflow

Asymmetric margins can form by multiphase exten-
sion (alternating periods of extension and tectonic
quiescence/cooling) as documented in Naliboff and
Buiter [2015]. It can also be the result of the selection
of one sub-basin out of the several ones that initially
form a wide rift as highlighted in the simulation with
weak lower crust and large **W** (Figure 3). In nature,
contrarily to simplistic simulations, the selection of
one of the sub-basins might be due to a structure
propagating out of the plane [e.g., Le Pourhiet et al.,
2017, 2018]. Here, we modify the kinematic parame-
ters through time to permit the selection of a specific
sub basin. The flexibility of our approach permits to
either let the physics play its role (pause and restart
stretching after some thermal relaxation), or to en-
force localisation when the physical reason for local-
isation is either imperfectly known (inheritance) or
<figure><img src="image_7.png" /></figure>

**Figure 4.** Influence of W kinematic parameter (initial width of the thermo-mechanically controlled zone) on the architecture of continental rift and symmetry of the rifted margin for an intermediate strength ($A_{scale} = 10^6$), a constant velocity $v_{xL} = v_{xR} = 0.4$ cm/yr and **zvel** equal to 40 km. Smaller **W** may favour higher level of necking that results in final asymmetry.

while the simulation with a central graben has re-
mained symmetric. The origin of this behaviour, ex-
cept for the initial forcing, is purely dynamic. Models
with initial central horsts develop two mantle necks
that compete until one horst evolves into an active
upwelling and inhibits the activity on the other horst.
This leads to the formation of an asymmetric mar-
gin. For W equal to 40 km, the left margin is more
extended (about 100 km) than its right counterpart
(about 60 km) and features two sub-basins instead
of one to the right. For W equal to 120 km, it is the
right margin which extends over 160 km with two
sub-basins. The left margin is smaller with 60 km in
length. It is difficult to control which side will break
in a completely dynamic simulation. The direction
of asymmetry indeed mostly relies on small numeri-
cal errors amplified by the non-linearities of the rhe-
ology. With our new kinematic module, we can stop
the simulation and pick-up the neck that was pre-
ferred in nature. Note that for values of W greater
than 200 km with this rheological configuration, the
rifted margins are symmetrical and 110 km long be-
cause the initial set up featured a low level of necking
that resulted in a single mantle neck at the centre of
the simulation.

The user and reader must keep in mind that de-
creasing **W** at constant $v_x$ lead to an increase of the

strain rate and that there are also some trade-offs
with the strength of the crust. Nonetheless, the rule
of thumb with the shallow/high and deep/low ini-
tial level of necking is robust at intermediate strength
and the impact of **W** is less important for both very
strong and very weak crusts as outlined in Figure 5.
A more comprehensive mapping of the parameter
space can be found in Supplementary Figure S2.

3.5. Workflow

Asymmetric margins can form by multiphase exten-
sion (alternating periods of extension and tectonic
quiescence/cooling) as documented in Naliboff and
Buiter [2015]. It can also be the result of the selection
of one sub-basin out of the several ones that initially
form a wide rift as highlighted in the simulation with
weak lower crust and large **W** (Figure 3). In nature,
contrarily to simplistic simulations, the selection of
one of the sub-basins might be due to a structure
propagating out of the plane [e.g., Le Pourhiet et al.,
2017, 2018]. Here, we modify the kinematic parame-
ters through time to permit the selection of a specific
sub basin. The flexibility of our approach permits to
either let the physics play its role (pause and restart
stretching after some thermal relaxation), or to en-
force localisation when the physical reason for local-
isation is either imperfectly known (inheritance) or
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Illustration of the existing yet limited trade-off between rheology and the choice of the width of rift zone. Weak and strong crust are not significantly affected by the choice of *W*, except for the age of sediments in the basins. For an intermediate rheology of the crust, the rate effect introduced by the choice of *W* is important in forming or not an initial mantle horst.</figcaption></figure>

coming out-of-plane (strike-slip or propagating ridge segment).

The method for the kinematic reconstruction of a rifted margin involves launching several models at each phase with a different range of velocity (*v*) and to choose the best model, i.e., the one that produces a crustal thickness (i.e., depth of the Moho) consistent with the data we try to adjust.

Among the launched models, we select the best fitting one and restart several new models from the time step at which we stopped. This operation is performed several times until the entire reconstruction of the rifted margin. The simulations are checkpointed after each phase of deformation, allowing us to modify the latter phase of deformation without re-computing the initial phase. The reconstruction approach of passive margins developed in this study involves determining the appropriate velocity (i.e., to find *W*, $X_0$, $v_{xL}$ and $v_{xR}$) at each phase of the simulation exactly like in the Jeanniot et al. [2016] approach.

A part of the calibration therefore consists in trials and errors improved by the results of the sen-

sitivity tests and the experience of the user with thermo-mechanical modelling. The other part relies on starting the reconstruction with a pre-existing reliable back stripping model of the sedimentary basins that will help constraining the reconstruction. The goal of the kinematic module is not to replace back-stripping tools, but to produce geometry and temperature fields in the lower crust, mantle and LAB that are compatible with the movement of the basement of the basins deduced from back-stripping.

**4. Benchmark: kinematic calibration of a synthetic 2D cross-section extracted from a 3D numerical model**

**4.1. Incentives to benchmark kinematic calibrations using complex numerical models**

In nature, the amount of available data, and particularly those for thermal calibration, is quite limited and will never be as complete as what is recovered from a numerical simulation. Moreover, all the data
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Illustration of the existing yet limited trade-off between rheology and the choice of the width of rift zone. Weak and strong crust are not significantly affected by the choice of *W*, except for the age of sediments in the basins. For an intermediate rheology of the crust, the rate effect introduced by the choice of *W* is important in forming or not an initial mantle horst.</figcaption></figure>

coming out-of-plane (strike-slip or propagating ridge segment).

The method for the kinematic reconstruction of a rifted margin involves launching several models at each phase with a different range of velocity (*v*) and to choose the best model, i.e., the one that produces a crustal thickness (i.e., depth of the Moho) consistent with the data we try to adjust.

Among the launched models, we select the best fitting one and restart several new models from the time step at which we stopped. This operation is performed several times until the entire reconstruction of the rifted margin. The simulations are checkpointed after each phase of deformation, allowing us to modify the latter phase of deformation without recomputing the initial phase. The reconstruction approach of passive margins developed in this study involves determining the appropriate velocity (i.e., to find *W*, $X_0$, $v_{xL}$ and $v_{xR}$) at each phase of the simulation exactly like in the Jeanniot et al. [2016] approach.

A part of the calibration therefore consists in trials and errors improved by the results of the sen-

sitivity tests and the experience of the user with thermo-mechanical modelling. The other part relies on starting the reconstruction with a pre-existing reliable back stripping model of the sedimentary basins that will help constraining the reconstruction. The goal of the kinematic module is not to replace back-stripping tools, but to produce geometry and temperature fields in the lower crust, mantle and LAB that are compatible with the movement of the basement of the basins deduced from back-stripping.

**4. Benchmark: kinematic calibration of a synthetic 2D cross-section extracted from a 3D numerical model**

**4.1. Incentives to benchmark kinematic calibrations using complex numerical models**

In nature, the amount of available data, and particularly those for thermal calibration, is quite limited and will never be as complete as what is recovered from a numerical simulation. Moreover, all the data
<figure><img src="image_9.png" /><figcaption>**Figure 6.** 2D cross-section bc extracted from a 3D thermomechanical numerical model with the determination of the reconstruction parameters. The section displays 4 discrete basins. For each phase/basin *i*, we determine $x_{bc}^i$ the abscissa of the depocenter, $w_{bc}^i$ the final width of the basin *i*, $d_{bc}^i$ the distance between each basin depocenter and $L_{bc}^i$ the crustal thickness at the depocenter of the basin *i*. Note that superscript *L* and *R* indicate respectively the left and right side of the basin 0 that was split at the moment of continental break-up.</figcaption></figure>

available for a case study are used for model calibration so that model prediction cannot really be tested unless a new borehole is drilled. Calibrating a kinematic model using the result of a complex fully self-consistent numerical simulation present the advantages of having an unlimited synthetic borehole data, or simply to access the depth of the LAB through time.

3D thermomechanical numerical simulation is currently the way to produce very complex physically consistent models of long-term deformation of lithosphere–asthenosphere systems. Moreover, as mentioned earlier, data are usually collected along 2D high resolution seismic profiles, while crustal and lithospheric deformation is by essence three-dimensional. In order to validate the new implemented kinematic module in pTatin2D, we have calibrated a 2D simulation to produce a model bt (blind test) that has a similar Moho and top basement geometry as the synthetic 2D cross-section bc (blind cross-section). Yet, the cross-section bc here is not a seismic profile but, instead, has been extracted from a 3D thermomechanical model of passive margin

formation published in Jourdon et al. [2020]. This model was chosen to have some significant out-of-plane deformation to make the exercise closer to typical case studies on natural cases using 2D high resolution seismic lines at the crustal scale. The first author, in charge of the calibration, was not informed from which 3D model the bc section that was given to him was extracted. Using a synthetic section permits to have access to the LAB a posteriori to validate that a thermo-mechanical simulation calibrated with the geometry of the crust to a geological cross-section might be predictive for the LAB and the mantle heat flow.

To perform the calibration, only the top-basement and the Moho geometry of the cross-section displayed in Figure 6 were used. The rest of the information available from the synthetic cross-section bc, such as the LAB evolution with time, is only used to evaluate the predictivity of the approach.

In other words, we utilize only a portion of the data, two interfaces at one time step, for calibration purposes, while retaining the remaining (temperature field, LAB, heat flow and cross-section evolving
<figure><img src="image_9.png" /><figcaption>**Figure 6.** 2D cross-section bc extracted from a 3D thermomechanical numerical model with the determination of the reconstruction parameters. The section displays 4 discrete basins. For each phase/basin *i*, we determine $x_{bc}^i$ the abscissa of the depocenter, $w_{bc}^i$ the final width of the basin *i*, $d_{bc}^i$ the distance between each basin depocenter and $L_{bc}^i$ the crustal thickness at the depocenter of the basin *i*. Note that superscript *L* and *R* indicate respectively the left and right side of the basin 0 that was split at the moment of continental break-up.</figcaption></figure>

available for a case study are used for model calibration so that model prediction cannot really be tested unless a new borehole is drilled. Calibrating a kinematic model using the result of a complex fully self-consistent numerical simulation present the advantages of having an unlimited synthetic borehole data, or simply to access the depth of the LAB through time.

3D thermomechanical numerical simulation is currently the way to produce very complex physically consistent models of long-term deformation of lithosphere–asthenosphere systems. Moreover, as mentioned earlier, data are usually collected along 2D high resolution seismic profiles, while crustal and lithospheric deformation is by essence three-dimensional. In order to validate the new implemented kinematic module in pTatin2D, we have calibrated a 2D simulation to produce a model bt (blind test) that has a similar Moho and top basement geometry as the synthetic 2D cross-section bc (blind cross-section). Yet, the cross-section bc here is not a seismic profile but, instead, has been extracted from a 3D thermomechanical model of passive margin

formation published in Jourdon et al. [2020]. This model was chosen to have some significant out-of-plane deformation to make the exercise closer to typical case studies on natural cases using 2D high resolution seismic lines at the crustal scale. The first author, in charge of the calibration, was not informed from which 3D model the bc section that was given to him was extracted. Using a synthetic section permits to have access to the LAB a posteriori to validate that a thermo-mechanical simulation calibrated with the geometry of the crust to a geological cross-section might be predictive for the LAB and the mantle heat flow.

To perform the calibration, only the top-basement and the Moho geometry of the cross-section displayed in Figure 6 were used. The rest of the information available from the synthetic cross-section bc, such as the LAB evolution with time, is only used to evaluate the predictivity of the approach.

In other words, we utilize only a portion of the data, two interfaces at one time step, for calibration purposes, while retaining the remaining (temperature field, LAB, heat flow and cross-section evolving
with time) for verification to ensure the accuracy of the calibration between synthetic data points. Although this approach does not guarantee a perfect calibration with natural data, it is satisfactory for validation of the calibration methods. If the calibration fails to produce accurate predictions between synthetic data points, it indicates that a similar approach based on real (and often even more fragmented) data is unlikely to be predictive.

## 4.2. Reconstruction parameters of the synthetic cross-section bc

The absence of sediments and therefore of stratigraphic markers in the synthetic cross-section makes the calibration step more complicated than that of a real cross-section calibrated by borehole data. A quick look at the architecture of the synthetic cross-section highlights an asymmetric morphology of the margins. The right margin can be characterized by four sub-basins or four structural units (Figure 6). The first one, at the far right (bc³), seems to be controlled by a nearly vertical structure linked to left-dipping normal faults which would be interpreted as having a strike-slip component in nature. The second one (bc²), to the left of the previous one, is bounded by two synthetic normal faults forming a classical graben system. The two last sub-basins (bc¹ and bc⁰) of this margin part show tilted blocks to the left and a hyper-extended necking structure. In contrast, the left margin side (Figure 6) is featured by a simpler margin morphology with tilted blocks to the right, forming a necking zone.

The initial thickness of the crust away from necking seems to be 40 km (Figure 6) so we can consider it as initial condition.

## 4.3. Calibration to top Basement and Moho impact on LAB and Tmax

In order to calibrate our simulation bt to the crustal geometry of the rifted margins extracted along the synthetic cross-section bc (Figure 6), we first compute the kinematic parameters using the reconstruction parameters listed in Table 2 and represented on Figure 6. As we only know the final geometry, we first measure the amount of extension *e* associated to each phase of rifting and, assuming that

<table><caption>**Table 2.** Reconstruction parameters for the basins from the synthetic cross-section presented in Figure 5</caption><thead><tr><th>Name</th><th>x</th><th>w</th><th>d</th><th>L</th><th>Tectonic phases</th></tr></thead><tbody><tr><td>bc<sup>0L</sup></td><td>595</td><td>20</td><td>40</td><td>5</td><td>4</td></tr><tr><td>bc<sup>0R</sup></td><td>370</td><td>20</td><td>40</td><td>5</td><td>4</td></tr><tr><td>bc<sup>1</sup></td><td>615</td><td>40</td><td>20</td><td>5</td><td>3</td></tr><tr><td>bc<sup>2</sup></td><td>740</td><td>30</td><td>125</td><td>20</td><td>1</td></tr><tr><td>bc<sup>3</sup></td><td>850</td><td>20</td><td>110</td><td>15</td><td>2</td></tr></tbody></table>

*x* is the localisation of the maximum depocenter of basins/sub-basins at the final stage, *w* is the width of the basins/sub-basins at the final stage, *d* is the spacing between basins/sub-basins at the depocenter, *L* is the crustal thickness at the depocenter, all given in km.

<table><caption>**Table 3.** Kinematic parameters of the calibration bt of the synthetic cross-section bc</caption><thead><tr><th>Phases</th><th>p<sup>0</sup></th><th>p<sup>1</sup></th><th>p<sup>2</sup></th><th>p<sup>3</sup></th><th>p<sup>4</sup></th><th>p<sup>5</sup></th></tr></thead><tbody><tr><td>v<sub>x</sub> (cm/yr)</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.227</td></tr><tr><td>t<sub>e</sub> (Myr)</td><td>2</td><td>8</td><td>14</td><td>20</td><td>24</td><td>80</td></tr><tr><td>t<sub>i</sub> (Myr)</td><td>0</td><td>2</td><td>8</td><td>14</td><td>20</td><td>24</td></tr><tr><td>X<sub>0</sub> (km)</td><td>600</td><td>600</td><td>690</td><td>485</td><td>465</td><td>465</td></tr><tr><td>W<sub>0</sub> (km)</td><td>150</td><td>50</td><td>30</td><td>20</td><td>20</td><td>150</td></tr><tr><td>e (km)</td><td>10</td><td>30</td><td>30</td><td>30</td><td>20</td><td>127</td></tr></tbody></table>

Note that the phases P⁰ and P⁵, respectively initiating and finalising the rifting are characterized by large W₀ so that mechanical equilibrium dominates over kinematic forcing during these phases.

$e = (v_{xL} - v_{xR})/(t_e - t_i)$, we modify the timing and the velocity using our stop-and-go approach to test different combinations. We found that the kinematic parameters, listed in Table 3, lead to a reasonably good adjustment of the Moho and top basement geometry of model bt to the synthetic cross-section bc presented in Figure 7. This combination is not unique because there are no erosion and sedimentation in the synthetic cross-section to date the different phase, and there are clearly no constraints on $t_e$ and $t_i$. The choice was solely based on the formation of structures at the correct location.

In order to reconstruct the synthetic cross-section (bc), the calibrated model bt needs 6 phases with
with time) for verification to ensure the accuracy of the calibration between synthetic data points. Although this approach does not guarantee a perfect calibration with natural data, it is satisfactory for validation of the calibration methods. If the calibration fails to produce accurate predictions between synthetic data points, it indicates that a similar approach based on real (and often even more fragmented) data is unlikely to be predictive.

## 4.2. Reconstruction parameters of the synthetic cross-section bc

The absence of sediments and therefore of stratigraphic markers in the synthetic cross-section makes the calibration step more complicated than that of a real cross-section calibrated by borehole data. A quick look at the architecture of the synthetic cross-section highlights an asymmetric morphology of the margins. The right margin can be characterized by four sub-basins or four structural units (Figure 6). The first one, at the far right (bc³), seems to be controlled by a nearly vertical structure linked to left-dipping normal faults which would be interpreted as having a strike-slip component in nature. The second one (bc²), to the left of the previous one, is bounded by two synthetic normal faults forming a classical graben system. The two last sub-basins (bc¹ and bc⁰) of this margin part show tilted blocks to the left and a hyper-extended necking structure. In contrast, the left margin side (Figure 6) is featured by a simpler margin morphology with tilted blocks to the right, forming a necking zone.

The initial thickness of the crust away from necking seems to be 40 km (Figure 6) so we can consider it as initial condition.

## 4.3. Calibration to top Basement and Moho impact on LAB and Tmax

In order to calibrate our simulation bt to the crustal geometry of the rifted margins extracted along the synthetic cross-section bc (Figure 6), we first compute the kinematic parameters using the reconstruction parameters listed in Table 2 and represented on Figure 6. As we only know the final geometry, we first measure the amount of extension *e* associated to each phase of rifting and, assuming that

<table><caption>**Table 2.** Reconstruction parameters for the basins from the synthetic cross-section presented in Figure 5</caption><thead><tr><th>Name</th><th>x</th><th>w</th><th>d</th><th>L</th><th>Tectonic phases</th></tr></thead><tbody><tr><td>bc<sup>0L</sup></td><td>595</td><td>20</td><td>40</td><td>5</td><td>4</td></tr><tr><td>bc<sup>0R</sup></td><td>370</td><td>20</td><td>40</td><td>5</td><td>4</td></tr><tr><td>bc<sup>1</sup></td><td>615</td><td>40</td><td>20</td><td>5</td><td>3</td></tr><tr><td>bc<sup>2</sup></td><td>740</td><td>30</td><td>125</td><td>20</td><td>1</td></tr><tr><td>bc<sup>3</sup></td><td>850</td><td>20</td><td>110</td><td>15</td><td>2</td></tr></tbody></table>

*x* is the localisation of the maximum depocenter of basins/sub-basins at the final stage, *w* is the width of the basins/sub-basins at the final stage, *d* is the spacing between basins/sub-basins at the depocenter, *L* is the crustal thickness at the depocenter, all given in km.

<table><caption>**Table 3.** Kinematic parameters of the calibration bt of the synthetic cross-section bc</caption><thead><tr><th>Phases</th><th>p<sup>0</sup></th><th>p<sup>1</sup></th><th>p<sup>2</sup></th><th>p<sup>3</sup></th><th>p<sup>4</sup></th><th>p<sup>5</sup></th></tr></thead><tbody><tr><td>v<sub>x</sub> (cm/yr)</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.227</td></tr><tr><td>t<sub>e</sub> (Myr)</td><td>2</td><td>8</td><td>14</td><td>20</td><td>24</td><td>80</td></tr><tr><td>t<sub>i</sub> (Myr)</td><td>0</td><td>2</td><td>8</td><td>14</td><td>20</td><td>24</td></tr><tr><td>X<sub>0</sub> (km)</td><td>600</td><td>600</td><td>690</td><td>485</td><td>465</td><td>465</td></tr><tr><td>W<sub>0</sub> (km)</td><td>150</td><td>50</td><td>30</td><td>20</td><td>20</td><td>150</td></tr><tr><td>e (km)</td><td>10</td><td>30</td><td>30</td><td>30</td><td>20</td><td>127</td></tr></tbody></table>

Note that the phases P⁰ and P⁵, respectively initiating and finalising the rifting are characterized by large W₀ so that mechanical equilibrium dominates over kinematic forcing during these phases.

$e = (v_{xL} - v_{xR})/(t_e - t_i)$, we modify the timing and the velocity using our stop-and-go approach to test different combinations. We found that the kinematic parameters, listed in Table 3, lead to a reasonably good adjustment of the Moho and top basement geometry of model bt to the synthetic cross-section bc presented in Figure 7. This combination is not unique because there are no erosion and sedimentation in the synthetic cross-section to date the different phase, and there are clearly no constraints on $t_e$ and $t_i$. The choice was solely based on the formation of structures at the correct location.

In order to reconstruct the synthetic cross-section (bc), the calibrated model bt needs 6 phases with
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Our favorite kinematically forced model bt of the cross-section bc. The brown transparent square delimits the zone of kinematic forcing at each phase.</figcaption></figure>

a total amount of extension about 247 km during 80 Ma. The amounts of extension, *e*, of each phase, are also gathered in Table 3.

In our calibration, two extra phases are added to the kinematic reconstruction ($P^0$ and $P^5$ in Table 3): (1) One is the initiation of the rifting with a large $W_0$ and central $X_0$ and (2) another one is set after the break-up with a large $W_0$ and with the same $X_0$ than the last phase to allow for a dynamic post-rift phase that accounts for realistic convection dynamics (Figure 7).

The phase $P^1$ is required to reconstruct the $bc^2$ structural unit in Figure 6. The next phase, $P^2$ permits to create $bc^3$ flower structure in Figure 6. The two last phases, $P^3$ and $P^4$, are needed to respectively build the structural unit $bc^1$ and to give their characteristic geometries to the $bc^{0R}$ and $bc^{0L}$ conjugate margins in Figure 6. We further discuss in part 4.5 the impact of the ordering of the phases.

## 4.4. Validation and limits for predicting mantle dynamics

Concerning predictivity of the calibration, we use the LAB geometry as a proxy for the thermal regime of the lithosphere. The comparison between the LAB of the synthetic cross-section and the LAB of our best bt calibration displays a small yet significant difference while comparing the shape of the LAB in the

centre of Figure 8. The difference in temperature on the lateral side is possibly due to different thicknesses of the domain in the 2D and 3D simulations. The “mushroom” aspect of the green bt curve compared to the “Gaussian” aspect of the bc curve is more significant. This “mushroom” shape of the LAB in our calibrated cross-section indeed materializes convection cells in the asthenosphere which are not present or less vigorous in the 3D model. In 3D, the drag acting on a buoyant cylinder is smaller than that acting on a buoyant sphere with the same diameter. As a result, the cylinders rise faster than the sphere causing more vigorous convection. In 2D simulations, all structures and temperature anomalies present in the cross-section are cylindrical (infinite in the third dimension). As a result, 2D plane strain simulations are known to overestimate the convective vigour of the flow. For a complete proof of this statement, the reader can refer to Weinberg and Podladchikov [1994].

Based on our understanding of mantle convection and 2D approximation, we therefore tested whether a reduction of the coefficient of thermal expansion from $\alpha = 3 \times 10^{-5} \text{ K}^{-1}$ to $\alpha = 1 \times 10^{-5} \text{ K}^{-1}$ (see Figure 8 caption and Appendix A) would improve the shape of the curve. This implies that when reconstructing a section with a large strike-slip component, it might be wise to use slightly smaller coefficients of thermal expansion to better predict the temperature and
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Our favorite kinematically forced model bt of the cross-section bc. The brown transparent square delimits the zone of kinematic forcing at each phase.</figcaption></figure>

a total amount of extension about 247 km during 80 Ma. The amounts of extension, *e*, of each phase, are also gathered in Table 3.

In our calibration, two extra phases are added to the kinematic reconstruction ($P^0$ and $P^5$ in Table 3): (1) One is the initiation of the rifting with a large $W_0$ and central $X_0$ and (2) another one is set after the break-up with a large $W_0$ and with the same $X_0$ than the last phase to allow for a dynamic post-rift phase that accounts for realistic convection dynamics (Figure 7).

The phase $P^1$ is required to reconstruct the $bc^2$ structural unit in Figure 6. The next phase, $P^2$ permits to create $bc^3$ flower structure in Figure 6. The two last phases, $P^3$ and $P^4$, are needed to respectively build the structural unit $bc^1$ and to give their characteristic geometries to the $bc^{0R}$ and $bc^{0L}$ conjugate margins in Figure 6. We further discuss in part 4.5 the impact of the ordering of the phases.

## 4.4. Validation and limits for predicting mantle dynamics

Concerning predictivity of the calibration, we use the LAB geometry as a proxy for the thermal regime of the lithosphere. The comparison between the LAB of the synthetic cross-section and the LAB of our best bt calibration displays a small yet significant difference while comparing the shape of the LAB in the

centre of Figure 8. The difference in temperature on the lateral side is possibly due to different thicknesses of the domain in the 2D and 3D simulations. The “mushroom” aspect of the green bt curve compared to the “Gaussian” aspect of the bc curve is more significant. This “mushroom” shape of the LAB in our calibrated cross-section indeed materializes convection cells in the asthenosphere which are not present or less vigorous in the 3D model. In 3D, the drag acting on a buoyant cylinder is smaller than that acting on a buoyant sphere with the same diameter. As a result, the cylinders rise faster than the sphere causing more vigorous convection. In 2D simulations, all structures and temperature anomalies present in the cross-section are cylindrical (infinite in the third dimension). As a result, 2D plane strain simulations are known to overestimate the convective vigour of the flow. For a complete proof of this statement, the reader can refer to Weinberg and Podladchikov [1994].

Based on our understanding of mantle convection and 2D approximation, we therefore tested whether a reduction of the coefficient of thermal expansion from $\alpha = 3 \times 10^{-5} \text{ K}^{-1}$ to $\alpha = 1 \times 10^{-5} \text{ K}^{-1}$ (see Figure 8 caption and Appendix A) would improve the shape of the curve. This implies that when reconstructing a section with a large strike-slip component, it might be wise to use slightly smaller coefficients of thermal expansion to better predict the temperature and
<figure><img src="image_11.png" /></figure>

**Figure 8.** Comparison between the 2D synthetic cross-section (bc) and different calibrations. bt is our preferred calibration without the LAB, bti is a calibration in which we inversed the timing between phase P² and P³; bta is a calibration with smaller coefficient of the thermal expansion $\alpha = 1 \times 10^{-5} \text{ K}^{-1}$. It was made to better adjust the fit to the LAB in the calibration.

mantle convection using 2D section calibrated with
crustal deformation.

## 4.5. Validation and limits for reconstruction

The kinematics retrieved during calibration is con-
sistent with the rifting propagation observed in the
a posteriori 3D model. Indeed, the extension starts
by forming a central rift which corresponds to the
structure bc² on the synthetic cross-section that is
abandoned thereafter. Then, a strike-slip fault shifts
the units to form the structure bc³ on the synthetic
cross-section. Finally, the continental break-up oc-
curs to the left of all these structures along with
basins identified as bcº and bc¹ on the synthetic
cross-section. Despite some resolution uncertain-
ties especially near the flower structure, the 2D cal-
ibration (bt) manages to reproduce this 3D problem
structure.

To assess whether the calibration is sensitive to
the order and the timing of the different phases, we
attempted a calibration bti in which P² and P³ are
switched in Table 3. When we compare this cali-
bration with our preferred calibration bt, as well as
with the synthetic cross-section bc in Figure 7, we
observe that the calibration bti fails in reproducing
the basin bc³ shown in Figure 5. This means that
the sequence in which we apply the stretching phase
not only influences the temperature but also the fi-
nal geometry of the model. This influence of the de-
formation path on the architecture primarily stems

from the thermo-mechanical feedback on the rheol-
ogy. This demonstrates that the thermo-mechanical
calibration approach can help mitigating uncertain-
ties inherent in a purely kinematic approach during
the back-stripping phase.

To conclude the kinematic module integrated to
the 2D thermomechanical code pTatin2D succeeds
in modelling complex 3D passive margins with a
good predictivity in term of lithospheric thermal
state. Moreover, this approach can be used to test
different reconstruction options as the final geome-
try is found to be path-dependent contrarily to the
purely kinematic approach. One should note that the
method is very costly and should only be used ad-
dress uncertainties between two scenarios when the
timing of the deformation, borehole data and clear
time correlations are missing.

Consequently, despite some limits, our new mod-
ule can be used to reconstruct some natural cases of
passive margins and to predict their thermal history
where calibration points are missing.

# 5. Application to Iberia–Newfoundland conjugate rifted continental margins

## 5.1. *Geological setting*

The Iberia–Newfoundland rifted continental mar-
gins are located in the North Atlantic Ocean (Fig-
ure 9a). The rifted margins are asymmetric, hyper-
extended [Sutra et al., 2013, Manatschal et al., 2007]
<figure><img src="image_11.png" /></figure>

**Figure 8.** Comparison between the 2D synthetic cross-section (bc) and different calibrations. bt is our preferred calibration without the LAB, bti is a calibration in which we inverted the timing between phase P² and P³; bta is a calibration with smaller coefficient of the thermal expansion $α = 1 × 10^{-5}$ K⁻¹. It was made to better adjust the fit to the LAB in the calibration.

mantle convection using 2D section calibrated with crustal deformation.

## 4.5. Validation and limits for reconstruction

The kinematics retrieved during calibration is consistent with the rifting propagation observed in the a posteriori 3D model. Indeed, the extension starts by forming a central rift which corresponds to the structure bc² on the synthetic cross-section that is abandoned thereafter. Then, a strike-slip fault shifts the units to form the structure bc³ on the synthetic cross-section. Finally, the continental break-up occurs to the left of all these structures along with basins identified as bc⁰ and bc¹ on the synthetic cross-section. Despite some resolution uncertainties especially near the flower structure, the 2D calibration (bt) manages to reproduce this 3D problem structure.

To assess whether the calibration is sensitive to the order and the timing of the different phases, we attempted a calibration bti in which P² and P³ are switched in Table 3. When we compare this calibration with our preferred calibration bt, as well as with the synthetic cross-section bc in Figure 7, we observe that the calibration bti fails in reproducing the basin bc³ shown in Figure 5. This means that the sequence in which we apply the stretching phase not only influences the temperature but also the final geometry of the model. This influence of the deformation path on the architecture primarily stems

from the thermo-mechanical feedback on the rheology. This demonstrates that the thermo-mechanical calibration approach can help mitigating uncertainties inherent in a purely kinematic approach during the back-stripping phase.

To conclude the kinematic module integrated to the 2D thermomechanical code pTatin2D succeeds in modelling complex 3D passive margins with a good predictivity in term of lithospheric thermal state. Moreover, this approach can be used to test different reconstruction options as the final geometry is found to be path-dependent contrarily to the purely kinematic approach. One should note that the method is very costly and should only be used address uncertainties between two scenarios when the timing of the deformation, borehole data and clear time correlations are missing.

Consequently, despite some limits, our new module can be used to reconstruct some natural cases of passive margins and to predict their thermal history where calibration points are missing.

# 5. Application to Iberia–Newfoundland conjugate rifted continental margins

## 5.1. Geological setting

The Iberia–Newfoundland rifted continental margins are located in the North Atlantic Ocean (Figure 9a). The rifted margins are asymmetric, hyperextended [Sutra et al., 2013, Manatschal et al., 2007]
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Bathymetry of the Iberia–Newfoundland rifted margins in the Atlantic North from ETOPO1 1 Arc-Minute Global Relief Model. Profiles locations, the ODP drill holes and the magnetic anomalies M0 and M3 are marked on the map modified from Jeanniot et al. [2016]. (b) Seismic reflection/refraction interpreted profile of North Iberia–Newfoundland conjugate rifted margins modified from Sutra et al. [2013]. The kinematics inputs (x, *w* and *L*) are compiled in Table 4.</figcaption></figure>

and are composed from east to west (Figure 9b) by the Galicia Interior Bank (GIB), the Galicia Bank (GB), the Deep Galicia Margin (DGM), the Newfoundland Basin (NB) and the Flemish Cap (FC). The formation of these conjugate margins is propagating from south to north [Nirrengarten et al., 2018] and results from the superposition of two consecutive rifts events [Péron-Pinvidic and Manatschal, 2009]: a first event dated as Late Triassic to Early Jurassic and a second, more important event dated as Late Jurassic to Early Cretaceous and leading to the breakup of the southern North Atlantic [Tucholke et al., 2007,

Tucholke and Sibuet, 2007]. During the first phase, the extension was distributed over a wide area and the crustal thinning was very small [Tucholke and Sibuet, 2007, Mohn et al., 2015]. To sum up, the calendar of these conjugate margins is as follows: Stretching (~161 Ma), necking (~145 Ma), mantle exhumation (~137 Ma) until break-up (~115–112 Ma). So, the duration of the whole rifting sequence is about ~50 Myr until break-up [Sutra et al., 2013, Eddy et al., 2017, Tucholke et al., 2007, Bronner et al., 2011, Tucholke and Sibuet, 2007, Péron-Pinvidic et al., 2007] and syn-rift sediment ages recovered at ODP Sites
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Bathymetry of the Iberia–Newfoundland rifted margins in the Atlantic North from ETOPO1 1 Arc-Minute Global Relief Model. Profiles locations, the ODP drill holes and the magnetic anomalies M0 and M3 are marked on the map modified from Jeanniot et al. [2016]. (b) Seismic reflection/refraction interpreted profile of North Iberia–Newfoundland conjugate rifted margins modified from Sutra et al. [2013]. The kinematics inputs (x, *w* and *L*) are compiled in Table 4.</figcaption></figure>

and are composed from east to west (Figure 9b) by the Galicia Interior Bank (GIB), the Galicia Bank (GB), the Deep Galicia Margin (DGM), the Newfoundland Basin (NB) and the Flemish Cap (FC). The formation of these conjugate margins is propagating from south to north [Nirrengarten et al., 2018] and results from the superposition of two consecutive rifts events [Péron-Pinvidic and Manatschal, 2009]: a first event dated as Late Triassic to Early Jurassic and a second, more important event dated as Late Jurassic to Early Cretaceous and leading to the breakup of the southern North Atlantic [Tucholke et al., 2007,

Tucholke and Sibuet, 2007]. During the first phase, the extension was distributed over a wide area and the crustal thinning was very small [Tucholke and Sibuet, 2007, Mohn et al., 2015]. To sum up, the calendar of these conjugate margins is as follows: Stretching (~161 Ma), necking (~145 Ma), mantle exhumation (~137 Ma) until break-up (~115–112 Ma). So, the duration of the whole rifting sequence is about ~50 Myr until break-up [Sutra et al., 2013, Eddy et al., 2017, Tucholke et al., 2007, Bronner et al., 2011, Tucholke and Sibuet, 2007, Péron-Pinvidic et al., 2007] and syn-rift sediment ages recovered at ODP Sites
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Our favorite kinematically forced model. Brown transparent square delimits the zone of kinematic forcing at each phase.</figcaption></figure>

638, 641 and 640 [Boillot and Winterer, 1988] permit to constrain the subsidence in the distal part.

Numerous modelling studies have been already carried out in order to better constrain the evolution and the peculiar architecture of the Iberia–Newfoundland conjugate margins [Nagel and Buck, 2004; Lavier and Manatschal, 2006; Pérez-Gussinyé et al., 2006; Burov, 2007; Mohn et al., 2015; Jeanniot et al., 2016; Brune et al., 2017]. Here, we focus on the north cross-section of Iberia–Newfoundland conjugate rifted margins which is documented as a magma-poor type margin [Manatschal et al., 2007] and therefore can be suitably modelled with our new module.

## 5.2. *Model calibration*

As previously noted, our new module is not the most efficient tool to test a large panel of back-stripping scenarios, but aims at computing a mantle scenario thermo-mechanically consistent with the crustal deformation. Following the same methodology as for the benchmark test, we have compartmentalized the Iberia–Newfoundland conjugate rifted margins in five tectonic units of final widths $w_{FC}$, $w_{NB}/w_{DGM}$, $w_{GB}$, $w_{GIBd}$ and $w_{GIBp}$. The Newfoundland Basin (NB) and Deep Galicia Margin (DGM) are considered

as a single structural unit because they represent a conjugate structure. We separated the Galicia Interior Bank (GIB) into two parts called the GIBd (deep) and GIBp (plateau).

Using the same methodology as for the benchmark test, we derived the values for reconstruction parameters *w*, *x*, *d*, *L*, and the duration of the modeled rifting phases for the Iberia–Newfoundland rifted margins case (Figure 9b). These parameters are summarized in Table 4. According to the methodology described in the previous part, we calibrated our model to the geophysical model of the conjugate rifted margins of the Iberia–Newfoundland (Figure 9). The results of our calibration are presented in Figure 10. In this showcase of our method on a natural example, contrarily to the benchmark, the sediment packages are accounted for by the calibration. The subsidence history of the forward kinematic modelling in the DGM unit is constrained by subsidence estimates available from Boillot et al. [1988] in the ODP Sites 638–639 (Figure 9 for localization and Figure 11b for data). This information also permits to calibrate the timing and duration of the different modelled rifting phases.

Lastly, the calibrated model presented in Figure 10 needs 8 phases listed in Table 5 in order to
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Our favorite kinematically forced model. Brown transparent square delimits the zone of kinematic forcing at each phase.</figcaption></figure>

638, 641 and 640 [Boillot and Winterer, 1988] permit to constrain the subsidence in the distal part.

Numerous modelling studies have been already carried out in order to better constrain the evolution and the peculiar architecture of the Iberia–Newfoundland conjugate margins [Nagel and Buck, 2004; Lavier and Manatschal, 2006; Pérez-Gussinyé et al., 2006; Burov, 2007; Mohn et al., 2015; Jeanniot et al., 2016; Brune et al., 2017]. Here, we focus on the north cross-section of Iberia–Newfoundland conjugate rifted margins which is documented as a magma-poor type margin [Manatschal et al., 2007] and therefore can be suitably modelled with our new module.

## 5.2. *Model calibration*

As previously noted, our new module is not the most efficient tool to test a large panel of back-stripping scenarios, but aims at computing a mantle scenario thermo-mechanically consistent with the crustal deformation. Following the same methodology as for the benchmark test, we have compartmentalized the Iberia–Newfoundland conjugate rifted margins in five tectonic units of final widths $w_{FC}$, $w_{NB} / w_{DGM}$, $w_{GB}$, $w_{GIBd}$ and $w_{GIBp}$. The Newfoundland Basin (NB) and Deep Galicia Margin (DGM) are considered

as a single structural unit because they represent a conjugate structure. We separated the Galicia Interior Bank (GIB) into two parts called the GIBd (deep) and GIBp (plateau).

Using the same methodology as for the benchmark test, we derived the values for reconstruction parameters *w*, *x*, *d*, *L*, and the duration of the modeled rifting phases for the Iberia–Newfoundland rifted margins case (Figure 9b). These parameters are summarized in Table 4. According to the methodology described in the previous part, we calibrated our model to the geophysical model of the conjugate rifted margins of the Iberia–Newfoundland (Figure 9). The results of our calibration are presented in Figure 10. In this showcase of our method on a natural example, contrarily to the benchmark, the sediment packages are accounted for by the calibration. The subsidence history of the forward kinematic modelling in the DGM unit is constrained by subsidence estimates available from Boillot et al. [1988] in the ODP Sites 638–639 (Figure 9 for localization and Figure 11b for data). This information also permits to calibrate the timing and duration of the different modelled rifting phases.

Lastly, the calibrated model presented in Figure 10 needs 8 phases listed in Table 5 in order to
<figure><img src="image_14.png" /></figure>

**Figure 11.** Comparison of kinematically forced simulation with the calibration data used for North Iberia–Newfoundland conjugate rifted continental margins. (a) Comparison of the architecture of passive margins (Moho and basement positions) between the best-fit simulation and data extracted from seismic reflection/refraction interpreted profile of North Iberia–Newfoundland conjugate rifted margins modified from Sutra et al. [2013]. (b) Subsidence model prediction near ODP Sites 638–639 (Figure 9a) for the best-fit simulation and comparison with subsidence estimates of ODP Sites 638–639 from Boillot et al. [1988] and model 3b in Jeanniot et al. [2016] for initial crustal thicknesses of 30, 35 and 40 km.

<table><caption>**Table 4.** Synthesis of different reconstruction parameters presented in Figure 9b</caption><thead><tr><th>Units/Basins</th><th>x</th><th>w</th><th>d</th><th>L</th></tr></thead><tbody><tr><td>GIBp</td><td>50</td><td>30</td><td>0</td><td>20</td></tr><tr><td>BIBd</td><td>110</td><td>30</td><td>60</td><td>10</td></tr><tr><td>GB</td><td>160</td><td>30</td><td>50</td><td>15</td></tr><tr><td>IB (NB+DGM)</td><td>230/240</td><td>50</td><td>70</td><td>10</td></tr><tr><td>FC</td><td>170</td><td>15</td><td>30</td><td>20</td></tr></tbody></table>

reconstruct the conjugate rifted margins of the Iberia–Newfoundland until break-up (Table 5). Some of them, with large changes in extension velocity, correspond to change in kinematics, others are needed to ensure that the deformation re-localizes at the right place. We have chosen a wide W of P⁰ to P¹ because Tucholke and Sibuet [2007], Mohn et al.

[2015] describe the first phase of extension as distributed over a wide area where the crustal thinning was very small. The first phase P⁰ is necessary to adjust the initial thickness of the crust and to concur to thermal subsidence during the period P¹ but the timing is not geologically significant. P² that starts at 161 Ma, represents the initiation of the second rifting event necessary to adjust the low subsidence rate observed in the subsidence estimates history (Figure 11b). The two next phases, P³ and P⁴, allow adjusting the calibration to the Galicia Interior Bank (Figure 9a) using a rift migration hypothesis different from those of Sutra et al. [2013] or Jeanniot et al. [2016] for the northern profile but compatible with the data. It allows separating GIB in two parts, the GIB deep and GIB plateau. The P⁵ and P⁶ are two short phases required to reconstruct the Galicia Bank and the Flemish Cap. The last kinematic phase P⁷ is needed to create the Newfoundland Basin (NB) and
<figure><img src="image_14.png" /></figure>

**Figure 11.** Comparison of kinematically forced simulation with the calibration data used for North Iberia–Newfoundland conjugate rifted continental margins. (a) Comparison of the architecture of passive margins (Moho and basement positions) between the best-fit simulation and data extracted from seismic reflection/refraction interpreted profile of North Iberia–Newfoundland conjugate rifted margins modified from Sutra et al. [2013]. (b) Subsidence model prediction near ODP Sites 638–639 (Figure 9a) for the best-fit simulation and comparison with subsidence estimates of ODP Sites 638–639 from Boillot et al. [1988] and model 3b in Jeanniot et al. [2016] for initial crustal thicknesses of 30, 35 and 40 km.

<table><caption>**Table 4.** Synthesis of different reconstruction parameters presented in Figure 9b</caption><thead><tr><th>Units/Basins</th><th>x</th><th>w</th><th>d</th><th>L</th></tr></thead><tbody><tr><td>GIBp</td><td>50</td><td>30</td><td>0</td><td>20</td></tr><tr><td>BIBd</td><td>110</td><td>30</td><td>60</td><td>10</td></tr><tr><td>GB</td><td>160</td><td>30</td><td>50</td><td>15</td></tr><tr><td>IB (NB+DGM)</td><td>230/240</td><td>50</td><td>70</td><td>10</td></tr><tr><td>FC</td><td>170</td><td>15</td><td>30</td><td>20</td></tr></tbody></table>

reconstruct the conjugate rifted margins of the Iberia–Newfoundland until break-up (Table 5). Some of them, with large changes in extension velocity, correspond to change in kinematics, others are needed to ensure that the deformation re-localizes at the right place. We have chosen a wide W of P⁰ to P¹ because Tucholke and Sibuet [2007], Mohn et al.

[2015] describe the first phase of extension as distributed over a wide area where the crustal thinning was very small. The first phase P⁰ is necessary to adjust the initial thickness of the crust and to concur to thermal subsidence during the period P¹ but the timing is not geologically significant. P² that starts at 161 Ma, represents the initiation of the second rifting event necessary to adjust the low subsidence rate observed in the subsidence estimates history (Figure 11b). The two next phases, P³ and P⁴, allow adjusting the calibration to the Galicia Interior Bank (Figure 9a) using a rift migration hypothesis different from those of Sutra et al. [2013] or Jeanniot et al. [2016] for the northern profile but compatible with the data. It allows separating GIB in two parts, the GIB deep and GIB plateau. The P⁵ and P⁶ are two short phases required to reconstruct the Galicia Bank and the Flemish Cap. The last kinematic phase P⁷ is needed to create the Newfoundland Basin (NB) and
<table><caption>**Table 5.** Different kinematic parameters for each phase leading to the calibrated model of the Iberia–Newfoundland conjugate rifted margins</caption><thead><tr><th>Phases</th><th>P<sup>0</sup></th><th>P<sup>1</sup></th><th>P<sup>2</sup></th><th>P<sup>3</sup></th><th>P<sup>4</sup></th><th>P<sup>5</sup></th><th>P<sup>6</sup></th><th>P<sup>7</sup></th><th>P<sup>8</sup></th></tr></thead><tbody><tr><td>vx (cm/yr)</td><td>0.8</td><td>0.01</td><td>0.10</td><td>0.25</td><td>0.38</td><td>0.30</td><td>0.40</td><td>0.20</td><td>0.2</td></tr><tr><td>te (Myrs)</td><td>2</td><td>19</td><td>40</td><td>44</td><td>47</td><td>49</td><td>50</td><td>80</td><td>85</td></tr><tr><td>(Ma)</td><td>178</td><td>161</td><td>140</td><td>136</td><td>133</td><td>131</td><td>130</td><td>100</td><td>95</td></tr><tr><td>ti (Myr)</td><td>0</td><td>2</td><td>19</td><td>40</td><td>44</td><td>47</td><td>49</td><td>50</td><td>80</td></tr><tr><td>(Ma)</td><td>180</td><td>178</td><td>161</td><td>140</td><td>136</td><td>133</td><td>131</td><td>130</td><td>100</td></tr><tr><td>X<sub>0</sub> (km)</td><td>705</td><td>705</td><td>705</td><td>860</td><td>800</td><td>750</td><td>650</td><td>680</td><td>680</td></tr><tr><td>W<sub>0</sub> (km)</td><td>200</td><td>200</td><td>200</td><td>30</td><td>30</td><td>30</td><td>15</td><td>50</td><td>200</td></tr><tr><td>E (km)</td><td>16</td><td>1.7</td><td>21</td><td>10</td><td>11.4</td><td>6</td><td>4</td><td>60</td><td>30</td></tr></tbody></table>

Deep Galicia Margin (DGM) conjugate margins. This chosen order agrees with the oceanward-younging age of sediments [Boillot and Winterer, 1988]. In the model, the crustal break-up occurs at 110 Ma, but the continental mantle breaks at around 90 Ma.

## 5.3. Discussion and comparison with previous models

The velocity ($vx$) imposed on the model varies from 0.8 to 0.01 cm/yr with a mean of 0.3 cm/yr all along the simulation until break-up (Table 5). These values are consistent with extensional velocities documented before the break-up in both Nirrengarten et al. [2018] that range from 2.25 cm/yr to 0.2 cm/yr and Sutra et al. [2013] that range from 0.83 cm/yr to 0.25 cm/yr for the Iberian part. However, our calibration significantly differs in terms of subsidence history and while it passes through the grey boxes representing the borehole data at ODP sites 638–639. Notable differences emerge when we compare the predictions in terms of subsidence between our calibration (Figure 11b, red curve) and those published by Jeanniot et al. [2016, model 3b], represented in Figure 11b (black curve). The main difference is the fast uplift and subsidence event associated to the re-localization of the synrift deformation some 100 km east from the early depocenter at around 50 Myr of model time (150 Ma). This large oscillation results from the large crustal strength at the site where the deformation needs to be re-localized in order to fit the crustal thickness. This might be regarded as an artefact of the method, or a wrong choice in the reconstruction of training mainly on the 3D benchmark presented earlier. Yet, this way of re-localizing

the deformation is consistent with the dynamics observed in 3D simulation of rift step-over [Le Pourhiet et al., 2017, Neuharth et al., 2021] or of single oblique rift propagators [Le Pourhiet et al., 2018, Jourdon et al., 2020], and is consistent also with the model of propagation proposed by Nirrengarten et al. [2018].

## 6. Conclusions

This new open-source forward kinematic thermomechanical modelling, associated with a “quick look” approach allows for the calibration and reconstruction of the thermal history of complex rifted margins using a crustal geophysical model and boreholes data. It performs as effectively in adjusting the data as simpler existing methods based solely on kinematics. However, the benchmark of our calibration of a 2D cross-section extracted from a complex 3D thermomechanical model shows that the mechanical feed-back limits the number of possible scenarios during the rifting stage. Our approach is predictive enough for quantitatively constraining the tectono-thermal evolution of the basins by calibrating both the crustal architecture and thermal history of the rifted margins, even when 2D simulations approximate 3D structural problems such as flower structures arising in 2D cross-sections offered by seismic lines. However, if the thermo-mechanical feedback may help to better choose between several possible reconstructions, the solution is non-unique and many flaws still need to be fixed in the future. On the technical side, adding a more quantitative measurement of the misfit of the model with the cross-section. Root means square minimisation of the
<table><caption>**Table 5.** Different kinematic parameters for each phase leading to the calibrated model of the Iberia–Newfoundland conjugate rifted margins</caption><thead><tr><th>Phases</th><th>P<sup>0</sup></th><th>P<sup>1</sup></th><th>P<sup>2</sup></th><th>P<sup>3</sup></th><th>P<sup>4</sup></th><th>P<sup>5</sup></th><th>P<sup>6</sup></th><th>P<sup>7</sup></th><th>P<sup>8</sup></th></tr></thead><tbody><tr><td>vx (cm/yr)</td><td>0.8</td><td>0.01</td><td>0.10</td><td>0.25</td><td>0.38</td><td>0.30</td><td>0.40</td><td>0.20</td><td>0.2</td></tr><tr><td>te (Myrs)</td><td>2</td><td>19</td><td>40</td><td>44</td><td>47</td><td>49</td><td>50</td><td>80</td><td>85</td></tr><tr><td>(Ma)</td><td>178</td><td>161</td><td>140</td><td>136</td><td>133</td><td>131</td><td>130</td><td>100</td><td>95</td></tr><tr><td>ti (Myr)</td><td>0</td><td>2</td><td>19</td><td>40</td><td>44</td><td>47</td><td>49</td><td>50</td><td>80</td></tr><tr><td>(Ma)</td><td>180</td><td>178</td><td>161</td><td>140</td><td>136</td><td>133</td><td>131</td><td>130</td><td>100</td></tr><tr><td>X<sub>0</sub> (km)</td><td>705</td><td>705</td><td>705</td><td>860</td><td>800</td><td>750</td><td>650</td><td>680</td><td>680</td></tr><tr><td>W<sub>0</sub> (km)</td><td>200</td><td>200</td><td>200</td><td>30</td><td>30</td><td>30</td><td>15</td><td>50</td><td>200</td></tr><tr><td>E (km)</td><td>16</td><td>1.7</td><td>21</td><td>10</td><td>11.4</td><td>6</td><td>4</td><td>60</td><td>30</td></tr></tbody></table>

Deep Galicia Margin (DGM) conjugate margins. This chosen order agrees with the oceanward-younging age of sediments [Boillot and Winterer, 1988]. In the model, the crustal break-up occurs at 110 Ma, but the continental mantle breaks at around 90 Ma.

### 5.3. Discussion and comparison with previous models

The velocity ($vx$) imposed on the model varies from 0.8 to 0.01 cm/yr with a mean of 0.3 cm/yr all along the simulation until break-up (Table 5). These values are consistent with extensional velocities documented before the break-up in both Nirrengarten et al. [2018] that range from 2.25 cm/yr to 0.2 cm/yr and Sutra et al. [2013] that range from 0.83 cm/yr to 0.25 cm/yr for the Iberian part. However, our calibration significantly differs in terms of subsidence history and while it passes through the grey boxes representing the borehole data at ODP sites 638–639. Notable differences emerge when we compare the predictions in terms of subsidence between our calibration (Figure 11b, red curve) and those published by Jeanniot et al. [2016, model 3b], represented in Figure 11b (black curve). The main difference is the fast uplift and subsidence event associated to the re-localization of the synrift deformation some 100 km east from the early depocenter at around 50 Myr of model time (150 Ma). This large oscillation results from the large crustal strength at the site where the deformation needs to be re-localized in order to fit the crustal thickness. This might be regarded as an artefact of the method, or a wrong choice in the reconstruction of training mainly on the 3D benchmark presented earlier. Yet, this way of re-localizing

the deformation is consistent with the dynamics observed in 3D simulation of rift step-over [Le Pourhiet et al., 2017, Neuharth et al., 2021] or of single oblique rift propagators [Le Pourhiet et al., 2018, Jourdon et al., 2020], and is consistent also with the model of propagation proposed by Nirrengarten et al. [2018].

## 6. Conclusions

This new open-source forward kinematic thermomechanical modelling, associated with a “quick look” approach allows for the calibration and reconstruction of the thermal history of complex rifted margins using a crustal geophysical model and boreholes data. It performs as effectively in adjusting the data as simpler existing methods based solely on kinematics. However, the benchmark of our calibration of a 2D cross-section extracted from a complex 3D thermomechanical model shows that the mechanical feed-back limits the number of possible scenarios during the rifting stage. Our approach is predictive enough for quantitatively constraining the tectono-thermal evolution of the basins by calibrating both the crustal architecture and thermal history of the rifted margins, even when 2D simulations approximate 3D structural problems such as flower structures arising in 2D cross-sections offered by seismic lines. However, if the thermo-mechanical feedback may help to better choose between several possible reconstructions, the solution is non-unique and many flaws still need to be fixed in the future. On the technical side, adding a more quantitative measurement of the misfit of the model with the cross-section. Root means square minimisation of the
distance between the mechanically modelled Moho and top-basement horizon and the geophysical model [e.g. Lovely et al., 2010] or volumetric strain field comparison [e.g. Durand-Riard et al., 2010] could help the automation of some parts of the calibration. We could test these approaches using additional synthetic cross-sections. Regarding the geological aspects, the most important bias and limitation of our approach is the impossibility to extract melt in order to model underplating and SDR at volcanic margins. This largely impairs the current possibilities of our software at simulating magmatic margins.

## Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The version of ptatin2D code that has been used to produce the results presented here is available publicly at https://bitbucket.org/ptatin/ptatin2d/src/restart/#0d4d8d1677edba59463041e934d9862c7824d921. All python scripts used to post-process the different are available at https://bitbucket.org/ptatin/ptatin2d/src/restart/ or as supplementary data.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Acknowledgements

We are most grateful to TOTAL SA, which funded this post-doctoral work. During the benchmark phase, the calibration was achieved as a “blind-test” by the first author, who did not know from which 3D simulation the third author extracted the 2D profile.

## Supplementary data

Supporting information for this article is available on the journal's website under https://doi.org/10.5802/crgeos.258 or from the author.

# Appendix A. Numerical method

We use the thermo-mechanical numerical code pTatin [May et al., 2014, 2015] in its 2D version [Jourdon et al., 2018b, Perron et al., 2021], where two new components are implemented, respectively a restart module and a kinematic module. This version of the open-source code is available publicly at https://bitbucket.org/ptatin/ptatin2d/src/restart/. The code relies on PETSc library [Balay et al., 2017] to solve conservation of momentum for an incompressible fluid flow described by its velocity vector **v** such as

$$ \nabla \cdot \boldsymbol{\sigma} = \rho \mathbf{g} \quad (A1) $$

$$ \nabla \cdot \mathbf{v} = 0 \quad (A2) $$

with the full stress tensor

$$ \boldsymbol{\sigma} = 2\eta\dot{\boldsymbol{\epsilon}} - P\mathbf{I} \quad (A3) $$

where P is the total pressure, I is the identity matrix, η the viscosity and

$$ \dot{\boldsymbol{\epsilon}} = \frac{1}{2}(\nabla\mathbf{v} + \nabla\mathbf{v}^T) \quad (A4) $$

is the strain rate tensor. pTatin2D uses the finite element method employing a mixed $Q_2 - P_{-1}$ space to discretize the velocity and pressure respectively.

This permits to model accurately the topography with a free surface.

To avoid deformation of the mesh, the lithologies are tracked with an ALE marker in cell approach [May et al., 2015]. Markers are also used to carry lithological properties. The density $ρ_0$ is one of them, but effective density $ρ$ also depends on temperature T and pressure P using Boussinesq approximation, where α and β are the coefficient of thermal expansion and adiabatic compressibility respectively:

$$ \rho = \rho_0(1 - \alpha(T - T_0) + \beta(P - P_0)) \quad (A5) $$

with $P_0$ and $T_0$ the reference temperature and pressure at which $ρ = ρ_0$.

Stress ($\boldsymbol{\sigma}$) and strain rate ($\dot{\boldsymbol{\epsilon}}$) are also computed on markers in two stages. In a first trial, the code evaluates the stress considering the fluid/rocks deform by dislocation creep

$$ \boldsymbol{\sigma} = A_{\text{scale}} A^{-\frac{1}{n}} \exp\left(\frac{Q}{nRT}\right) (\dot{\boldsymbol{\epsilon}}^{\text{II}})^{\frac{1-n}{n}} \dot{\boldsymbol{\epsilon}} \quad (A6) $$
distance between the mechanically modelled Moho and top-basement horizon and the geophysical model [e.g. Lovely et al., 2010] or volumetric strain field comparison [e.g. Durand-Riard et al., 2010] could help the automation of some parts of the calibration. We could test these approaches using additional synthetic cross-sections. Regarding the geological aspects, the most important bias and limitation of our approach is the impossibility to extract melt in order to model underplating and SDR at volcanic margins. This largely impairs the current possibilities of our software at simulating magmatic margins.

## Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The version of ptatin2D code that has been used to produce the results presented here is available publicly at https://bitbucket.org/ptatin/ptatin2d/src/restart/#0d4d8d1677edba59463041e934d9862c7824d921. All python scripts used to post-process the different are available at https://bitbucket.org/ptatin/ptatin2d/src/restart/ or as supplementary data.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Acknowledgements

We are most grateful to TOTAL SA, which funded this post-doctoral work. During the benchmark phase, the calibration was achieved as a “blind-test” by the first author, who did not know from which 3D simulation the third author extracted the 2D profile.

## Supplementary data

Supporting information for this article is available on the journal's website under https://doi.org/10.5802/crgeos.258 or from the author.

# Appendix A. Numerical method

We use the thermo-mechanical numerical code pTatin [May et al., 2014, 2015] in its 2D version [Jourdon et al., 2018b, Perron et al., 2021], where two new components are implemented, respectively a restart module and a kinematic module. This version of the open-source code is available publicly at https://bitbucket.org/ptatin/ptatin2d/src/restart/. The code relies on PETSc library [Balay et al., 2017] to solve conservation of momentum for an incompressible fluid flow described by its velocity vector **v** such as

$$ \nabla \cdot \boldsymbol{\sigma} = \rho \mathbf{g} \qquad (A1) $$

$$ \nabla \cdot \mathbf{v} = 0 \qquad (A2) $$

with the full stress tensor

$$ \boldsymbol{\sigma} = 2\eta\dot{\boldsymbol{\epsilon}} - P\mathbf{I} \qquad (A3) $$

where $P$ is the total pressure, $\mathbf{I}$ is the identity matrix, $\eta$ the viscosity and

$$ \dot{\boldsymbol{\epsilon}} = \frac{1}{2}(\nabla\mathbf{v} + \nabla\mathbf{v}^T) \qquad (A4) $$

is the strain rate tensor. pTatin2D uses the finite element method employing a mixed $Q_2 - P_{-1}$ space to discretize the velocity and pressure respectively.

This permits to model accurately the topography with a free surface.

To avoid deformation of the mesh, the lithologies are tracked with an ALE marker in cell approach [May et al., 2015]. Markers are also used to carry lithological properties. The density $\rho_0$ is one of them, but effective density $\rho$ also depends on temperature $T$ and pressure $P$ using Boussinesq approximation, where $\alpha$ and $\beta$ are the coefficient of thermal expansion and adiabatic compressibility respectively:

$$ \rho = \rho_0(1 - \alpha(T - T_0) + \beta(P - P_0)) \qquad (A5) $$

with $P_0$ and $T_0$ the reference temperature and pressure at which $\rho = \rho_0$.

Stress ($\boldsymbol{\sigma}$) and strain rate ($\dot{\boldsymbol{\epsilon}}$) are also computed on markers in two stages. In a first trial, the code evaluates the stress considering the fluid/rocks deform by dislocation creep

$$ \boldsymbol{\sigma} = A_{\text{scale}} A^{-\frac{1}{n}} \exp\left(\frac{Q}{nRT}\right) (\dot{\boldsymbol{\epsilon}}^{\text{II}})^{\frac{1-n}{n}} \dot{\boldsymbol{\epsilon}} \qquad (A6) $$
and therefore, its effective viscosity depends on temperature $T$, lithology ($A, n, Q$ see Table 1) and the second invariant of strain rate ($\dot{\epsilon}^{II}$). However, when this trial viscous stress exceeds brittle frictional strength

$$ \sigma^b = \sin\phi + C \cos\phi \quad (A7) $$

or maximum plastic strength $\sigma^p$ fixed at 450 MPa [stress/viscosity limiter; Watremez et al., 2013], the effective viscosity is adjusted to

$$ \eta = \frac{\sigma^\gamma}{2\dot{\epsilon}^{\text{II}}} $$

in order to keep the stress on the yield cap defined as $\sigma^\gamma = \min(\sigma^b, \sigma^p)$.

Conservation of momentum is coupled with conservation of heat

$$ \frac{\partial T}{\partial t} = \nabla \cdot \kappa \nabla T + \nu \nabla T + \frac{H}{Cp}. $$

The heat diffusivity $\kappa$, heat production $H$ and heat capacity $Cp$ do not vary for the different simulations.

Sediment transport is simulated using advection diffusion of the topography in 1D

$$ \frac{\partial h}{\partial t} = v_y + \frac{\partial}{\partial x} \left( \kappa_e \frac{\partial h}{\partial x} \right) + S $$

with a source term $S$, which permits out-of-plane sediments inflow and outflow. Details about the implementation and tracking of the stratigraphy may be found in Jourdon et al. [2018a].

## References

Ammann, N., Liao, J., Gerya, T., and Ball, P. (2018). Oblique continental rifting and long transform fault formation based on 3D thermomechanical numerical modeling. *Tectonophysics*, 746, 106–120.

Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L. D., Eijkhout, V., Gropp, W., Kaushik, D., Knepley, M., May, D., McInnes, L. C., Munson, T., Rupp, K., Sanan, P., Smith, B., Zampini, S., Zhang, H., and Zhang, H. (2017). *PETSc Users Manual Revision 3.8*. Argonne National Lab. (ANL), Argonne, IL (United States).

Balázs, A., Matenco, L., Vogt, K., Cloetingh, S., and Gerya, T. (2018). Extensional polarity change in continental rifts: inferences from 3-D numerical modeling and observations. *J. Geophys. Res. Solid Earth*, 123, 8073–8094.

Blaich, O. A., Faleide, J. I., and Tsikalas, F. (2011). Crustal breakup and continent-ocean transition at South Atlantic conjugate margins. *J. Geophys. Res.*, 116, article no. B01402.

Boillot, G., Girardeau, J., and Kornprobst, J. (1988). Rifting of the Galicia margin: crustal thinning and emplacement of mantle rocks on the seafloor. In Boillot, G., Winterer, E. L., et al., editors, *Proceedings of the ODP, Scientific Results*, volume 103, pages 741–756. Ocean Drilling Program, College Station, TX.

Boillot, G. and Winterer, E. L. (1988). Drilling on the Galicia margin: retrospect and prospect. In Boillot, G., Winterer, E. L., et al., editors, *Proceedings of the ODP, Scientific Results*, volume 103, pages 809–828. Ocean Drilling Program, College Station, TX.

Bronner, A., Sauter, D., Manatschal, G., Péron-Pinvidic, G., and Munschy, M. (2011). Magmatic breakup as an explanation for magnetic anomalies at magma-poor rifted margins. *Nat. Geosci.*, 4, 549–553.

Brune, S. and Autin, J. (2013). The rift to break-up evolution of the Gulf of Aden: Insights from 3D numerical lithospheric-scale modelling. *Tectonophysics*, 607, 65–79.

Brune, S., Heine, C., Clift, P. D., and Pérez-Gussinyé, M. (2017). Rifted margin architecture and crustal rheology: Reviewing Iberia-Newfoundland, Central South Atlantic, and South China Sea. *Mar. Pet. Geol.*, 79, 257–281.

Brune, S., Heine, C., Pérez-Gussinyé, M., and Sobolev, S. V. (2014). Rift migration explains continental margin asymmetry and crustal hyper-extension. *Nat. Commun.*, 5, article no. 4014.

Buck, W. R. (1991). Modes of continental lithospheric extension. *J. Geophys. Res. Solid Earth*, 96, 20161–20178.

Burov, E. (2007). The role of gravitational instabilities, density structure and extension rate in the evolution of continental margins. *Geol. Soc. Lond. Spec. Publ.*, 282, 139–156.

Carter, N. L. and Tsenn, M. C. (1987). Flow properties of continental lithosphere. *Tectonophysics*, 136, 27–63.

Chenin, P. and Beaumont, C. (2013). Influence of offset weak zones on the development of rift basins: activation and abandonment during continental extension and breakup. *J. Geophys. Res. Solid Earth*, 118, 1698–1720.
and therefore, its effective viscosity depends on temperature $T$, lithology ($A, n, Q$ see Table 1) and the second invariant of strain rate ($\dot{\epsilon}^{II}$). However, when this trial viscous stress exceeds brittle frictional strength

$$ \sigma^b = \sin\phi + C \cos\phi \quad (A7) $$

or maximum plastic strength $\sigma^p$ fixed at 450 MPa [stress/viscosity limiter; Watremez et al., 2013], the effective viscosity is adjusted to

$$ \eta = \frac{\sigma^\gamma}{2\dot{\epsilon}^{\text{II}}} $$

in order to keep the stress on the yield cap defined as $\sigma^\gamma = \min(\sigma^b, \sigma^p)$.

Conservation of momentum is coupled with conservation of heat

$$ \frac{\partial T}{\partial t} = \nabla \cdot \kappa \nabla T + \nu \nabla T + \frac{H}{C p} $$

The heat diffusivity $\kappa$, heat production $H$ and heat capacity $Cp$ do not vary for the different simulations.

Sediment transport is simulated using advection diffusion of the topography in 1D

$$ \frac{\partial h}{\partial t} = v_y + \frac{\partial}{\partial x} \left( \kappa_e \frac{\partial h}{\partial x} \right) + S $$

with a source term $S$, which permits out-of-plane sediments inflow and outflow. Details about the implementation and tracking of the stratigraphy may be found in Jourdon et al. [2018a].

## References

Ammann, N., Liao, J., Gerya, T., and Ball, P. (2018). Oblique continental rifting and long transform fault formation based on 3D thermomechanical numerical modeling. *Tectonophysics*, 746, 106–120.

Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L. D., Eijkhout, V., Gropp, W., Kaushik, D., Knepley, M., May, D., McInnes, L. C., Munson, T., Rupp, K., Sanan, P., Smith, B., Zampini, S., Zhang, H., and Zhang, H. (2017). *PETSc Users Manual Revision 3.8*. Argonne National Lab. (ANL), Argonne, IL (United States).

Balázs, A., Matenco, L., Vogt, K., Cloetingh, S., and Gerya, T. (2018). Extensional polarity change in continental rifts: inferences from 3-D numerical modeling and observations. *J. Geophys. Res. Solid Earth*, 123, 8073–8094.

Blaich, O. A., Faleide, J. I., and Tsikalas, F. (2011). Crustal breakup and continent-ocean transition at South Atlantic conjugate margins. *J. Geophys. Res.*, 116, article no. B01402.

Boillot, G., Girardeau, J., and Kornprobst, J. (1988). Rifting of the Galicia margin: crustal thinning and emplacement of mantle rocks on the seafloor. In Boillot, G., Winterer, E. L., et al., editors, *Proceedings of the ODP, Scientific Results*, volume 103, pages 741–756. Ocean Drilling Program, College Station, TX.

Boillot, G. and Winterer, E. L. (1988). Drilling on the Galicia margin: retrospect and prospect. In Boillot, G., Winterer, E. L., et al., editors, *Proceedings of the ODP, Scientific Results*, volume 103, pages 809–828. Ocean Drilling Program, College Station, TX.

Bronner, A., Sauter, D., Manatschal, G., Péron-Pinvidic, G., and Munschy, M. (2011). Magmatic breakup as an explanation for magnetic anomalies at magma-poor rifted margins. *Nat. Geosci.*, 4, 549–553.

Brune, S. and Autin, J. (2013). The rift to break-up evolution of the Gulf of Aden: Insights from 3D numerical lithospheric-scale modelling. *Tectonophysics*, 607, 65–79.

Brune, S., Heine, C., Clift, P. D., and Pérez-Gussinyé, M. (2017). Rifted margin architecture and crustal rheology: Reviewing Iberia-Newfoundland, Central South Atlantic, and South China Sea. *Mar. Pet. Geol.*, 79, 257–281.

Brune, S., Heine, C., Pérez-Gussinyé, M., and Sobolev, S. V. (2014). Rift migration explains continental margin asymmetry and crustal hyper-extension. *Nat. Commun.*, 5, article no. 4014.

Buck, W. R. (1991). Modes of continental lithospheric extension. *J. Geophys. Res. Solid Earth*, 96, 20161–20178.

Burov, E. (2007). The role of gravitational instabilities, density structure and extension rate in the evolution of continental margins. *Geol. Soc. Lond. Spec. Publ.*, 282, 139–156.

Carter, N. L. and Tsenn, M. C. (1987). Flow properties of continental lithosphere. *Tectonophysics*, 136, 27–63.

Chenin, P. and Beaumont, C. (2013). Influence of offset weak zones on the development of rift basins: activation and abandonment during continental extension and breakup. *J. Geophys. Res. Solid Earth*, 118, 1698–1720.
Chenin, P., Jammes, S., Lavier, L. L., Manatschal, G., Picazo, S., Müntener, O., Karner, G. D., Figueredo, P. H., and Johnson, C. (2019). Impact of mafic underplating and mantle depletion on subsequent rifting: a numerical modeling study. *Tectonics*, 38, 2185–2207.

Crosby, A., White, N., Edwards, G., and Shillington, D. J. (2008). Evolution of the Newfoundland–Iberia conjugate rifted margins. *Earth Planet. Sci. Lett.*, 273, 214–226.

Culling, W. E. H. (1965). Theory of erosion on soil-covered slopes. *J. Geol.*, 73, 230–254.

Dias, A. E. S., Lavier, L. L., and Hayman, N. W. (2015). Conjugate rifted margins width and asymmetry: the interplay between lithospheric strength and thermomechanical processes. *J. Geophys. Res. Solid Earth*, 120, 8672–8700.

Durand-Riard, P., Caumon, G., and Muron, P. (2010). Balanced restoration of geological volumes with relaxed meshing constraints. *Comput. Geosci.*, 36(4), 441–452.

Duretz, T., Petri, B., Mohn, G., Schmalholz, S. M., Schenker, F. L., and Müntener, O. (2016). The importance of structural softening for the evolution and architecture of passive margins. *Sci. Rep.*, 6, article no. 38704.

Eddy, M. P., Jagoutz, O., and Ibañez-Mejia, M. (2017). Timing of initial seafloor spreading in the Newfoundland-Iberia rift. *Geology*, 45, 527–530.

Goetze, C. and Evans, B. (1979). Stress and temperature in the bending lithosphere as constrained by experimental rock mechanics. *Geophys. J. R. Astron. Soc.*, 59, 463–478.

Gorczyk, W., Willner, A. P., Gerya, T. V., Connolly, J. A. D., and Burg, J.-P. (2007). Physical controls of magmatic productivity at Pacific-type convergent margins: numerical modelling. *Phys. Earth Planet. Inter.*, 163, 209–232.

Gouiza, M. and Naliboff, J. (2021). Rheological inheritance controls the formation of segmented rifted margins in cratonic lithosphere. *Nat. Commun.*, 12, article no. 4653.

Gueydan, F. and Précigout, J. (2014). Modes of continental rifting as a function of ductile strain localization in the lithospheric mantle. *Tectonophysics*, 612, 18–25.

Heine, C. and Brune, S. (2014). Oblique rifting of the equatorial Atlantic: why there is no Saharan Atlantic Ocean. *Geology*, 42, 211–214.

Huet, B., Le Pourhiet, L., Labrousse, L., Burov, E., and Jolivet, L. (2011a). Post-orogenic extension and metamorphic core complexes in a heterogeneous crust: the role of crustal layering inherited from collision. Application to the Cyclades ( Aegean domain). *Geophys. J. Int.*, 184(2), 611–625.

Huet, B., Le Pourhiet, L., Labrousse, L., Burov, E. B., and Jolivet, L. (2011b). Formation of metamorphic core complex in inherited wedges: a thermomechanical modelling study. *Earth Planet. Sci. Lett.*, 309(3–4), 249–257.

Huismans, R. S. and Beaumont, C. (2014). Rifted continental margins: the case for depth-dependent extension. *Earth Planet. Sci. Lett.*, 407, 148–162.

Jeanniot, L., Kusznir, N., Mohn, G., Manatschal, G., and Cowie, L. (2016). Constraining lithosphere deformation modes during continental breakup for the Iberia–Newfoundland conjugate rifted margins. *Tectonophysics*, 680, 28–49.

Jourdon, A., Kergaravat, C., Duclaux, G., and Huguen, C. (2021). Looking beyond kinematics: 3D thermomechanical modelling reveals the dynamics of transform margins. *Solid Earth Discuss.*, 12(5), 1211–1232.

Jourdon, A., Le Pourhiet, L., Mouthereau, F., and May, D. (2020). Modes of propagation of continental breakup and associated oblique rift structures. *J. Geophys. Res. Solid Earth*, 125, article no. e2020JB019906.

Jourdon, A., Le Pourhiet, L., Petit, C., and Rolland, Y. (2018a). Impact of range-parallel sediment transport on 2D thermo-mechanical models of mountain belts: Application to the Kyrgyz Tien Shan. *Terra Nova*, 30, 279–288.

Jourdon, A., Le Pourhiet, L., Petit, C., and Rolland, Y. (2018b). The deep structure and reactivation of the Kyrgyz Tien Shan: modelling the past to better constrain the present. *Tectonophysics*, 746, 530–548. ISSN 0040-1951.

Lavier, L. L. and Manatschal, G. (2006). A mechanism to thin the continental lithosphere at magma-poor margins. *Nature*, 440, 324–328.

Le Pourhiet, L., Burov, E., and Moretti, I. (2004). Rifting through a stack of inhomogeneous thrusts (the dipping pie concept). *Tectonics*, 23, article no. TC4005.

Le Pourhiet, L., Chamot-Rooke, N., Delescluse, M., May, D. A., Watremez, L., and Pubellier, M. (2018). Continental break-up of the South China Sea
Chenin, P., Jammes, S., Lavier, L. L., Manatschal, G., Picazo, S., Müntener, O., Karner, G. D., Figueredo, P. H., and Johnson, C. (2019). Impact of mafic underplating and mantle depletion on subsequent rifting: a numerical modeling study. *Tectonics*, 38, 2185–2207.

Crosby, A., White, N., Edwards, G., and Shillington, D. J. (2008). Evolution of the Newfoundland–Iberia conjugate rifted margins. *Earth Planet. Sci. Lett.*, 273, 214–226.

Culling, W. E. H. (1965). Theory of erosion on soil-covered slopes. *J. Geol.*, 73, 230–254.

Dias, A. E. S., Lavier, L. L., and Hayman, N. W. (2015). Conjugate rifted margins width and asymmetry: the interplay between lithospheric strength and thermomechanical processes. *J. Geophys. Res. Solid Earth*, 120, 8672–8700.

Durand-Riard, P., Caumon, G., and Muron, P. (2010). Balanced restoration of geological volumes with relaxed meshing constraints. *Comput. Geosci.*, 36(4), 441–452.

Duretz, T., Petri, B., Mohn, G., Schmalholz, S. M., Schenker, F. L., and Müntener, O. (2016). The importance of structural softening for the evolution and architecture of passive margins. *Sci. Rep.*, 6, article no. 38704.

Eddy, M. P., Jagoutz, O., and Ibañez-Mejia, M. (2017). Timing of initial seafloor spreading in the Newfoundland-Iberia rift. *Geology*, 45, 527–530.

Goetze, C. and Evans, B. (1979). Stress and temperature in the bending lithosphere as constrained by experimental rock mechanics. *Geophys. J. R. Astron. Soc.*, 59, 463–478.

Gorczyk, W., Willner, A. P., Gerya, T. V., Connolly, J. A. D., and Burg, J.-P. (2007). Physical controls of magmatic productivity at Pacific-type convergent margins: numerical modelling. *Phys. Earth Planet. Inter.*, 163, 209–232.

Gouiza, M. and Naliboff, J. (2021). Rheological inheritance controls the formation of segmented rifted margins in cratonic lithosphere. *Nat. Commun.*, 12, article no. 4653.

Gueydan, F. and Précigout, J. (2014). Modes of continental rifting as a function of ductile strain localization in the lithospheric mantle. *Tectonophysics*, 612, 18–25.

Heine, C. and Brune, S. (2014). Oblique rifting of the equatorial Atlantic: why there is no Saharan Atlantic Ocean. *Geology*, 42, 211–214.

Huet, B., Le Pourhiet, L., Labrousse, L., Burov, E., and Jolivet, L. (2011a). Post-orogenic extension and metamorphic core complexes in a heterogeneous crust: the role of crustal layering inherited from collision. Application to the Cyclades ( Aegean domain). *Geophys. J. Int.*, 184(2), 611–625.

Huet, B., Le Pourhiet, L., Labrousse, L., Burov, E. B., and Jolivet, L. (2011b). Formation of metamorphic core complex in inherited wedges: a thermomechanical modelling study. *Earth Planet. Sci. Lett.*, 309(3–4), 249–257.

Huismans, R. S. and Beaumont, C. (2014). Rifted continental margins: the case for depth-dependent extension. *Earth Planet. Sci. Lett.*, 407, 148–162.

Jeanniot, L., Kusznir, N., Mohn, G., Manatschal, G., and Cowie, L. (2016). Constraining lithosphere deformation modes during continental breakup for the Iberia–Newfoundland conjugate rifted margins. *Tectonophysics*, 680, 28–49.

Jourdon, A., Kergaravat, C., Duclaux, G., and Huguen, C. (2021). Looking beyond kinematics: 3D thermomechanical modelling reveals the dynamics of transform margins. *Solid Earth Discuss.*, 12(5), 1211–1232.

Jourdon, A., Le Pourhiet, L., Mouthereau, F., and May, D. (2020). Modes of propagation of continental breakup and associated oblique rift structures. *J. Geophys. Res. Solid Earth*, 125, article no. e2020JB019906.

Jourdon, A., Le Pourhiet, L., Petit, C., and Rolland, Y. (2018a). Impact of range-parallel sediment transport on 2D thermo-mechanical models of mountain belts: Application to the Kyrgyz Tien Shan. *Terra Nova*, 30, 279–288.

Jourdon, A., Le Pourhiet, L., Petit, C., and Rolland, Y. (2018b). The deep structure and reactivation of the Kyrgyz Tien Shan: modelling the past to better constrain the present. *Tectonophysics*, 746, 530–548. ISSN 0040-1951.

Lavier, L. L. and Manatschal, G. (2006). A mechanism to thin the continental lithosphere at magma-poor margins. *Nature*, 440, 324–328.

Le Pourhiet, L., Burov, E., and Moretti, I. (2004). Rifting through a stack of inhomogeneous thrusts (the dipping pie concept). *Tectonics*, 23, article no. TC4005.

Le Pourhiet, L., Chamot-Rooke, N., Delescluse, M., May, D. A., Watremez, L., and Pubellier, M. (2018). Continental break-up of the South China Sea
stalled by far-field compression. *Nat. Geosci.*, 11, 605–609.

Le Pourhiet, L., May, D. A., Huille, L., Watremez, L., and Leroy, S. (2017). A genetic link between transform and hyper-extended margins. *Earth Planet. Sci. Lett.*, 465, 184–192.

Liao, J. and Gerya, T. (2015). From continental rifting to seafloor spreading: insight from 3D thermomechanical modeling. *Gondwana Res.*, 28(4), 1329–1343.

Lovely, P., Zahasky, C., and Pollard, D. D. (2010). Fold geometry at Sheep Mountain anticline, Wyoming, constructed using airborne laser swath mapping data, outcrop-scale geologic mapping, and numerical interpolation. *J. Geophys. Res.*, 115, article no. B12414.

Manatschal, G., Lavier, L., and Chenin, P. (2015). The role of inheritance in structuring hyperextended rift systems: some considerations based on observations and numerical modeling. *Gondwana Res.*, 27, 140–164.

Manatschal, G., Müntener, O., Lavier, L. L., Minshull, T. A., and Péron-Pinvidic, G. (2007). Observations from the Alpine Tethys and Iberia–Newfoundland margins pertinent to the interpretation of continental breakup. *Geol. Soc. Lond. Spec. Publ.*, 282, 291–324.

May, D. A., Brown, J., and Le Pourhiet, L. (2015). A scalable, matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow. *Comput. Methods Appl. Mech. Eng.*, 290, 496–523.

May, D. A., Brown, J., and Pourhiet, L. L. (2014). pTatin3D: high-performance methods for long-term lithospheric dynamics. In *SC'14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*, pages 274–284. IEEE, New York.

McKenzie, D. (1978). Some remarks on the development of sedimentary basins. *Earth Planet. Sci. Lett.*, 40, 25–32.

Mohn, G., Karner, G. D., Manatschal, G., and Johnson, C. A. (2015). Structural and stratigraphic evolution of the Iberia–Newfoundland hyper-extended rifted margin: a quantitative modelling approach. *Geol. Soc. Lond. Spec. Publ.*, 413, 53–89.

Nagel, T. J. and Buck, W. R. (2004). Symmetric alternative to asymmetric rifting models. *Geology*, 32, 937–940.

Naliboff, J. and Buiter, S. J. H. (2015). Rift reactivation and migration during multiphase extension. *Earth Planet. Sci. Lett.*, 421, 58–67.

Neuharth, D., Brune, S., Glerum, A., Heine, C., and Welford, J. K. (2021). Formation of continental microplates through rift linkage: numerical modelling and its application to the Flemish Cap and Sao Paulo Plateau. *Geochem. Geophys. Geosyst.*, 22(4), article no. e2020GC009615.

Nirrengarten, M., Manatschal, G., Tugend, J., Kusznir, N., and Sauter, D. (2018). Kinematic evolution of the Southern North Atlantic: implications for the formation of hyperextended rift systems. *Tectonics*, 37, 89–118.

Pérez-Gussinyé, M., Morgan, J. P., Reston, T. J., and Ranero, C. R. (2006). The rift to drift transition at non-volcanic margins: Insights from numerical modelling. *Earth Planet. Sci. Lett.*, 244, 458–473.

Péron-Pinvidic, G. and Manatschal, G. (2009). The final rifting evolution at deep magma-poor passive margins from Iberia-Newfoundland: a new point of view. *Int. J. Earth Sci.*, 98, 1581–1597.

Péron-Pinvidic, G., Manatschal, G., Minshull, T. A., and Sawyer, D. S. (2007). Tectonosedimentary evolution of the deep Iberia-Newfoundland margins: Evidence for a complex breakup history. *Tectonics*, 26, article no. TC2011.

Perron, P., Pourhiet, L. L., Guiraud, M., Vennin, E., Moretti, I., Portier, É., and Konaté, M. (2021). Control of inherited accreted lithospheric heterogeneity on the architecture and the low, long-term subsidence rate of intracratonic basins. *Bull. Soc. Géol. Fr. - Earth Sci. Bull.*, 192, article no. 15.

Petri, B., Duretz, T., Mohn, G., Schmalholz, S. M., Karner, G. D., and Müntener, O. (2019). Thinning mechanisms of heterogeneous continental lithosphere. *Earth Planet. Sci. Lett.*, 512, 147–162.

Ranalli, G. and Murphy, D. C. (1987). Rheological stratification of the lithosphere. *Tectonophysics*, 132, 281–295.

Sapin, F., Ringenbach, J.-C., and Clerc, C. (2021). Rifted margins classification and forcing parameters. *Sci. Rep.*, 11, 1–17.

Seton, M., Müller, R. D., Zahirovic, S., Gaina, C., Torsvik, T., Shephard, G., Talsma, A., Gurnis, M., Turner, M., Maus, S., and Chandler, M. (2012). Global continental and ocean basin reconstructions since 200 Ma. *Earth-Sci. Rev.*, 113, 212–270.
stalled by far-field compression. *Nat. Geosci.*, 11, 605–609.

Le Pourhiet, L., May, D. A., Huille, L., Watremez, L., and Leroy, S. (2017). A genetic link between transform and hyper-extended margins. *Earth Planet. Sci. Lett.*, 465, 184–192.

Liao, J. and Gerya, T. (2015). From continental rifting to seafloor spreading: insight from 3D thermomechanical modeling. *Gondwana Res.*, 28(4), 1329–1343.

Lovely, P., Zahasky, C., and Pollard, D. D. (2010). Fold geometry at Sheep Mountain anticline, Wyoming, constructed using airborne laser swath mapping data, outcrop- scale geologic mapping, and numerical interpolation. *J. Geophys. Res.*, 115, article no. B12414.

Manatschal, G., Lavier, L., and Chenin, P. (2015). The role of inheritance in structuring hyperextended rift systems: some considerations based on observations and numerical modeling. *Gondwana Res.*, 27, 140–164.

Manatschal, G., Müntener, O., Lavier, L. L., Minshull, T. A., and Péron-Pinvidic, G. (2007). Observations from the Alpine Tethys and Iberia–Newfoundland margins pertinent to the interpretation of continental breakup. *Geol. Soc. Lond. Spec. Publ.*, 282, 291–324.

May, D. A., Brown, J., and Le Pourhiet, L. (2015). A scalable, matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow. *Comput. Methods Appl. Mech. Eng.*, 290, 496–523.

May, D. A., Brown, J., and Pourhiet, L. L. (2014). pTatin3D: high-performance methods for long-term lithospheric dynamics. In *SC'14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*, pages 274–284. IEEE, New York.

McKenzie, D. (1978). Some remarks on the development of sedimentary basins. *Earth Planet. Sci. Lett.*, 40, 25–32.

Mohn, G., Karner, G. D., Manatschal, G., and Johnson, C. A. (2015). Structural and stratigraphic evolution of the Iberia–Newfoundland hyper-extended rifted margin: a quantitative modelling approach. *Geol. Soc. Lond. Spec. Publ.*, 413, 53–89.

Nagel, T. J. and Buck, W. R. (2004). Symmetric alternative to asymmetric rifting models. *Geology*, 32, 937–940.

Naliboff, J. and Buiter, S. J. H. (2015). Rift reactivation and migration during multiphase extension. *Earth Planet. Sci. Lett.*, 421, 58–67.

Neuharth, D., Brune, S., Glerum, A., Heine, C., and Welford, J. K. (2021). Formation of continental microplates through rift linkage: numerical modelling and its application to the Flemish Cap and Sao Paulo Plateau. *Geochem. Geophys. Geosyst.*, 22(4), article no. e2020GC009615.

Nirrengarten, M., Manatschal, G., Tugend, J., Kusznir, N., and Sauter, D. (2018). Kinematic evolution of the Southern North Atlantic: implications for the formation of hyperextended rift systems. *Tectonics*, 37, 89–118.

Pérez-Gussinyé, M., Morgan, J. P., Reston, T. J., and Ranero, C. R. (2006). The rift to drift transition at non-volcanic margins: Insights from numerical modelling. *Earth Planet. Sci. Lett.*, 244, 458–473.

Péron-Pinvidic, G. and Manatschal, G. (2009). The final rifting evolution at deep magma-poor passive margins from Iberia-Newfoundland: a new point of view. *Int. J. Earth Sci.*, 98, 1581–1597.

Péron-Pinvidic, G., Manatschal, G., Minshull, T. A., and Sawyer, D. S. (2007). Tectonosedimentary evolution of the deep Iberia-Newfoundland margins: Evidence for a complex breakup history. *Tectonics*, 26, article no. TC2011.

Perron, P., Pourhiet, L. L., Guiraud, M., Vennin, E., Moretti, I., Portier, É., and Konaté, M. (2021). Control of inherited accreted lithospheric heterogeneity on the architecture and the low, long-term subsidence rate of intracratonic basins. *Bull. Soc. Géol. Fr. - Earth Sci. Bull.*, 192, article no. 15.

Petri, B., Duretz, T., Mohn, G., Schmalholz, S. M., Karner, G. D., and Müntener, O. (2019). Thinning mechanisms of heterogeneous continental lithosphere. *Earth Planet. Sci. Lett.*, 512, 147–162.

Ranalli, G. and Murphy, D. C. (1987). Rheological stratification of the lithosphere. *Tectonophysics*, 132, 281–295.

Sapin, F., Ringenbach, J.-C., and Clerc, C. (2021). Rifted margins classification and forcing parameters. *Sci. Rep.*, 11, 1–17.

Seton, M., Müller, R. D., Zahirovic, S., Gaina, C., Torsvik, T., Shephard, G., Talsma, A., Gurnis, M., Turner, M., Maus, S., and Chandler, M. (2012). Global continental and ocean basin reconstructions since 200 Ma. *Earth-Sci. Rev.*, 113, 212–270.
Spadini, G. and Podladchikov, Y. (1996). Spacing of consecutive normal faulting in the lithosphere: A dynamic model for rift axis jumping (Tyrrhenian Sea). *Earth Planet. Sci. Lett.*, 144(1–2), 21–34.

Sutra, E., Manatschal, G., Mohn, G., and Unternehr, P. (2013). Quantification and restoration of extensional deformation along the Western Iberia and Newfoundland rifted margins. *Geochem. Geophys. Geosyst.*, 14, 2575–2597.

Taylor, B., Goodliffe, A. M., and Martinez, F. (1999). How continents break up: insights from Papua New Guinea. *J. Geophys. Res. Solid Earth*, 104, 7497–7512.

Tetreault, J. L. and Buiter, S. J. H. (2018). The influence of extension rate and crustal rheology on the evolution of passive margins from rifting to break-up. *Tectonophysics*, 746, 155–172.

Tucholke, B. E., Sawyer, D. S., and Sibuet, J.-C. (2007). Breakup of the Newfoundland–Iberia rift. *Geol. Soc. Lond. Spec. Publ.*, 282, 9–46.

Tucholke, B. E. and Sibuet, J.-C. (2007). Leg 210 synthesis: tectonic, magmatic, and sedimentary

evolution of the Newfoundland-Iberia rift. In Tucholke, B. E., Sibuet, J.-C., and Klaus, A., editors, *Proceedings of the ODP, Scientific Results*, volume 210, pages 1–56. Ocean Drilling Program, College Station, TX.

Vine, F. J. and Matthews, D. H. (1963). Magnetic anomalies over oceanic ridges. *Nature*, 199, 947–949.

Watremez, L., Burov, E., d'Acremont, E., Leroy, S., Huet, B., Le Pourhiet, L., and Bellahsen, N. (2013). Buoyancy and localizing properties of continental mantle lithosphere: insights from thermomechanical models of the eastern Gulf of Aden. *Geochem. Geophys. Geosyst.*, 14, 2800–2817.

Weinberg, R. F. and Podladchikov, Y. (1994). Diapiric ascent of magmas through power law crust and mantle. *J. Geophys. Res.: Solid Earth*, 99(B5), 9543–9559.

Wenker, S. and Beaumont, C. (2018). Effects of lateral strength contrasts and inherited heterogeneities on necking and rifting of continents. *Tectono-physics*, 746, 46–63.
Spadini, G. and Podladchikov, Y. (1996). Spacing of consecutive normal faulting in the lithosphere: A dynamic model for rift axis jumping (Tyrrhenian Sea). *Earth Planet. Sci. Lett.*, 144(1–2), 21–34.

Sutra, E., Manatschal, G., Mohn, G., and Unternehr, P. (2013). Quantification and restoration of extensional deformation along the Western Iberia and Newfoundland rifted margins. *Geochem. Geophys. Geosyst.*, 14, 2575–2597.

Taylor, B., Goodliffe, A. M., and Martinez, F. (1999). How continents break up: insights from Papua New Guinea. *J. Geophys. Res. Solid Earth*, 104, 7497–7512.

Tetreault, J. L. and Buiter, S. J. H. (2018). The influence of extension rate and crustal rheology on the evolution of passive margins from rifting to break-up. *Tectonophysics*, 746, 155–172.

Tucholke, B. E., Sawyer, D. S., and Sibuet, J.-C. (2007). Breakup of the Newfoundland–Iberia rift. *Geol. Soc. Lond. Spec. Publ.*, 282, 9–46.

Tucholke, B. E. and Sibuet, J.-C. (2007). Leg 210 synthesis: tectonic, magmatic, and sedimentary

evolution of the Newfoundland-Iberia rift. In Tucholke, B. E., Sibuet, J.-C., and Klaus, A., editors, *Proceedings of the ODP, Scientific Results*, volume 210, pages 1–56. Ocean Drilling Program, College Station, TX.

Vine, F. J. and Matthews, D. H. (1963). Magnetic anomalies over oceanic ridges. *Nature*, 199, 947–949.

Watremez, L., Burov, E., d'Acremont, E., Leroy, S., Huet, B., Le Pourhiet, L., and Bellahsen, N. (2013). Buoyancy and localizing properties of continental mantle lithosphere: insights from thermomechanical models of the eastern Gulf of Aden. *Geochem. Geophys. Geosyst.*, 14, 2800–2817.

Weinberg, R. F. and Podladchikov, Y. (1994). Diapiric ascent of magmas through power law crust and mantle. *J. Geophys. Res.: Solid Earth*, 99(B5), 9543–9559.

Wenker, S. and Beaumont, C. (2018). Effects of lateral strength contrasts and inherited heterogeneities on necking and rifting of continents. *Tectono-physics*, 746, 46–63.