# Crustal-Scale Signatures of Steady-State Thermal Inheritance: Insights from the South China Sea

Laetitia Le Pourhiet, Manuel Pubellier, Anthony Jourdon, Thomas Francois

► To cite this version:

Laetitia Le Pourhiet, Manuel Pubellier, Anthony Jourdon, Thomas Francois. Crustal-Scale Signatures of Steady-State Thermal Inheritance: Insights from the South China Sea. Geodynamica, In press, not per reviewed. hal-05322955

HAL Id: hal-05322955

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

Submitted on 21 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

# Crustal-Scale Signatures of Steady-State Thermal

## Inheritance: Insights from the South China Sea

Laetitia Le Pourhiet ID*<fnref n="1" />,<fnref n="2" />, Manuel Pubellier ID<fnref n="3" />, Anthony Jourdon ID<fnref n="1" />,<fnref n="4" />, Thomas Francois ID<fnref n="1" />

## Abstract

Long-lived lateral variations in radiogenic heat production create persistent thermal heterogeneities that shape continental lithosphere over geological timescales. We introduce a steady-state concept of thermal inheritance, linking these variations to crustal-scale strain localization and tectonic architecture. Using numerical models, we explore both crustal- and lithospheric-scale consequences of heterogeneous heat production. A key finding is that lateral variations in heat production leave a distinct crustal-scale tectonic signature, controlling patterns of strain localization. The South China Sea serves as a proof-of-concept: the segmented, oblique extension observed there aligns with zones of mechanically weaker crust, reflecting the underlying inherited thermal heterogeneity. These results highlight that crustal-scale tectonic features can emerge from steady-state thermal conditions, independently of transient anomalies. They provide a quantitative framework linking inherited thermal structure to observable deformation patterns. More broadly, our study suggests that laterally heterogeneous heat production offers a physically motivated alternative to traditional exponential-decay models, better capturing the spatial complexity and persistence of lithospheric thermal structure. By emphasizing the crustal imprint of thermal inheritance, we demonstrate that radiogenic heat variations are a fundamental control on strain localization and tectonic segmentation. This approach opens a new perspective on how long-lived thermal heterogeneities shape continental deformation and the architecture of lithospheric structures over hundreds of millions of years.

<footnote n="1">¹Institut des Sciences de la Terre de Paris, Sorbonne Université, Paris, France</footnote>

<footnote n="2">²Institut Universitaire de France</footnote>

<footnote n="3">³Laboratoire de Géologie, Ecole Normale Supérieure, PSL, Paris, France</footnote>

<footnote n="4">⁴GeoAzur, Université de cote d'Azur, Nice, France</footnote>

<footnote>* laetitia.le_pourhiet@sorbonne-universite.fr</footnote>
# 1 Structural vs. Thermal Inheritance: Distinct Controls on Lithospheric Strain Localization

Strain localization in layered or foliated media occurs through shear instabilities such as necking, folding, and buckling, operating across scales from millimetric to lithospheric (e.g., Schmalholz and Mancktelow, 2016; Gerbault et al., 1999). This mechanism is sometimes referred to as structural softening (e.g. Duretz et al., 2016). At outcrop scale, layering reflects sedimentary processes or pre-existing metamorphic fabrics, whereas at lithospheric scale, mechanical layering is largely governed by vertical thermal gradients and crustal thickness, which control the flow capacity of the lower crust and the mechanical coupling between crust and mantle (Burov and Diament, 1995).

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

**FIGURE 1** – Plutonic belts and rift orientations in a. the European Cenozoic Rift System (geological information based on Maierová et al. (2016)) and b. Southeast Asia (geological information based on Pubellier and Morley (2014)). Colored areas indicate zones of elevated radiogenic heat production associated with post-orogenic processes or former volcanic arcs. Light yellow shows the location and orientation of rift structures. In both cases, inherited thermal heterogeneity is oblique to the main extension direction, represented by thick white arrows, motivating our investigation of how lateral variations in heat production influence rift localization, orientation, and their interaction with propagating extension. Both basemaps were produced using GeoMapApp (www.geomapapp.org). Ryan et al. (2009).

Shear instabilities require perturbations to the mechanical layering; perfectly horizontal layers remain stable. Linear stability analyses (Smith, 1977) demonstrate that the fastest-growing instability wavelength is set by the contrast in mechanical strength, layer thickness, and perturbation amplitude. While idealized white-noise perturbations predict a dominant wavelength, finite-amplitude perturbations arising from structural or thermal inheritance can accelerate growth at specific wavelengths and influence the resulting deformation pattern.

At crustal scales, structural inheritance provides pre-existing zones of weakness, including weak or dipping tectonic units (Le Pourhiet et al., 2004; Jammes and Lavier, 2018; Zhou et al., submitted), dis-
45

46

47

continuous structures (Petri et al., 2019), or contacts between lithologies of contrasting rheology (Huet et al., 2011). These features are often transient: they localize strain during rifting but may be erased or reworked by subsequent tectonic events.

48

49

50

51

52

53

54

55

56

57

58

59

60

61

At lithospheric scale, tectonic inheritance is expressed as lithospheric scars, i.e., weak suture zones (e.g., Heron et al., 2016; Balázs et al., 2017, 2018; Jourdon et al., 2018; Heron et al., 2019). These structures are generally assumed to be long-lived, owing to factors such as preserved grain size, water content, or compositional contrasts inherited from past tectonic events. Beyond discrete scars, variations in lithospheric thickness are often introduced using different initial geotherms, commonly by adjusting the depth of the 1300 °C isotherm for different terranes. In two-dimensional models, such variations can completely dominate strain localization at lithospheric scale (e.g., Burov et al., 2007), whereas in three dimensions they influence the orientation and geometry of newly formed faults or sedimentary bodies and the distribution of strain within the lithosphere (e.g., Autin et al., 2013; Brune et al., 2017). These thickness variations are effectively a form of thermal inheritance, although not always explicitly named as such. In the literature, structural inheritance is generally considered long-lived at all scales (Bercovici and Ricard, 2014; Manatschal et al., 2015), while thermal inheritance is often treated as lithospheric-scale and transient, controlled by thermal age (Vacherat et al., 2014; Manatschal et al., 2015; Oravecz et al., 2024, 2025).

62

63

64

65

66

67

68

69

Here, we expand the concept of thermal inheritance to include a steady-state component: lateral variations in crustal radiogenic heat production. Unlike transient geotherms used in most previous studies, these lateral variations persist over plate-tectonic timescales, providing long-lived heterogeneities that influence lithospheric rheology and strain localization (Gouiza and Naliboff, 2021; Perron et al., 2021). By incorporating radiogenic heat as a steady-state component, we bridge the gap between transient thermal effects and structural inheritance: structural inheritance provides geometrically fixed weaknesses, while steady-state thermal inheritance provides long-lived rheological variations guiding the growth of shear instabilities across mobile belts.

70

71

72

73

74

75

76

77

Continental lithospheres often display lateral variations in crustal composition and radiogenic heat production, reflecting the presence of plutonic belts, metamorphic units, or former volcanic arcs. Figure 1 illustrates the distribution of plutonic belts relative to rift structures in Southeast Asia (based on Pubellier and Morley (2014)) and in western Europe (based on Maierová et al. (2016)). In both regions, inherited thermal heterogeneity, expressed as zones of elevated radiogenic heat production, is oblique to the direction of extension. These configurations motivated us to test whether such long-lived thermal perturbations influence rift localization, orientation, and propagation, and how they interact with other geodynamic controls such as slab-pull-driven extension.
The South China Sea (SCS) region, in particular, provides an exceptionally well-documented example of upper crust intruded by numerous Triassic to Late Cretaceous plutons (Yan et al., 2010; Huang et al., 2017; Cao et al., 2018), concentrated primarily along the NW margin (Savva et al., 2014), forming a belt of elevated radiogenic heat production. In this study, we use this regional example to motivate lateral heat-production patterns in numerical models, where plutonic belts are represented as zones of enhanced radiogenic heat in the upper crust, thereby providing a steady-state component of thermal inheritance.

## 2 Geotherm of continent 102: from 1D to 2D and the concept of Thermal Inheritance

<figure><img src="image_3.png" /><figcaption>**FIGURE 2** – Illustration of 1D continental geotherms highlighting the effect of crustal radiogenic heat production. Pink: upper crust, blue: lower crust, green: mantle lithosphere. (a) Variation of LAB depth with thermal age A; $\delta T_{tt}$ represents the transient deviation from the steady-state solution. Geotherm are colored by age (b) Variation of LAB depth with lateral heat production at steady state, showing the contribution of crustal radiogenic heating ($\Delta T_M$) to Moho temperature. Geotherm are colored by $\Delta T_M$. The vertical double dashed lines corresponds respectively to 95% of $T_{LAB}$ and $T_{LAB}$. The depth of the LAB on the lithospheric column is represented for 95% of $T_{LAB}$.</figcaption></figure>

### 2.1 1D Continental Geotherms

Classically, long-term tectonic modeling of the continental lithosphere uses the 1D analytical solution of the heat equation to compute temperature as a function of depth down to the lithosphere-asthenosphere boundary (LAB), typically taken at 1300 °C (McKenzie, 1978). Following Burov and Diament (1995), a plate of thickness $h_l = 100–250$ km is considered, with thermal age A of the last thermo-tectonic event defining a non-dimensional time:

$$ \tau = \kappa \frac{A}{h_l^2}, \qquad (1) $$

used to compute the transient plate-cooling component:
$$ \delta T_{tt} = \frac{1}{2} \sum_{n=1}^{\infty} \sin \left( n\pi \frac{z}{h_l} \right) \frac{(-1)^n}{n} e^{-(n\pi)^2 \tau}, \qquad (2) $$

so that the total transient geotherm is

$$ T = T_{ss} + \delta T_{tt}. \qquad (3) $$

Here, $T_{ss}$ is the steady-state geotherm as in Turcotte and Schubert (2002), assuming either a basal heat flux $q_m$ or an LAB depth $z_{LAB}$, linked by

$$ z_{LAB} = \frac{T_{LAB} - \Delta T_M - T_s}{q_m/k}, \qquad (4) $$

where $T_{LAB}$ is conventionally 1300°C, $T_s$ is surface temperature, $k$ the thermal conductivity, and $\Delta T_M$ the contribution of crustal radiogenic heat production to Moho temperature:

$$ T_M = \frac{\int_{z_M}^{z_s} \rho H dy}{k}, \qquad (5) $$

with $H$ the volumetric heat production, $\rho$ the density, $y_M$ the Moho depth, and $y_s$ the surface. $H(y)$ may be constant (parabolic solution) or decay exponentially with depth (exponential solution) (Turcotte and Schubert, 2002); the integral defines the total Moho contribution. Steady-state geotherms often assign different $\Delta T_M$ values to lithospheric blocks, effectively “stitching” independent 1D solutions and neglecting lateral heat diffusion following McKenzie (1978).

Figure 2 shows two 1D geotherms. Panel a illustrates the transient solution, where $\delta T_{tt}$ decays with age, converging to steady state after 100–150 Myr for $h_l \approx 100$ km. Panel b shows the steady-state case: lateral variations in radiogenic heating ($\Delta T_M$) modify Moho temperature and LAB depth, producing long-lived lithospheric strength contrasts. Our work focuses on this first-order effect of steady-state thermal inheritance as a control on continental rift evolution.

## 2.2 Effects of Lateral Heat Conduction

The robust linear correlation observed between surface heat flow and surface heat production is a fundamental feature of continental crust. Lachenbruch (Lachenbruch, 1968) proposed that this relationship could be explained with a 1D solution if heat production decreases exponentially with depth:
$$H(z) = H_0 e^{-z/h_r}, \qquad (6)$$

where $H_0$ is surface heat production and $h_r$, the decay depth. He suggested that differential erosion exposes deeper levels with lower production, leading to reduced surface heat flow in the Sierra Nevada region. His model relied on a limited dataset from four batholith locations. More recent data (Brady et al., 2006), however, show that radiogenic heat production initially increases down to 5 km before declining, indicating that the simple exponential decay assumption is inconsistent with observations. The classical emphasis on Lachenbruch's model in textbooks (e.g., Turcotte and Schubert, 2002) originally written in 1982 may therefore be overstated given these newer findings.

As an alternative, Jaupart (1983) proposed that lateral heat conduction—considered in models with more than one dimension—can naturally generate a near-linear correlation between surface heat flow and surface heat production. Lateral conduction acts as a low-pass filter on temperature variations, damping short-wavelength anomalies exponentially with distance from the source. For a perturbation of wavelength $λ$, the characteristic decay length is

$$d_c = \frac{\lambda}{2\pi}, \qquad (7)$$

so that anomalies with size smaller than $3d_c$ are strongly attenuated at depth. For example, the 20 km perturbation size case, which results in a sinusoidal wavelength $λ = 40$ km (Figure 3a), is drastically smoothed by diffusion, whereas wider perturbations (50–100 km) can produce enduring Moho temperature anomalies of 5–50°C.

This mechanism explains why the observed linear correlation is robust without invoking a specific vertical decay of heat production. Considering 2D or 3D lateral diffusion thus provides a more realistic description of the geotherm than stitching together independent 1D solutions, which neglect horizontal heat transfer, and underpins the study of lateral heat production variations presented in this paper.

## 2.3 Thermo-mechanical coupling: a 2D proof-of-concept

To illustrate how steady-state thermal inheritance organizes strain, we perform simple 2D thermo-mechanical experiments. The domain is 1000 km wide and 250 km deep, with a 40 km homogeneous crust overly-ing dry olivine mantle. Two crustal rheologies are tested: Strong and Weak highlighting the effect of crust-mantle coupling (Figure 3b). Flow parameters are listed in Table 1 and method in found in the annex A

Lateral perturbations of radiogenic heat production are applied in the upper 20 km of the crust ($4 \times H_0$).
<figure><img src="image_4.png" /></figure>

**FIGURE 3** – Effect of lateral variations in crustal radiogenic heat production on lithospheric temperature and strain in 2D models. (a) Surface heat flow (top) and Moho temperature (bottom) for a homogeneous crust (solid black), elevated heat production (dashed black), and perturbation size of 20, 50, 100 km. The 20 km perturbation size case illustrates strong attenuation by lateral diffusion. (b) Initial viscosity structure for strong (mafic) and weak (felsic) crust columns, highlighting lateral rheological variations induced by temperature anomalies. (c,d) Strain rate after 2 Myr for strong and weak crust, respectively, showing that lateral thermal anomalies control strain localization patterns and intensity. Narrow perturbations are strongly attenuated at depth, while wider perturbations significantly influence lithospheric rheology, demonstrating the role of lateral heat production as a steady-state component of thermal inheritance.

Although limited to the upper crust, these perturbations generate temperature anomalies throughout the crust and lithosphere, modulating viscosity and organizing strain localization. The domain is extended at 4 cm/yr, and Gaussian random noise in the plastic damage parameter provides a mechanical seed for localization. Perturbation size of 50 km and 100 km are tested; at depth, these rectangular anomalies produce quasi-sinusoidal Moho temperature variations of 5°C and 20°C, respectively with wavelength $\lambda$ of 100 km and 200 km.

After 2 Myr, strain rate fields and isotherms (300°C, 600°C, 900°C) show that even modest lateral variations in heat production efficiently organize strain (Figure 3c–d). In Strong lithosphere, initial rift mode can switch from narrow to wide, while in Weak lithosphere, wide rifts form more rapidly and systematically within thermal highs. Larger-scale anomalies amplify strain localization, with fewer but more pronounced shear bands accommodating deformation.
These results demonstrate that steady-state thermal inheritance in the upper crust is sufficient to modulate rift initiation and early localization in 2D, providing a robust proof-of-concept for more realistic 3D experiments. In the following section, we extend these ideas to three dimensions, applying steady-state lateral variations of crustal heat production to a South China Sea-inspired geometry, where crustal heterogeneities, oblique extension, and propagating rifts interact to organize strain localization and fault patterns.

## 3 Tectono-thermal context of multi-propagator rifting in the South China Sea

<figure><img src="image_5.png" /><figcaption>**FIGURE 4** – Map illustrating the northern South China Sea margin, including E–W tear basins (yellow arrows) and the distribution of Triassic to Cretaceous plutons.</figcaption></figure>

### 3.1 Crustal context and pluton distribution

The South China Sea margins have been studied extensively due to the large number of seismic lines revealing variable crustal composition, including granitoids, metamorphic blocks, and clear structures of crustal boudinage (Pichot et al., 2014).

The continental margin of the South China Sea corresponds to the southern part of the continental
South China Block (Holloway, 1982a; Larsen et al., 2018; Taylor, 2025). Since the Triassic, this margin occupied the upper plate of a subduction zone. As a result, the crust has been injected by rounded plutons (Yan et al., 2010) from the Triassic until the Late Cretaceous (Huang et al., 2017; Cao et al., 2018). Subsequently, a block largely devoid of plutons, except for isolated Miocene intrusions and inferred to be the former Luconia continental block or a large oceanic plateau (Fyhn et al., 2010; Pubellier and Sautter, 2022; Huang et al., 2017), collided with the South China Block. From that moment, the future SCS margins remained in the upper plate of a subduction zone.

## 3.2 Rift geometry and propagators

The overall geometry of the SCS is rhomboidal, shaped by the construction of the continental shelf, mostly on the northern margin, and by large deltas along the NW coast of Borneo from the Early Miocene. Most of the oceanic basin exhibits a V-shaped morphology pointing southwest, reflecting the action of a propagator active from 24 Ma to 15.5 Ma (Sibuet et al., 2016; Chang et al., 2022). Earlier studies indicated that spreading initiated around 33 Ma in a N-S direction, with E-W magnetic anomalies predominantly in the Xisha Trough (Taylor and Hayes, 1980, 1983; Briais et al., 1993).

Along the northern margin, the continent-ocean transition displays a series of E-W tear basins forming a saw-tooth trace, with the Xisha Trough as the largest and northernmost feature (yellow arrows in Figure 4). These tear basins are analogous to small basins forming at the tip of a propagator during continental rifting in numerical simulations (Le Pourhiet et al., 2018), potentially representing relics of early westward extension visible in the continent-ocean transition (Pubellier et al., 2017). They are, however, more en-echelon than a simple propagator model, oblique to the extension direction, and terminating with the E-W Penyu basin at the Peninsular Malaysia margin. In this sense, they resemble en-echelon basins formed along localizing transform boundaries between offset spreading segments (Le Pourhiet et al., 2017; Ammann et al., 2018).

## 3.3 The role of plutons and lithospheric strength

A striking observation is that all N-S extensional notches vanish when they encounter a large band of Triassic to Cretaceous plutons (Savva et al., 2014). Some studies suggest these plutons are mechanically strong and can halt continental breakup propagation (Ding and Li, 2016; Qing et al., 2024), while others demonstrate that the lithosphere must remain weak to achieve distributed deformation along a rift segment (Gouiza and Naliboff, 2021). Both 3D modeling and simulations of continental rifting at propagator tips support this requirement for lithospheric weakness (Le Pourhiet et al., 2018; Jourdon et al., 2020).

Although structural inheritance may locally influence strain (Zhou et al., submitted; Li et al., 2024), in this work we focus on how lateral variations in crustal heat production associated with plutons might
<table><caption>TABLE 1 – Parameters used in 2D and 3D numerical models, Mantle is dry olivine (Brace and Kohlstedt, 1980), Strong crust is diorite (Tsenn and Carter, 1987) Weak crust is based on Quartz (Gleason and Tullis, 1995)</caption><thead><tr><th>Param.</th><th>Units</th><th>Eq.</th><th>Weak crust</th><th>Strong crust</th><th>Mantle</th></tr></thead><tbody><tr><td>A</td><td>MPa<sup>-n</sup>.s<sup>-1</sup></td><td>(16)</td><td>6.7 × 10<sup>-6</sup></td><td>0.063</td><td>2400</td></tr><tr><td>Q</td><td>kJ.mol<sup>-1</sup></td><td>(16)</td><td>156</td><td>275</td><td>540</td></tr><tr><td>n</td><td></td><td>(16)</td><td>2.4</td><td>3</td><td>3.5</td></tr><tr><td>V</td><td>m<sup>3</sup>.mol<sup>-1</sup></td><td>(16)</td><td>0</td><td>3.8 × 10<sup>-5</sup></td><td>8 × 10<sup>-6</sup></td></tr><tr><td>φ<sub>0</sub></td><td>°</td><td>(15)</td><td>30</td><td>30</td><td>30</td></tr><tr><td>φ<sub>1</sub></td><td>°</td><td>(15)</td><td>5</td><td>5</td><td>5</td></tr><tr><td>C<sub>0</sub></td><td>MPa</td><td>(15)</td><td>20</td><td>20</td><td>20</td></tr><tr><td>C<sub>1</sub></td><td>MPa</td><td>(15)</td><td>5</td><td>5</td><td>5</td></tr><tr><td>ρ<sub>0</sub></td><td>kg.m<sup>-3</sup></td><td>(11)</td><td>2850</td><td>2850</td><td>3300</td></tr><tr><td>α</td><td>K<sup>-1</sup></td><td>(11)</td><td>3 × 10<sup>-5</sup></td><td>3 × 10<sup>-5</sup></td><td>3 × 10<sup>-5</sup></td></tr><tr><td>β</td><td>Pa<sup>-1</sup></td><td>(11)</td><td>10<sup>-11</sup></td><td>10<sup>-11</sup></td><td>10<sup>-11</sup></td></tr><tr><td>H<sub>0</sub></td><td>μW.m<sup>-3</sup></td><td>(12)</td><td>0.365</td><td>0.365</td><td>0</td></tr><tr><td>κ</td><td>m<sup>2</sup>.s<sup>-1</sup></td><td>(12)</td><td>10<sup>-6</sup></td><td>10<sup>-6</sup></td><td>10<sup>-6</sup></td></tr><tr><td>ε<sub>0</sub> (2D)</td><td></td><td>(15)</td><td>0</td><td>0</td><td>0</td></tr><tr><td>ε<sub>1</sub> (2D)</td><td></td><td>(15)</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td>ε<sub>0</sub> (3D)</td><td></td><td>(15)</td><td>0.02</td><td>0.02</td><td>0.02</td></tr><tr><td>ε<sub>1</sub> (3D)</td><td></td><td>(15)</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td>ḣ (3D)</td><td>s<sup>-1</sup></td><td>(14)</td><td>10<sup>-15</sup></td><td>10<sup>-15</sup></td><td>10<sup>-15</sup></td></tr></tbody></table>

control continental rifting at both crustal and lithospheric scales.

## 3.4 3D South China Sea setup: thermal and mechanical configuration

The 3D simulations are designed to capture the effect of lateral variations in heat production on rift localization along the South China Sea margins. The model domain inherits the thermal concept from Section 2, extended to three dimensions. The geotherm is imposed as a thermal ribbon embedded within a lithosphere of variable thickness and crustal composition, representing the distribution of plutons and continental lithosphere (Figure 5a).

The simulations are based on the setup of Le Pourhiet et al. (2018). The Dirichlet boundary conditions impose lateral extension and a small amount of compression designed to slow down the propagation of rifts, while a two-dimensional Gaussian plastic seed is located on the opposite side of the modeling box (Figure 5b). A white noise in the initial plastic damage parameter is imposed everywhere and in all simulations.

As in 2D, we test a strong and a weak crust scenario, but we also extend to two other scenarios. The stratified scenario, based on Duretz et al. (2016), alternates horizontal layers of strong and weak crustal rheology and favors structural softening. The Post-Orogenic Lower Crust scenario (POLC) corresponds to an inverted rheological profile due to underthrusting of upper crustal material beneath an ophiolite or lower crustal units (Huet et al., 2011). Figure 5c presents the initial viscosity structure along a vertical cross section *yz* located at *x* = 400 km. Combining the thermal profile with four different crustal
213

214

215

216

217

organizations highlights how the thermal ribbon imposes lateral variations in lithospheric strength. It
shows clearly that the thermal ribbon affects the viscosity at different wavelengths: very short in the
upper crust, intermediate in the lower crust, and large in the mantle lithosphere. It is clear that differ-
ent crustal compositions will probably be more or less sensitive to the different wavelengths. The POLC
crust, in particular, seems to be quite insensitive to the short wavelength in the upper crust.

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

**FIGURE 5** – 3D South China Sea simulation setup. (a) Maps of the 300°C 600°C and 900°C isotherm depth with isocontours of upper crustal radiogenic heat production (thermal ribbon) that illustrate the 3D structure of the geotherm and its relationship with heat production. (b) Map view illustrating Dirichlet boundary conditions (arrows) and the location of propagator seeds (white outlines). (c) Vertically exaggerated cross section of initial viscosity for four lithospheric columns, including the yield strength envelope (taken outside the thermal ribbon). Dark red line indicates the location of the thermal ribbon.

218

219

220

This configuration allows us to explore how steady-state thermal inheritance interacts with rift prop-
agation, oblique extension, and crustal rheology to control strain localization and could lead to the de-
velopment of wide rift and en-echelon rift structures.
221

# 4 Results of 3D South China Sea Simulations

222

## 4.1 Effect of seeding

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

**FIGURE 6** – Maps of strain rate in grey scale overlapping the heat production after 4 Myr for weak crust. Columns corresponds to different distribution of heat production, lines to different distribution of initial plastic strain to seed the localisation. Black color for strain rate can be interpreted as faults/localised shear bands, gray-ish tones correspond to zone that deforms rather in a ductile manner. Zones that are in the transparent part deforms at rates that would not be detectable by GPS.

223

224

225

226

227

228

229

230

231

232

In 3D, just as in 2D, we need to seed the system to obtain strain localisation. Here, we explore three seeding possibilities. The first one, homogeneous isotropic seeding, consists of white noise in the plastic damage parameter. This does not impose any wavelength or anisotropy in the mechanical properties and is imposed in all simulations. The second is a two-dimensional Gaussian seed located on one side of the modelling box. It allows us to study the propagation of localisation across the lithosphere in the direction perpendicular to extension (Allken et al., 2012; Le Pourhiet et al., 2018; Jourdon et al., 2020). The third is a strip of random Gaussian-shaped anomalies of heat production that we refer to as the thermal ribbon. This ribbon is oriented at an angle to the extension direction in order to introduce obliquity. Figure 6 displays the distribution of heat production in 3D with shades of strain rate highlighting the deformation pattern after 4 Myr for these 3 scenarios and an additional scenario which combines the
thermal ribbon with the propagator. All simulations are shown for the weak lithosphere scenario. For the first two types of seeding, left column in Figure 6, heat production is set to $4H_0$, the highest value used in 2D, corresponding to wide continental rifting according to the 1D classification of Buck (1991).

The deformation patterns obtained under these seeding strategies are shown in Figure 6. Panel a illustrates the reference case with homogeneous plastic damage: deformation is distributed across numerous small grabens that develop orthogonal to the extension direction. In panel b, the addition of a single propagator focuses deformation at the propagator front, producing a few localised normal faults, while the surrounding areas remain relatively undeformed. In panel c, the presence of the thermal ribbon leads to a contrasting strain pattern: deformation outside the ribbon localises into narrow shear zones, whereas inside the ribbon strain remains distributed, with multiple small shear zones arranged both orthogonal and oblique to the extension direction. Such oblique and en-échelon structures are consistent with the type of fault linkage and obliquity described in the South China Sea by Savva et al. (2014). Finally, in panel d, the combined effect of a propagator and the thermal ribbon results in two competing strain corridors: one diffuse and conjugate, the other dominated by en-échelon faults contained within the ribbon, cross-cut by a single oblique fault. Outside the ribbon, deformation localises in a large, narrow graben.

Taken together, these four cases highlight the competition between kinematically imposed localisation and thermal inheritance. Kinematic forcing by a propagator tends to confine localisation to its tips, whereas in the absence of such forcing deformation remains distributed. In contrast, the presence of a thermal ribbon produces a bimodal strain pattern, with wide en-échelon distributed rifting inside the ribbon and narrow, strongly localised rifts outside. When both are combined, the propagator controls the initiation of deformation, but the ribbon provides an asymmetry that was not visible in the model with only a propagator.

## 4.2 Time evolution of a model without propagator

The narrow rift structures visible outside the ribbon at 4 Myr in Figure 6c and d are not random. To understand their origin, it is instructive to examine the temporal evolution of the weak crust model with a thermal ribbon but no propagator (Figure 7).

At 1 Myr, deformation inside the ribbon is distributed across numerous small faults, while outside the ribbon two narrow rifts form near the edges of the model domain. These are located where the ribbon comes closest to the free-slip sidewalls. Their orientation perpendicular to the boundary reflects the stress conditions imposed by the absence of shear stress on the free-slip walls: in effect, the ribbon corners act as secondary propagators into the stronger lithosphere outside the ribbon.
<figure><img src="image_8.png" /><figcaption>**FIGURE 7** – Temporal evolution of localisation in map view for the weak crust with a thermal ribbon and no propagator. Strain rate (grey scale) is overlain on the lithology with projected heat production.</figcaption></figure>

By 2 Myr, a new rift develops further outside the ribbon, along the $z_{max}$ boundary. Its position does not correspond to the geometric centre of the model but rather to a location where the 300°C and 600°C isotherms are strongly perturbed (see Figure 5a). This new rift begins to interact with the earlier rift on the opposite side of the ribbon. As a result, the first-generation rift zones lose activity, as is evident by 4 Myr in Figure 6c.

At 8 Myr, the two active lithospheric necks still do not connect as a transform fault, unlike in the models of Le Pourhiet et al. (2017) for similar distances of rift interaction. Instead, they delimit a rotating micro-block of continental lithosphere, comparable to the scenario described by Neuharth et al. (2021). In our models, the block is smaller than in that study, likely due to the overall weaker lithosphere used here.

Although this scenario does not directly apply to the South China Sea, it provides two important insights. First, in the absence of a propagator, deformation localises within the ribbon and then propagates into the surrounding stronger lithosphere, leading to a gradient from distributed to highly localised deformation. This transition is similar to the rheological inheritance described by Gouiza and Naliboff (2021). Second, free-slip boundaries act as strong attractors for localised structures within the ribbon, influencing where rifts nucleate and propagate.

## 4.3 Fault pattern and rheological stratification on Time evolution of models with propagators

At 4 Myr (Fig. 8), three of the four simulations with propagators—strong, stratified, and weak crust—have abandoned the initial lithospheric neck at the centre of the cross-section and evolved into large-offset
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**FIGURE 8** – Localisation as a function of time and rheology the crust in map view. Rendering of strain rate overlapping the type of rocks with a projection of the heat production.</figcaption></figure>

proto-transform systems, while the POLC model shows a distinct outcome. To understand these differences, we now turn to the early stages (1–2 Myr) and examine the 3D blocks of strain localisation at the propagator tip (Fig. 9). These snapshots illustrate how rheology controls the depth of necking, the
<figure><img src="image_11.png" /></figure>

**FIGURE 9** – Zooms on 3D blocks highlighting the mode of strain localisation as a function of crustal rheology. Each row corresponds to one rheological profile introduced in Figure 5c (Weak, Strong, Stratified, and POLC), and the two columns show snapshots at 1 Myr and 2 Myr. Strain rate is shown in grey scale with transparency: black tones correspond to faults or localised shear bands, whitish tones to more distributed ductile deformation, and transparent zones to deformation rates too low to be detected by GPS. The blocks also display the type of rocks using the colour code given at the bottom of the figure. Heat production is projected in map view at the surface but only where it exceeds $5H_0$, thus highlighting plutons.

287 geometry of faulting, and the influence of plutons on localisation.

288 In the strong crust case, the first structures visible at 1 Myr are already lithospheric in scale. Faulting roots in wide X-shaped mantle necks, producing normal faults that cut plutons from the outside rather than through their cores. By 4 Myr, the system develops into a transtensional corridor bounded by two oblique normal faults, with bookshelf-style strike-slip faults accommodating shearing inside the corridor.

290 The stratified crust also shows deep X-rooted mantle necks, but additional décollement levels within the crust generate antithetic normal faults and a more complex surface fault network. Strike-slip localisation simplifies to one dominant fault after 2 Myr, contrasting with the distributed bookshelf faulting in
296

the strong case.

297

In the weak crust case, plutons dominate the initial response. Grabens nucleate inside the high-heat production domains (in red) and propagate laterally as they accumulate displacement. By 2 Myr, these crustal faults link downwards to mantle necks. At 4 Myr, the system reorganises onto a single oblique strike-slip fault that cuts across the earlier pluton-centred grabens, marking the late recoupling of crust and mantle.

298

299

300

301

302

The POLC crust behaves differently. Localisation at 1–2 Myr is governed by the strong brittle upper crust, which produces evenly spaced grabens. The weak green lower crust flows in the footwall of the main detachment, forming metamorphic core complexes. At 4 Myr, plutons exert little influence, and deformation is partitioned into a rift and a strike-slip fault only loosely related to the thermal ribbon.

303

304

305

306

Taken together, Figs. 8 and 9 highlight a systematic control of crustal rheology on propagator dynamics. Strong and stratified crusts are mantle-driven, with localisation rooted at lithospheric depth, while weak and POLC crusts are crust-driven, with localisation imposed from above and only later recoupling to the mantle.

307

308

309

310

# 5 Discussion

311

## 5.1 Comparison with previous models

312

An important first remark is the distinct signature of the POLC model, which bypasses the influence of lateral heat production entirely, demonstrating that the effect of a thermal ribbon may be annihilated by a specific, though not very common, inverted rheological profile. It should be noted that the model used here included 25 km of strong crust, originally described as an ophiolite in Huet et al. (2011). Applying a thermal ribbon to an ophiolite-type rheology is therefore not geologically very relevant.

313

314

315

316

317

The simulations reported by Le Pourhiet et al. (2018) and Jourdon et al. (2020), with respective extension rates of 2 cm/yr and 1 cm/yr, showed that the application of a small shortening rate could delay the propagation of continental break-up for up to 10 Myr. By contrast, at an extension rate of 4 cm/yr, we find that it is not possible to distribute the stretching long enough to form a wide continental rift at the tip of a propagator, even in weak crust (Fig. 6), because of vertical advection of heat beneath the notch. The thermal ribbon allows strain to be distributed for 4 to 5 Myr in the best cases, thanks to the formation of en-échelon oceanic basins that share the extension. However, the ribbon systematically introduces a second lithospheric-scale propagator on the other side of the domain, reducing the problem to one of interacting propagators in a proto-transform passive margin (e.g. Allken et al., 2012; Le Pourhiet et al., 2017; Ammann et al., 2018; Neuharth et al., 2021).

318

319

320

321

322

323

324

325

326
In the absence of a seeded propagator, deformation reorganizes into a rotating micro-continent (Fig. 7) as the rift that develops within the ribbon propagates toward the free-slip boundary. This type of dynamics was also obtained in the simulations of Neuharth et al. (2021). Their models used two seeds propagating synchronously, and they found a “sweet spot” in lithospheric strength versus seed spacing that produced a similar outcome. Here, localized rift segments emerge from the hot thermal ribbon at approximately the same time, just as in Gouiza and Naliboff (2021), but because the weak and warm zone is oblique, the segments are offset, likely causing the rotation of the rift. Further investigation is needed to distinguish the relative influence of the applied shortening from that of the ribbon’s obliquity, but such questions are beyond the scope of this study, which focuses on the crustal and lithospheric tectonic signatures of the thermal ribbon.

## 5.2 Comparision with the SCS

All model scenarios, except POLC, develop secondary propagators in the southwest corner of the domain. This outcome is not completely incompatible with the existence and location of the Penyu Basin (Pe, Fig. 4), south of the Indochinese Peninsula just outside the thermal ribbon, but under our boundary conditions the models produce excessive lithospheric thinning compared to natural observations. As we pointed out, this limitation arises from (1) the free-slip boundary condition on the Vietnam side, which allows the southwest basin to act as a secondary spreading center, and (2) an excess of extension applied to the simulation. For this reason, the models are best compared with the early rifting stage of the SCS, before this propagator dominates the dynamics. A more realistic setup would involve transferring extension along a strike-slip boundary located, depending on the chosen tectonic hypothesis, either to the south along the proto-SCS fault system (Taylor and Hayes, 1983; Holloway, 1982b) or to the north along the Red River Fault (Briais et al., 1993).

In these models, the presence of a thermal ribbon in the weak crust scenario systematically organizes early rift development. Within the ribbon, deformation is distributed across numerous small en-échelon grabens, with successive basins developing along the ribbon's axis. Outside the ribbon, strain localizes into narrower rifts with higher subsidence rates. This dual behavior reproduces the observed north-to-south evolution of rifting in the SCS, where early large basins formed to the north and en-échelon basins appeared further south in the plutonic domain, particularly highlighted by the contour of the COT (Savva et al., 2014; Pubellier et al., 2017). These include the Xisha Trough, the Paracel, and the Phu Khanh basins, highlighted by arrows in Fig. 4. Based on reconstructions, Chang et al. (2022) proposed that they were en-échelon pull-apart basins that later coalesced during the second stage of SCS opening. Here, we interpret them instead as distributed proto-transform margins circumscribed to the thermal ribbon. The two interpretations are not inconsistent, and only simulations with alternative boundary conditions
360

could provide more conclusive insight.

361

As previously stated, the en-échelon basins form systematically within the ribbon, but they are longer lived—i.e., better developed—in the Weak Crust simulation (Fig. 8). Strong or stratified crust models are more sensitive to lithospheric-scale effects, and as a result, strain localizes very early into oblique, strike-slip-dominated structures. Although these oblique features align geometrically with the later oblique rifting (secondary rifting) phase of the SCS, they open unrealistically as pure wrench zones rather than transtensional basins. Nevertheless, the plutons of the SCS appear in the bathymetry as rounded entities with blurred seismic signals, considered as the pre-rift basement (Franke et al., 2014). This suggests that these bodies were only weakly dissected by rifting—a characteristic best reproduced by the Strong and Stratified experiments, though still present to a lesser extent in the Weak Crust case (Fig. 9, top). Because the wrenching is in part controlled by the artificial secondary propagator, which does not exist in nature, we conclude that our modeling does not provide a strong argument for a weak crust. Models with more appropriate boundary conditions that suppress this artifact are needed to properly test crustal rheology.

362

363

364

365

366

367

368

369

370

371

372

373

## 5.3 Thermal inheritance

374

The debate on inheritance in the geoscientific community arises from the difference between surface observations and subsurface physical parameters. Geologists observe upper crustal deformation imprinted in rocks, whereas geophysical data reveal deeper lithospheric properties that are inaccessible to direct observation. At former convergent plate boundaries, orogens record a history of subduction and magmatism before, during, and after collision. Belts of arc-related magmatism, syn-orogenic melting, and post-collisional magmatic flare-ups related to lithospheric delamination are systematic features of orogenies.

375

376

377

378

379

380

381

Subsequent extension may be influenced by the thermo-tectonic age of the lithosphere, especially during orogenic collapse. However, lateral variations in radiogenic heat production, represented in our experiments as thermal ribbons, cannot be neglected because they decay on much longer timescales than tectonic thermal anomalies. Moreover, even in a warm lithospheric environment like the South China Sea, our models demonstrate that these variations can significantly influence tectonics at the crustal scale.

382

383

384

385

386

387

The next steps are therefore: (1) to better quantify how these anomalies interact with mantle convection, rather than assuming a constant mantle heat flow as we did in the introduction, and (2) to ensure that more precise 3D geomodels of crustal heat-production anomalies become available to the community.

388

389

390
# 6 Conclusions

This study therefore expands the concept of thermal inheritance beyond transient geotherms to include a steady-state component: lateral variations in crustal radiogenic heat production. Whereas previous models of lithospheric evolution largely relied on thermal age, we show that lateral thermal heterogeneity, rooted in plutonic belts and arc-related magmatism, provides a long-lived control on rift initiation and early strain localization.

Our 2D proof-of-concept demonstrates that even modest lateral variations in crustal heat production efficiently organize strain: within thermal highs, rifts develop as distributed and wide systems, while outside they remain narrow and strongly localized. Extending this framework to 3D models inspired by the South China Sea highlights how steady-state thermal inheritance interacts with rheology and obliquity to produce contrasting fault patterns, including the systematic development of en-échelon basins within plutonic domains. These results support the interpretation that plutonic belts, through their elevated radiogenic heat, influence the geometry and distribution of rift structures at both crustal and lithospheric scales.

Our 3D numerical experiments provide new insights into how lateral variations in crustal radiogenic heat production—steady-state thermal inheritance—interact with crustal rheology and initial heterogeneities to control the initiation and localization of continental rifts, with the South China Sea as a natural test case.

Overall, our experiments underline the importance of considering lateral variations in crustal heat production when modeling continental rifting. Thermal inheritance acts as a persistent control on lithospheric rheology, guiding strain localization across crustal and lithospheric scales. When combined with realistic rheological assumptions and structural inheritance, it provides a predictive framework to explain the distribution of early rifts and the development of en-échelon basins in regions such as the SCS, where the extension established on top of a former orogenic mobile belt.

Together, these findings suggest that steady-state thermal inheritance is a fundamental ingredient of early continental rifting, not only in stable cratonic lithosphere but also in mobile belts where plutonic bodies thin the lithosphere and leave behind a crustal radiogenic signature. Recognizing this long-lived thermal control provides a framework for linking crustal composition to rift architecture and complements the more transient effects of thermal age and structural inheritance in governing how continental lithosphere deforms.
421

Acknowledgements

422

423

424

425

426

427

TF acknowledge salary support by the French National Research Agency (ANR), project ANR-19-CE05-0043 GERESFAULT. A.J. acknowledge salary support from the European Commission via the HORIZON-EUROHPC-JU-2021-COE-01 program, under Grant Agreement No 101093038 (ChEESE-2P: Center of Excellence for Exascale in Solid Earth). The simulation were perform on the OSU ECCE TERRA Rheolith Platform in 2D and using GENCI allocation A0160415068 allocated to LLP on Irene Skylake and Rome partition in 3D. We used OpenAI's ChatGPT (GPT-5, 2025) to assist with language editing.

428

Author contributions

429

430

431

432

List the contribution of each author following Contributor Role Taxonomy (CRediT): https://credit.niso.org/: LLP: Conceptualization, Investigation, Data curation, Funding acquisition, Visualization, Methodology, writing original and review; MP: Conceptualization, Investigation, writing original draft and review; AJ: Conceptualization, Software, Methodology, Writing - review and editing; TF: Conceptualization

433

Data availability

434

435

436

437

Results of the simulations are available at doi:10.5281/zenodo.17005985. The codes are freely available on bitbucket and ghub as listed in A.3 and in the metadata of the results repository. (As we do not want to publish before we are sure reviewer are satisfied with data availability, For review please use: zenodo link for review)

438

Competing interests

439

Authors declare no competing interests.

440

References

441

442

Allken, V., R. S. Huismans, and C. Thieulot (2012), Factors controlling the mode of rift interaction in brittle-ductile coupled systems: A 3D numerical study, *Geochemistry, Geophysics, Geosystems*, 13(5).

443

444

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

445

446

Autin, J., N. Bellahsen, S. Leroy, L. Husson, M.-O. Beslier, and E. d'Acremont (2013), The role of structural inheritance in oblique rifting: Insights from analogue models and application to the Gulf of Aden, *Tectonophysics*, 607, 51–64.

447

Balay, S., and et al. (2024), PETSc Users Manual, Tech. Rep. ANL-21/39, Argonne National Laboratory, revision 3.21.

448

Balay, S., W. D. Gropp, L. C. McInnes, and B. F. Smith (1997), Efficient Management of Parallelism in Object Oriented
449

450

Numerical Software Libraries, in *Modern Software Tools in Scientific Computing*, edited by E. Arge, A. M. Bruaset, and H. P. Langtangen, pp. 163–202, Birkhäuser Press.

451

452

453

Balay, S., S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, and H. Zhang (2015), PETSc Web page, http://www.mcs.anl.gov/petsc.

454

455

456

Balázs, A., E. Burov, L. Matenco, K. Vogt, T. Francois, and S. Cloetingh (2017), Symmetry during the syn-and post-rift evolution of extensional back-arc basins: The role of inherited orogenic structures, *Earth and Planetary Science Letters*, **462**, 86–98.

457

458

459

Balázs, A., L. Matenco, K. Vogt, S. Cloetingh, and T. Gerya (2018), Extensional polarity change in continental rifts: Inferences from 3-D numerical modeling and observations, *Journal of Geophysical Research: Solid Earth*, **123**(9), 8073–8094.

460

Bercovici, D., and Y. Ricard (2014), Plate tectonics, damage and inheritance, *Nature*, **508**(7497), 513–516.

461

462

Brace, W. F., and D. Kohlstedt (1980), Limits on lithospheric stress imposed by laboratory experiments, *Journal of Geophysical Research: Solid Earth*, **85**(B11), 6248–6252.

463

464

Brady, R. J., M. N. Ducea, S. B. Kidder, and J. B. Saleeby (2006), The distribution of radiogenic heat production as a function of depth in the Sierra Nevada Batholith, California, *Lithos*, **86**(3-4), 229–244.

465

466

467

Briais, A., P. Patriat, and P. Tapponnier (1993), Updated interpretation of magnetic anomalies and seafloor spreading stages in the South China Sea: Implications for the Tertiary tectonics of Southeast Asia, *Journal of Geophysical Research: Solid Earth*, **98**(B4), 6299–6328.

468

469

470

Brune, S., G. Corti, and G. Ranalli (2017), Controls of inherited lithospheric heterogeneity on rift linkage: Numerical and analog models of interaction between the Kenyan and Ethiopian rifts across the Turkana depression, *Tectonics*, **36**(9), 1767–1786.

471

472

Buck, W. R. (1991), Modes of continental lithospheric extension, *Journal of Geophysical Research: Solid Earth*, **96**(B12), 20,161–20,178.

473

474

Burov, E., L. Guillou-Frottier, E. d'Acremont, L. Le Pourhiet, and S. Cloetingh (2007), Plume head–lithosphere interactions near intra-continental plate boundaries, *Tectonophysics*, **434**(1-4), 15–38.

475

476

Burov, E. B., and M. Diament (1995), The effective elastic thickness (T e) of continental lithosphere: What does it really mean?, *Journal of Geophysical Research: Solid Earth*, **100**(B3), 3905–3927.

477

478

479

Cao, J., X. Yang, J. Du, Q. Wu, H. Kong, H. Li, Q. Wan, X. Xi, Y. Gong, and H. Zhao (2018), Formation and geodynamic implication of the Early Yanshanian granites associated with W-Sn mineralization in the Nanling range, South China: An overview, *International Geology Review*, **60**(11-14), 1744–1771.
480

481

482

Chang, S.-P., M. Pubellier, M. Delescluse, Y. Qiu, M. Nirrengarten, G. Mohn, N. Chamot-Rooke, and Y. Liang (2022), Crustal architecture and evolution of the southwestern South China Sea: Implications to continental breakup, *Marine and Petroleum Geology*, 136, 105,450.

483

484

Ding, W., and J. Li (2016), Propagated rifting in the Southwest Sub-basin, South China Sea: Insights from analogue modelling, *Journal of Geodynamics*, 100, 71–86.

485

486

Duretz, T., B. Petri, G. Mohn, S. Schmalholz, F. Schenker, and O. Müntener (2016), The importance of structural softening for the evolution and architecture of passive margins, *Scientific reports*, 6(1), 38,704.

487

488

Franke, D., D. Savva, M. Pubellier, S. Steuer, B. Mouly, J.-L. Auxietre, F. Meresse, and N. Chamot-Rooke (2014), The final rifting evolution in the South China Sea, *Marine and Petroleum Geology*, 58, 704–720.

489

490

491

Fyhn, M. B., S. A. Pedersen, L. O. Boldreel, L. H. Nielsen, P. F. Green, P. T. Dien, L. T. Huyen, and D. Frei (2010), Palaeocene–early Eocene inversion of the Phuquoc–Kampot Som Basin: SE Asian deformation associated with the suturing of Luconia, *Journal of the Geological Society*, 167(2), 281–295.

492

493

Gerbault, M., E. B. Burov, A. N. Poliakov, and M. Daignieres (1999), Do faults trigger folding in the lithosphere?, *Geophysical Research Letters*, 26(2), 271–274.

494

495

Gleason, G. C., and J. Tullis (1995), A flow law for dislocation creep of quartz aggregates determined with the molten salt cell, *Tectonophysics*, 247(1-4), 1–23.

496

497

Gouiza, M., and J. Naliboff (2021), Rheological inheritance controls the formation of segmented rifted margins in cratonic lithosphere, *Nature communications*, 12(1), 4653.

498

499

Heron, P. J., R. N. Pysklywec, and R. Stephenson (2016), Identifying mantle lithosphere inheritance in controlling intraplate orogenesis, *Journal of Geophysical Research: Solid Earth*, 121(9), 6966–6987.

500

501

502

Heron, P. J., A. L. Peace, K. J. W. McCaffrey, J. K. Welford, R. Wilson, J. van Hunen, and R. N. Pysklywec (2019), Segmentation of Rifts Through Structural Inheritance: Creation of the Davis Strait, *Tectonics*, 38(7), 2411–2430, doi: http://doi.org/https://doi.org/10.1029/2019TC005578.

503

504

Holloway, N. (1982a), North Palawan block, Philippines—Its relation to Asian mainland and role in evolution of South China Sea, *AAPG Bulletin*, 66(9), 1355–1383.

505

506

Holloway, N. (1982b), North Palawan block, Philippines—Its relation to Asian mainland and role in evolution of South China Sea, *AAPG Bulletin*, 66(9), 1355–1383.

507

508

509

Huang, C.-D., T.-Y. Lee, C.-H. Lo, S.-L. Chung, J.-C. Wu, C.-L. Tien, M.-W. Yeh, S.-C. Chen, Y.-L. Chan, and C.-Y. Hu (2017), Structural inversion in the northern South China Sea continental margin and its tectonic implications., *Terrestrial, Atmospheric & Oceanic Sciences*, 28(6).

510

Huet, B., L. Le Pourhiet, L. Labrousse, E. Burov, and L. Jolivet (2011), Post-orogenic extension and metamorphic
511

512

core complexes in a heterogeneous crust: the role of crustal layering inherited from collision. Application to the Cyclades ( Aegean domain), *Geophysical Journal International*, **184**(2), 611–625.

513

514

Jammes, S., and L. Lavier (2018), Effect of contrasting strength from inherited crustal fabrics on the development of rifting margins, *Geosphere*, **15**(2), 407–422.

515

516

Jaupart, C. (1983), Horizontal heat transfer due to radioactivity contrasts: causes and consequences of the linear heat flow relation, *Geophysical Journal International*, **75**(2), 411–435.

517

518

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

519

520

Jourdon, A., L. Le Pourhiet, F. Mouthereau, and D. May (2020), Modes of propagation of continental breakup and associated oblique rift structures, *Journal of Geophysical Research: Solid Earth*, **125**(9), e2020JB019,906.

521

522

Lachenbruch, A. H. (1968), Preliminary geothermal model of the Sierra Nevada, *Journal of Geophysical Research*, **73**(22), 6977–6989.

523

524

525

Larsen, H. C., G. Mohn, M. Nirrengarten, Z. Sun, J. Stock, Z. Jian, A. Klaus, C. Alvarez-Zarikian, J. Boaga, S. Bowden, et al. (2018), Rapid transition from continental breakup to igneous oceanic crust in the South China Sea, *Nature Geoscience*, **11**(10), 782–789.

526

527

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

528

529

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

530

531

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

532

533

534

Li, K., S. Brune, Z. Erdős, D. Neuharth, G. Mohn, and A. Glerum (2024), From orogeny to rifting: The role of inherited structures during the formation of the South China Sea, *Journal of Geophysical Research: Solid Earth*, **129**(12), e2024JB029,006.

535

536

537

Maierová, P., K. Schulmann, O. Lexa, S. Guillot, P. Štípská, V. Janoušek, and O. Čadek (2016), European Variscan orogenic evolution as an analogue of Tibetan-Himalayan orogen: Insights from petrology and numerical modeling, *Tectonics*, **35**(7), 1760–1780.

538

539

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

540

541

May, D. A., J. Brown, and L. L. Pourhiet (2014), pTatin3D: Efficient parallel finite-element modeling of lithospheric dynamics, in *Proceedings of the International Conference for High Performance Computing, Networking, Storage*
and Analysis (SC14), pp. 805-816, IEEE, doi: http://doi.org/10.1109/SC.2014.70.

May, D. A., J. Brown, and L. L. Pourhiet (2015), pTatin3D: a finite-element toolkit for efficient long-term lithospheric dynamics, *Computers & Mathematics with Applications*, 70, 413-430, doi: http://doi.org/10.1016/j.camwa.2015.06.031.

McKenzie, D. (1978), Some remarks on the development of sedimentary basins, *Earth and Planetary science letters*, 40(1), 25-32.

Neuharth, D., S. Brune, A. Glerum, C. Heine, and J. K. Welford (2021), Formation of continental microplates through rift linkage: Numerical modeling and its application to the Flemish Cap and Sao Paulo Plateau, *Geochemistry, Geophysics, Geosystems*, 22(4), e2020GC009,615.

Oravecz, É., A. Balázs, T. Gerya, D. A. May, and L. Fodor (2024), Competing effects of crustal shortening, thermal inheritance, and surface processes explain subsidence anomalies in inverted rift basins, *Geology*, 52(6), 447-452.

Oravecz, É., T. Gerya, and A. Balázs (2025), The location of compression-induced subduction initiation controlled by structural versus thermal inheritance, *Communications Earth & Environment*, 6(1), 652.

Perron, P., L. Le Pourhiet, M. Guiraud, E. Vennin, I. Moretti, É. Portier, and M. Konaté (2021), Control of inherited accreted lithospheric heterogeneity on the architecture and the low, long-term subsidence rate of intracratonic basins, *Bulletin de la Société Géologique de France*, 192(1).

Petri, B., T. Duretz, G. Mohn, S. M. Schmalholz, G. D. Karner, and O. Müntener (2019), Thinning mechanisms of heterogeneous continental lithosphere, *Earth and Planetary Science Letters*, 512, 147-162.

Pichot, T., M. Delescluse, N. Chamot-Rooke, M. Pubellier, Y. Qiu, F. Meresse, G. Sun, D. Savva, K. P. Wong, L. Wartmez, et al. (2014), Deep crustal structure of the conjugate margins of the SW South China Sea from wide-angle refraction seismic data, *Marine and Petroleum Geology*, 58, 627-643.

Pubellier, M., and C. Morley (2014), The basins of Sundaland (SE Asia): Evolution and boundary conditions, *Marine and Petroleum Geology*, 58, 555-578.

Pubellier, M., and B. Sautter (2022), Key structural elements around the East Vietnam sea (south China sea) and implications on reconstructions: Towards a clarification, *Vietnam Journal of Marine Science and Technology*, 22(3), 217-229.

Pubellier, M., D. Savva, M. Aurelio, and F. Sapin (2017), *Structural map of the South China Sea at 1:3 000 000.*, 1st ed., © CCGM-CGMW Commission for the Geological Map of the World, Paris, doi: http://doi.org/10.14682/2017STRUCTUSCS.

Qing, J., J. Liao, and S. Brune (2024), Rift propagation interacting with pre-existing microcontinental blocks, *Journal of Geophysical Research: Solid Earth*, 129(3), e2023JB028,109.
573

Ryan, W. B., S. M. Carbotte, J. O. Coplan, S. O'Hara, A. Melkonian, R. Arko, R. A. Weissel, V. Ferrini, A. Goodwillie, F. Nitsche, et al. (2009), Global multi-resolution topography synthesis, *Geochemistry, Geophysics, Geosystems*, 10(3).

574

575

576

Savva, D., M. Pubellier, D. Franke, N. Chamot-Rooke, F. Meresse, S. Steuer, and J.-L. Auxietre (2014), Different expressions of rifting on the South China Sea margins, *Marine and Petroleum Geology*, 58, 579-598.

577

578

Schmalholz, S. M., and N. S. Mancktelow (2016), Folding and necking across the scales: a review of theoretical and experimental results and their applications, *Solid Earth*, 7(5), 1417-1465.

579

580

Sibuet, J.-C., Y.-C. Yeh, and C.-S. Lee (2016), Geodynamics of the south China sea, *Tectonophysics*, 692, 98-119.

581

Smith, R. B. (1977), Formation of folds, boudinage, and mullions in non-Newtonian materials, *Geological Society of America Bulletin*, 88(2), 312-320.

582

583

Taylor, B. (2025), Reappraisal of the Continental Rifting and Seafloor Spreading That Formed the South China Sea, *Geosciences*, 15(4), 152.

584

585

Taylor, B., and D. E. Hayes (1980), *The Tectonic Evolution of the South China Basin*, pp. 89-104, American Geophysical Union (AGU), doi: http://doi.org/10.1029/GM023p0089.

586

587

Taylor, B., and D. E. Hayes (1983), Origin and history of the South China Sea basin, *Geophysical Monograph Series*, 27, 23-56.

588

589

Tsenn, M. C., and N. L. Carter (1987), Upper limits of power law creep of rocks, *Tectonophysics*, 136(1-2), 1-26.

590

Turcotte, D. L., and G. Schubert (2002), *Geodynamics*, Cambridge university press.

591

Vacherat, A., F. Mouthereau, R. Pik, M. Bernet, C. Gautheron, E. Masini, L. Le Pourhiet, B. Tibari, and A. Lahfid (2014), Thermal imprint of rift-related processes in orogens as recorded in the Pyrenees, *Earth and Planetary Science Letters*, 408, 296-306.

592

593

594

Yan, Q., X. Shi, J. Liu, K. Wang, and W. Bu (2010), Petrology and geochemistry of Mesozoic granitic rocks from the Nansha micro-block, the South China Sea: Constraints on the basement nature, *Journal of Asian Earth Sciences*, 37(2), 130-139.

595

596

597

Zhou, F., L. Le Pourhiet, M. Pubellier, and M. Delescluse (submitted), Impact of structural inheritance and mantle potential temperature on wide asymmetric rifts, *Earth and Planetary Science Letters*.

598
599

# A Numerical modelling method

600

## A.1 Governing equations

601

602

603

We simulate the long-term deformation of the lithosphere by solving the conservation of momentum, mass, and thermal energy for an incompressible, temperature-dependent, visco-plastic medium. The governing equations are:

$$ \nabla \cdot (\underline{\underline{\tau}}(\vec{u}, p, T) - p\underline{\underline{I}}) + \rho(p, T)\vec{g} = \vec{0}, \qquad (8) $$

$$ \nabla \cdot \vec{u} = 0, \qquad (9) $$

604

605

606

607

where $\vec{u}$ is velocity, $p$ is pressure, $T$ is temperature, $\rho$ is density, and $\vec{g}$ is gravitational acceleration. The stress tensor $\underline{\underline{\sigma}}$ has been decomposed into deviatoric and isotropic components, $\underline{\underline{\sigma}} = \underline{\underline{\tau}} - p\underline{\underline{I}}$, with $\underline{\underline{I}}$ the identity tensor. The deviatoric stress is related to the strain-rate tensor through a temperature- and pressure-dependent viscosity,

$$ \underline{\underline{\tau}} = 2\eta\dot{\underline{\underline{\epsilon}}}, \quad \dot{\underline{\underline{\epsilon}}} := \frac{1}{2}(\nabla\vec{u} + \nabla\vec{u}^T). \qquad (10) $$

608

Density variations are described using the Boussinesq approximation,

$$ \rho(p, T) = \rho_0 (1 - \alpha(T - T_0) + \beta(p - p_0)), \qquad (11) $$

609

with $\rho_0$ the reference density, $\alpha$ the thermal expansion coefficient, and $\beta$ the compressibility.

610

Thermal evolution is governed by the energy equation,

$$ \rho_0 C_p \left( \frac{\partial T}{\partial t} + \vec{u} \cdot \nabla T \right) - \nabla \cdot (k \nabla T) = H_0 + H_s, \qquad (12) $$

611

612

where $C_p$ is the heat capacity, $k$ the thermal conductivity, $H_0$ internal heating (e.g., radiogenic), and $H_s = \underline{\underline{\tau}} : \dot{\underline{\underline{\epsilon}}}$ the shear heating term.

613

The initial geotherm is obtained by solving the steady-state heat equation with internal heating,

$$ \nabla \cdot (k \nabla T) = H_0, $$

614

and adjusting the asthenospheric conductivity ($k = 70 \text{ W.m}^{-1}.\text{K}^{-1}$) to position the 1300°C isotherm at
615 105 km depth.

## A.2 Rheological model

617 Lithosphere and asthenosphere deformation is described by a visco-plastic rheology. Plastic (brittle) 618 failure is captured using a Drucker-Prager yield criterion,

619 where C is cohesion and φ the friction angle. The corresponding effective plastic viscosity is

620 To allow for fault weakening, we define a damage-like variable θ such that

621 and the friction angle evolves as

622 with φ₀ and φ₁ the initial and final friction angles, and ε₀, ε₁ the onset and completion of weakening. 623 Cohesion C follows an equivalent softening law.

624 Viscous deformation is represented by a dislocation creep law of Arrhenius type,

625 where A is the pre-exponential factor, n the stress exponent, Q the activation energy, V the activation 626 volume, and R the gas constant.

627 Finally, the effective viscosity is taken as the minimum of the plastic and viscous contributions,

$$ \eta = \min(\eta_p, \eta_v). $$

$$ \sigma_y(p) = C \cos \phi + p \sin \phi, $$

$$ \eta_p = \frac{\sigma_y}{\dot{\varepsilon}^{II}}, \quad \dot{\varepsilon}^{II} := \sqrt{\frac{1}{2} \dot{\varepsilon} : \dot{\varepsilon}}. \qquad (13) $$

$$ \dot{\theta} = \dot{\varepsilon}_p - \dot{h}, \qquad (14) $$

$$ \phi(\theta) = \max \left( \phi_0 - \frac{\theta - \varepsilon_0}{\varepsilon_1 - \varepsilon_0} (\phi_0 - \phi_1), \phi_0 \right), \qquad (15) $$

$$ \eta_v = A^{-\frac{1}{n}} \dot{\varepsilon}^{II} n^{-\frac{1}{n}} \exp\left(\frac{Q+pV}{nRT}\right), \qquad (16) $$
628

### A.3 Numerical discretization and solver

629

630

631

632

633

634

635

636

637

638

All simulations are performed using the finite-element code pTatin (May et al., 2015, 2014) in 2D (https://bit-bucket.org/ptatin/ptatin2d/src/bath2-sepide/) and 3D (https://github.com/laeticialp/ptatin-gene), which discretizes the governing equations on unstructured quadrilateral (2D) or hexahedral (3D) meshes. The weak form of the momentum-mass system is implemented via mixed finite elements (Q2/P1), with second-order interpolation for velocity and linear interpolation for pressure and temperature to ensure stability. Nonlinearities from visco-plastic rheology and thermal coupling are solved using quasi-Newton iterations, while the thermal and damage equations are integrated in time using an implicit (backward Euler) scheme. The energy equation is solved with finite elements in 2D and finite volumes in 3D. Linear algebra and solver operations rely on PETSc (Balay et al., 1997, 2015; Balay and et al., 2024), with scalable multigrid preconditioners ensuring efficient convergence in high-resolution 3D models.