Rift linkage and inheritance determine collisional mountain belt evolution

Domain: earth_science
Published: 2025-12-04
Editorial outcome: minor_revision

Abstract
Many mountain belts—such as the Pyrenees, European Alps, Greater Caucasus, or Atlas—form through inversion of pre-collisional extensional basins. These orogens exhibit three-dimensional complexity, with along-strike variations in topography, orientation, and deformation patterns. Yet, the relationship between these characteristics and the inherited extensional architecture remains enigmatic. Here, we use 3D geodynamic models coupled to a landscape evolution model to investigate how pre-collisional rift-linkage controls rift inversion and mountain belt evolution. Our results show that initial mountain belt structure reflects the inherited basin geometry, while later growth is governed by subduction polarity. This polarity depends on the magnitude of basin offset and the presence of pre-existing structural weaknesses. Comparison with natural examples suggests that along-strike variations observed in the Greater Caucasus, Atlas, and Pyrenees can be explained by the inversion of segmented and offset rift basins. Building on these insights, we propose a simple diagnostic framework that establishes a direct link between topography and deep lithospheric structures, showing how extensional inheritance influences mountain building on Earth.

Introduction
Mountain belts forming through continent-continent collision are one of the fundamental expressions of plate tectonics and are shaped by the interplay of tectonics, surface processes, and climate1–4. As part of the Wilson cycle5, many collisional mountain belts form at the locations of previous rift systems that are inverted during orogenesis and now constitute a significant part of the orogen. Examples in a continent-continent collision setting are the Pyrenees6,7, the Atlas mountains8,9, the Greater Caucasus10, the European Alps11, and the Bowen basin in Australia12. Inversion of extensional basins is also prominent in ocean-continent collision systems, as for instance observed along-strike in the Andes13–15. Looking in plan view at collisional mountain belts on Earth, we observe that mountain belts are not linear features, but typically contain steps with variations in local orogen orientation (1), and branching into several mountain belt domains separated by lowlands (2) (Fig. 1a, b). Looking in plan-view at active continental rifts like the East African rift system, we observe that also rifts are not linear features (Fig. 1c). Rather, extensional systems typically show offset basins and grabens that are often linked and connected through transfer zones, creating discrete steps (i), and branching of one rift arm into several domains (ii). A reconstruction of the pre-extensional template in the Pyrenees also exhibits linked and offset rift basins with discrete steps (i) and several rift branches (ii), which were reactivated and inverted during orogenesis (Fig. 1d). The observation of an ubiquitous extensional history pre-dating mountain building, and the observations that extensional systems and mountain belts show similar 3D complexities poses the question: How does the 3D architecture of pre-collisional extension influence mountain building?

Geodynamic modelling studies investigating the dynamics of inversion-tectonics during orogenesis found that extensional inheritance can have a strong influence on the structural evolution of mountain belts16–23. Most studies investigated the influence of extensional inheritance in a 2D cross section and found that weak inherited extensional structures are preferentially reactivated during collision and determine the structural style in a mountain belt16–21, 23. 3D studies of inversion tectonics primarily focused on the brittle part of the lithosphere and showed how inversion of extensional basins influences the evolution of fold-thrust belts, see review in ref. 22. These predominantly analogue modelling studies often exclude mantle dynamics and are inherently constrained by material properties, limiting their ability to reproduce inherited extensional structures. Hence, the importance of a 3D rift setting with stepping, offset, and linked rift basins on inversion tectonics remains unknown when including a consistent upper mantle-scale perspective.

Here, we use 3D upper-mantle-scale geodynamic models coupled to a surface processes model to investigate how rift linkage influences subsequent inversion during continental collision and mountain building. We choose a relatively simple model setup with two offset extensional basins that are subsequently inverted during continent-continent collision. The aim of this article is (a) to provide key characteristics of 3D inversion tectonics during mountain building, that can be used as a simple template to explain mountain belts on Earth, (b) to quantify the underlying physical principles that guide model evolution, and (c) to apply model results to well-studied mountain belts, specifically the Greater Caucasus, Atlas, and Pyrenees.

Results
Modelling results
We developed a coupling between the 3D thermo-mechanical geodynamic model pTatin3D24,25 and the surface processes model FastScape26,27, to investigate the effects of extensional inheritance on mountain building and the associated surface processes response. The tectonic model is laterally uniform and filled with materials corresponding to a typical layered continental lithosphere on Earth (see “Methods” section). The surface processes model solves for fluvial erosion, hillslope creep and sediment deposition. We do not vary the efficiency of fluvial erosion, but choose constant parameters for all models. The mountain belts produced by the chosen set of parameters can be classified as tectonics-dominated mountain belts with Beaumont number Bm=3.5, i.e., mountain belts that do not reach flux steady state, which is the dominant mountain belt type on Earth4. See Methods section for computation of Bm.

Model deformation is localised and initiated using two weak seeds. To investigate the influence of pre-collisional offset rift basins on mountain building, we vary the offset between the weak seeds and the shape of the seeds. Each model is extended for 12 Myr with a rate of 1 cm/yr in the east-west direction, followed by a reversal of the velocity boundary conditions and shortening with the same rate. This period of extension leads to full crustal break-up; the applied rate of extension/shortening represents typical plate velocities on Earth. Models run for 63 Myr but were stopped when the evolving orogen reached the eastern or western model side. This period and rate of shortening create orogens reaching intermediate size, slightly larger than, for instance, the European Alps. We first present the reference model without any offset between the seeds (M1), followed by a model with a very large offset of 400 km (M2), and a model with an intermediate offset of 200 km (M3). Finally, we present one model with an intermediate offset between the weak seeds (200 km) but with a preferential orientation of elongated weak seeds that mimic a pre-existing weak fabric (M4). Supplementary model output showing deformational regimes and deviatoric stresses is presented in Supplementary Figs. S1–S4.

Model 1 (Fig. 2) has simple and uniform weak seeds without offset between them. During extension, deformation first localises at the weak seeds and form two grabens that quickly propagate and link to form one cylindrical rift. At the end of the extension phase, the continuous rift system shows full crustal break-up with more offset in the weak seed domains than in the model centre (Fig. 2a). The symmetric rift faults are reactivated during the inversion stage and form a small juvenile mountain belt with pieces of (sub-) lithospheric mantle trapped in the inverted domain (Fig. 2b, cross sections). Shortening and topography are highest in the central domain, which experienced the least amount of extension. During the ensuing mountain building stage (Fig. 2c), the model grows by subduction of the western lithospheric mantle, and formation of thick-skinned crustal thrust sheets originating from the western, pro-side of the orogen. This mature inversion orogen is translated towards the pro-side of the orogen and keeps the cylindrical structure inherited from the rifting phase. The landscape is dominated by rivers flowing perpendicular to the orogen strike. Throughout the manuscript, we refer to ‘juvenile’ and ‘mature’ rifts or orogens. Juvenile rifts contain unconnected sub-basins, unlike mature rifts, where these are linked (Fig. 2a). Juvenile orogens lack clear subduction polarity and are only composed of inverted extensional structures (Fig. 2b), whereas mature orogens exhibit distinct lithospheric subduction and additional crustal shortening (Fig. 2c).

Model 2 (Fig. 3) is characterised by the same initial and boundary conditions as Model 1, but with a large offset between the weak seeds of 400 km. Initially, extension creates two isolated rift basins at the weak seeds. After around 6 Myr, the eastern rift basin becomes aborted, and the western basin propagates to the southern boundary of the model and turns into the main locus of extension (see Supplementary Movie 2). At the end of the extension phase, the western rift shows full crustal break-up near the northern boundary, and less extension near the southern boundary (Fig. 3a). During the inversion phase, the western rift is the main locus of deformation, but the eastern, aborted rift also gets inverted and forms a small and low mountain belt. Again, small pieces of lower crust and lithospheric mantle are trapped in the inverted domain (Fig. 3b). With ongoing convergence, the western, main mountain belt behaves similarly to Model 1: It grows by one-sided subduction and formation of thick-skinned thrust sheets originating primarily from the pro-side with translation of the orogen towards the west. The growing orogen is cylindrical with river flow perpendicular to the orogen boundary, and mimics the pre-existing extensional structure (Fig. 3c).

Model 3 (Fig. 4) is characterised by the same boundary and initial conditions as Model 1, with the difference of an intermediate offset of the weak seeds of 200 km. Extension first localises at the two weak seeds, forming two offset rift basins. With ongoing extension, the two rift basins propagate and link after 6 Myr of model evolution. Subsequently, they form an offset, linked and continuous rift basin with enhanced extension in the initial rift basins compared to the model centre (Fig. 4a). The ‘inner faults’ of the offset rift basin, i.e., the eastern fault in the northern basin and the western fault in the southern basin record more offset and therefore also slightly higher rift-flank topography (Fig. 4a). During the inversion stage, the extensional structures are reactivated and form a small and continuous mountain belt that has the same plan-view shape as the linked rift basin. The primary orientation of the drainage networks is perpendicular to the local strike of the orogen, so that the first-order river orientation is changing along-strike (Fig. 4b). Continued shortening of the orogen leads to the development of opposite-polarity subduction, and linkage of the two inner faults of the former rift basin. During growth, the orogen increases in size by thick-skinned thrusting, with thrust sheets forming on the respective pro-sides of the orogen. This mechanism leads to orogen growth and translation towards the west in the southern part of the model, and growth and translation towards the east in the northern part, reversing the inherited rift structure. Interaction of the two subducting slabs is characterised by several large pieces of lower crust and lithospheric mantle sheared off from the opposite-polarity subducting lithospheres, and creating a topographic low in the transfer zone (Fig. 4c).

Model 4 (Fig. 5) is characterised by the same boundary and initial conditions, and intermediate seed offset as Model 3, but contains two elongated and angled weak seeds that dip in the same direction. This model setup mimics a pre-existing weak fabric with homogeneous along-strike orientation. During extension the model evolves similarly to Model 3, with the difference of larger margin asymmetry related to the shape of the weak seeds (Fig. 5a). Model evolution is also similar to M3 during the inversion stage, forming a continuous mountain belt with the same plan-view shape as the pre-existing extensional basins and river-flow perpendicular to the local strike of the orogen (Fig. 5b). With ongoing convergence, the model develops same-polarity subduction. Same-polarity subduction is facilitated by the linkage of the eastern basin bounding faults of the northern and southern sub-basins that are related to the pre-existing weak zones. The orogen grows by thick-skinned thrusting of thrust sheets forming on the pro-side of the orogen, inducing orogen-translation towards the pro-side. To first order, the plan-view rift structure is maintained during orogenic growth, and the orogen continuously exhibits river flow perpendicular to the local orogen boundary orientation. As a second-order structural feature, we observe that a piece of lower crust and lithospheric mantle is sheared off in the southern part of the linkage zone (Fig. 5c, cross section II) and becomes trapped at shallow depths.

Supplementary Models 1–3 (SM1-3) show sensitivities to a very small seed offset (SM1), to an intermediate seed offset lower than in M3 (SM2), and to opposite polarity in dip of the weak seeds (SM3). SM1, with a very small offset between simple weak seeds of 20 km, shows similar behaviour to Model 1, with same-polarity subduction leading to first-order to a conservation of the offset basin structure during inversion and mountain growth (Supplementary Fig. S5). SM2, with an intermediate offset of 100 km between simple weak seeds, exhibits very similar behaviour to M3, with opposite polarity subduction as a result of inner fault-linkage (Supplementary Fig. S6). Additional models with similar offsets all develop into opposite polarity subduction. These results show that opposite polarity subduction in intermediate-offset simple seed models is a systematic and repeatable model result. SM3, with the same initial and boundary conditions as Model 4 but with an opposite polarity in dip of the weak seeds, shows a very similar evolution as Model 3: Opposite polarity subduction develops through linkage of the inner faults that are also related to the inherited weaknesses (Supplementary Fig. S7).

Model analysis
Our models show that rift basins tend to link when the lateral offset between them is small to moderate, but fail to link at large offsets, which is consistent with previous work28–31. To explain this behaviour, we present a simple force-minimisation analysis that quantifies how rift linkage depends on offset magnitude. See supplementary information for detailed information on the definition of the different forces involved, including characteristic values and simplifying assumptions. We assume that deformation follows the path of least resistance, minimising total mechanical work, and that rifting starts in the southern weak zone (segment P1-P2, Fig. 6a). From this initial segment, deformation can follow two possible paths: (a) it links to a second weak zone and utilises the segmented path D0+Dlink, requiring a force Flink, or (b) it bypasses the linkage and progresses directly to the model boundary along path D, requiring a force FNoLink. Linkage occurs only when Flink < FNoLink. Since the direct path of (b) is shorter, linkage and hence path (a) will only be utilised if the energetic gain G from deforming along the weak zone exceeds the penalty P of following a longer path. Applying these force considerations to Model 2, we can calculate Flink ≈ 6.4 × 1018 N > FNoLink ≈ 5.2 × 1018 N, which explains why the basins do not link. On the contrary, in Model 3 with small offset Flink ≈ 4.4 × 1018 N < FNoLink ≈ 5.2 × 1018 N, which explains the observed basin linkage. For our set of parameters the maximum offset for rift linkage to occur, i.e., Flink = FNoLink, is reached at 292 km offset.

To summarise, model behaviour can be successfully quantified by assuming that deformation follows the path of least resistance, such that the integral of the depth-integrated lithospheric strength along any given path is minimised. This implies that, in nature, weak zones in closest proximity are more likely to link during rifting. While this outcome may appear intuitive, the force-minimisation analysis offers a physically grounded explanation and predictive framework for understanding rift-basin linkage.

Using force minimisation principles, we can also quantify the evolution of subduction polarity in our models and its implications for orogenic systems on Earth. Consider the inversion of an intermediate-offset rift basin, as shown in Models 3 and 4 and illustrated in Fig. 6b. To first order, rifting produces a conjugate set of shear zones that dip toward the basin centre. During subsequent inversion, these shear zones are reactivated. As the lithosphere is restored to its normal thickness, the resulting mountain belt retains this conjugate shear zone geometry (represented in green and red in Fig. 6b). At this stage, deformation may proceed via three potential fault linkage pathways: 1. Inner fault linkage—connecting points PI1 and PI2 via segment LII. 2. Same-dip fault linkage—connecting PI1 and PI2 via LIO (or equivalently PO1 and PO2). 3. Outer fault linkage—connecting PO1 and PO2 via LOO. Minimising deformational work, it is expected that deformation utilises the shortest path, so that the inner faults link and LII is deformed (option 1). This linkage creates opposite polarity subduction related to the different dip-angles of the linked faults and can also be observed in Model 3. Linking faults with the same dip direction via LIO, option 2, requires greater deformational work because of the longer path, incurring a mechanical penalty P. For this pathway to be favoured, it must yield a mechanical advantage, or gain G, such that G > P. One potential source of gain is a difference in fault strength. For example, if the west-dipping shear zone (of length Lweak in Fig. 6b) is mechanically weaker than its east-dipping conjugate, it may be preferentially reactivated during shortening. This behaviour is observed in Model 4, wherein version proceeds via reactivation of the west-dipping shear zones. Here, the presence of pre-existing mechanical weaknesses provides sufficient gain to offset the additional work required by the longer deformation path. Finally, linkage of the two outer faults via LOO, option 3, also leads to opposite polarity subduction, but involves the longest deformation path and highest mechanical cost, making the least energetically favourable option.

These force minimisation principles also explain same-polarity subduction in cases of a straight rift system without significant offset, as in Models 1 and 2. In such configurations, LIO > LII, so the system preferentially reactivates one set of same-dipping faults (option 2), establishing a subduction interface with uniform polarity. Whether the left- or right-dipping faults are favoured depends on small perturbations during model evolution, such as localised erosion or variations in fault offset. The transition of the energetically most favourable linkage from same- to opposite-polarity subduction is governed by the magnitude of basin offset O. Once O > 0.5 dO, where dO is the surface distance between the conjugate and inverted rift-bounding shear zones (approximately 60 km in our models), opposite polarity subduction becomes the favourable option. Accordingly, any linked basin with an offset greater than 30 km will favour opposite polarity subduction. This transition explains why supplementary Model 2 (20 km offset and simple weak seeds) exhibits same-polarity subduction, whereas Model 3 (200 km offset and simple weak seeds) evolves into opposite-polarity subduction.

To summarise, the default behaviour for inversion of linked basins is same-polarity subduction. Once basin offset is larger than a few tens of km, i.e., LIO > LII, the default model behaviour is opposite polarity subduction and linkage of the inner faults. If, however, there are inherited weaknesses with a preferential orientation that provide a gain over the penalty of longer linkage G > P, same-polarity subduction may still prevail.

Comparison to the Greater Caucasus, Atlas, and Pyrenees
We next compare our modelling results with three mountain belts on Earth that are well known for their pre-collisional extensional history: the Greater Caucasus, the Atlas, and the Pyrenees.

The Greater Caucasus is an actively growing orogen forming between the Black and Caspian Seas as a result of the Arabia-Eurasia collision (Fig. 7a). The orogen is primarily a result of inversion of the Greater Caucasus basin, a continental extensional back-arc basin that resulted from subduction of the Neotethys below Eurasia and separated the Transcaucasian-Lesser Caucasus (South Caucasus) from the Scythian Platform (North Caucasus)32. Back-arc rifting occurred during the Jurassic, which formed a hyper-extended continental rift system, and was followed by a period of tectonic quiescence and basin filling10,33,34. The onset of inversion is debated but happened during late Eocene-Oligocene times (~35 Ma–30 Ma)33–35 and is likely related to the closure of the Neotethys and initiation of the Arabia-Eurasia collision36. The orogen is underlain by a N-dipping slab that is possibly tearing off in the W37,38. The slab length of 100–200 km indicates a similar amount of shortening of the mountain belt. Thickening and shortening are primarily located at the southern, pro-side, of the orogen35,39. The orogen also exhibits deformation on the retro-side and several thin-skinned foreland fold-thrust belts detaching in shallow weak horizons33,35,39. Following the drainage divide of the mountain belt, we see that the orogen has an overall NW-SE trend and is compartmentalised into longer NW-SE trending domains (B) that are offset by 30–40 km and connected by shorter W-E-trending and NNW-SSE-trending segments (T; Fig. 7a). This segmentation can also be observed in the Main Central Thrust (MCT), which is inferred to constitute the boundary between the Pro- and Retro-sides10. The orogen has a continuous and high topography with maximum elevations >5 km. The fluvial network follows, to first order, the changing trends in strike of the mountain belt, with drainage basins that are predominantly perpendicular to the local orientation of the drainage divide. The Greater Caucasus with its continuous high topography, high amounts of shortening, drainage network orientation, extensional history and changes in orogen-strike fits very well with a mature inversion orogen that exhibits the same subduction-polarity along strike, as modelled here in Model 4 (Fig. 7a). This characterisation implies that the orogen-segmentation reflects the configuration of the pre-collisional extensional basins, which we interpret as several extensional sub-basins (B) that are connected by transfer zones (T; Fig. 7a). Our modelling and force balance analysis indicates that same-polarity subduction orogens are favoured by weak-zones with a preferential orientation that provided deformational gain (Model 4, Figs. 5, 6b). A reconstruction of the extensional structure describes several N-dipping detachment levels10, which are likely reactivated during orogenesis and provide the weaknesses necessary to develop the observed same-polarity subduction orogen.

A well-known example of widespread basin inversion during continental orogenesis is the Atlas. The Atlas extends over 2000 km from the Moroccan coast through Algeria to Tunisia and the coastal Tell-mountains, and also branches out into the Rif-mountain through the Middle Atlas (Fig. 7b). Triassic and Jurassic rifting related to the opening of the Atlantic formed a linked intra-continental rift system that reactivated pre-existing Hercynian structures40. Shortening and inversion of the filled continental rift basins started in the Eocene with the main phase in the Oligocene and Miocene9,41. Shortening reactivated extensional structures, forming a mountain belt that is located on top of the pre-collisional basins9,41,42. Inversion was very mild with the highest amount of shortening of at most 35 km in the Central High Atlas, which potentially exhibits incipient subduction of lower crust and lithospheric mantle toward the north (Fig. 7b)8,9,42. Low amounts of shortening and exhumation explain why upper crustal rocks are in some locations at lower crustal levels inside the mountain belt compared to the forelands (Fig. 7b). Pre-collisional rifting was less pronounced in the Western High Atlas compared to the Central High Atlas, which explains exhumation of pre-collisional Palaeozoic rocks in the west43. Low topography in the eastern part of the Atlas indicates very low amounts of shortening and inversion in this domain. A relatively thin crustal root combined with a thinned mantle lithosphere suggests that some topography is supported by mantle processes9,44. Looking at the large-scale trends of the Atlas, we find that the western branch, the so-called High Atlas, trends SW-NE, with a change in strike at the Eastern High Atlas (EHA) towards W-E, which changes again after ~300 km to a SW-NE strike. The northern branch termed Middle Atlas (MA) has yet another orientation with a more northerly direction (Fig. 7b). The different large-scale domains exhibit heterogeneous orogen heights, with the High Atlas being a domain of continuous high topography. In the High Atlas, we observe a stepping drainage divide with river flow perpendicular to the local orientation of the divide, i.e., perpendicular to the local strike of the orogen. The mild inversion of extensional basins with low amounts of shortening and incipient subduction indicates that the Atlas is a juvenile inversion orogen where our models show that subduction polarity does not yet play a significant role, as shown in Models 3 and 4 after ~20 Myr of model evolution (Figs. 4b, 5b). The large-scale stepping of the mountain belt indicates that orogenesis reactivated a linked system of extensional basins. Different amounts of shortening and inversion, and potentially different amounts of mantle contribution, likely explain the variations in topography along the entire strike of the Atlas. Focussing only on the High Atlas, we interpret the stepping of the mountain belt strike with associated changes in fluvial drainage orientation as resembling the previous small-scale rift structure, with local SW-NE-trending basins with more E-W-trending linkage zones. The Atlas is a fitting example of a juvenile inversion orogen where subduction polarity has no primary influence, and highlights that basin linkage and associated changes in mountain belt structure during inversion occur on different spatial scales, from small-scale steps (10s of km) to large-scale branching (100s of km).

A classical example of basin-inversion during orogenesis are the Pyrenees. This double-vergent orogen formed as a result of the Iberia-Eurasia collision between Late-Cretaceous (~85 Ma) and Miocene times (~20 Ma)6,45,46. Collision inverted a laterally non-uniform extensional system which resulted in different characteristics along strike (Fig. 1d). In the E, between France and Spain, the Pyrenees sensu stricto (s.s.) formed by continent-continent collision that inverted the Pyrenean extensional system and also propagated into the un-thinned Iberian pro-side (Figs. 1d, 7c)6,7,47. Collision was accommodated by N-ward subduction of the Iberian lower crust and lithospheric mantle and crustal thickening, with shortening values of at most 90–165 km, decreasing from E to W6,7,46. The structural style is strongly influenced by pre- and syn-collisional weak salt layers that facilitated a pro-wedge thin-skinned foreland fold thrust-belt, and induced antiformal stacking of thick-skinned thrust sheets in the orogen core (Fig. 7c, section C-C’)6,7,18,48. Exhumed fragments of mantle rocks on the N-flank (retro-wedge) of the orogen, are found embedded in metamorphosed Mesozoic sedimentary cover between inverted crustal extensional blocks49,50. A large positive gravity anomaly in the W Pyrenees s.s., the Labourd anomaly, is also attributed to higher-density mantle rocks, and has been interpreted based on seismic tomography, as a continuous rise of the lithospheric mantle related to inversion of the hyperextended passive margin (Fig. 7c, (1))51. The decrease in shortening from E to W resulted also in a lower topographic relief and shallower exhumational levels in the W Pyrenees (Fig. 7c). Further to the W, the Pyrenees s.s. transition into the Basque Pyrenees, which grew by inversion of the Basque-Cantabrian Basin and northward subduction of Iberia (Fig. 1d). Shortening was lower than in the Pyrenees s.s. and involved primarily inversion of the extensional structures45,52. A very high magnetic and gravimetric anomaly in the vicinity of Bilbao can be explained by mantle and lower crustal rocks near the surface, underlying the sedimentary layers and continuous with the lithospheric mantle (Fig. 7c, (2), section A-A’)52. The reconstruction of the pre-collisional extensional basins indicates that the Pyrenean and Basque-Cantabrian basins were linked and offset by at least 100 km, while the Parentis basin further N did not link to the Pyrenean basin (Fig. 1d)45. Inversion of the extensional system initially produced limited shortening in the Parentis Basin, with deformation primarily localised in the Pyrenean-Basque Cantabrian basins45, consistent with the model evolution as shown in M2. Despite the intermediate basin offset between the Pyrenean and Basque-Cantabrian basins, which favours opposite polarity subduction, the orogen exhibits the same subduction polarity along strike. Same-polarity subduction with intermediate basin offset requires pre-existing weaknesses (M4, Fig. 6b), which are likely provided by a pervasive N-dipping Variscan pre-extensional fabric6,45. The basin offset can nowadays be seen in the kink in the transition between the Pyrenees s.s. and Basque Pyrenees (Fig. 7c, (3)). Higher amounts of shortening and pro-side deformation leading to southward translation of the Pyrenees s.s. likely counteracted the basin offset in the E and explains why the Basque Pyrenees and Pyrenees s.s. roughly form a continuous mountain belt with the basin offset only visible in the kink of the transition zone. The exhumed fragments of mantle rocks on the N-(retro-)side can also be observed in all of our inversion models. Model 4 furthermore shows that in the case of intermediate offset and same-polarity subduction conditions, a larger piece of lower crust and lithospheric mantle can be sheared off the retro-plate and become trapped at shallow depth on the retro-side. This mechanism potentially explains the shallow lower lithosphere in the Basque Pyrenees near Bilbao (Fig. 7c, (1)). We speculate that the location of the Labourd anomaly (Fig. 7c, (2)), situated at the W end of the pre-collisional Pyrenean basin, is similarly related to deformation of the basin linkage zone. The Pyrenees are thus a prime example of a mature same-polarity subduction orogen that also exhibits several secondary inversion features and highlights the complexity involved in the inversion of linked extensional basins.

Methods
Geodynamic model
We use the 3D ALE finite element model pTatin3D24,25 to solve for incompressible momentum (1), mass (2), and thermal energy conservation (3) in a 3D cartesian box:

∇·(2ηϵ̇) - ∇P = ρg, (1)

∇·v = 0, (2)

cp ρ (∂T/∂t + v·∇T) = ∇·(k∇T) + H + vz αg Tρ, (3)

where η is the material viscosity, ϵ̇ is the strain rate tensor, P the pressure, ρ the material density, g the gravitational acceleration vector, v=(vx, vy, vz) is velocity, cp is the specific heat capacity, T is the temperature, t is time, k is the heat conductivity, H the radiogenic heat production rate. Note that the subscript z indicates the vertical direction. Density varies as a function of temperature

ρ(T) = ρ0(1 - α(T - T0)), (4)

with the reference density ρ0 corresponding to the density at temperature T0. We do not include changes in density as a function of pressure.

Model materials deform either by frictional-plastic deformation or non-linear temperature-dependent viscous creep. Frictional-plasticity is modelled using a Drucker-Prager yield criterion,

τY = Ccos(ϕ) + Psin(ϕ), (5)

where τY is the deviatoric yield stress, ϕ is the angle of internal friction, and C is material cohesion. Material fails by frictional-plastic deformation when

σ0II > τY, (6)

where σ0II is the square root of the second invariant of the deviatoric stress tensor σ0II = 2ηϵ̇. We track material damage (ϵplastic) via the second invariant of the accumulated plastic strain. Strain weakening is accounted for by a linear decrease in the angle of friction from 15° to 2° over the pre-defined strain interval 0.5 < ϵplastic < 0.7. Strain weakening, as parameterised here, accounts for the combined effects of permanent material damage occurring during faulting. We do not apply regularisation methods to eliminate the mesh-dependence of the strain weakened material. Tests with different strain weakening thresholds did not show a first-order effect on model results.

Viscous deformation is modelled using a T- and P- dependent power-law rheology

σ0visc = fA-1/n (ϵ̇II)1/n exp(Q + PV / nRT), (7)

where σ0visc is the square root of the second invariant of the deviatoric stress, f is the scaling factor used to model different geological scenarios, A is the pre-exponential factor, n is the power-law exponent, Q the activation energy, V the activation volume, and R the universal gas constant.

Surface processes model
We use the landscape evolution model FastScape to model landscape evolution above sea level26,27. FastScape solves for stream power law erosion (K-term in eq. (8)), hillslope creep (Kd-term), and in-stream sediment deposition (G-term):

∂h/∂t = U - Kf AmSn + Kd ∇2h + G ∂h/∂t dA, (8)

where h is elevation, t is time, U is uplift, Kf is the stream power law coefficient termed fluvial erodibility, A is drainage area, S is slope, m, n are the stream power law exponents, Kd is the hillslope diffusion coefficient, and G is the depositional coefficient. Rivers are always connected to one of the four model boundaries through bridging of local minima, ensuring hydraulic connectivity. We use the lake-filling algorithm to fill local minima. The sea level is set at -500 m with respect to the model sides. Sediments transported to the shoreline are deposited with a mass-conserving aggradation scheme, filling the marine extensional and foreland basins with sediments from the bottom. m, n and G are relatively well known, with G ≈ 156, m/n varying from 0.3 to 0.557, and n varying from 1–358. Here, we use the typical values n=1, m=0.4, G=14. The fluvial erodibility spans a wider range as it incorporates variations related to e.g., rock strength, climate, or vegetation, with values ranging from at least 1×10−6 m0.2/yr to 1×10−4 m0.2/yr, given m/n=0.4, n=157. The efficiency of erosional versus tectonic flux can be expressed by the non-dimensional Beaumont number Bm4, which, given m, n as used here, simplifies to:

Bm = (ρc g zdec v) / (0.5 0.2 Kf W0.8 min Fint), (9)

where zdec = 25 km is the depth to the main crustal decoupling horizon, ρc = 2750 kg/m3 is the average crustal density, vc = 1 cm/yr is the convergence rate, Wmin = 80 km is a typical length-scale of minimum orogen width, and Fint ≈ 2.7 × 1012 N/m is approximately equivalent to the maximum integrated crustal strength of undeformed foreland crust above zdec during mountain building as measured in the model (see supplementary material for deviatoric stress plots). We chose Kf = 0.5 × 10−5 m0.2/yr for all our models. This set of parameters creates strength-limited mountain belts with Beaumont number Bm=3.5, a value typical for mountain belts on Earth that are not exposed to wet and seasonal climatic conditions. We use a constant hillslope diffusion coefficient of Kd = 1 × 10−2 m2/yr59.

The landscape evolution model is two-way coupled to the geodynamic model. FastScape has the same amount of nodes as the mechanical Q2-nodes of pTatin3D’s finite element mesh. Hence, FastScape truly sits on top of the geodynamic model. After each mechanical timestep, the landscape evolution model is advected according to the computed velocity field at the surface of the geodynamic model. Changes in topography are then computed using eq. (8) and the associated lake-filling and aggradational algorithms for deposition. The new topography is given back to the geodynamic model, and timestepping continues. We note that we limit the FastScape timestep to 1000 yr so that several FastScape timesteps are normally executed per mechanical timestep, which is dynamically computed based on a Courant-Friedrichs-Lewy condition and typically in the order of 9000 yr. Deposited sediments have the same properties as the upper crust, and we track the stratigraphy of deposited sediments in the geodynamic model at 5 Myr intervals.

Modelling setup
The geodynamic model is laterally uniform and vertically stratified consisting of upper and middle crust to 25 km depth, 10 km of lower crust up to the Moho at 35 km, lithospheric mantle down to a depth of 120 km, and sub-lithospheric upper mantle down to the lower model boundary (Fig. 9). The model has dimensions of 1200 km x 600 km x 600 km and is subdivided into 256x128x128 finite elements, resulting in a horizontal resolution of 4.7 km per Q2-element. Vertically, the mesh is non-uniform, with a vertical resolution of 2.0 km in the upper 150 km of lithosphere and 8.6 km in the lower 450 km of the model. We base our model materials on a few well-established laboratory-derived flow laws and use the scaling factor f to account for different geological scenarios and uncertainties when extrapolating laboratory-derived flow laws to nature. Specifially, we base model materials on the ‘wet’ quartz flow law from ref. 60 for sediments and upper and middle crust, the ‘dry’ Maryland diabase flow law from ref. 61 for lower crust, and the ‘wet’ olivine flow law from ref. 62 for mantle material. The lithospheric mantle is scaled by a factor f=10 to account for compositionally stronger, e.g., melt-depleted, conditions. The crustal materials based on the ‘wet’ quartz flow law are scaled by a factor f=0.5. The latter leads to mountain belts with at most 8 km in height, similar to the maximum topography as observed on Earth. Material densities are chosen to be simple and represent typical values of the Phanerozoic lithosphere (Fig. 9 and Supplementary Table 1). The lithospheric mantle is depleted by 15 kg m−3 with respect to the underlying sub-lithospheric mantle, which is typical for Phanerozoic lithosphere63.

The initial temperature distribution in the model domain is laterally uniform and at conductive steady state, with 550 °C at the base of the crust and 1330 °C at the base of the lithosphere, resulting in a surface heat flow of ~53 mW m−2 and a mantle mantle heat flow ~21 mW m−2, which is typical for Phanerozoic lithosphere64. k is linearly increased from k = 2.25 W m−1 K−1 to k = 52.0 W m−1 K−1 between 1330 °C and 1340 °C, in order to maintain the heat flow at the base of the lithosphere, prevent the model from cooling, and resembling mantle convection at high Nusselt number65. This approach leads to an adiabatic gradient of 0.4 °C/km in the sub-lithospheric upper mantle. All the side boundaries are insulated, the bottom boundaries are described by Dirichlet boundary conditions with 0 °C at the top and 1522 °C at the bottom boundary.

To model continental rifting and orogenesis, we apply velocity boundary conditions on the E and W faces of the model; material inflow or outflow on the sides is compensated by outflow/inflow at the bottom boundary. We apply extensional boundary conditions of 0.5 cm/yr at each side, for 12 Myr and then apply shortening with 0.5 cm/yr at each side. The transition from extension to shortening occurs at a linear rate of velocity change between 11.5 Myr and 12.5 Myr of model evolution. The side boundaries have vertical roller conditions, the bottom boundary is characterised by free slip, and the top boundary is a true free surface subject to erosion and deposition by the surface processes model.

Many mountain belts on Earth exhibit to first order a decoupling between crust and mantle at ~20–25 km depth, one-sided subduction of the lithospheric mantle, and mountain building related to crustal thickening6,38,66–68. The chosen rheological and thermal setup (see Fig. 9) allows for a decoupled evolution of crust and mantle and represents this first-order observation. It’s very typical that mountain belts not only exhibit thick-skinned deformation of crustal-scale thrust sheets, but also thin-skinned thrusting related to weak and shallow décollement horizons in the upper crust, as observed amongst others in the Pyrenees, Alps, or Caucasus. Thin-skinned thrusting changes the structural style, e.g., ref. 18, but is deliberately not part of the model, as the numerical resolution is not high enough to properly resolve this deformational style.

Model Limitations and extended parameter space
Many pre-collisional extensional basins exhibit a long phase of several 10s to 100 Myr of post-rift cooling before basin inversion. We do not include a phase of post-rift quiescence, but seamlessly transition from extension to shortening in our models. Buiter et al. 16 showed with 2D geodynamic models that varying amounts of post-rift cooling have a very limited effect on the structures developing during inversion, because the weak fault zones formed during extension retain their weakness during post-rift cooling and are therefore reactivated during inversion. We also tested the influence of post-rift cooling for one of our models and found that it has no first-order influence. Furthermore, the force minimisation principles and analysis presented above also hold for geological situations of prolonged post-rift subsidence. We therefore expect that the first-order features observed in our modelling are also applicable to mountain belts exhibiting a phase of post-rift cooling. Further studies should nonetheless investigate the effect of post-rift cooling in detail, including an investigation of transient weakening effects and fault healing mechanisms.

Many mountain belts, like for instance the Alps or Himalaya-Tibet, exhibit a phase of oceanic subduction before continent-continent collision. We do not apply extensional boundary conditions in our models beyond crustal break-up and therefore do not form a large oceanic domain. Future research should focus on this aspect and also include melting processes and the formation of new oceanic lithosphere.

We only apply frictional-plastic softening as a weakening mechanism in our models. Although the reactivation of weak frictional-plastic shear zones is characteristic of inversion orogens - such as those examined in this study - we note that additional permanent or transient weakening processes may play a role during extension and orogenesis, like grain-size variations, shear heating, or metamorphic phase changes. Future modelling efforts should explore these additional mechanisms. The model results and force balance analysis presented here can serve as a comparative baseline for such investigations.

We only use a uniform surface process efficiency along-strike in the orogen. A strong lateral gradient in the efficiency of erosion, e.g., related to strong climatic gradients or variability in rock erodibility, would induce a higher widening rate in areas with lower erosional efficiency, and vice versa4. In this case, the along-strike differences in surface processes would potentially overprint the mountain belt structure inherited from a previous extensional event, and the previous basin structure is likely only visible in the transfer zone, i.e., in a laterally confined area. A similar effect can be obtained by varying the amount of convergence along strike, as discussed for the Pyrenees. Furthermore, along-strike differences in the strength of the colliding plates would have a similar effect. Higher plate strength would lead to a higher and narrower mountain belt, and vice versa4.

We investigated the simplified case of inverting two offset rift basins. In natural systems, however, multiple offset basins are likely to occur. Based on our model results and force-balance analysis, we predict that inverting such systems would be behave as if composed of several basin pairs similar to those examined here, potentially involving multiple subduction polarity flips in the absence of dominant inherited weaknesses guiding deformation. Future research should explore this topic.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The paper presents new, state-of-the-art numerical scenarios that combine a 3D thermomechanical model with a landscape evolution model to simulate the influence of rifting geometry on the evolution of a collisional mountain belt. The authors explore scenarios with different offsets between two rift branches, demonstrating that the orogenic structure is controlled by the rifting geometry. Additionally, they present an analytical approach that offers a convincing physical explanation for the numerical results, highlighting the conditions under which rift branches with different offsets may connect.\nI recommend publication of the article after the following points are addressed:\nIn scenario M2, which involves a large offset of the weak seed (400 km), the simulation evolves into one-sided subduction with an essentially two-dimensional (cylindrical) geometry. As described by the authors, this occurs because the eastern rift branch is aborted, and the western basin propagates southward, resulting in a linear rift geometry. However, this result likely depends on the size and rheology of the weak seeds: if both branches contain large, sufficiently weak seeds, rifting might develop on both the western and eastern sides, potentially leading to a successful rifting scenario with a large offset. It is probable that the small size of the weak seed contributed to the failure of one branch, making the resulting scenario redundant with scenario M1.\nIn scenario M3, which uses an intermediate offset of 200 km between weak seeds, the simulation produces an interesting subduction pattern with opposite polarity, resulting in a complex, three-dimensional orogenic system (Figure 4c). My question concerns the robustness of this result: is opposite-polarity subduction a systematic outcome for intermediate-offset (~200 km) scenarios with simple weak seeds, or is it a fortuitous result? Given that the weak seeds in these scenarios are symmetric in East–West cross-sections, and that 2D thermomechanical models with similar seeds do not allow for a priori prediction of subduction polarity (since small perturbations in the thermal and compositional structure of the lithosphere can change the polarity) the observed pattern may not be systematic. Additionally, the numerical setup explores only a single offset per scenario. What would be expected if multiple branches with intermediate offsets (~200 km) were present in the domain? As the authors convincingly argue, polarity is sensitive to inherited geometries of the weak zones. Therefore, I suggest that the offset plays only a secondary (or negligible) role in controlling subduction polarity.\nMinor points:\nThe authors use cardinal points to indicate features in the numerical scenarios. This is a helpful approach, especially for 3D models with complex geometries, as it facilitates understanding of the numerical results. However, since the scenarios are generic and not intended to represent any specific orogenic system on Earth, it is important to define the direction of the initial extension. This would improve the reader's comprehension of the results. I suggest adding a clarification in lines 79–81, such as: \"We assumed the direction of extension to be East–West.\" Although this is mentioned in the Methods section and indicated in Figures 2–5, including it in the main text would enhance clarity.\nContinuing on this topic, the description of the Model 3 results is unclear. The authors state in lines 119–120 that \"the right fault in the northern basin and the left fault in the southern basin record more offset and therefore also slightly higher rift-flank topography.\" However, the combination of cardinal points with terms such as \"right\" and \"left\" is confusing. Please consider replacing these with \"eastern\" and \"western\" for clarity.\nModel names appear in both lowercase (e.g., line 109) and uppercase (e.g., lines 221–222). Please standardize their formatting throughout the manuscript.\nReference numbers in lines 461–469 should be placed within square brackets [ ] to avoid confusion with exponents.", "reviewer_2": "Review of Wolfe, S.G., Huismans, R.S., Munoz, J.A., May, D.A.\nRift Linkage and inheritance determine collisional mountain belt evlution\nFor Nature communications\nThis paper presents the results of a series of sophisticated 3D thermomechanical simulations of lithospheric extension followed by shortening. The objective is to consider how offset or linked rift segments affect subsequent orogen growth. The thermomechanical software is, in addition, coupled with a 3D surface processes modelling package that simulates erosion from growing relief and deposition of sediments in developing basins. 3D models simulate orthogonal, segmented rifting that is allowed to develop to the point of continental breakup followed by convergence to juvenile (or partial) inversion, mature (or complete) inversion and beyond. Extensional and compressional stresses are parallel. The simulations consider the effect of offset distance between rift segments on orogen evolution. The authors conclude that this factor influences initial growth of the mountain belt being strongly impacted by the offset distance of rift segments. Later orogen growth is governed by subduction polarity, which itself is mainly dictated by the asymmetry of the main conjugate rift faults. The model input parameters and uncertainties are described at the end of the paper. Simulation results are compared to three mountain chains, the Causases, Atlas and Pyrenees, which are in turn interpreted in terms of rift geometry inheritance.\nThis is an excellent paper from one of the leading thermomechanical modelling groups in geosciences. The methodology is sound and well established with the innovative combination of two softwares to provide a more complete 3D earth system simulation. It is well written and illustrated with impressive videos of model development. The authors state clearly the uncertainties and limitations of their models. The conclusions and comparisons with natural examples are valid, even if some queries need addressing. The results will provoke much discussion and interest amongst a wide readership interested in Wilson cycle lithospheric processes. It is also a nice demonstration of cutting edge multidimensional thermomechanical modelling combined with surface processes modelling. I feel that the surface processes aspects of the modelling do not get the attention that perhaps they merit. However, this may perhaps be the subject of a later paper. I recommend this paper for publication once the points below and other minor comments on the pdf are addressed.\nMary Ford,\nUniversité de Lorraine, Nancy.\nLines 53_55: Some analogue modelling studies do actually address these issues. What is missing from the majority of analogue models, it seems to me, is a full crustal/lithospheric treatment of these issues.\n79: Can you explain why this rate and period were chosen? This appears to thin the lithosphere to the point of mantle exhumation and breakup?\n81: State how long the shortening phases were and what the rate of shortening was. The duration of shortening appears to be variable (roughly 44 Myrs to 63 M Yrs) but always a lot longer than 12 Ma. Give the reasoning used for the duration and rate of the shortening ? Attaining a set amount of shortening for example ?\n95: southern?\n96: southern? Why use left and right? Why not stick to the cardinal directions of the compass as on Figures 2-5 of models and elsewhere in text?\n96: What do you mean by \"Mature\"? Fully or completely inverted? Maybe explain what you mean here and/or give reference. Same point for Juvenile inversion.\n119 – western fault\n125-6 : Oppositely dipping subduction - how do these link ? Can you give some detail here of how this works in the models?\n129 : not clear here. Reverting is not the appropriate word as it means to return to a former state. Inverting or switching would perhaps be clearer.\n178 : But you say you have already assumed this in line 166.\n273-275 : Incomplete sentences... punctuation problem?\n288 : EHA is not defined. Also this sentence is not clear and needs a rewrite.\n295 : Compare to what stage of these models ? You could define this in the videos\n304 : Its not clear how do the Atlas mountain belt highlights this point that basin linkage and associated changes in mountain belt structure occur at different scales during inversion.\n310 : The Mauleon basin is a small (but important) depocentre within the western Pyrenean system (SW of Pau). It is misleading and inappropriate here to refer to the whole Pyrenean rift (ss) by this name.\n310 ‘rifted margin’. Is the Pyrenean depocentre a rift or a rifted margin? So far you have been talking about rifts. If the latter then which margin are your referring to?\n310 : Suggest propagated would be clearer than ‘progressed’.\n317 : Suggest rewording : « Exhumed fragments of mantle rocks on the N-flank (retrowedge) of the orogen, are found embedded in metamorphosed Mesozoic sedimentary cover between inverted crustal extensional blocks ». I suggest citing not just Teixell et al. (2018) but also, for example, Angrand et al. (2020, J. Geol Soc London) and/or Asti et al. (2022 ; Earth Science Reviews)\n320 : In model of Wang et al. the rise of the mantle lithosphere occurred during hyper-extension. This configuration was essentially preserved due to very low convergence in this area of the westernmost north Pyrenees.\n393 : I dont think this explains the presence of large high density bodies. The offset between the Basque Pyrenees and the main Pyrenees ss can certainly be explained by inversion of a linked offset rift basin. But the presence of high density bodies is surely related to exhumation of deep lithosphere during hyper-extension and their subsequent incorporation into the orogen. In other words their presence is related primarily to the initial style and amount of extension in the rift system...\n409 : and the amount of shortening?\n411-414 : Is there any natural mountain belt that demonstrate this behaviour?\n481 : This section and Figure 9 would be better in my opinion early in the paper to allow the reader to understand the model results. For example, the colours codes for results figures are given on figure 9.\nFigures\nFigure 1. I have located the actual Mauleon basin on figure 1d. As this is a well known and defined local basin, its name is not appropriate for the whole Pyrenean rift system.\nFigures 2 to 5. It would be helpful to put the compass points on at least the first model image Need a key or labels for colours on the cross sections.\nCan you separate the three cross sections as complex details are not clear in the areas where they overlap ? Same comment for all cross sections.\nGiven that the 3D map views are labelled a, b, c, it would probably be better not to use A, B, C for the cross sections ?", "reviewer_3": "The manuscript titled “Rift linkage and inheritance determine collisional mountain belt evolution“ by Wolf and co-authors presents 3D numerical models of lithosphere extension directly followed by convergence. The main aims of the modelling are to better understand linkage of rift segments during extension and the impact of the rift geometry on the subsequent convergence, particularly the subduction polarity. The applied numerical model consists of a lithosphere deformation model coupled to a landscape evolution model. The study includes a very interesting application of force analysis to explain linkage of rift segments and the subduction polarity during subsequent convergence and rift inversion. Furthermore, the study discusses the application of the model results to the Greater Caucasus, Atlas and Pyrenees. There are still many open questions concerning the impact of 3D features on rifting and rift inversion and the manuscript, hence, presents a very interesting study on the application of deterministic mathematical models to better understand these geodynamic processes.\nThe presented results demonstrate the importance of the pre-orogenic 3D rift geometry on the subsequent orogeny related to the inversion of the rift basins. The authors apply a force analysis to explain the 3D model results which is very useful for interpreting and understanding the numerical results. I find the study of interest for a wide readership in Geosciences and, hence, find the study suitable for publication. I have several comments, with the main aim to make the model results more transparent, which the authors could consider.\nThe applied lithosphere deformation model represents an end-member scenario in which the deformation of the crust and lithospheric mantle is dominated by brittle-frictional processes. The maximal friction angles of both crust and mantle lithosphere are only 15 degrees (providing a deeper brittle-ductile transition compared to a friction angle of 30 degree) and only strain weakening due to a reduction of the friction angle with strain is considered which furthermore generates a permanent weakening. No other weakening mechanisms are considered, such as due to grain size variation or shear heating, and no transient weakening mechanisms are considered which would allow for healing and re-strengthening of shear zones and faults that were active during rifting. Since weakening mechanisms other than friction-weakening and healing most likely occur in nature, it is not clear how representative and robust (with respect to model variations by including other weakening mechanisms and healing) the presented results are. Ideally, also another end-member model would be presented with different weakening mechanisms and non-permanent weakening to test how dependent the first order results are on the particular model assumptions. I can imagine that such additional models are outside the scope of the study, but I maybe the authors can at least explain the particularities of their model and discuss the potential impact of different weakening mechanisms and healing on their results.\nFor the force analysis, it would be useful to also know the values of the depth-integrated strengths, for example F_weak and F_strong (in N/m ; discussed in the supplementary information) and how they vary laterally because these values could then be compared to other estimates for these strengths, such as plate driving forces (e.g. slab pull, ridge push or force balance at continental plateaus). It would be also useful to have figures that show the cross section of the models with the corresponding deviatoric stress field to be able to see the stress magnitudes and ideally the brittle-to-ductile transition in the models. Similarly, it would be useful to have additional figures of the cross sections that indicate with a colormap which regions deform in the brittle-plastic regime and which ones in the viscous regime. This can further help to compare the model results with field observations or earthquake observations. Such figures could be put in the supplementary material. Currently, there are no figures that show any stress fields, although they are the fundamental values used in calculating the forces, and it is not clear how the brittle-frictional and viscous deformations are distributed during orogeny.\nSome minor comments and questions:\nLine 75 : Please provide the formula for the Bm number. Also, please provide the values used in the simulations; Bm > 1 is vague.\nDoes the landscape evolution model has a significant impact on the collisional model evolution or would a simpler erosion model yield similar results?\n2.2.1. : As mentioned, it would be useful to know also the values of the vertically integrated strengths to compare them with plate driving forces (in N/m).\n2.2.2. : The model assumes significant permanent weakening (minimum friction angle of 2 degree in both crust and mantle) and no healing. Hence, it is expected that in the presented models the permanent weak zones developed during rifting control the subsequent collision stage.\nEquation (3) : It seems that shear heating due to dissipative deformation (friction and viscous shearing) is not considered? If associated shear heating by dissipative deformation is not considered in equation (3), then energy is not balanced in the model.\nEquation (4) : The model is 600 km deep. Why density is not also a function of pressure?\nLine 439 : Strain weakening by reducing the friction angle is usually mesh-dependent if no regularization method is applied. How is the mesh-dependence handled in the numerical model?\nLine 466 : Please provide the formula for the Bm number.\nLine 505 to 509. When the lithosphere is subducted into the sub-lithospheric upper mantle, how can this parameterization of k then produce a correct temperature field around the subducting lithosphere? This parameterization may work well during the extension phase, but how can it be justified for the subduction process?\nLine 534 to 537 : This limited effect is controlled by the model assumptions because only a frictional strain weakening mechanism is considered which causes permanent weakening. No transient weakening mechanisms are considered which would allow healing and re-strengthening rift-related faults and shear zones. Therefore, the limited impact of a cooling phase is model related, but not related to natural processes which can involve healing mechanisms.\nBest regards,\nStefan Schmalholz"}, "author_response": "Response to Reviewers\nSebastian G. Wolf, Ritske S. Huismans, Josep Anton Muñoz, Dave A. May\nThis document concerns our revision of the paper entitled Rift linkage and inheritance\ndetermine collisional mountain belt evolution. We thank all three reviewers for very positive,\nconstructive, and useful comments that helped improve our manuscript. We accommodated\nall changes suggested by the reviewers. Below follows a point-by-point reply to the reviewer\ncomments. Comments are in black font with serif; our answer is written in blue font without\nserif. Lighter blue text without serif highlights changes to the manuscript text.\nReviewer #1\nThe paper presents new, state-of-the-art numerical scenarios that combine a 3D\nthermomechanical model with a landscape evolution model to simulate the influence of\nrifting geometry on the evolution of a collisional mountain belt. The authors explore\nscenarios with different offsets between two rift branches, demonstrating that the orogenic\nstructure is controlled by the rifting geometry. Additionally, they present an analytical\napproach that offers a convincing physical explanation for the numerical results, highlighting\nthe conditions under which rift branches with different offsets may connect.\nI recommend publication of the article after the following points are addressed:\nWe thank the reviewer for very positive, useful and constructive comments that helped\nimprove our manuscript.\nComment 1: In scenario M2, which involves a large offset of the weak seed (400 km), the\nsimulation evolves into one-sided subduction with an essentially two-dimensional\n(cylindrical) geometry. As described by the authors, this occurs because the eastern rift\nbranch is aborted, and the western basin propagates southward, resulting in a linear rift\ngeometry. However, this result likely depends on the size and rheology of the weak seeds: if\nboth branches contain large, sufficiently weak seeds, rifting might develop on both the\nwestern and eastern sides, potentially leading to a successful rifting scenario with a large\noffset. It is probable that the small size of the weak seed contributed to the failure of one\nbranch, making the resulting scenario redundant with scenario M1.\nWe agree with the comment by the reviewer. As shown by the models and the force\nbalance analysis, a linked rift forms if the force required to deform the connection\nbetween the two o>set rift basins (D’; Fig. 6a) plus the force required to deform the\nbasin with length (Dseed) that is being linked to is smaller than the force required to link\nto the model side (Db). If there would be weaker and possibly larger seeds, that reduce\nthe integrated strength in the o>set basins to be smaller than F_weak, then the force\nrequired to deform Dseed would be smaller and linkage could occur with larger basin\no>sets. The comment from the reviewer is therefore fully consistent with our analysis\nand modelling results. To clarify, we added the following sentence to the supplementary\nmaterial of the manuscript:\nbasin o:sets. A smaller F_weak could be related to very weak inherited weaknesses,\nand could be modelled with weaker and larger initial weak seeds.\nWe note though that we prefer to use a simple model setup with relatively small weak\nseeds that are fully strain weakened. This setup prevents over-constraining the models\nand provides the biggest amount of “freedom” for the models to develop. We do not\nbelieve that running more models with artificially even weaker and bigger seeds\nprovides more insight beyond the produced modelling results.\nComment 2a: In scenario M3, which uses an intermediate offset of 200 km between weak\nseeds, the simulation produces an interesting subduction pattern with opposite polarity,\nresulting in a complex, three-dimensional orogenic system (Figure 4c). My question concerns\nthe robustness of this result: is opposite-polarity subduction a systematic outcome for\nintermediate-offset (~200 km) scenarios with simple weak seeds, or is it a fortuitous result?\nGiven that the weak seeds in these scenarios are symmetric in East–West cross-sections, and\nthat 2D thermomechanical models with similar seeds do not allow for a priori prediction of\nsubduction polarity (since small perturbations in the thermal and compositional structure of\nthe lithosphere can change the polarity) the observed pattern may not be systematic.\nWe thank the reviewer for this comment, that we deliberately split into two for simplicity\n(Comment 2a and Comment 2b below).\nYes, this scenario is systematic and not fortuitous. We tried di>erent amounts of\nextension and di>erent amounts of o>set, but o>sets > ~30 km consistently lead to\nopposite polarity subduction. Only if there is an existing directional weakness that\nprovides a gain so that the longer same-polarity deformation path is taken, the system\nwould develop into same-polarity subduction (M4). The force balance analysis\nquantifies this behaviour and provides an explanation why the polarity flip is happening.\nThe nature of the opposite polarity subduction is one of the interesting key outcomes of\nthis study. To highlight that this result is systematic, we added another model with o>set\n100 km (SM2) to the supplementary material, that also develops opposite polarity\nsubduction, as expected.\nWe added the following sentences to the results section:\nSM2, with an intermediate o:set of 100 km between simple weak seeds exhibits very\nsimilar behaviour to M3, with opposite polarity subduction as a result of inner fault-\nlinkage. Additional models with similar o:sets all develop into opposite polarity\nsubduction. These results show that opposite polarity subduction in intermediate-o:set\nsimple seed models is a systematic and repeatable model result.\nWe refer the reviewer to section 2.2.2, the force balance analysis, that gives a physical\nexplanation for subduction polarity showing that the model results are not only\nconsistent but also quantifiable.\nAnd yes, we also made the same experience that a simple weak seed does not allow for\na priori prediction of subduction polarity in 2D models; these models behave like a\nstraight rift system without o>set (Model 1). However, as our study shows, when\nconsidering small o>sets and 3D linkage, subduction polarity is a priori predictable. We\nthink that this part is a very nice result of this study as it provides, to our knowledge for\n2\nthe first time, a quantifiable understanding for subduction polarity development\nfollowing rift basin inversion. On a side note, from a modellers perspective who has\nbeen frustrated with the a priori unpredictability of 2D models, this result is especially\nsatisfying as it provides physical reasoning. We believe that this topic does not need\nadditional textual changes as it is covered by section 2.2.2.\nComment 2b: Additionally, the numerical setup explores only a single offset per scenario.\nWhat would be expected if multiple branches with intermediate offsets (~200 km) were\npresent in the domain? As the authors convincingly argue, polarity is sensitive to inherited\ngeometries of the weak zones. Therefore, I suggest that the offset plays only a secondary (or\nnegligible) role in controlling subduction polarity.\nThe reviewer raises an interesting question about what would happen with several o>set\nweaknesses. In this case we can only speculate and use simple logic based on the force\nbalance analysis. Deformation will always choose the path of least resistance.\nAssuming that there are no other inherited weak zone geometries that govern\nsubduction polarity, we would expect these systems to behave as if composed of\nseveral basin pairs similar to those examined here. We added the following text to the\nlimitations section of this manuscript:\n“We investigated the simplified case of inverting two o:set rift basins. In natural\nsystems, however, multiple o:set basins are likely to occur. Based on our model results\nand force-balance analysis, we predict that inverting such systems would behave as if\ncomposed of several basin pairs similar to those examined here, potentially involving\nmultiple subduction polarity flips in the absence of dominant inherited weaknesses\nguiding deformation. Future research should explore this topic.\nWe agree that polarity is sensitive to inherited geometries of the weak zones, as shown\nby the models and analysis. The physical intuition of the reviewer might be correct that\nin mountain belts on Earth inherited weaknesses play a dominant role in determining\nsubduction polarity – otherwise we would likely see a subduction polarity switch in\nmost orogens. Our study quantifies this process and shows when subduction polarity\nswitches, how structures and topography evolve in response to the polarity of\nsubduction, and how weak inherited weakness need to be to change subduction\npolarity. We think this is new and a useful contribution. To highlight the point made by\nthe reviewer, we added the following sentence to the second to last paragraph of the\nmanuscript:\nComparison of model results with natural systems shows that opposite-polarity\nsubduction is rather the exception than the norm, indicating that inherited directional\nweaknesses play an important role in determining subduction polarity in inversion\norogens on Earth.\nMinor points:\nThe authors use cardinal points to indicate features in the numerical scenarios. This is a\nhelpful approach, especially for 3D models with complex geometries, as it facilitates\nunderstanding of the numerical results. However, since the scenarios are generic and not\nintended to represent any specific orogenic system on Earth, it is important to define the\n3\ndirection of the initial extension. This would improve the reader's comprehension of the\nresults. I suggest adding a clarification in lines 79–81, such as: \"We assumed the direction of\nextension to be East–West.\" Although this is mentioned in the Methods section and indicated\nin Figures 2–5, including it in the main text would enhance clarity.\nWe thank the reviewer for this comment and agree. We update the sentence in lines 79-\n81 as follows (addition in bold font):\nEach model is extended for 12 Myr with a rate of 1 cm/yr in East-West direction,\nfollowed by a reversal of the velocity boundary conditions and shortening with the same\nrate.\nContinuing on this topic, the description of the Model 3 results is unclear. The authors state in\nlines 119–120 that \"the right fault in the northern basin and the left fault in the southern basin\nrecord more offset and therefore also slightly higher rift-flank topography.\" However, the\ncombination of cardinal points with terms such as \"right\" and \"left\" is confusing. Please\nconsider replacing these with \"eastern\" and \"western\" for clarity.\nWe agree and changed all \"right\" and \"left\" to “west(ern)” and “east(ern)”.\nModel names appear in both lowercase (e.g., line 109) and uppercase (e.g., lines 221–222).\nPlease standardize their formatting throughout the manuscript.\nWe thank the reviewer for very careful reading! We agree and changed the inconsistent\ncases. Now the models are named either Model 1 or M1 etc.\nReference numbers in lines 461–469 should be placed within square brackets [ ] to avoid\nconfusion with exponents.\nWe agree and wrapped these citations with square brackets.\nReviewer #2:\nThis paper presents the results of a series of sophisticated 3D thermomechanical simulations\nof lithospheric extension followed by shortening. The objective is to consider how offset or\nlinked rift segments affect subsequent orogen growth. The thermomechanical software is, in\naddition, coupled with a 3D surface processes modelling package that simulates erosion from\ngrowing relief and deposition of sediments in developing basins. 3D models simulate\northogonal, segmented rifting that is allowed to develop to the point of continental breakup\nfollowed by convergence to juvenile (or partial) inversion, mature (or complete) inversion\nand beyond. Extensional and compressional stresses are parallel. The simulations consider\nthe effect of offset distance between rift segments on orogen evolution. The authors conclude\nthat this factor influences initial growth of the mountain belt being strongly impacted by the\noffset distance of rift segments. Later orogen growth is governed by subduction polarity,\nwhich itself is mainly dictated by the asymmetry of the main conjugate rift faults. The model\ninput parameters and uncertainties are described at the end of the paper. Simulation results\nare compared to three mountain chains, the Causases, Atlas and Pyrenees, which are in turn\ninterpreted in terms of rift geometry inheritance.\nThis is an excellent paper from one of the leading thermomechanical modelling groups in\ngeosciences. The methodology is sound and well established with the innovative combination\n4\nof two softwares to provide a more complete 3D earth system simulation. It is well written\nand illustrated with impressive videos of model development. The authors state clearly the\nuncertainties and limitations of their models. The conclusions and comparisons with natural\nexamples are valid, even if some queries need addressing. The results will provoke much\ndiscussion and interest amongst a wide readership interested in Wilson cycle lithospheric\nprocesses. It is also a nice demonstration of cutting edge multidimensional thermomechanical\nmodelling combined with surface processes modelling. I feel that the surface processes\naspects of the modelling do not get the attention that perhaps they merit. However, this may\nperhaps be the subject of a later paper. I recommend this paper for publication once the points\nbelow and other minor comments on the pdf are addressed.\nMary Ford,\nUniversité de Lorraine, Nancy.\nWe thank the reviewer Mary Ford for a very useful and constructive review and very\nprecise comments that helped improve the manuscript.\nConcerning the surface processes aspects: We agree that the surface processes part of\nthe article is very exciting and should be followed up with future work.\nWe thank the reviewer for providing comments directly in the manuscript PDF and for\nalso exporting these comments in the list below. Some minor textual comments that\nrequired only very small changes to the text are not listed below. We have implemented\nall of these minor textual suggestions. In case we did not follow them, they are listed\nbelow.\nLines 53_55: Some analogue modelling studies do actually address these issues. What is\nmissing from the majority of analogue models, it seems to me, is a full crustal/lithospheric\ntreatment of these issues.\nWe agree with the reviewer. We already write that analogue models only focus on\ncrustal scale in the sentences before. To be even clearer and specific, we changed the\nrespective sentence to:\nHence, the importance of a 3D rift setting with stepping, o:set, and linked rift basins on\ninversion tectonics remains unknown when including a consistent upper mantle-scale\nperspective.\n79: Can you explain why this rate and period were chosen? This appears to thin the\nlithosphere to the point of mantle exhumation and breakup?\nWe thank the reviewer for this comment. The rate was chosen because it represents\ntypical plate velocities on Earth; and yes, the period was chosen to create full crustal\nbreak-up. We tested shorter periods as well, which had no significant e>ect. We added\nthe following sentence:\nThis period of extension leads to full crustal break-up; the applied rate of\nextension/shortening represents typical plate velocities on Earth.\n81: State how long the shortening phases were and what the rate of shortening was. The\nduration of shortening appears to be variable (roughly 44 Myrs to 63 M Yrs) but always a lot\n5\nlonger than 12 Ma. Give the reasoning used for the duration and rate of the shortening ?\nAttaining a set amount of shortening for example ?\nThe rate was again chosen to represent typical values of plate movement on Earth; the\nperiod was chosen to create at most an intermediate-size orogen larger than the\nEuropean Alps. The rate is accounted for by the previous comment; we added the\nfollowing sentence, also addressing model duration:\nModels run for 63 Myr but were stopped when the evolving orogen reached the Eastern\nor Western model side. This period and rate of shortening creates orogens reaching\nintermediate size, slightly larger than, for instance, the European Alps.\n95: southern?w\nWe agree and changed to western\n96: southern? Why use left and right? Why not stick to the cardinal directions of the compass\nas on Figures 2-5 of models and elsewhere in text?\nWe agree and changed all “left” and “right” to cardinal directions throughout the\nmanuscript.\n96: What do you mean by \"Mature\"? Fully or completely inverted? Maybe explain what you\nmean here and/or give reference. Same point for Juvenile inversion.\nWe thank the reviewer for detailed reading, agree, and added the following explanation\nto the end of the first results paragraph:\nThroughout the manuscript, we refer to ‘juvenile’ and ‘mature’ rifts or orogens. Juvenile\nrifts contain unconnected sub-basins, unlike mature rifts where these are linked (Fig.\n2a). Juvenile orogens lack clear subduction polarity and are only composed of inverted\nextensional structures (Fig. 2b), whereas mature orogens exhibit distinct lithospheric\nsubduction and additional crustal shortening (Fig. 2c).\n119 – western fault\nWe agree and changed accordingly.\n125-6 : Oppositely dipping subduction - how do these link ? Can you give some detail here of\nhow this works in the models?\nWe observe that the linkage between the oppositely dipping slabs is accommodated by\nsuccessive separation of discrete pieces of lower crust and lithospheric mantle in the\ntransitional domain stemming alternatingly from both slabs (see Fig. 4c). We describe\nthese structures and the associated process two sentences later. To clarify, we changed\nthe respective text to:\nInteraction of the two subducting slabs is characterised by several large pieces of lower\ncrust and lithospheric mantle sheared o: from the opposite-polarity subducting\nlithospheres, and creating a topographic low in the transfer zone (Fig. 4c).\n6\n129 : not clear here. Reverting is not the appropriate word as it means to return to a former\nstate. Inverting or switching would perhaps be clearer.\nWe agree. We think that reversing is the right term, as the plan-view shape is reversed.\nWe changed accordingly.\n178 : But you say you have already assumed this in line 166.\nWe agree and changed to:\nTo summarize, model behaviour can be successfully quantified by assuming that\ndeformation follows the path of least resistance, such that the integral of the depth-\nintegrated lithospheric strength along any given path is minimised.\nYes, we agree. We changed this sentence to:\nThese force minimisation principles also explain same-polarity subduction in cases of a\nstraight rift system without significant o:set, as in Models 1 and 2.\n227: provides.\nWe think provide without ‘s’ is correct, as provide refers to inherited weaknesses\n273-275 : Incomplete sentences... punctuation problem?\nWe agree and removed the subsentence:\n, and uplift continuing today.\n288 : EHA is not defined. Also this sentence is not clear and needs a rewrite.\nIt was defined in the sentence before. However, we rewrote the respective sentences to:\nLooking at the large-scale trends of the Atlas, we find that the Western branch, the so-\ncalled High Atlas, trends SW-NE, with a change in strike at the Eastern High Atlas (EHA)\ntowards W-E, which changes again after ∼300 km to a SW-NE strike. The northern\nbranch termed Middle Atlas (MA) has yet another orientation with a more Northerly\ndirection (Fig. 7b).\n295 : Compare to what stage of these models ? You could define this in the videos\nWe agree, the sentence now reads as follows:\nThe mild inversion of extensional basins with low amounts of shortening and incipient\nsubduction indicate that the Atlas is a juvenile inversion orogen where our models show\nthat subduction polarity does not yet play a significant role, as shown in Models 3 and 4\nafter ~20 Myr of model evolution (Fig. 4b and Fig. 5b )\nWe think that this is clear definition. We did not label this stage in the videos so as to not\ncreate the wrong impression that a specific model represents a specific orogen.\n7\n304 : Its not clear how do the Atlas mountain belt highlights this point that basin linkage and\nassociated changes in mountain belt structure occur at different scales during inversion.\nWe thank the reviewer for detailed reading and specified to:\nThe Atlas is a fitting example of a juvenile inversion orogen where subduction polarity\nhas no primary influence, and highlights that basin linkage and associated changes in\nmountain belt structure during inversion occur on di:erent spatial scales, from small-\nscale steps (10s of km) to large-scale branching (100s of km).\n310 : The Mauleon basin is a small (but important) depocentre within the western Pyrenean\nsystem (SW of Pau). It is misleading and inappropriate here to refer to the whole Pyrenean\nrift (ss) by this name.\nWe agree and changed the sentence to:\nIn the E, between France and Spain, the Pyrenees sensu stricto (s.s.) formed by\ncontinent-continent collision that inverted the Pyrenean extensional system and also\npropagated into the un-thinned Iberian pro-side (Fig. 1d, Fig. 7c)6,7,50.\n310 ‘rifted margin’. Is the Pyrenean depocentre a rift or a rifted margin? So far you have been\ntalking about rifts. If the latter then which margin are your referring to?\nWe agree and changed to extensional system, see sentence above\n310 : Suggest propagated would be clearer than ‘progressed’.\nWe agree and changed accordingly, see sentence above.\nWe agree and changed accordingly.\n317 : Suggest rewording : “Exhumed fragments of mantle rocks on the N-flank (retrowedge)\nof the orogen, are found embedded in metamorphosed Mesozoic sedimentary cover between\ninverted crustal extensional blocks”. I suggest citing not just Teixell et al. (2018) but also, for\nexample, Angrand et al. (2020, J. Geol Soc London) and/or Asti et al. (2022 ; Earth Science\nReviews)\nWe thank the reviewer for a detailed suggestion for rewording and used the proposed\nsentence. We also added a citation to Angrand et al. 2022 J Geol Soc London. We are\nalready at the limit of possible citations and therefore only chose to cite the field study\nover the review paper.\n320 : In model of Wang et al. the rise of the mantle lithosphere occurred during hyper-\nextension. This configuration was essentially preserved due to very low convergence in this\narea of the westernmost north Pyrenees.\nYes, we agree. We are aware of this interpretation, which is not in conflict with the\nsentence. See also answer below to line 393.\n8\n350: inversion?\nWe prefer to retain the word deformation instead of inversion.\n353: no s in highlights\nWe thank for careful reading, but think that the 3rd person s should be retained:\n...the same-polarity orogen…. highlights the complexity\n392: Given that the greater Causasus and Pyrenees are both interpreted as mature same\npolarity orogens, it would be very relevant to discuss the differences between the two orogens\nin terms of the principal aspects of your models (amount of initial extension? amount of\nshortening? orientation of linkage zones and branches of rifts?)\nWe agree that this is interesting. There are many di>erences between both systems.\nOne of the key di>erences is that in the Pyrenees lateral di>erences in orogen structure\nrelated to rift inheritance are overprinted by di>erences in the amount of shortening, as\ndescribed in the text. We expanded on this sentence as follows:\nThe Pyrenees are a mature same-polarity orogen, in which lateral di:erences in orogen\nstructure related to rift inheritance are overprinted by di:erences in the amount of\ncrustal shortening, which is likely one of the major di:erences when compared to the\nGreater Caucasus, another mature same polarity orogen.\nWe agree that it is very interesting to discuss further di>erences, also compared to the\nAtlas, but this discussion is beyond the scope of this study (which is at its length limit),\nand also beyond the scope of this paragraph, which provides a short summary of the\ndi>erent orogenic systems that we compare to and extends to other orogenic systems.\n393 : I dont think this explains the presence of large high density bodies. The offset between\nthe Basque Pyrenees and the main Pyrenees ss can certainly be explained by inversion of a\nlinked offset rift basin. But the presence of high density bodies is surely related to\nexhumation of deep lithosphere during hyper-extension and their subsequent incorporation\ninto the orogen. In other words their presence is related primarily to the initial style and\namount of extension in the rift system...\nWe thank the reviewer for this comment concerning the two large high density bodies.\nWe agree that other studies have argued that these high-density bodies are related to\nexhumation of deep lithosphere during extension that was subsequently incorporated\ninto the mountain belt at the inversion stage. However, following this interpretation, we\nwould expect to find these types of large high density bodies along-strike the whole\nPyrenees and not only in two specific locations. Yet, we don’t. Interestingly, these\nlocations coincide with rift linkage zones. Our models show that inverting an o>set\nlinked rift basin can produce these structures in the linkage zones. We think that this is\nan intriguing result as it provides a very simple alternative explanation. We though also\nbelieve that this feature is not the main point of this article and deserves follow-up\ninvestigation. To soften our statement, we changed it to:\nObserved large high-density bodies at shallow depth in the Pyrenees can furthermore\npotentially be explained by the inversion of a linked o:set rift basin.\n9\n409 : and the amount of shortening?\nWe agree and added the proposed phrase.\n411-414 : Is there any natural mountain belt that demonstrate this behaviour?\nYes, the Hindu-Kush Pamir region is potentially such a mountain belt. The last two\nsentences of the previous paragraph discuss this, here repeated for clarity:\nSubduction polarity reversal has also been proposed for mountain building in the Hindu\nKush - Pamir region57. The distribution and propagation-direction of shortening57, as well\nas lower topography in the transition zone between both opposite-polarity slabs fits well\nwith the characteristics of an opposite polarity subduction orogen as shown in Fig. 8b.\n481 : This section and Figure 9 would be better in my opinion early in the paper to allow the\nreader to understand the model results. For example, the colours codes for results figures are\ngiven on figure 9.\nWe appreciate the suggestion by the reviewer to move this model setup section into the\nmain body of the article. However, we think that adding this section up front makes the\narticle unnecessarily more technical. Therefore, we would like to retain the section and\nmodel setup at the end of the article. However, to make sure that the reader is properly\npointed to the Methods section, we added to the first paragraph in section 2.1:\n… See Methods section for a detailed model setup description.\nWe also added model materials legends to the model results figures, which will help for\nclarity.\nFigures\nFigure 1. I have located the actual Mauleon basin on figure 1d. As this is a well known and\ndefined local basin, its name is not appropriate for the whole Pyrenean rift system.\nWe agree and changed the name to Pyrenean basin, as proposed by the reviewer in the\ntext.\nFigures 2 to 5. It would be helpful to put the compass points on at least the first model image\nNeed a key or labels for colours on the cross sections.\nCan you separate the three cross sections as complex details are not clear in the areas where\nthey overlap ? Same comment for all cross sections.\nWe thank the reviewer for detailed comments on the figures. We added a compass-\narrow to the first landscape plot. We also want to note that the cardinal directions are\nalso displayed on the model setup cartoons in the top right corner of the result figures.\nWe added a legend for the material colours to all plots. We also moved the cross\nsections a bit. Note that the cross sections and their distance represent the real model\ndimensions. We want to retain this feature as it is most true to the results. We think that\nthe model figures in addition to the animations show the model development very well\nnow, including details.\nC for the cross sections ?\nWe agree and changed to I, II, III\nReviewer #3:\nThe manuscript titled “Rift linkage and inheritance determine collisional mountain belt\nevolution“ by Wolf and co-authors presents 3D numerical models of lithosphere extension\ndirectly followed by convergence. The main aims of the modelling are to better understand\nlinkage of rift segments during extension and the impact of the rift geometry on the\nsubsequent convergence, particularly the subduction polarity. The applied numerical model\nconsists of a lithosphere deformation model coupled to a landscape evolution model. The\nstudy includes a very interesting application of force analysis to explain linkage of rift\nsegments and the subduction polarity during subsequent convergence and rift inversion.\nFurthermore, the study discusses the application of the model results to the Greater Caucasus,\nAtlas and Pyrenees. There are still many open questions concerning the impact of 3D features\non rifting and rift inversion and the manuscript, hence, presents a very interesting study on\nthe application of deterministic mathematical models to better understand these geodynamic\nprocesses. The presented results demonstrate the importance of the pre-orogenic 3D rift\ngeometry on the subsequent orogeny related to the inversion of the rift basins. The authors\napply a force analysis to explain the 3D model results which is very useful for interpreting\nand understanding the numerical results. I find the study of interest for a wide readership in\nGeosciences and, hence, find the study suitable for publication. I have several comments,\nwith the main aim to make the model results more transparent, which the authors could\nconsider.\nWe thank the reviewer, Stefan Schmalholz, for a positive and very helpful review. We\naccommodated the suggestions from the reviewer to make the model results more\ntransparent, as outlined below.\nThe applied lithosphere deformation model represents an end-member scenario in which the\ndeformation of the crust and lithospheric mantle is dominated by brittle-frictional processes.\nThe maximal friction angles of both crust and mantle lithosphere are only 15 degrees\n(providing a deeper brittle-ductile transition compared to a friction angle of 30 degree) and\nonly strain weakening due to a reduction of the friction angle with strain is considered which\nfurthermore generates a permanent weakening. No other weakening mechanisms are\nconsidered, such as due to grain size variation or shear heating, and no transient weakening\nmechanisms are considered which would allow for healing and re-strengthening of shear\nzones and faults that were active during rifting. Since weakening mechanisms other than\nfriction-weakening and healing most likely occur in nature, it is not clear how representative\nand robust (with respect to model variations by including other weakening mechanisms and\nhealing) the presented results are. Ideally, also another end-member model would be\npresented with different weakening mechanisms and non-permanent weakening to test how\ndependent the first order results are on the particular model assumptions. I can imagine that\nsuch additional models are outside the scope of the study, but I maybe the authors can at least\nexplain the particularities of their model and discuss the potential impact of different\nweakening mechanisms and healing on their results.\nweakening as source of weakness. We only focussed on this mechanism, because in\nfield studies, we can observe that brittle-frictional fault zones that formed during basin\nformation are reactivated during orogenesis, hence likely stay weak even many Myr after\nextension. Prime-examples are the Pyrenees, inverted Salta rift basin in the Andes, or\nthe Atlas, but really these structures are similarly also observed in all inversion orogens.\nHowever, we fully agree with the reviewer that there could be other mechanisms\ncreating permanent or also transient weakening, like grain-size variations, shear\nheating, or metamorphic phase changes creating a dominant metamorphic fabric.\nIndeed, the directional weakness needed to create same-polarity subduction in the\nPyrenees might be related to a pre-existing Hercynian metamorphic fabric, as described\nin the text. We believe that running simple models and properly quantifying them, as\ndone in the study, is a good approach and provides transferability of results. We believe\nthat running additional models with di>erent types of weakening is beyond the scope of\nthis study. A proper exploration of all di>erent types of weakening, including (partial)\nhealing would be work for several additional articles. The models presented in this study\ncan serve as a comparative baseline for such additional studies. To account for this\ndiscussion, we added the following paragraph to the model limitations section:\nWe only apply frictional-plastic softening as a weakening mechanism in our models.\nAlthough the reactivation of weak frictional-plastic shear zones is characteristic of\ninversion orogens - such as those examined in this study - we note that additional\npermanent or transient weakening processes may play a role during extension and\norogenesis, like grain-size variations, shear heating, or metamorphic phase changes.\nFuture modelling e:orts should explore these additional mechanisms. The model\nresults and force balance analysis presented here can serve as a comparative baseline\nfor such investigations.\nFor the force analysis, it would be useful to also know the values of the depth-integrated\nstrengths, for example F_weak and F_strong (in N/m ; discussed in the supplementary\ninformation) and how they vary laterally because these values could then be compared to\nother estimates for these strengths, such as plate driving forces (e.g. slab pull, ridge push or\nforce balance at continental plateaus). It would be also useful to have figures that show the\ncross section of the models with the corresponding deviatoric stress field to be able to see the\nstress magnitudes and ideally the brittle-to-ductile transition in the models. Similarly, it\nwould be useful to have additional figures of the cross sections that indicate with a colormap\nwhich regions deform in the brittle-plastic regime and which ones in the viscous regime. This\ncan further help to compare the model results with field observations or earthquake\nobservations. Such figures could be put in the supplementary material. Currently, there are no\nfigures that show any stress fields, although they are the fundamental values used in\ncalculating the forces, and it is not clear how the brittle-frictional and viscous deformations\nare distributed during orogeny.\nWe thank the reviewer for this very useful comment. We now report the required values.\nHere a few notes about these forces. We considered computing F_weak and F_strong\nand particularly their integration along strike of the modelled orogen directly from the\ndeviatoric stress fields of the models. However, we could not find a good criterion to do\nthis computation automatically, for instance based on a strain rate-cuto>. Rather, we\nfound that computing F_weak and F_strong analytically and integrating along simple\n12\ndistances as shown in the block diagrams (Fig. 6) gives a very good approximation of the\nmodel behaviour. For the computations we used a constant strain rate of 1e-14 1/s that\nis typical for the deformation rates in our models, assumed lithostatic pressure, and\nassumed a steady-state thermal field corresponding to the initial setup of the model.\nUsing these simplifications, we found that the transition from same-polarity subduction\n(seed o>set 0 – 30 km) to opposite polarity subduction (seed o>set > 30km) can be\npredicted very accurately by the analytical solution when compared to the models.\nSimilarly, the predicted transition from basin-linkage to no linkage during extension (at\n292 km) also fits very well with what we observe in the models. This tells us that despite\nthese simplifications, the analysis works well. Assuming lithostatic pressure is a\nsimplification, as it neglects tectonic overpressure stemming from gravitational\npotential energy variations. We tested the e>ect of overpressure through simple\nmultiplication of the lithostatic pressure. We found that it had no first-order influence\non the analytical results, but would of course increase F_weak and F_strong. For\ninstance, a factor 2 on the lithostatic pressure only changed the basin linkage transition\nfrom 292 km to 290 km, which is secondary. We believe that presenting a simple\nanalytical solution that is easy to understand and successfully predicts model\nbehaviour is useful and better than a very complex computation that has no gain in\npredictive power. We added the following text to the supplementary material:\nAssuming lithostatic pressure conditions, a constant typical strain rate of 1e-14 s−1, and\na steady-state thermal field corresponding to the initial thermal field of the models gives\nF_weak = 2.8e12 N/m and F_strong =12.1e12 N/m. These assumptions neglect tectonic\noverpressure, e.g. related to gravitational potential energy variability, which can be\nsignificant [Schmalholz et al., 2014]. Sensitivity tests including substantial overpressure\nshowed no first-order e:ects on the resulting relationships for basin linkage or\nsubduction polarity, but would increase F_weak and F_strong. These tests suggests that\nthe analytical results are robust within the framework of our assumptions, and show\nthat the reported values of F_weak and F_strong are minimum values.\nThe citation to [Schmalholz et al., 2014] was also added to the main manuscript.\nWe agree with the reviewer concerning the additional plots and added model plots of\nthe deviatoric stress fields as well as plots of the deformational regimes. These plots\nare shown in the supplementary material as figures S1-S4. We think that these are nice\nand useful plots. We note that the fields were computed during post-processing, with\none average value per cell. In pTatin3D, properties are computed at each Lagrangian\nmarker particle and then averaged over a cell. Because of this averaging, we expect very\nsmall di>erences between the fields presented here and the fields used for\ncomputation in pTatin3D. The latter three sentences are also stated in the\nsupplementary material. We also added the following sentence to the Results section:\nSupplementary model output showing deformational regimes and deviatoric stresses is\npresented in Supplementary Figures S1 - S4.\nSome minor comments and questions:\nLine 75 : Please provide the formula for the Bm number. Also, please provide the values used\nin the simulations; Bm > 1 is vague.\nnumber of variables so that we moved the computation to the Methods section and\nadded a sentence to line 75 referring to the methods section. See also answer below.\nDoes the landscape evolution model has a significant impact on the collisional model\nevolution or would a simpler erosion model yield similar results?\nThis is an interesting, yet rather hypothetical discussion, because we did not try with\nother landscape evolution models. Bm = 3.5 implies that the model is “tectonics\ndominated”, i.e. the orogen is widening and in height limited by the lithospheric strength\nand not the e>iciency of erosion, which seems to be the norm on Earth. Therefore, we\nspeculate that the observed subduction pattern and first order evolution of mountain\nbuilding is likely also happening in models without surface processes or with a simpler\nerosion model, e.g. approximated by di>usion, as long as tectonics is “dominant” over\nerosion. However, the sedimentation pattern and diagnostic landscape evolution\npattern with river valleys being a marker of deformation requires a landscape evolution\nmodel similar to the one used here. In that sense, having a landscape evolution that\n“properly” parametrizes erosion and sediment transport is important and very useful.\nAlthough interesting, we don’t think this discussion needs to be in the manuscript, as it\nis tangential to the work and rather hypothetical.\n2.2.1. : As mentioned, it would be useful to know also the values of the vertically integrated\nstrengths to compare them with plate driving forces (in N/m).\nYes, we agree and report these in the supplementary material, see also answer above.\nWe added the following to this paragraph to highlight the values:\nSee supplementary information for detailed information on the definition of the di:erent\nforces involved, including characteristic values.\n2.2.2. : The model assumes significant permanent weakening (minimum friction angle of 2\ndegree in both crust and mantle) and no healing. Hence, it is expected that in the presented\nmodels the permanent weak zones developed during rifting control the subsequent collision\nstage.\nWe agree with the reviewer. We added a discussion of healing in the model limitations\nsection, see answer above. Specifically, we also added the following sentence to the\nfirst paragraph:\nFurther studies should nonetheless investigate the e:ect of post-rift cooling in detail,\nincluding an investigation of transient weakening e:ects including healing mechanisms.\nEquation (3) : It seems that shear heating due to dissipative deformation (friction and viscous\nshearing) is not considered? If associated shear heating by dissipative deformation is not\nconsidered in equation (3), then energy is not balanced in the model.\nYes, we agree. This shortcoming is discussed in the model limitations section. We did\ntest the importance of shear heating in very similar 2D models, and could not find a first\norder e>ect on model results. We therefore expect no first order changes to the model\n14\nresults when including shear heating in the models presented here, but would do so in\nthe future.\nEquation (4) : The model is 600 km deep. Why density is not also a function of pressure?\nWe opted to not include density as function of pressure, because the deformational\ndomain of interest is the upper 200 km, where changes as function of temperature and\ncomposition are much more important than changes as function of pressure, which are\nsecondary. We agree though that this could be added. We added the following\nsentence:\nWe do not include changes in density as function of pressure.\nLine 439 : Strain weakening by reducing the friction angle is usually mesh-dependent if no\nregularization method is applied. How is the mesh-dependence handled in the numerical\nmodel?\nWe agree with the reviewer that this is a shortcoming of the simplified strain-weakening\nused. We note that this parametrisation is nonetheless commonly used throughout the\ncommunity and reported in the literature. We tested di>erent strain weakening\nthresholds and found that they did not have a first-order impact on model results. We\nadded the following sentence:\nWe do not apply regularization methods to eliminate mesh-dependence of the strain-\nweakening algorithm. Tests with di:erent strain weakening thresholds did not show a\nfirst-order e:ect on model results.\nLine 466 : Please provide the formula for the Bm number.\nWe agree, and added the formula of the Bm number as well as an explanation of its\nvariables.\nLine 505 to 509. When the lithosphere is subducted into the sub-lithospheric upper mantle,\nhow can this parameterization of k then produce a correct temperature field around the\nsubducting lithosphere? This parameterization may work well during the extension phase, but\nhow can it be justified for the subduction process?\nWe thank the reviewer for this comment. As written in the text, we use this\nparametrisation to maintain the heat flow at the base of the lithosphere, prevent the\nmodel from cooling, and resembling mantle convection at high Nusselt number. We\nbelieve that this a very useful parametrisation, as it mimics a natural process and keeps\nconstant thermal conditions in the lithosphere. We agree that when modelling whole-\nmantle subduction, then this parametrisation is not needed. However, our models do\nnot exhibit deep lithosphere subduction. Rather, the lithosphere subducts at most to\ndepths of 200 – 250 km, and only reaches this depth during the late stages of model\nevolution. Therefore, we opted to use this parametrisation for this study.\nLine 534 to 537 : This limited effect is controlled by the model assumptions because only a\nfrictional strain weakening mechanism is considered which causes permanent weakening. No\ntransient weakening mechanisms are considered which would allow healing and re-\n15\nstrengthening rift-related faults and shear zones. Therefore, the limited impact of a cooling\nphase is model related, but not related to natural processes which can involve healing\nmechanisms.\nWe agree and changed the last sentence of this paragraph to:\nFurther studies should nonetheless investigate the e:ect of post-rift cooling in detail,\nincluding an investigation of transient weakening e:ects including healing mechanisms.\nWe note though that in orogenic systems that had significant amounts of post-extension\nquiescience, like the Atlas, we still observe that extensional fault zones are\npreferentially reactivated during shortening. This indicates that they do not completely\nheal. Exploring this in detail is though beyond the content and scope of this article.\nBest regards,\nStefan Schmalholz\nSchmalholz, S. M., S. Medvedev, S. M. Lechmann, and Y. Podladchikov (2014),\nRelationship between tectonic overpressure, deviatoric stress, driving force, isostasy\nand gravitational potential energy, Geophys J Int, 197(2), 680-696,\ndoi:10.1093/gji/ggu040.\nSebastian G. Wolf, Ritske S. Huismans, Josep Anton Muñoz, Dave A. May\nThis document concerns our final revision of the paper entitled Rift linkage and inheritance\ndetermine collisional mountain belt evolution. We thank all three reviewers again for very\npositive, constructive, and useful comments that helped improve our manuscript and are\nhappy that the reviewers considered our revisions suitable. We accommodated the few\nremaining changes suggested by the reviewers. Below follows a point-by-point reply to the\nreviewer comments. Comments are in black font with serif; our answer is written in blue font\nwithout serif. Lighter blue text without serif highlights changes to the manuscript text.\nReviewer #1\nThe authors have adequately addressed all of my comments and concerns in their revised\nmanuscript. I find the current version suitable for publication and recommend acceptance in\nits present form.\nWe thank the reviewer for useful comments and are pleased to hear that the reviewer\nconsiders her/his comments and concerns adequately addressed.\nReviewer #2\nI have read through the revised manuscript and the reply of authors to reviewers. The authors\nhave carefully addressed and discussed the comments and revisions suggested by reviewers.\nThe paper is now much improved in clarity and represents an important and illuminating\ncontribution to our understanding of rift inversion and orogenesis at a lithospheric scale.\nOne minor correction needs to be made on figure 1d - change Mauleon Basin to Pyrenean\nBasin as revised in the text.\nI congratulate the authors on this very nice work.\nI recommend publication\nMary Ford, Université de Lorraine, Nancy, France.\nWe thank the reviewer, Mary Ford, for detailed reading and useful comments, and are\npleased to read that our revisions are considered suitable. We agree with the reviewer\nand changed the basin name on figure 1d to Pyrenean Basin as suggested.\nReviewer #3:\nThe authors have addressed my comments and clarified the text with respect to model\nassumptions. They also added new panels in the figures in the supplementary material\nshowing the deviatoric stress which I find very useful. Only in case the authors will make\n1\nanother revision (which is not necessary), I have few minor comments to the manuscript.\nI also have two comments to their rebuttal.\nWe thank the reviewer, Stefan Schmalholz for useful constructive feedback and are\npleased to read that he considers further revisions not necessary. We nonetheless\naccommodated the minor comments by the reviewer, see below.\nMinor comments\nLine 377: “…if their offset is not too large…“; maybe the authors could quantify what is “not\ntoo large“ with a dimensionless number, for example, if the ratio of offset to crustal thickness\nis less than approximately…\nWe thank the reviewer for this comment. We quantified “what is not too large” in section\n2.2 under Continental rifting: To link or not to link. There we show that linkage depends\non distance between the weak zones and their relative lithospheric strength. It’s thus\nnot only one factor that can be easily described by a non-dimensional number. We\nagree though that this sentence can be clarified and changed it to the following\nsentence:\n(i) extensional basins with closest proximity and lowest lithospheric strength are most\nlikely to link during rifting\nLine 451: If H is only the radiogenic heat production rate, please specify this.\nYes, it is. We followed the reviewer suggestion and added the word radiogenic to the\nrespective sentence.\nLine 510: For completeness, the authors could also specify the value or criteria for the\nmechanical time step.\nWe agree and added the following sub-sentence:\n…, which is dynamically computed based on a Courant–Friedrichs–Lewy condition and\ntypically in the order of 9000 yr.\nComments to rebuttal and new text in supplementary material\nLithostatic pressure and force balance analysis: It is still not clear to me why the authors\nassume lithostatic pressure and a constant strain rate for the force calculation. Why don’t the\nauthors not use the values of deviatoric stresses calculated in their 3D thermo-mechanical\nmodel to calculate F_int?\nWe thank the reviewer for this comment, which we split in two for simplicity. As written\nin the previous reply, we tried to compute F_weak and F_strong directly from the model.\nThis is straight forward for F_strong, as it can be computed with a 1D depth-integrated\nstrength profile. It is, however, not a simple 1D profile for F_weak, which would need to\nfollow a 3D path along a shear zone. Doing this numerically is possible, but comes with\nadditional complexity, e.g. which criteria to use in choosing the path along which to\nintegrate. Therefore, we opted to use yield-strength envelopes. We then tested whether\n2\nthis simple approach can reproduce the observed deviatoric stresses when we simply\nscale the lithostatic pressure by a factor, and find that our approach is a useful first-\norder approximation. Furthermore, as written in the previous reply, the resulting values\nthat we get for instance for the maximum oSset for linkage of oSset rift basins (292 km)\nare not sensitive to the overpressure-factor used. Our calculations are also successful\nin predicting model behaviour so that we consider them useful. We think that this point\nis already discussed in the supplementary material. To point the reader more clearly to\nthe supplementary material, we amended the following sentence to the manuscript\n(addition in bold font):\nSee supplementary information for detailed information on the definition of the diHerent\nforces involved, including characteristic values and simplifying assumptions.\nConcerning deviations from lithostatic pressure: figures S1 to S4 all indicate deviatoric\nstresses of several hundreds of MPa up to 1 GPa according to the colobar. Since tectonic\noverpressure magnitudes are similar to deviatoric stress magnitudes, overpressure should be\nsignificant in the models. Therefore, if there are significant deviatoric stresses, then the\nassumption of lithostatic pressure is not justified. In the rebuttal, the authors say that tectonic\noverpressure stems from gravitational potential energy (GPE) variations. However, this\napplies only to situations without far-field deformation. In a lithosphere under\ncompression/shortening, there can be significant tectonic overpressure also without GPE\nvariations if deviatoric stresses are large, as is the case in the presented models. For the yield\nstrength envelope for F_strong, shown in figure 6, there is significant overpressure, for\nexample, in the brittle regions of the mantle lithosphere. This has been shown, for example,\nin the study of Petrini & Podladchikov, Journal of metamorphic Geology, 2000.\nWe thank the reviewer and fully agree with his assessment that there are significant\noverpressures. As described above, we assume lithostatic pressure and use a scaling\nfactor to approximate overpressure conditions. We think that this is a simple and useful\napproach, and discussed in the supplementary material. We agree that overpressure\nalso arises from lithosphere being under compression/shortening as for instance\ndescribed by Petrini & Podladchikov, 2000. Therefore, we wrote in the supplementary\nmaterial:\nThese assumptions neglect tectonic overpressure, e.g. related to gravitational potential\nenergy variability, which can be significant [Schmalholz et al., 2014]\nThe “e.g.” should indicate that there are other “sources” of overpressure, without going\ninto a longer discussion. We nevertheless think that it is useful to mention far-field\ndeformation eSects and that overpressure is observed in our models, and changed the\nsentence above (additions in bold font):\nThese assumptions neglect tectonic overpressure, e.g. related to gravitational potential\nenergy variability or far-field deformation, which can be significant [Petrini and\nPodladchikov, 2000; Schmalholz et al., 2014] and is also observed in our models\n(Figs. S1 - S4).\nResponse concerning line 439: “We note that this parametrisation is nonetheless commonly\nused throughout the community and reported in the literature.“. This is not a convincing\nreply, because the fact that others use a certain method does not necessarily mean that the\n3\nmethod is correct. I am also not sure which community is meant. The mesh-dependence can\naffect the calculated load-bearing capacity of a strain-softening material and hence can affect\nthe force calculation. But since the authors do anyway not use their 3D model results,\ninvolving strain softening, to calculate the forces, but a simplified 1D yield strength envelope,\ntheir force analysis is not affected by the mesh-dependence.\nWe thank the reviewer for this comment. As the reviewer writes, it is not aSecting the\nforce analysis. With the community we mean other research groups that also\ninvestigate upper mantle-scale deformation with non-linear rheologies and large strain.\nMany of them, using also other codes than the one used here, use strain-dependent\nweakening of the frictional parameters, which is the same approach.\nThe overall results of the force balance analysis make sense, but the accuracy and numerical\nrobustness of the calculated values of F_weak and F_strong (in units of N/m) and the\nresulting forces (in units N) is not clear to me.\nWe thank the reviewer for this comment. As discussed in the supplementary material,\nthe presented values are minimum estimates because of the assumption of lithostatic\npressure. We think that the additional supplementary model plots give a good overview\nof the involved deviatoric stresses.\nBest regards,\nStefan Schmalholz\nReferences:\nPetrini, K., and Y. Podladchikov (2000), Lithospheric pressure-depth relationship in\ncompressive regions of thickened crust, J Metamorph Geol, 18(1), 67-77.\nSchmalholz, S. M., S. Medvedev, S. M. Lechmann, and Y. Podladchikov (2014),\nRelationship between tectonic overpressure, deviatoric stress, driving force, isostasy\nand gravitational potential energy, Geophys J Int, 197(2), 680-696,\ndoi:10.1093/gji/ggu040.\n4"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have adequately addressed all of my comments and concerns in their revised manuscript. I find the current version suitable for publication and recommend acceptance in its present form.", "reviewer_2": "I have read through the revised manuscript and the reply of authors to reviewers. The authors have carefully addressed and discussed the comments and revisions suggested by reviewers. The paper is now much improved in clarity and represents an important and illuminating contribution to our understanding of rift inversion and orogenesis at a lithospheric scale.\nOne minor correction needs to be made on figure 1d - change Mauleon Basin to Pyrenean Basin as revised in the text.\nI congratulate the authors on this very nice work.\nI recommend publication\nMary Ford, université de Lorraine, Nancy, France.", "reviewer_3": "The authors have addressed my comments and clarified the text with respect to model assumptions. They also added new panels in the figures in the supplementary material showing the deviatoric stress which I find very useful. Only in case the authors will make another revision (which is not necessary), I have few minor comments to the manuscript. I also have two comments to their rebuttal.\nMinor comments\nLine 377: “…if their offset is not too large…“; maybe the authors could quantify what is “not too large“ with a dimensionless number, for example, if the ratio of offset to crustal thickness is less than approximately…\nLine 451: If H is only the radiogenic heat production rate, please specify this.\nLine 510: For completeness, the authors could also specify the value or criteria for the mechanical time step.\nComments to rebuttal and new text in supplementary material\nLithostatic pressure and force balance analysis: It is still not clear to me why the authors assume lithostatic pressure and a constant strain rate for the force calculation. Why don’t the authors not use the values of deviatoric stresses calculated in their 3D thermo-mechanical model to calculate F_int? Concerning deviations from lithostatic pressure: figures S1 to S4 all indicate deviatoric stresses of several hundreds of MPa up to 1 GPa according to the colobar. Since tectonic overpressure magnitudes are similar to deviatoric stress magnitudes, overpressure should be significant in the models. Therefore, if there are significant deviatoric stresses, then the assumption of lithostatic pressure is not justified. In the rebuttal, the authors say that tectonic overpressure stems from gravitational potential energy (GPE) variations. However, this applies only to situations without far-field deformation. In a lithosphere under compression/shortening, there can be significant tectonic overpressure also without GPE variations if deviatoric stresses are large, as is the case in the presented models. For the yield strength envelope for F_strong, shown in figure 6, there is significant overpressure, for example, in the brittle regions of the mantle lithosphere. This has been shown, for example, in the study of Petrini & Podladchikov, Journal of metamorphic Geology, 2000.\nResponse concerning line 439: “We note that this parametrisation is nonetheless commonly used throughout the community and reported in the literature.“. This is not a convincing reply, because the fact that others use a certain method does not necessarily mean that the method is correct. I am also not sure which community is meant. The mesh-dependence can affect the calculated load-bearing capacity of a strain-softening material and hence can affect the force calculation. But since the authors do anyway not use their 3D model results, involving strain softening, to calculate the forces, but a simplified 1D yield strength envelope, their force analysis is not affected by the mesh-dependence.\nThe overall results of the force balance analysis make sense, but the accuracy and numerical robustness of the calculated values of F_weak and F_strong (in units of N/m) and the resulting forces (in units N) is not clear to me.\nBest regards,\nStefan Schmalholz"}, "author_response": "Response to Reviewers for final revision\nSebastian G. Wolf, Ritske S. Huismans, Josep Anton Muñoz, Dave A. May\nThis document concerns our final revision of the paper entitled Rift linkage and inheritance\ndetermine collisional mountain belt evolution. We thank all three reviewers again for very\npositive, constructive, and useful comments that helped improve our manuscript and are\nhappy that the reviewers considered our revisions suitable. We accommodated the few\nremaining changes suggested by the reviewers. Below follows a point-by-point reply to the\nreviewer comments. Comments are in black font with serif; our answer is written in blue font\nwithout serif. Lighter blue text without serif highlights changes to the manuscript text.\nReviewer #1\nThe authors have adequately addressed all of my comments and concerns in their revised\nmanuscript. I find the current version suitable for publication and recommend acceptance in\nits present form.\nWe thank the reviewer for useful comments and are pleased to hear that the reviewer\nconsiders her/his comments and concerns adequately addressed.\nReviewer #2\nI have read through the revised manuscript and the reply of authors to reviewers. The authors\nhave carefully addressed and discussed the comments and revisions suggested by reviewers.\nThe paper is now much improved in clarity and represents an important and illuminating\ncontribution to our understanding of rift inversion and orogenesis at a lithospheric scale.\nOne minor correction needs to be made on figure 1d - change Mauleon Basin to Pyrenean\nBasin as revised in the text.\nI congratulate the authors on this very nice work.\nI recommend publication\nMary Ford, Université de Lorraine, Nancy, France.\nWe thank the reviewer, Mary Ford, for detailed reading and useful comments, and are\npleased to read that our revisions are considered suitable. We agree with the reviewer\nand changed the basin name on figure 1d to Pyrenean Basin as suggested.\nReviewer #3:\nThe authors have addressed my comments and clarified the text with respect to model\nassumptions. They also added new panels in the figures in the supplementary material\nshowing the deviatoric stress which I find very useful. Only in case the authors will make\n1\nanother revision (which is not necessary), I have few minor comments to the manuscript.\nI also have two comments to their rebuttal.\nWe thank the reviewer, Stefan Schmalholz for useful constructive feedback and are\npleased to read that he considers further revisions not necessary. We nonetheless\naccommodated the minor comments by the reviewer, see below.\nMinor comments\nLine 377: “…if their offset is not too large…“; maybe the authors could quantify what is “not\ntoo large“ with a dimensionless number, for example, if the ratio of offset to crustal thickness\nis less than approximately…\nWe thank the reviewer for this comment. We quantified “what is not too large” in section\n2.2 under Continental rifting: To link or not to link. There we show that linkage depends\non distance between the weak zones and their relative lithospheric strength. It’s thus\nnot only one factor that can be easily described by a non-dimensional number. We\nagree though that this sentence can be clarified and changed it to the following\nsentence:\n(i) extensional basins with closest proximity and lowest lithospheric strength are most\nlikely to link during rifting\nLine 451: If H is only the radiogenic heat production rate, please specify this.\nYes, it is. We followed the reviewer suggestion and added the word radiogenic to the\nrespective sentence.\nLine 510: For completeness, the authors could also specify the value or criteria for the\nmechanical time step.\nWe agree and added the following sub-sentence:\n…, which is dynamically computed based on a Courant–Friedrichs–Lewy condition and\ntypically in the order of 9000 yr.\nComments to rebuttal and new text in supplementary material\nLithostatic pressure and force balance analysis: It is still not clear to me why the authors\nassume lithostatic pressure and a constant strain rate for the force calculation. Why don’t the\nauthors not use the values of deviatoric stresses calculated in their 3D thermo-mechanical\nmodel to calculate F_int?\nWe thank the reviewer for this comment, which we split in two for simplicity. As written\nin the previous reply, we tried to compute F_weak and F_strong directly from the model.\nThis is straight forward for F_strong, as it can be computed with a 1D depth-integrated\nstrength profile. It is, however, not a simple 1D profile for F_weak, which would need to\nfollow a 3D path along a shear zone. Doing this numerically is possible, but comes with\nadditional complexity, e.g. which criteria to use in choosing the path along which to\nintegrate. Therefore, we opted to use yield-strength envelopes. We then tested whether\n2\nthis simple approach can reproduce the observed deviatoric stresses when we simply\nscale the lithostatic pressure by a factor, and find that our approach is a useful first-\norder approximation. Furthermore, as written in the previous reply, the resulting values\nthat we get for instance for the maximum oSset for linkage of oSset rift basins (292 km)\nare not sensitive to the overpressure-factor used. Our calculations are also successful\nin predicting model behaviour so that we consider them useful. We think that this point\nis already discussed in the supplementary material. To point the reader more clearly to\nthe supplementary material, we amended the following sentence to the manuscript\n(addition in bold font):\nSee supplementary information for detailed information on the definition of the diHerent\nforces involved, including characteristic values and simplifying assumptions.\nConcerning deviations from lithostatic pressure: figures S1 to S4 all indicate deviatoric\nstresses of several hundreds of MPa up to 1 GPa according to the colobar. Since tectonic\noverpressure magnitudes are similar to deviatoric stress magnitudes, overpressure should be\nsignificant in the models. Therefore, if there are significant deviatoric stresses, then the\nassumption of lithostatic pressure is not justified. In the rebuttal, the authors say that tectonic\noverpressure stems from gravitational potential energy (GPE) variations. However, this\napplies only to situations without far-field deformation. In a lithosphere under\ncompression/shortening, there can be significant tectonic overpressure also without GPE\nvariations if deviatoric stresses are large, as is the case in the presented models. For the yield\nstrength envelope for F_strong, shown in figure 6, there is significant overpressure, for\nexample, in the brittle regions of the mantle lithosphere. This has been shown, for example,\nin the study of Petrini & Podladchikov, Journal of metamorphic Geology, 2000.\nWe thank the reviewer and fully agree with his assessment that there are significant\noverpressures. As described above, we assume lithostatic pressure and use a scaling\nfactor to approximate overpressure conditions. We think that this is a simple and useful\napproach, and discussed in the supplementary material. We agree that overpressure\nalso arises from lithosphere being under compression/shortening as for instance\ndescribed by Petrini & Podladchikov, 2000. Therefore, we wrote in the supplementary\nmaterial:\nThese assumptions neglect tectonic overpressure, e.g. related to gravitational potential\nenergy variability, which can be significant [Schmalholz et al., 2014]\nThe “e.g.” should indicate that there are other “sources” of overpressure, without going\ninto a longer discussion. We nevertheless think that it is useful to mention far-field\ndeformation eSects and that overpressure is observed in our models, and changed the\nsentence above (additions in bold font):\nThese assumptions neglect tectonic overpressure, e.g. related to gravitational potential\nenergy variability or far-field deformation, which can be significant [Petrini and\nPodladchikov, 2000; Schmalholz et al., 2014] and is also observed in our models\n(Figs. S1 - S4).\nResponse concerning line 439: “We note that this parametrisation is nonetheless commonly\nused throughout the community and reported in the literature.“. This is not a convincing\nreply, because the fact that others use a certain method does not necessarily mean that the\n3\nmethod is correct. I am also not sure which community is meant. The mesh-dependence can\naffect the calculated load-bearing capacity of a strain-softening material and hence can affect\nthe force calculation. But since the authors do anyway not use their 3D model results,\ninvolving strain softening, to calculate the forces, but a simplified 1D yield strength envelope,\ntheir force analysis is not affected by the mesh-dependence.\nWe thank the reviewer for this comment. As the reviewer writes, it is not aSecting the\nforce analysis. With the community we mean other research groups that also\ninvestigate upper mantle-scale deformation with non-linear rheologies and large strain.\nMany of them, using also other codes than the one used here, use strain-dependent\nweakening of the frictional parameters, which is the same approach.\nThe overall results of the force balance analysis make sense, but the accuracy and numerical\nrobustness of the calculated values of F_weak and F_strong (in units of N/m) and the\nresulting forces (in units N) is not clear to me.\nWe thank the reviewer for this comment. As discussed in the supplementary material,\nthe presented values are minimum estimates because of the assumption of lithostatic\npressure. We think that the additional supplementary model plots give a good overview\nof the involved deviatoric stresses.\nBest regards,\nStefan Schmalholz\nReferences:\nPetrini, K., and Y. Podladchikov (2000), Lithospheric pressure-depth relationship in\ncompressive regions of thickened crust, J Metamorph Geol, 18(1), 67-77.\nSchmalholz, S. M., S. Medvedev, S. M. Lechmann, and Y. Podladchikov (2014),\nRelationship between tectonic overpressure, deviatoric stress, driving force, isostasy\nand gravitational potential energy, Geophys J Int, 197(2), 680-696,\ndoi:10.1093/gji/ggu040.\n4"}]

Decision letter