Anatomy of a post-subduction collision

Domain: biology
Published: 2026-03-26
Editorial outcome: major_revision

Abstract
The spatiotemporal interplay of long-lasting post-orogenic deformations make continental
collision zones among Earth's most enigmatic systems. Here, we employ three-dimensional
thermomechanical modeling to decode mantle dynamics of the Arabian-Eurasian collision—an
archetype of post-subduction tectonics. Our key findings reveal that plumelet-plate interactions
drive deformations both within and at the margins of convergent plates, forming modern
kinematics, surface tectonics, and plate boundary configurations. We document previously
unrecognized segmentation of the subducted Neotethyan slabs (Bitlis and Zagros) accompanied
by upper-plate tearing, which fundamentally modifies the seismotectonic stress accumulation
along the Arabian-Turkish-Iranian boundary. The convective support from the plumelet beneath
the former Tethyan magmatic arc drives drip-like lithospheric removal under the southern
Georgian highland, providing a regional-scale example of arc-to-intraplate deformation
transformation. Our results offer unified framework for understanding how upper mantle processes
control surface deformation in post-subduction systems dominated by plumelet dynamics.

Introduction
Plate tectonic collision zones can experience prolonged crustal and upper mantle
deformation that defines post-subduction tectonics in continental orogenic settings. Geodynamic
interpretations of these processes, such as slab break-off (e.g., Eastern Anatolia; Memiş et al.,
2020), mantle indentation (e.g., Alps; Schwartz et al., 2024), relamination (e.g., Zagros; Tesauro
et al., 2024) and dripping of continental and/or oceanic lithospheric mantle and mafic lower crust
(e.g., Central Anatolia; Andersen et al., 2024; Colorado Plateau; Levander et al., 2011; Andean
system; Curie et al., 2015; Göğüş et al., 2022; Beck & Zandt, 2002; Carrapa et al., 2025), polarity
reversal of plate consumption (Pysklywec, 2001) accompanied by crustal shortening and surface
uplift have been studied globally (van Hinsbergen et al., 2024 and references therein). However,
an important unresolved aspect of post-subduction tectonics lies in deciphering the multiscale
anatomy of upper mantle dynamics and its surface manifestations. This persistent knowledge gap
stems from the spatiotemporal overlap of collisional processes, which produce complex, often
ambiguous magmatic, tectonic, and geophysical fingerprints.

In this work, we focus on the Arabian-Eurasian collision region (Fig. 1) as a plate collision
archetype where a multitude of post-subduction mantle dynamics are actively shaping the surface
geology and tectonic configuration. Based on seismic imaging (Piromallo & Morelli 2003) and
our geodynamic model of the upper mantle structures, we resolve the three-dimensional (3D)
mantle circulation patterns with large- and small-scale convective cells, and their dynamic
interactions with overlying plates. Different from previous studies, including numerical
experiments in this part of the collision belt (Fig. 1b, e.g., Faccenna et al., 2014; Uluocak et al.,
2021; 2024 and references therein), our 3D thermomechanical model, extending east to the
Caspian Sea, is sensitive to both large (e.g., plumelet, Uluocak et al., 2024) and regional-scale
changes, which may correspond to specific zones with thinned/missing lithospheric mantle (e.g.,
the proposed delamination-modified regions of eastern Anatolia and western Greater Caucasus),
high plateaus (e.g., Turkish-Georgian-Armenian-TGA, the East Anatolian and the NW Iranian
plateaus), deep basins (e.g., Mesopotamian-Zagros, Kura, eastern Black and western Caspian seas
basins) and plate boundaries (e.g., Bitlis-Zagros suture and fold-thrust belt) (Fig. 1).

The Arabian-Eurasian collision zone involves a complex interplay of crustal tectonic
processes and lithospheric variations. Widespread Neogene-Quaternary volcano-sedimentary
rocks occur along the collisional front, young volcanic centers, and high plateaus (~2 km, TGA
and North Iranian Plateau) surrounded by suture and fold-thrust belts define major physiographic
units of the region (Fig. 1). Syn and post-collisional magmatic rocks with arc-type geochemical
compositions extend from the Eastern Pontide and the Lesser Caucasus arcs to the south, along the
Alborz magmatic belt (Adamia et al., 2017, Sugden et al., 2019, Fig. 1a). The northeastern Arabian
margin is dominated by multi-detachment fold-thrust systems (i.e., the Zagros-Sanandaj-Sirjan
basal fold and thrust belt), deforming a thick (~10–12 km) sedimentary cover and underlying crust
(Verges et al., 2011; Fig. 1b). North of the intramountain Cenozoic foreland basins of the
Transcaucasia (e.g., Terek, Rioni, Kartli and Kura basins), the Greater Caucasus orogen
accommodates the NW-SE striking Caucasus Main Thrust (e.g., Ismail-Zadeh et al., 2020 and
references therein). Within this regional context, a transform fault, the Aras (Araks) Fault (Fig.
1b), separates a region undergoing mostly strike-slip deformation to the west from an area
experiencing shortening to the east (Fig. 2), implying a segmentation of collisional tectonics
between the eastern Anatolian, eastern Greater Caucasus, and western Iranian regions (Rolland,
2017 and references therein). However, the active mantle root of this complex collisional system
remains poorly quantified.

Sub-crustal seismicity (≥ ~40 km depth) is predominantly recorded in the central and eastern
Greater Caucasus, Kura Basin, and south-central Caspian Sea basins (Mumladze et al., 2015;
Rolland 2017; Shumlyanskaya & Burmin, 2022, Fig. 1b). These events are commonly attributed
to the NE-dipping lithospheric root of the Kura Basin and/or the remnant subducting oceanic slab
beneath the eastern Greater Caucasus (Zor, 2008; Koulakov et al., 2012; Mumladze et al., 2015;
Gunnels et al., 2021; Shiva et al., 2025 and references therein). Fast seismic velocity anomalies in
the large and regional-scale seismic studies image the thick lithospheric root (~ 200 km depth) of
the Kura Basin and the eastern part of the Lesser Caucasus (e.g., Piromallo & Morelli 2003; Zor,
2008; Koulakov et al., 2012). Similarly, fast seismic velocities beneath the Terek (~ 4.7 km
sedimentary thickness) and Caspian basins reflect dense crustal and lithospheric roots underlying
thick sedimentary sequences (Adamia et al., 2011; Gunnels et al., 2021). Comparable high-velocity
anomalies at shallow depths (≤ ~ 150 km) are observed beneath the Black Sea, Mesopotamian-
Zagros foreland basins, and Zagros fold-thrust belt, consistent with high Pn seismic wave velocity
perturbations in these regions (Zor, 2008; Lu et al., 2017; Kounoudis et al., 2020; Uluocak et al.,
2021; 2024).

In contrast to high amplitudes of guided seismic waves (Pn and Sn) that travel in the
lithospheric mantle, crustal seismicity and slow seismic velocities correlate well with zones of high
Sn and Pn attenuation, indicating thinned or absent lithosphere (≤ 50–90 km) beneath the East
Anatolian, TGA, the North Iranian plateaus and the West Greater Caucasus (Koulakov et al., 2012;
Mumladze et al., 2015; Lu et al., 2017; Kaviani et al., 2022; Uluocak et al., 2021). Concerning
lithospheric mantle structures with numerous lines of geophysical and petrographic evidence,
including a large-scale low-velocity zone extending from Afar to the Greater Caucasus (≤ ~300
km depth, Fig. 1c, e.g., Civiero et al., 2023), the latest numerical model by Uluocak et al. (2024)
implies a lithospheric channel with SW-NE mantle flows (Ershov & Nikishin, 2004) –– a plumelet
(i.e., upper-mantle plume migration without significant tail and a mushroom head; Uluocak et al.,
2024)––here.

The lithospheric heterogeneity documented beneath the region (Uluocak et al., 2021) can be
attributed to Tethyan plate orogeny, as evidenced by numerous observations. The northern branch
of the slab subducted beneath Eurasia created arc and back-arc volcanism, mostly dating from the
Jurassic to Eocene (e.g., Şengör & Yılmaz, 1981; Adamia et al., 2011; Rolland et al. 2012; Rolland,
2017 and references therein). This part of the slab broke off around the Early Eocene, triggering
crustal shortening in the Greater and the Lesser Caucasus by the Oligocene (Avdeev & Niemi,
2011; Cowgill et al., 2016). Concurrently, the southern branch (Bitlis-Zagros slab) initiated
northward subduction, and the collisional front started to form in response to the Eocene-
Oligocene Zagros slab roll-back with the mantle window (slab gap) opening between the Zagros
and the Bitlis slabs (Rabiee et al., 2020 and references therein). Subduction jumped northward (~
Eocene, e.g., Rolland et al., 2012), and the closure of the oceanic basin persisted due to the ongoing
collision between the Anatolian-Southern Armenian Block and Pontides active margin (Rolland
et al., 2012; Shiva et al., 2025). This phase was alternatively characterized by the initiation of
south-dipping subduction of the Caucasus Basin along the arc (e.g., Rolland 2017). Terminal
closure along the Bitlis-Zagros suture (Late Oligocene–Middle Miocene, Fig. 2a) may have led to
slab detachment beneath eastern Anatolia at ~25-10 Ma (Keskin, 2007; Şengör et al., 2008;
Rabayrol et al., 2019 and references therein , Fig. 2b). Ongoing convergence reactivated Pliocene
uplift, folding in the Kura Basin with crustal shortening in the southern Greater Caucasus,
intensified exhumation across the Greater Caucasus, and strike slip deformation to the south (e.g.,
North and East Anatolian faults, Figs. 1a, 1b, & 2c) (Adamia et al. 2011; 2017; Avdeev & Niemi,
2011; Cowgill et al. 2016; Rolland, 2017; Ismail- Zadeh et al., 2020; Molin et al., 2023). Consistent
with numerous studies, such as river network analyses (Molin et al., 2023), the initial break-off of
the southern Neotethys slab beneath the Bitlis-Zagros suture zone could have triggered mantle
upwelling, propagating from SW to NE beneath the Arabian Plate - a plumelet (Uluocak et al.,
2024, Fig. 2d). Among various slab deformation models, removal mechanisms involving
spatiotemporal overlapping processes (Uluocak et al., 2024) can be primarily characterized by: i)
wholesale or progressive Bitlis-Zagros slab delamination along the Arabian collisional front in the
south (Keskin 2007; Memiş et al., 2020 and references therein); and ii) delamination beneath the
Transcaucasus (e.g., Shiva et al., 2025 and references therein) accompanied by northward plate
tearing toward the Greater Caucasus (e.g., Mumladze et al., 2015). Although 3D heterogeneities
beneath the region are still debated, magmatic records provide evidence of mantle-crust
interactions as a consequence of a discontinuous/semi-episodic closure of the Neotethys Ocean
(e.g., Rolland et al., 2012).

Post-collisional volcanism initiated (around in the early–mid-Miocene) on the Erzurum-Kars
Plateau (locations shown by T in Fig. 1d) and migrated south and eastward, towards the Bitlis-
Zagros Suture and the Armenian border, respectively, consistent with radiometric ages (e.g., Lin
et al., 2020; age compilation in Memis et al., 2020). The late Miocene-Quaternary calc-alkaline to
alkaline magmatism is documented across the collisional front with ages decreasing westward and
subduction signatures weakening southward, towards the Bitlis-Zagros suture (Keskin, 2007 and
references therein). Based on chemical and petrological analyses of widespread alkaline lavas
active since the Late Cenozoic in the Turkish-Iranian Plateau, it has been argued that the volcanism
may be linked to a low seismic velocity zone/body extending to depths of ~300 km
(Soltanmohammadi et al., 2018). The western and central Greater Caucasus host extensive
Neogene–Quaternary magmatism (e.g., Elbrus, Kazbek; Adamia et al., 2017; Forte et al., 2022),
contrasting with the eastern segment, and highlighting potential along-strike variations in mantle
dynamics. Quaternary volcanism in the Lesser Caucasus (e.g., including Kazbek and Georgian
volcanic rocks with varying degrees of subduction components, Bewick et al., 2020) show
compositional variability and distinct melting conditions, indicating small-scale convective
upwellings and lateral contrast in lithospheric depths (Neill et al., 2015; Sugden et al., 2019). Yet,
the regional-scale interactions of long-wavelength asthenospheric processes—including potential
plumelet activity (Uluocak et al., 2024, Fig. 2d)—remain poorly constrained in such settings.
Overall, the Arabia-Eurasia collision domain exhibits a complex interplay of multistage
tectonics with diverse volcanism, plateau uplift and active fault systems - all linked to Neotethyan
slab dynamics involving advance, retreat, and eventual break-off (Adamia et al., 2011; Avdeev &
Niemi, 2011; Cowgill et al., 2016; Rabiee et al., 2020; Forte et al., 2022; Cavazza et al., 2024)
with growing evidence suggesting key contributions from mantle convective flows (Uluocak et
al., 2021; 2024).

Here, our comprehensive geology-geophysics-geodynamic analyses may help to illuminate
the tectonic anatomy of this complex post-subduction system. By coupling mantle flow models
with thermal and isostatic constraints, we reveal not only a SW-NE oriented plumelet from the
Arabian foreland injecting into the subduction system (Uluocak et al., 2024), but also 3D
thermomechanical erosion of the upper plates and the active Neotethys slab deformation along the
southern collisional front. Integrated data further expose distinct west-to-east differentiation in
deformation patterns across the Greater Caucasus, which modify tectonics at the northern margin
of the collision. Our results establish constraints on drip-like removal beneath the post-orogenic
magmatic arc under convective forces. Compared to the well-studied seismotectonics beneath the
Greater Caucasus Range, the poorly defined tectonic implications of crustal and subcrustal
earthquakes along the Arabian collision front between the Arabian-Anatolian-Iranian plate
boundary are tried to clarify in this study. Crucially, our data-model integration uncovers a
previously unrecognized complexity in post-subduction tectonics, providing a transformative
framework for understanding active continental collisions worldwide.

Results
We determine the regional mantle circulation and lithosphere deformation with a 3D
instantaneous thermomechanical model based on seismically defined upper mantle structures (Fig.
3) by following the conventional modeling approach discussed in detailed by Uluocak et al. (2021;
2024) in the Arabian-Eurasian collisional region. The results define the anatomy of the structures
and dynamics that interact in a complex way in this active post-subduction system of the Arabian-
Eurasian collision zone.

Mantle dynamics in the region are dominated by a focused rapid SW-NE directed inflow
from the Arabian foreland to the eastern Greater Caucasus—a plumelet— penetrating through the
complex morphology of the post-collision environment (Fig. 3). In the southeast, the upwelling
mantle beneath the Urumieh-Dokhtar magmatic belt in the Iranian Plateau migrates along-strike
of the Zagros mountains, reaching the stable Scythian Plate. The surface response to convective
flows (i.e., dynamic topography, e.g., Flament et al. 2013 and references therein) is dominated by
a positive dynamic support on topography (up to ~800 m) across the Turkish-Georgian-Armenian
volcanic zone, East Anatolian Plateau, including Lake Van rift province, North Iranian Plateau
and the western Greater Caucasus (Fig. 4a). Negative dynamic topography (down to -1.6 km)
occurs in basins, e.g., the Kura, Terek, western Caspian and eastern Black Sea. In contrast to
highlands (e.g., the East Anatolian and the Armenian plateaus), basins with underlying cold mantle
lithosphere (e.g., blue zones with ≤ 1000 K; Fig. 3b) and low upper mantle seismic wave
attenuations (Kaviani et al., 2022) show a similar high shear stress pattern (>10 MPa) in modeling
results (Fig. 4b). The 3D flow model reveals lithospheric-scale zones of prolonged high strain rates
(e.g., > 2x10-15 1/s at the depth of 100 km; Fig. 4c), delineating a SW–NE trending deformation
belt along the northern Lake Van intraplate rift zone and the SE-NW oriented Sanandaj-Sirjan and
Zagros fold-thrust belts (Fig. 3c). High strain rate anomalies also emerge beneath the North
Arabian lithosphere, particularly surrounding the deep NW Zagros, western Caspian and Kura
basins (Fig. 4c).

The flow modelling indicates regional deformation zones in the central and eastern Greater
Caucasus extending to the stable Scythian Plate, marked by visible lateral contrast of low shear
stresses in the crust (≤ 1 MPa; Fig. 4b) with high strain rates beneath the crust (Fig. 4c). The low
shear stress and strain rates with the positive dynamic topography manifest in the TGA volcanic
highlands and the West Greater Caucasus is in concordance with hot upper mantle anomalies (≥
1500 K; Fig. 4d).

To elucidate mantle controls on crustal dynamics, we integrate seismically-derived thermal
state with crustal isostasy analysis using various independent datasets (Fig. 5). We show P- and S-
wave velocity ratios, which can give insights into melting conditions in the case of upwelling
mantle at 100 km depth and for the crust of the region (e.g., Zhu et al., 2015; Zhu, 2018). The
highest subcrustal velocity ratios (≥1.9) occur along magmatic arcs across the region. Slightly
lower values (~1.8 -1.9; Fig. 5a) are observed in the North Arabian Platform and the East Anatolian
Plateau. In terms of crustal velocity ratios (green contour; Fig. 5a), compared to global crustal
averages of ≥1.77 (Lin et al., 2020), the highest average values are spatially correlated with the
highest mantle velocity ratios in magmatic arcs and the West Greater Caucasus with an interesting
exception of the Central Greater Caucasus region. This result demonstrates a direct mantle-crust
linkage in concordance with crustal and subcrustal thermal anomalies beneath young volcanic
centers (Fig. 5a) and our 3D temperature model (Figs. 3b & 4d).

We also provide a comparison of regional velocity ratios, crustal thickness, and topography,
along with overall data from other areas (Fig. 5b, Gvirtzman et al., 2016). Although there are local
differences arising from different data sets, the changes in the South Caspian, Zagros, and Black
Sea basins and dynamically uplifted continents can be identified by cold (i.e., yellow zones with
lower elevation) and hot features (i.e., red zones with ~ 40-45 km crustal thickness), respectively,
in Fig. 5b. Vp/Vs values in the Armenian Plateau, for instance, fall in the zone (A; Fig. 5b) that is
classified by a crustal thickness of ≤ ~36 km (Lin et al., 2020), yet has a high elevation (≥ 1.5 km).
These results agree well with positive dynamic topography, estimated by our thermomechanical
model (Figs. 4a & 6a).

Consistent with the hot crustal thermal regime (Fig. 5), regional mantle structures are defined
by relatively shallow, arc-parallel relatively cold anomalies (< 1500 K) beneath the Transcaucasian
zone (including the eastern Pontides to the Lesser Caucasus, Fig. 4d). Our temperature model
reveals a distinct Pontide-cold anomaly buoyant at shallow mantle depths (~80-250 km; Figs. 4d
& 6b) that exhibits two key characteristics: i) it is decoupled from the overlying crust beneath the
volcanic plateau (Dzhavakheti Province of South Georgia); ii) it is spatially correlated with the
surface expression of young volcanoes and the adjacent Armenian magmatic highland (Figs. 3a &
5a).

In the active Arabian collisional front, the results reveal an interaction between the plumelet
and cold upper mantle structures, namely a mantle anomaly beneath the thick northern Arabian
lithosphere (~250 km depths; Fig. 6b). This anomaly differs from the NW Zagros depression area
(P2; Figs. 6d-e) in that it is decoupled from the upper plate, implying fully detached cold material-
termed as a part of a detached slab or delaminated Bitlis fragments- sinking beneath the North
Arabian Platform. Deep basins, such as Zagros and Kura basins in Fig. 6d and 6e, are characterized
by their cold roots and negative dynamic support on surface. We define an active slab breaking off
beneath the East Greater Caucasus with a cold thick root of stable Scythian Plate at the northern
boundary of the collision (P2 in Fig. 6e).

The temperature model (Fig. 6e) depicts the southeast subducted, actively tearing Zagros
Slab along the southern collisional margin (Fig 7c). We define the fully and partially detached
Bitlis slab in detail, along a cross-section in Fig. 6 and as a 3D temperature model in Fig. 7.
Notably, our findings provide compelling evidence for a plumelet-induced thermomechanical
erosion with the 3D mantle thermal structure (Fig. 7b) revealing a distinctive spatial signature
within the detached/sinking Bitlis Slab. Results demonstrate the regional thermomechanical
erosion beneath the crust (~40 km) of the Arabian foreland in the form of a plumelet-fed low-
viscosity zone (red star in Figs. 7).

Methods
We employed state-of-the-art numerical modeling using a well-established framework
applied to diverse tectonic settings (e.g., Central Anatolia, SE Carpathians; Uluocak et al., 2024
and references therein). We used seismic data (Piromallo & Morelli, 2003) to constrain the
geometry of lithospheric structures and mantle features for the model (Uluocak et al., 2024:
Uluocak, 2024a). Since estimations of instantaneous thermomechanical models are data-
dependent, sensitivity analyses for the modeling approach used here were performed and
documented in previous work on the collisional zone (Uluocak et al., 2021, 2024). Our model
leverages ASPECT (Bangerth et al., 2022 and references therein) to solve the governing equations
of the conservation of mass:
𝛻.𝐯 = 0 (1)
where 𝐯 is the velocity (m/s), momentum:
∇.σ +ρg = 0 (2)
ij
𝜎 is the stress tensor; g is gravitational acceleration (9.8 m/s2), ρ is density (kg/m3) and energy:
𝑖𝑗
ρC ( ∂T k∇2T+ρH (3)
p
∂t+ 𝐯.∇T)=
C is specific heat capacity (J/kg/K), T is the temperature (K), t is time(s), k is thermal conductivity
p
(W/m/K), H is the rate of internal heat production per unit mass (W/m3) incorporating the visco-
plastic rheology for an incompressible medium (Bangerth et al., 2022 and references therein). The
effective viscosity (η ) in the viscous stress is defined for the power law creep based on the
eff
temperature and the square root of the second invariant of the strain rate tensor (ε̇);
1 −1 (1−n) RQ+PV
η eff (ε̇,T) = βBnε̇ n exp ( ) (4)
2 PnRT
where n is the non-Newtonian viscosity exponent, Q is the activation energy, V is activation
volume, R is the universal gas constant (8.31 J/mol/K), B is the viscosity pre-factor and β is the
scaling factor, P is the total pressure (Bangerth et al., 2022, Table 1). The deviatoric stress; σ′ =
ij
min {σ ;σ } in relation to the deformation of the material was defined based on the plastic yield
y v
stress (σ = Psinφ+Ccosφ, i.e., Drucker-Prager yields criterion for the brittle failure, where φ
y
is the internal angle of friction, C is cohesion) and the shear stress tensor (σ = 2η ε̇). The
v eff
temperature model (e.g., Figs. 3b, 4d & 5e) was obtained by using seismic data with a density
scaling approach based on the thermal-expansion expression; ρ(T)= ρ (1- α ΔT) (reference
0
densities, ρ , are listed in Table 1; a coefficient of thermal expansivity, α, 2x10-5 1/K was used)
0
with the sampling intervals of 20 km in each direction for the investigated region (Uluocak et al.,
2024). The time stepping in model was carried out using a small Courant–Friedrichs–Lewy
number of 0.1, yielding a timestep of up to ~1x105 years and the relatively small‐velocity errors
(i.e., RMS ≤ 2 cm/year) with the effective computation time (Uluocak et al., 2021). Dynamic
support in the topography- dynamic topography (e.g., Faccenna et al., 2014)- is defined along the
free surface based on the vertical stress exerted on the lithosphere as a response to normal fluid
stresses in the underlying mantle. Any possible responses of elasticity, mass transfer, and
radioactive heating were not considered in the incompressible solution domain. Since different
numerical models with various temperature configurations and rheological features were tested
(Uluocak et al., 2021; 2024), in this study we followed a convection modeling approach with inputs
that have already been confirmed/validated by considering large sets of geophysical data and
analyses documented in the Arabian-Eurasian collisional region. For a more detailed discussion of
the modeling methodology and parameters, the reader should refer to Uluocak et al., (2024);
Uluocak (2024b) and Bangerth et al. (2022).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I read this manuscript with great interest as it presents a novel perspective on the Arabia-Eurasia collision zone, drawing upon a synthesis of geological and geophysical data, as well as thermo-mechanical modelling. This new perspective is supported by the data, and I recommend that it be considered for publication in Nature Communications. However, I believe that some points still need to be clarified before it can be accepted.\nFirstly, it seems that the authors are pitting their model against a 'traditional model', but the latter is not clearly explained in the present paper. In my view, there are several theories regarding the evolution of the Caucasus region. Some previous papers suggest a significant role for delamination (cf. Grosjean et al., 2022), a process which is not clearly discussed or taken into account in this paper. Some papers have already proposed the contribution of the deep asthenosphere to the magmatism of the area (e.g. Soltanmohammadi et al., 2018), thus excluding any purely 'collisional' component. Previous models did not strictly fall into a single category, and I agree that the horizontal mantle flow or 'plumelet' model proposed here is novel, as previous models assumed that the mantle flowed mainly vertically.\nOn the other hand, in the 'Results' section, the authors immediately propose a model without first presenting the results in a cautious or analytical manner. For example, at lines 197–199, they claim that 'Mantle dynamics are dominated by a focused, rapid, SW–NE-directed inflow from the Arabian foreland to the eastern Greater Caucasus—a plumelet—penetrating through the complex morphology of the post-collision environment'. While I believe this interpretation is entirely plausible, it is crucial to present results more factually and demonstrate how they support the plumelet model.\nFurther points should be clarified:\nWhat role do vertical structures on the lithosphere scale play in propagating the plume? The inception of horizontal flow within a structure featuring vertical slabs requires zones of weakness inherited from the subduction-early collision history (step faults etc).\n- The so-called 'Bitlis slab' may also be the result of delamination, because the Bitlis subduction zone is characterised by HP metamorphism at 74–71 Ma (e.g. Rolland et al., 2012), which likely resulted in sinking at very deep levels of the mantle (deeper than 300 km) and should have been N-dipping (no volcanic arc developed towards the south). This tectonic zone is generally considered to be complex and composed of numerous tectonic blocks and sutures. Therefore, the 'closure of Neotethys' should not be considered a progressive process, but rather a discontinuous 'step-by-step' process in which basins were successively closed. This explains why the timing of 'collision' is disputed in this region.\nOverall, the text is well written, but some sections could be improved (I have underlined some examples in the text). Please avoid self-satisfaction (e.g. 'Our comprehensive geological, geophysical and geodynamic analyses illuminate the tectonic anatomy...').\nThe figures are well drafted and informative. However, when dealing with seismicity (Fig. 6), it is difficult to understand the clustering of seismicity. It might be better to present a representative vertical section showing the seismicity in relation to the crustal levels.\nI also think the title does not fully reflect the content, as it should emphasise the 'plumelet' hypothesis.\nI hope these comments are useful. With best regards.\nrefs\nGrosjean, M., Moritz, R., Rezeau, H., Hovakimyan, S., Ulianov, A., Chiaradia, M., & Melkonyan, R. (2022). Arabia-Eurasia convergence and collision control on Cenozoic juvenile K-rich magmatism in the South Armenian block, Lesser Caucasus.\nEarth-Science Reviews, 226, 103949.\nSoltanmohammadi, A., Grégoire, M., Rabinowicz, M., Gerbault, M., Ceuleneer, G., Rahgoshay, M., ... & Benoit, M. (2018). Transport of volatile-rich melt from the mantle transition zone via compaction pockets: implications for mantle metasomatism and the origin of alkaline lavas in the Turkish–Iranian plateau. Journal of Petrology, 59(12), 2273-2310.\nRolland, Y., Perincek, D., Kaymakci, N., Sosson, M., Barrier, E., & Avagyan, A. (2012). Evidence for∼ 80–75 Ma subduction jump during Anatolide–Tauride–Armenian block accretion and∼ 48 Ma Arabia–Eurasia collision in lesser caucasus–East Anatolia. Journal of Geodynamics, 56, 76-85.", "reviewer_2": "# Anatomy of a Post-Subduction Collision\n**Authors:** Ebru Şengül Uluocak, Russell N. Pysklywec, Claudio Faccenna & Taylor Schildgen\nThis paper addresses the problem of understanding mantle dynamics in continental collision zones where the driving oceanic slab has detached or broken off from the overriding plate. These post-collisional systems are characterized by prolonged crustal and upper mantle deformation that can persist long after the initial convergence and slab consumption has ceased. The Arabian-Eurasian collision zone exemplifies such complexity, where the Neotethyan oceanic slabs have undergone various stages of subduction, break-off, and foundering, leaving behind a fragmented lithospheric architecture.\nUnderstanding what drives continued deformation in these post-subduction environments remains a fundamental challenge in tectonics, as traditional plate-driving forces are absent yet significant surface deformation, volcanism, and seismicity persist.\n## Overview\nThe manuscript is well-written and organised and it is of significance because it showcases how mantle dynamics strongly influences strain rate, stress, topography, and mantle flow of post-collisional settings. It also demonstrates how tectonic activity post–dates the interruption of subduction due to the effect of mantle dynamics and provides a framework for the integration of geodynamic modelling with geophysical observables.\nThe approach is based on correlating various geophysical observables - seismic velocities, topography, crustal thickness, and seismicity patterns - with their 3D thermomechanical model results. The authors systematically demonstrate how regions of high Vp/Vs ratios correspond to magmatic arcs and elevated topography, how cold mantle anomalies correlate with thick lithospheric roots and deep basins, and how modeled strain rates align with observed seismicity distributions. It is these (positive) correlations that form the empirical foundation for validating their model predictions and serve as the primary evidence supporting their interpretations.\nAnother central concept is the \"plumelet\" which is defined by the authors as the migration of upper-mantle plume material into the region. It may be drawn in by flow in the region or by its own buoyancy. The plumelet is proposed as a flow of hot mantle material from the Arabian region to the eastern Greater Caucasus. They propose that, from their 3D thermomechanical modeling, this lateral mantle flow *drives* regional deformation, influences slab break-off processes, *causes* thermomechanical erosion of overlying plates, and creates the observed dynamic topography patterns across the Arabian-Eurasian collision zone. There are a number of other areas where subduction influences mantle flow - for example the paper by Miller & Becker (https://doi.org/10.1130/G34959.1) argues that uplift of the Atlas is enhanced by flow of plume material around the foundering Alboran slab. Or indeed the paper by Lanari et al (2023) which has some author overlap with this paper. Are the other notable examples that would help to establish this plumelet hypothesis and generalise it ?\nA fundamental question raised by this interpretation is whether the plumelets **are** actively driving the observed tectonic signals, or whether they represent more passive structures that are being drawn into the region by the process of slab breaking off and foundering. The distinction between active versus passive roles has significant implications for understanding the causal relationships in post-collisional tectonics and the broader applicability of the plumelet concept to other collision zones. Is the plumelet genuinely causing these effects ? Can it be determined convincingly through the models ?\nThe Tectonic context and motivation section is well-written, but is hard to follow due to the absence of a figure that provides context for the tectonic evolution of the study area. In Figure 1, the main magmatic units are shown; however, it would be beneficial to provide information about the composition of these volcanic units. Therefore, a figure showing the tectonic evolution, the location of magmatic units and their chemistry is needed. This is very important for readers who are not familiar with the geology of the study area.\nIn general, figures are very pretty but not always easy to follow. Figure 5 is quite straightforward, but Figure 6a is not. Figure 6b is difficult to see compared to figure 6c. We know the 3D renderings are a lot of work to produce but in static form they are less useful than the cross sections (sorry).\nFigure 6 is helpful because it integrates all the observables and the results from the numerical modelling. However, the mantle structure of the study area is quite complex and the 3D geometry of the plumelet propagation channel is not fully clear. We suggest you could draw your interpretation of this geometry in Figure 6b, and this will help the reader to better visualize how it correlates with variations in stress, strain rate, dynamic topography and Vp/Vs ratios.\n## Detailed Feedback\n- Pn and Sn are quite well known but not universally, it would be useful to add a short explanation\n- Lines 86-89: In the captions of Figure 1d, are the red triangles showing the magmatic units or active volcanoes?\n- LInes 149-150: It is worth mentioning which observations this refer to (geological, geophysical, geodynamic) ?\n- LInes 204-205: Figure 2 would benefit from adding a label with the location of the Scythian plate.\n- Lines 221-222: differences of less than 1MPa do not seem that “high”\n- Line 228: Is there any reason for Figure 2 to display the north at an orientation different from all other figures?. Although you include the labels of the main tectonic domains, it is confusing.\n- Lines 249-251: you mentioned “crustal velocity ratios”; however, you didn’t mention the crustal depth or depths considered in the contours.\n- Lines 253-254: “corroborates” is not the most suitable word for your explanation\n- Lines 283-284: What does “Pontide cold” mean?\n- Lines 333-336: The vectors in Figure 6 can be improved. This would facilitate the comparison with the observed GPS velocity field. Also, Figure 6 shows a zoomed region, but the extent of this zoomed region is not clear. It would be useful to include an inset with a regional map showing the location of this region.\n- Lines 274-275: What is the purpose of the vertical line at ~35 km of crustal thickness?. Is it related to the “A” zone?. This has to be clarified in the captions of Fig. 4b.\n- Lines 296-298: Interpretation should go into the discussion\n- Lines 309-310: Does “shaded” refer to the areas filled with red?. This has to be clarified.\n- Line 369 - 'indicating upper mantle-driven isostatic instabilities'. This is quite a strange phrase because isostasy implies a static, buoyancy equilibrium whereas this paragraph refers to mantle-driving of topography (which implies non-isostatic support)\n- Lines 403-406 Can you elaborate on the observed changes in mantle composition?. You mentioned this in the “Tectonic context and motivation” section, but this is an opportunity to link the specific changes to your modelling results.\n- Lines 416-420 Shortening this sentence would improve clarity\n- Line 469 - advection / diffusion equation is scrambled. It would be helpful to just write these out with equation numbers as it is really difficult to read these equations in-line.\n## Recommendation\nThis paper covers a novel and interesting topic and is potentially publishable in Nature Communications. The writing is currently somewhat breathless and that seems to get in the way of readability. Similarly, the figures are elegant but they are cluttered in a way that gets in the way of understanding. We are happy to recommend publication upon addressing comments", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "Review of “Anatomy of a post-subduction collision” by Uluocak et al.\nThe manuscript by Uluocak and co-authors presents an interesting study on the relationship between mantle thermal structure and various observables from the Arabian–Eurasian collision zone. Overall, the paper is well written and addresses a topic of broad relevance. The potential implications are significant. However, I believe the current modeling approach requires substantial strengthening before the conclusions can be fully supported.\nMajor comments\n1) Modeling strategy\nTo my understanding, only a single simulation is presented. This limits the robustness of the conclusions. Sensitivity tests exploring variations in the initial thermal structure (inferred from geophysical constraints) are essential. Similarly, the role of mantle rheology, crustal architecture, inherited sutures, and other weak zones should either be explicitly investigated or at least discussed in detail.\n2) Model setup and methodology\nThe modeling approach is not sufficiently detailed in its current form. Model dimensions and resolution should be specified.\nEven in instantaneous flow models, runs should be evolved long enough to achieve isostatic equilibration. What timescales were used?\nInitial conditions must be clearly introduced: what is the starting thermal and density structure?\nA table summarizing all material properties should be provided, ideally in the appendix.\n3) Figures and presentation\nFigure 2: The contours are unclear. What variable is being contoured? This needs clarification.\nFigure 3A: The method used to calculate dynamic topography must be described.\nA direct comparison between model outputs and observations is missing. For example, modeled surface velocities should be plotted against GPS data in a single figure. Such comparisons are crucial to convincingly demonstrate the model’s validity.\n4) Data availability.\nIf I recall correctly, Nature Communications requires that the data necessary to reproduce the models and generate the figures be made publicly available, for example via Zenodo.\nThe manuscript presents a substantial body of work and demonstrates the authors’ strong understanding of the Arabian–Eurasian collision zone. However, the modeling component requires further work and sensitivity testing to ensure the robustness of the conclusions. Addressing these points would significantly strengthen the study and enhance its impact.\nLine per line comments:\nL22. I am not quite sure to know/understand what is a “plumelet” at that stage (I am not on top of this litterature). As this is likely a key concept of the paper you may want to clarify this for a general audience, right from the start. Is that a small plume? Is it super-adiabatic? Is its buoyancy chemically controlled?\nL54. “and the geodynamic model of the upper mantle and lithosphere structures”. Is that based on previous work and reconstructions? You may need a reference here.\n§ L51-L65. It is not clear how the model is built from geophysical and geodynamic data at that stage.\nL68-L69. “By coupling mantle flow models with thermal 68 and isostatic constraints” this partially answer my previous comment. Maybe this should be moved in previous section.\nL73-L74 . “Our results for the first time establish new constraints”, remove “for the first time “, as you have “new constraints” in the same sentence, and also because you use “for the first time” L78.\nFig. 1c. Shouldn’t the caption be delta_Vs (or delta S-wave and not just S-wave) as you display % variation? I am not a geophysicist, but S-wave should be in km/s. If you have % it likely means that this is a variation from a reference model.\nL147-148. “The northern branch of the slab subducted beneath Eurasia by creating arc and back-arc volcanism”. Something is wrong in this sentence. I suppose you meant “by creating” -> “created”?\nL165. I don’t see the Erzurum-Kars plateau “T” on figure 1D.\nL192-L195. This should go in method and be expanded (and partly moved to appendix too). How did you integrate the data? Did you invert solely for temperature and you that to define the thermal structure of the model?\nWhat are the densities of your material? How do you choose them? Are they computed from phase equilibrium modelling? How sensitive is you modeling results when propagating error of the inversion? Are you only imposing the initial temperature structure? What are the assumptions in terms of rheologies and material properties for the different crustal and lithospheric units?\nDo you define fault, shear zones and weak zones? How are they distributed? Although it does not look like you impose such structures, how important can they be in controlling the stress-strain state of the region? I would suspect that they are rather important.\nYou also need to provide a table with all used material parameters.\nMoreover, I understand you perform instantaneous flow models. Do you account for elasticity? How is isostatic equilibrium achieved? How many timesteps are performed to reach elastic and isostatic equilibrium?\nWhat is your model size and resolution?\nL197-L199. At that stage I am still not quite sure to understand what is a plumelet. It seems to me that this is simply a buoyancy contrasts driving upward motion. An upwell?\nFig. 2. What are the iso-contours you are showing on the figure?\nFig. 3a. How do you compute the dynamic topography? This is quite critical to explain how this is computed in your model as this is not trivial especially for non-modelers.\nFigs 3b,c Shear stress and Strain rate. Are those invariants? Second invariant the shear stress and strain rates? Otherwise they should be tensors…\nFig. 3d. 3D temperature model with a slice. Are you plotting an iso-surface of the temperature? + a slice? Moreover, there are no yellow dashed lines in Fig. 3d. It needs be to moved in caption 3c.\nFig. 5. Is that realistic to have temperature reaching 2000°C at 300 km depth? Also it seems that large regions of very low temperature are predicted at great depth ~1000°C. Without contouring I have trouble to exactly see the temperature structure, but it seems rather extreme.\nL279-L284. Where is the cold anomaly on the figure 3d. I don’t see it. Actually, the figure shows really well the trench parallel hot anomaly due to slab break off.\nL288-L290. The “Sinking Britis Slab” of figure 5b?\nL298. The temperature model should be presented way earlier in my opinion as this is the basis of your modeling! Variation in the temperature model should also be explored to test how sensitive the modelling outcome is.\nL301-L305. First, this slab sinking and related upwell and potential erosion are very close to what seems the boundary of the model. How is this constraining the results? Or is it because of how the figure is clipped?\nL328-L330. This is rather difficult to see both modelled velocities and GPS velocities are not super-imposed on a figure. Having a figure showing how well the model and the data fit would really help to convince the reader.\nL337-L340. This is also difficult to visualize this based on figure 6. Would be better to have a figure super-imposing seismicity and modeled second invariant of the strain rate for instance.\nFig. 6c Given that you are using non-linear T dependent rheology, the discretized viscosity colorbar is an odd choice\nL349-L355. How can the southeastern propagation of the deformation be discussed in the light of an instantaneous flow model? I may be misunderstanding here, but from my point of view this can only be explored when performing forward models.\nL390-L392. It seems like a bit of a circular reasoning here. As the small scale anomaly are the outcome of using the tomography model as input for the model.\nL435- Where is the plumelet coming from? What drives its ascent?"}, "author_response": "Response Letter to Editor\n---------------------------------\nWe sincerely thank you and Reviewers for their insightful reviews of the manuscript. We agree\nthat these additional clarifications are helpful in improving the explanation of our approach. A\nnumber of changes have been done in the revised manuscript. In particular, we have made the\nfollowing primary modifications:\n1) The Introduction has been revised to include a new Figure 2 that more clearly illustrates\nthe tectonic evolution of the region as suggested by Reviewers 1 and 2. Corresponding\nreferences for this conceptual model have been added in the revised manuscript.\n2) We have modified figures (revised Fig 1d, Fig 3, Fig. 7) to show variations, symbols and\nresults more clearly as suggested by the Reviewers. Related captions have been\ncorrected/modified to explain figures and our findings clearly.\n3) In response to Reviewer comments, we have added new references to the manuscript. To\naccommodate these within the reference limit, we removed one reference on local-scale\nchanges (Oyan et al., 2023) and some modeling references (Uluocak et al., 2016; 2019),\nwhose context is covered by other cited works.\n4) As suggested by Reviewer 4, the Method section has been expanded with additional\ninformation, references and Table to clarify material parameters, modeling time and space\nresolutions. We think that the revised manuscript will make this modeling strategy more\nexplicit for a broader audience.\n5) We clarified Data Availability section in the response letter. The full suite of input\nparameters is detailed in the following link (https://doi.org/10.5281/zenodo.12568287) that\nhas already been shared publicly by Uluocak, et al. (2024, e2024GC011639) and Uluocak\n(2024, https://doi.org/10.5281/zenodo.12568287). All codes used in the numerical\ncalculations are also available in open-sourced ASPECT directory\n(https://aspect.geodynamics.org/). Seismic data used in modeling is also available in\nPiromallo & Morelli (2003), as stated in the text (e.g., Introduction, Method). Further, we\nstate that ESU, the corresponding author, will gladly assist anyone who needs help\ncompiling these open sources. The core objective of this study is to utilize a\nproven/confirmed numerical modeling methodology to investigate complex processes in\ncollisional settings by considering standard space and reference limitations of the journal.\nWe think this strategy can allow us to focus on elucidating deformation mechanisms and to\neffectively communicate these new findings to not only modelers but also a broad\nreadership of Nature Communications.\n6) We thank Reviewers for also bringing our attention to minor corrections, we have addressed\nthese and typos in the text to ensure they are fixed in the revised manuscript as shown in\nthe tracked-changes file. We believe the revised manuscript presents a much clearer and\nmore cogent contribution representing our original new research.\nSincerely yours.\nOn behalf of the Authors\nDr. Ebru Şengül Uluocak (sengul@gfz.de, ebrusengul@comu.edu.tr, ebrusengul@gmail.com)\n(Corresponding Author)\nREVIEWER COMMENTS AND RESPONSES\nReviewer #1 (Remarks to the Author):\n-I read this manuscript with great interest as it presents a novel perspective on the Arabia-\nEurasia collision zone, drawing upon a synthesis of geological and geophysical data, as well\nas thermo-mechanical modelling. This new perspective is supported by the data, and I\nrecommend that it be considered for publication in Nature Communications. However, I\nbelieve that some points still need to be clarified before it can be accepted.\nResponse: We thank the reviewer for their careful review and comments. We have modified\nthe manuscript with additional explanations, figures, and references. Our point-by-point\nresponses are listed below. Also please see the Revised tracked file that highlights tracked\nchanges and revised file for numbers of the lines.\n-Firstly, it seems that the authors are pitting their model against a 'traditional model', but the\nlatter is not clearly explained in the present paper. In my view, there are several theories\nregarding the evolution of the Caucasus region. Some previous papers suggest a significant\nrole for delamination (cf. Grosjean et al., 2022), a process which is not clearly discussed or\ntaken into account in this paper. Some papers have already proposed the contribution of the\ndeep asthenosphere to the magmatism of the area (e.g. Soltanmohammadi et al., 2018), thus\nexcluding any purely 'collisional' component.\nPrevious models did not strictly fall into a single category, and I agree that the horizontal\nmantle flow or 'plumelet' model proposed here is novel, as previous models assumed that the\nmantle flowed mainly vertically.\nResponse: Thanks to the reviewer for raising this point. We agree that the Arabian-Eurasian\ncollisional zone has been the subject of numerous studies with diverse datasets aiming to\ncharacterize processes both prior to and during its post-subduction phases. Given the journal's\nword and reference limits, a significant challenge for our manuscript was to concisely\nsynthesize this large body of previous work. For this reason, in the first section, we focus\nprimarily on neotectonics and present-day observations with brief information on the tectonic\nevolution of the region (with appropriate references), while we discuss previous models in the\nlast section by also citing Grosjean et al. (2022).\nContrary to models that attribute deformation in the region solely to lateral mantle flow,\nboundary forces, or small-scale vertical convective flows linking to the slab dynamics (e.g.,\ndelamination/detachment/break-off; Keskin, 2007; Göğüş, &Pysklywec, 2008; Bewick et al.,\n2022), we propose a more complex framework for collisional settings. Namely, we posit that\ndeformation mechanisms involve the spatial and temporal interplay of both vertical mantle\nimpulses and lateral forces owing to plumelet propagation. Yet, we agree with the reviewer\nthat it is an oversimplification to categorize all prior works under a single, unifying deformation\nhypothesis. It was not our intention to present them as such, and we thank the reviewer for\ngiving us the opportunity to clarify this point in our revised manuscript. Therefore, we have\nimproved the text with additional explanations, references and a new Figure 2 (please see the\nfollowing comment) to properly acknowledge diversity of existing models. We believe that\nFigure 2 and the revised manuscript (we thank for your contributions) will help readers to better\nunderstand the tectonic evolution of the region, including the Caucasus segment of the\ncollision. In the Introduction, we’ve cited Faccenna et al., (2013), Rolland (2017), Rolland et\nal (2012), Grosjean et al. (2022) (but not Avagyan et al., 2010 suggested by Rev.1 because of\nspace limitations) to develop our arguments-- the important role of mantle dynamics in plate\nconsumption, magmatism and partial melting in the region. These revisions can be followed\nbelow and in the tracked changes file (both in Introduction and Discussion):\nIntroduction, Tectonic context and motivation: “The lithospheric heterogeneity documented\nbeneath the region (Uluocak et al., 2021) can be attributed to Tethyan plate orogeny, as\nevidenced by numerous observations. The northern branch of the slab subducted beneath\nEurasia by creating arc and back-arc volcanism, mostly dating from the Jurassic to Eocene\n(e.g., Şengör & Yılmaz, 1981; Adamia et al., 2011; Rolland et al. 2012; Rolland, 2017 and\nreferences therein). This part of the slab broke off around the Early Eocene, triggering crustal\nshortening in the Greater and the Lesser Caucasus by the Oligocene (Avdeev & Niemi, 2011;\nCowgill et al., 2016). Concurrently, the southern branch (Bitlis-Zagros slab) initiated\nnorthward subduction, and the collisional front started to form in response to the Eocene-\nOligocene Zagros slab roll-back with the mantle window (slab gap) opening between the\nZagros and the Bitlis slabs (Rabiee et al., 2020 and references therein). Subduction jumped\nnorthward (~ Eocene, e.g., Rolland et al., 2012), and the closure of the oceanic basin persisted\ndue to the ongoing collision between the Anatolian-Southern Armenian Block and Pontides\nactive margin (Rolland et al., 2012; Shiva et al., 2025). This phase was alternatively\ncharacterized by the initiation of south-dipping subduction of the Caucasus Basin along the arc\n(e.g., Rolland 2017). Terminal closure along the Bitlis-Zagros suture (Late Oligocene–Middle\nMiocene, Fig. 2a) may have led to slab detachment beneath eastern Anatolia at ~25-10 Ma\n(Keskin, 2007; Şengör et al., 2008; Rabayrol et al., 2019 and references therein, Fig. 2b).\nOngoing convergence reactivated Pliocene uplift, folding in the Kura Basin with crustal\nshortening in the southern Greater Caucasus, intensified exhumation across the Greater\nCaucasus, and strike slip deformation to the south (e.g., North and East Anatolian faults, Figs.\n1a, 1b, & 2c) (Adamia et al. 2011; 2017; Avdeev & Niemi, 2011; Cowgill et al. 2016; Rolland,\n2017; Ismail- Zadeh et al., 2020; Molin et al., 2023). Consistent with numerous studies,\nsuch\nas river network analyses (Molin et al., 2023), the initial break-off of the southern Neotethys\nslab beneath the Bitlis-Zagros suture zone could have triggered mantle upwelling, propagating\nfrom SW to NE beneath the Arabian Plate --a plumelet (Uluocak et al., 2024, Fig. 2d). Among\nvarious slab deformation models, removal mechanisms involving spatiotemporal overlapping\nprocesses (Uluocak et al., 2024) can be primarily characterized by: i) wholesale or progressive\nBitlis-Zagros slab delamination along the Arabian collisional front in the south (Keskin 2007;\nMemiş et al., 2020 and references therein); and ii) delamination beneath the Transcaucasus\n(e.g., Shiva et al., 2025 and references therein) accompanied by northward plate tearing toward\nthe Greater Caucasus (e.g., Mumladze et al., 2015). Although 3D heterogeneities beneath the\nregion are still debated, magmatic records provide evidence of mantle-crust interactions as a\nconsequence of a discontinuous/semi-episodic closure of the Neotethys Ocean (e.g., Rolland\net al., 2012).”\n-The Tectonic context and motivation section is well-written, but is hard to follow due to the\nabsence of a figure that provides context for the tectonic evolution of the study area.\nResponse: The revised text has been enhanced by the inclusion of an additional figure -\nFigure 2- showing the generalized tectonic evolution of the region (as also recommended by\nRev.2).\n-On the other hand, in the 'Results' section, the authors immediately propose a model without\nfirst presenting the results in a cautious or analytical manner. For example, at lines 197–199,\nthey claim that 'Mantle dynamics are dominated by a focused, rapid, SW–NE-directed inflow\nfrom the Arabian foreland to the eastern Greater Caucasus—a plumelet—penetrating through\nthe complex morphology of the post-collision environment'. While I believe this interpretation\nis entirely plausible, it is crucial to present results more factually and demonstrate how they\nsupport the plumelet model.\nResponse: Upper mantle flows propagate from SW to NE in the Arabian-Eurasian collision\nzone, a region encompassing much of the western part of our study area. A similar flow was\nintroduced as a plumelet in previous work (Uluocak et al., 2024) that synthesizes/reviews\ndeformation models and observations both supporting and contesting the plumelet hypothesis\n(e.g., Ershov & Nikishin, 2004; Doi: 10.4236/ijg.2014.55047). While we previously cited this\nstudy in our tectonic context with Figure 1c (and new figure; Figure 2d), we have now\nexpanded our discussion in the revised Introduction and Discussion sections to further explain\nplumelet dynamics. Namely, although the plumelet is not a new concept (as stated in the text\nclearly), we claim that impacts of a plumelet in collisional settings, such as plumelet-induced\nplate erosion, low viscosity zone and its seismogenic implications are novel. We have modified\nrevised text to make the narrative clearer:\nIntroduction: “By coupling mantle flow models with thermal and isostatic constraints, we\nreveal not only a SW-NE oriented plumelet from the Arabian foreland injecting into the\nsubduction system (Uluocak et al., 2024), but also 3D thermomechanical erosion of the upper\nplates and the active Neotethys slab deformation along the southern collisional front.\n...Consistent with numerous studies, such as river network analyses (Molin et al., 2023),\nthe\ninitial break-off of the southern Neotethys slab beneath the Bitlis-Zagros suture zone could\nhave triggered mantle upwelling, propagating from SW to NE beneath the Arabian Plate--a\nplumelet (Uluocak et al., 2024).”\nFurther points should be clarified:\n-What role do vertical structures on the lithosphere scale play in propagating the plume? The\ninception of horizontal flow within a structure featuring vertical slabs requires zones of\nweakness inherited from the subduction-early collision history (step faults etc).\nResponse: We aim to define the complex morphology of the plumelet's propagation channel\nand its interaction with deep basins, high plateaus, and uppermost mantle cold anomalies (e.g.,\nthe Pontide-cold). We show thermomechanical erosion caused by the plumelet (e.g., the\nactively sinking cold material beneath the Arabian collisional front), by linking these processes\nto the region's subduction and early collision history. We further identify deep-seated inherited\nplate boundaries (e.g., the Aras Fault) as the inter-plate deformation between converging\nblocks. Because of limitations of the journal, further discussion will have to be reserved for\nsubsequent work. Thank you for your comment and contributions.\n- The so-called 'Bitlis slab' may also be the result of delamination, because the Bitlis\nsubduction zone is characterised by HP metamorphism at 74–71 Ma (e.g. Rolland et al.,\n2012), which likely resulted in sinking at very deep levels of the mantle (deeper than 300 km)\nand should have been N-dipping (no volcanic arc developed towards the south). This tectonic\nzone is generally considered to be complex and composed of numerous tectonic blocks and\nsutures. Therefore, the 'closure of Neotethys' should not be considered a progressive process,\nbut rather a discontinuous 'step-by-step' process in which basins were successively closed.\nThis explains why the timing of 'collision' is disputed in this region.\n-Rev. 1: maybe this has nothing to see with a Bitlis slab (more likely N-dipping...), rather\ndelamination?\nResponse: We thank the Reviewer for these comments on the Bitlis feature. We agree that\nuppermost mantle structures beneath the region are quite ambiguous and need to be clarified\nthrough high resolution observations and further evolutionary numerical models. In general,\nthese anomalies (both Pontide-cold and Bitlis slab anomalies as termed in the manuscript)\nare\ninterpreted as “Bitlis slab fragments” in observations and previous deformation models (e.g.,\nKeskin 2007; Portner et al., 2018; Memis et al., 2020). In other words, relatively small-\nscale/shallow (not deeper than ~300 km)/cold bodies beneath the collision zone, extending\nfrom the northern Arabian Platform to the eastern Pontides, are interpreted as a remnant part\nof fully detached Bitlis slab fragments. However, to open up the interpretation in the revised\nmanuscript we now describe the southern anomaly as “detached slab or delaminated Bitlis\nfragments”. We also interpret two groups of anomalies spatially separated (beneath the arc and\nthe Arabian Platform) within the complex Cenozoic tectonic concepts of the region: i) arc-\nalong deformation pattern (drip-like deformation) in the north (beneath the Pontides-\nTranscaucasus arc); ii) fully and partly detached (active sinking) Bitlis slab beneath the thick\nArabian lithospheric mantle connected with the eastward propagation of the plate tearing.\nThis\nnew concept may create a link between active eastern Zagros subduction zone along the\nArabian collisional front/Bitlis-Zagros suture zone. Observations, including magmatism and\nseismic tomography datasets, indicate a gap between Zagros subduction and its eastern part\nalong the Bitlis-Zagros suture zone -i.e., Bitlis Slab- in response to postulated Eocene-\nOligocene Zagros slab roll-back with a mantle window (slab gap) opening between the\nZagros\nand the Bitlis slabs (e.g., Rabiee et al., 2020), which we believe supports our interpretation\nhere. Following the Reviewer's advice, we revised the text to clearly express tectonic evolution\nand our findings:\nResults: “This anomaly differs from the NW Zagros depression area (P2; Figs. 6d-e) in that it\nis decoupled from the upper plate, implying fully detached cold material- termed as a part of\na detached slab or delaminated Bitlis fragments - sinking beneath the North Arabian\nPlatform.”\nTectonic context and motivation: “Magmatic records provide evidence of mantle-crust\ninteractions as a consequence of a discontinuous/semi-episodic closure of the Neotethys\nOcean (e.g., Rolland et al., 2012).”\nDiscussion: “These results reconcile previously disparate observations - from shallow (< 300\nkm), fast seismic velocities (at ~250 km depth; Koulakov et al., 2012) to overlying slow seismic\nanomalies (e.g., Koulakov et al., 2012; Civiero et al., 2023; the plumelet Uluocak et al., 2024)-\ninto a unified geodynamic framework validated by different numerical models (Uluocak et al.,\n2021) based on recent seismic tomography data (Kounoudis et al., 2020).”\n-Overall, the text is well written, but some sections could be improved (I have underlined some\nexamples in the text). Please avoid self-satisfaction (e.g. 'Our comprehensive geological,\ngeophysical and geodynamic analyses illuminate the tectonic anatomy...').\nResponse: Thank you for your detailed review. We have edited the text again, trying to\ncorrect such instances—including the underlined examples and the quoted passage (please\nsee the tracked changes file).\n-The figures are well drafted and informative. However, when dealing with seismicity (Fig.\n6), it is difficult to understand the clustering of seismicity. It might be better to present a\nrepresentative vertical section showing the seismicity in relation to the crustal levels.\nResponse: We agree with the reviewer regarding the complexity of earthquake distributions,\nalthough we already show seismicity more clearly in Figure 1b. In former Figure 6c, we only\naimed to focus on subcrustal seismic events in a specific part of the study area. Following the\nreviewer's suggestion, in the revised manuscript we have added a new figure -Figure 7d- that\nspecifically shows distribution of both crustal and subcrustal earthquakes (Magnitude of ≥ 2.5,\nUSGS, 1900-2025) as a longitudinal profile for the zone between 34oN-35oN (south of profile-\nP3). Revised Figure 7 (former Fig. 6) has also been improved by an inset to clarify directions\n(east, west, north, & south, depth) used in revised Fig. 7b-c. We think that these additional\nfigures, as well as Figure 1b, clearly image how the seismically active region in the east differs\nfrom the west of the collisional front, which is subject to thermomechanical erosion.\n-I also think the title does not fully reflect the content, as it should emphasise the 'plumelet'\nhypothesis.\nResponse: We thank the reviewer and carefully considered their suggestion, but respectfully\nwe choose to retain the original title. We wish to emphasize that our study focuses on the\nanatomy of post-subduction collisional settings with complex deformation patterns.\n-Not clear: “...the plumelet may reworks the subducted lithosphere actively converging with\nEurasia”\n-the plumelet supposedly causes break-off? “is that really considered \"collsion related\ndeformation\"?\n-is there just one traditional model (then what is it?) or a multiple pannel of proposed models?\n“In the overall tectonic anatomy of the complex plate boundary, during its migration from deep mantle source(s) in the south, the plumelet reworks the subducted lithosphere actively converging with Eurasia. As well as causing collision-related deformation (e.g., slab breaking-off), the plumelet induces thermomechanical erosion of the foundering slab at the active collisional margin. The results discover a more intricate multi-scale behavior to the post-collisional tectonic boundary than the traditional model...”\nResponse: The multiple previous deformation models (including numerical experiments)\nmostly focus on the Neotethyan slab dynamics involving advance, retreat, and eventual break-\noff of the slab to explain complex post-subduction tectonics of the region. Here, we aim to\nclarify forces caused by mantle dynamics in the form of relatively large- (plumelet) and small-\nscale convective cells (plumelet-induced upwellings) and provide a link between tectonic\nevolution and present-day deformation patterns. In addition to the thermal and isostatic\nconstraints shown in the text, our results are supported by numerous observations, such as\nintense magmatism with mantle-sourced fingerprints and high plateau evolution (pls see Fig. 2\nin Uluocak et al., 2024, below). All indicate that convective forces may have been active since\nthe Pliocene and could have formed post-subduction tectonics in the Arabian-Eurasian\ncollisional zone. We have rephrased the passage quoted by the Reviewer to clarify:\nDiscussion: “In the tectonic anatomy of the complex plate boundary, during its migration from\ndeep mantle source(s) in the south, the plumelet may rework the subducted lithosphere actively\nconverging with Eurasia. That is, the plumelet plays an active role in plate consumption, as\ndemonstrated by thermomechanical erosion in the overlying plate and sinking slab along the\nactive collision margin and by drip-type removal beneath the volcanic arc. The results reveal a\nmore intricate multi-scale behavior of the post-collisional tectonic boundary than the\ntraditional models (e.g., Keskin, 2007; Memiş et al., 2020) based solely on slab dynamics and\nboundary forces.”\n-what is a foundering slab?\nResponse: We used this term referring to slab detachment processes in a broad sense, but the\nphrase has been deleted to avoid possibly being misleading.\n-Plate tectonic collision zones can experience prolonged crustal and upper mantle\ndeformation that defines often(marked by Rev) post-subduction tectonics in continental\norogenic settings.\nResponse: Removed.\n-I hope these comments are useful. With best regards.\nResponse: We sincerely thank you for your careful review and constructive comments.\n-refs\nGrosjean, M., Moritz, R., Rezeau, H., Hovakimyan, S., Ulianov, A., Chiaradia, M., &\nMelkonyan, R. (2022). Arabia-Eurasia convergence and collision control on Cenozoic\njuvenile K-rich magmatism in the South Armenian block, Lesser Caucasus. Earth-Science\nReviews, 226, 103949.\nSoltanmohammadi, A., Grégoire, M., Rabinowicz, M., Gerbault, M., Ceuleneer, G.,\nRahgoshay, M., ... & Benoit, M. (2018). Transport of volatile-rich melt from the mantle\ntransition zone via compaction pockets: implications for mantle metasomatism and the origin\nof alkaline lavas in the Turkish–Iranian plateau. Journal of Petrology, 59(12), 2273-2310.\nRolland, Y., Perincek, D., Kaymakci, N., Sosson, M., Barrier, E., & Avagyan, A. (2012).\nEvidence for∼ 80–75 Ma subduction jump during Anatolide–Tauride–Armenian block\naccretion and∼ 48 Ma Arabia–Eurasia collision in lesser caucasus–East Anatolia. Journal of\nGeodynamics, 56, 76-85.\nAlizadeh Noudeh, S. Yann, R., Magali, R., Delphine, B., Philippe, M., Arthur, I., ... &\nRahgoshay, M. (2025). Geochronological, isotopic and petrogenetic investigations of\nCenozoic Volcanic rocks in the Talysh Massif, NW Iran: Insights for the Eocene magmatic\nflare-up. Lithos, 496, 107954.\nResponse: The revised text has been modified by adding the suggested references, we thank\nthe Reviewer for their contribution.\nReviewer #2 (Remarks to the Author):\n# Anatomy of a Post-Subduction Collision\n**Authors:** Ebru Şengül Uluocak, Russell N. Pysklywec, Claudio Faccenna & Taylor\nSchildgen\nThis paper addresses the problem of understanding mantle dynamics in continental collision\nzones where the driving oceanic slab has detached or broken off from the overriding plate.\nThese post-collisional systems are characterized by prolonged crustal and upper mantle\ndeformation that can persist long after the initial convergence and slab consumption has\nceased. The Arabian-Eurasian collision zone exemplifies such complexity, where the\nNeotethyan oceanic slabs have undergone various stages of subduction, break-off, and\nfoundering, leaving behind a fragmented lithospheric architecture. Understanding what drives\ncontinued deformation in these post-subduction environments remains a fundamental\nchallenge in tectonics, as traditional plate-driving forces are absent yet significant surface\ndeformation, volcanism, and seismicity persist.\n## Overview\nThe manuscript is well-written and organised and it is of significance because it showcases\nhow mantle dynamics strongly influences strain rate, stress, topography, and mantle flow of\npost-collisional settings. It also demonstrates how tectonic activity post–dates the interruption\nof subduction due to the effect of mantle dynamics and provides a framework for the\nintegration of geodynamic modelling with geophysical observables.\nThe approach is based on correlating various geophysical observables - seismic velocities,\ntopography, crustal thickness, and seismicity patterns - with their 3D thermomechanical\nmodel results. The authors systematically demonstrate how regions of high Vp/Vs ratios\ncorrespond to magmatic arcs and elevated topography, how cold mantle anomalies correlate\nwith thick lithospheric roots and deep basins, and how modeled strain rates align with\nobserved seismicity distributions. It is these (positive) correlations that form the empirical\nfoundation for validating their model predictions and serve as the primary evidence\nsupporting their interpretations.\nResponse: We sincerely thank Reviewer 2 for their careful review.\n-Another central concept is the \"plumelet\" which is defined by the authors as the migration of\nupper-mantle plume material into the region. It may be drawn in by flow in the region or\nby its own buoyancy. The plumelet is proposed as a flow of hot mantle material from the\nArabian region to the eastern Greater Caucasus. They propose that, from their 3D\nthermomechanical modeling, this lateral mantle flow *drives* regional deformation,\ninfluences slab break-off processes, *causes* thermomechanical erosion of overlying plates,\nand creates the observed dynamic topography patterns across the Arabian-Eurasian collision\nzone. There are a number of other areas where subduction influences mantle flow - for\nexample the paper by Miller & Becker (https://doi.org/10.1130/G34959.1) argues that uplift\nof the Atlas is enhanced by flow of plume material around the foundering Alboran slab. Or\nindeed the paper by Lanari et al (2023) which has some author overlap with this paper. Are\nthe other notable examples that would help to establish this plumelet hypothesis and\ngeneralise it ?\nResponse: We appreciate this constructive feedback. Unfortunately, due to the journal's\nreference and space limitations, we were unable to cite many other articles in the manuscript.\nHowever, we have rephrased this section with suggested reference (Miller & Becker) by\nremoving Lanari et al., 2023 (due to the journal's reference limitation)- that may help build\nand generalize the plumelet argument.\n-A fundamental question raised by this interpretation is whether the plumelets **are**\nactively driving the observed tectonic signals, or whether they represent more passive\nstructures that are being drawn into the region by the process of slab breaking off and\nfoundering. The distinction between active versus passive roles has significant implications\nfor understanding the causal relationships in post-collisional tectonics and the broader\napplicability of the plumelet concept to other collision zones. Is the plumelet genuinely\ncausing these effects? Can it be determined convincingly through the models?\nResponse: Based on our analyses, we believe that mantle dynamics play an active role in\ndriving tectonics in the region. As well as the thermal and isostatic constraints provided in the\nmanuscript, it can be concluded from previous studies (e.g., Fig. 2 by Uluocak et al., 2024 and\nreferences) that convective forces may have been active at least since the Pliocene and formed\npost-subduction tectonics in the Arabian-Eurasian collision zone. The link between present-\nday observations and post-subduction related deformation patterns imaged in our work here\nprovides support that mantle flow is playing an active role in the tectonics of the region. This\napproach gives a new and integrated assessment for the collisional settings. We revised the text\nto make clearer our findings and conclusion as follows (Discussion):\nDiscussion: “In the tectonic anatomy of the complex plate boundary, during its migration from\ndeep mantle source(s) in the south, the plumelet reworks the subducted lithosphere actively\nconverging with Eurasia. That is, the plumelet plays an active role in plate consumption, as\ndemonstrated by thermomechanical erosion in the overlying plate and the sinking slab along\nthe active collision margin, and by drip-type removal beneath the volcanic arc. The results\nreveal a more intricate multi-scale behavior of the post-collisional tectonic boundary than the\ntraditional models (e.g., Keskin, 2007; Memiş et al., 2020) based solely on slab dynamics and\nboundary forces.”\n-In Figure 1, the main magmatic units are shown; however, it would be beneficial to provide\ninformation about the composition of these volcanic units. Therefore, a figure showing the\ntectonic evolution, the location of magmatic units and their chemistry is needed. This is very\nimportant for readers who are not familiar with the geology of the study area.\nResponse: The revised manuscript has been improved with an additional figure (new Figure\n2) and references to clarify the tectonic evolution of the region. Because of the journal’s\nlimitations, providing another figure with discussion on magmatism (chemistry, ages and even\nsamples’ locations) seems not feasible/possible. Yet, following the reviewer’s suggestions,\nwe\nhave improved the text with additional references and explanation:\nTectonic context and motivation: “Post-collisional volcanism initiated (around in the early–\nmid-Miocene) on the Erzurum-Kars Plateau (locations shown by “T” in Fig. 1d) and migrated\nsouth and eastward, towards the Bitlis-Zagros Suture and the Armenian border, respectively,\nconsistent with radiometric ages (e.g., Lin et al., 2020; age compilation in Memis et al., 2020).\nThe late Miocene-Quaternary calc-alkaline to alkaline magmatism is documented across the\ncollisional front with ages decreasing westward and subduction signatures weakening\nsouthward, towards the Bitlis-Zagros suture (Keskin, 2007 and references therein). Based on\nchemical and petrological analyses of widespread alkaline lavas active since the Late Cenozoic\nin the Turkish-Iranian Plateau, it has been argued that the volcanism may be linked to a low\nseismic velocity zone/body extending to depths of ~300 km (Soltanmohammadi et al., 2018).”\n-In general, figures are very pretty but not always easy to follow. Figure 5 is quite\nstraightforward, but Figure 6a is not. Figure 6b is difficult to see compared to figure 6c.\nWe\nknow the 3D renderings are a lot of work to produce but in static form they are less useful\nthan the cross sections (sorry).\nResponse: We understand the complexity of the figures raised by the Reviewer. For the revised\nmanuscript we have improved figures (Figures 1d, 3, 5, 6, 7) with different perspectives/map\nthat may help to clarify our results here. Namely, we have improved revised Figure 3 (previous\nFig. 2) with additional figure (Figure 3a) to clarify locations of 3D anomalies and profile in\nrevised Figure 3b and revised Figure 7. We also improved Figure 7 (former Fig.6) with an\nadditional figure and an inset. Figure 7d images crustal and subcrustal earthquake distributions\nalong the profile P3 for a zone between 34oN-35oN (“The longitudinal distribution of\nearthquakes between 34oN-35oN (USGS, 1900-2025). Red stars indicate thermomechanically\neroded zones.”). An inset figure shows cartesian coordinates to make the 3D orientations\nclearer.\n-Figure 6 is helpful because it integrates all the observables and the results from the numerical\nmodelling. However, the mantle structure of the study area is quite complex and the 3D\ngeometry of the plumelet propagation channel is not fully clear. We suggest you could draw\nyour interpretation of this geometry in Figure 6b, and this will help the reader to better visualize\nhow it correlates with variations in stress, strain rate, dynamic topography and Vp/Vs ratios.\n- Lines 333-336: The vectors in Figure 6 can be improved. This would facilitate the\ncomparison with the observed GPS velocity field. Also, Figure 6 shows a zoomed region, but\nthe extent of this zoomed region is not clear. It would be useful to include an inset with a\nregional map showing the location of this region.\nResponse: We agree with the reviewer that cross-sections projecting results onto a 2D plane\nare easier to follow. However, they can also be oversimplified and potentially misleading.\nThus, we have included both 2D and 3D images to accurately represent the complex\ndeformation in the collisional zone. In fact, we aim to make this complexity visible in terms of\ndeformational patterns in the collisional setting.\nIn the text, as well as 3D views (revised Figure 3, 4 and 7) in the manuscript, we show plumelet\npropagation and geographic locations in the revised Figure 3 that includes cross-section of P1\nalso shown in revised Figures 6 and 7. Namely, in response to Reviewer’s concern regarding\nthe link between geographic locations and anomalies, the revised Figure 3 has been improved\nwith a new map (revised Figure 3a) showing locations of the profile (P1) and various features\n(names, arc, plateaus). By following the suggestions, we have also improved revised Figure 7\n(previous Fig. 6) by adding new figure (revised Fig. 7d) to establish the link between our\nresults, the locations of deformations, and seismicity. Furthermore, in order to avoid additional\ncomplexity of Figure 7, an inset has been added in revised Figure 7 to clarify the coordinate\nsystem and orientations of 3D images in Figures 7a-c. Captions and text have been revised\naccordingly. Due to the journal’s space limitations, we address Figure 1 to compare GPS\norientation with our results. We thank the reviewer for this comment, which has helped us\nimprove the manuscript.\n## Detailed Feedback\n- Pn and Sn are quite well known but not universally, it would be useful to add a short\nexplanation\nResponse: This part has been improved as follows:\nTectonic context and motivation: “Comparable high-velocity anomalies at shallow depths (≤\n~ 150 km) are observed beneath the Black Sea, Mesopotamian-Zagros foreland basins, and\nZagros fold-thrust belt, consistent with high Pn seismic wave velocity perturbations in these\nregions …\nIn contrast to high amplitudes of guided seismic waves (Pn and Sn) that travel in the\nlithospheric mantle, crustal seismicity and slow seismic velocities correlate well with zones of\nhigh Sn and Pn attenuation, indicating thinned or absent lithosphere (≤ 50–90 km) beneath the\nEast Anatolian, TGA, the North Iranian plateaus and …”\n- Lines 86-89: In the captions of Figure 1d, are the red triangles showing the magmatic units\nor active volcanoes?\nResponse: Volcanoes. The caption has been revised to clarify.\n- LInes 149-150: It is worth mentioning which observations this refer to (geological,\ngeophysical, geodynamic) ?\n“The northern 148 branch of the slab subducted beneath Eurasia by creating arc and back-arc volcanism,\nmostly 149 dating from the Early Cretaceous to Eocene (e.g., Şengör & Yılmaz, 1981; Adamia et al., 2011).“\nResponse: The tectonic context and motivation section has been revised with an additional\nfigure (Fig.2), explanations, and related references to clarify evolution history and\nobservations. This revised description addresses geochronological, isotopic and petrogenetic\ninvestigations of volcanics of the region (in an abbreviated form owing to the journal’s space\nlimitations).\nTectonic context and motivation: “The northern branch of the slab subducted beneath Eurasia\nby creating arc and back-arc volcanism, mostly dating from the Jurassic to Eocene (e.g., Şengör\n& Yılmaz, 1981; Adamia et al., 2011; Rolland et al. 2012; Rolland, 2017 and references\nintherein)”\n- LInes 204-205: Figure 2 would benefit from adding a label with the location of the Scythian\nplate.\n- Line 228: Is there any reason for Figure 2 to display the north at an orientation different\nfrom all other figures?. Although you include the labels of the main tectonic domains, it is\nconfusing.\nResponse: In the revised manuscript, Figure 3 (previous Fig. 2) has been modified for clarity.\nParticularly in this figure, we aim to show plumelet propagation with convection vectors clearly\nwith a view that is different from the revised Figures 6 & 7. We added a new topographic map\n(revised Figure 3a) showing the profile location (P1) with relevant place names. The location\nof the Scythian Plate is now included in the revised Figs. 3a and b. Also, the figure captions\nhave been revised in accordance with the Reviewer's comments.\n- Lines 221-222: differences of less than 1MPa do not seem that “high”\nResponse: This statement has been modified in the revised text:\nResults: “The flow modelling indicates regional deformation zones in the central and eastern\nGreater Caucasus extending to the stable Scythian Plate, marked by visible lateral contrast of\nlow shear stresses in the crust (≤ 1 MPa; Fig. 4b)…”\n- Lines 249-251: you mentioned “crustal velocity ratios”; however, you didn’t mention the\ncrustal depth or depths considered in the contours.\nResponse: To avoid complexity in revised Figure 5 (previous Fig. 4), only average values of\nthe Vp/Vs ratios were obtained at depths of 20 km and 40 km (Zhu, 2018). These depths have\nbeen added to the revised caption for Figure 5.\nFigure 5. “… Average of Vp/Vs ratios for the crust (at 20 km and 40 km depths) and Vp/Vs\nratios for the depth of 100 km (Zhu, 2018).”\n- Lines 274-275: What is the purpose of the vertical line at ~35 km of crustal thickness?. Is it\nrelated to the “A” zone?. This has to be clarified in the captions of Fig. 4b.\nResponse: That is right, it is related to south Armenian highland; the caption of figure has\nbeen corrected to clarify this.\nFigure 5b: “A” stands for the Armenian volcanic highland with relatively thin crustal thickness\n(gray dashed line) and is interpreted as concerning partially melted lower crustal features (Lin\net al., 2020).”\n- Lines 253-254: “corroborates” is not the most suitable word for your explanation\nResponse: Thank you for the correction, it has been rephrased:\nResults: “This result demonstrates a direct mantle-crust linkage in concordance with crustal\nand subcrustal thermal anomalies beneath young volcanic centers (Fig. 5a) and our 3D\ntemperature model (Figs. 3b & 4d).”\n- Lines 283-284: What does “Pontide cold” mean?\nResponse: We use this term to indicate relatively shallow localized cold bodies by following\nthe same terminology in Uluocak et al. (2019).\n- Lines 296-298: Interpretation should go into the discussion\n“ We interpret an active slab breaking off beneath the East Greater Caucasus with a cold\nthick root of stable Scythian Plate at the northern boundary of the collision (P2 in Fig. 5e).”\nResponse: The statement has been rephrased as follows:\nResults: “We define an active slab breaking off beneath the East Greater Caucasus with a\ncold thick root of stable Scythian Plate at the northern boundary of the collision (P2 in Fig.\n6e).”\n- Lines 309-310: Does “shaded” refer to the areas filled with red?. This has to be clarified.\nResponse: Thank you for the correction, it has been done:\nFigure 6. “… observed topography (gray) along profile..”\n- Line 369 - 'indicating upper mantle-driven isostatic instabilities'. This is quite a strange\nphrase because isostasy implies a static, buoyancy equilibrium whereas this paragraph refers\nto mantle-driving of topography (which implies non-isostatic support)\n“Our integrated analysis (Fig. 4) demonstrates a systematic correlation between high average Vp/Vs ratios (≥\n1.8) at elevated regions (≥ 1 km) and relatively thin crusts (≤ 40 km), indicating upper mantle-driven isostatic\ninstabilities. The Turkish-Georgian-Armenian (TGA) volcanic provinces, East Anatolian Plateau, western\nGreater Caucasus, Talesh (Azerbaijan), and Urumieh-Dokhtar magmatic zone are defined by isostatically\nunder-compensated topography with thermally heated crustal and uppermost mantle structures. “\nResponse: Agreed--these variations show non-isostatic support from mantle; to clarify the\nstatement ‘isostatic instabilities’ has been rephrased as `non-isostatic support` in the revised\nmanuscript.\nDiscussion: “Our integrated analysis (Fig. 5) demonstrates a systematic correlation between\nhigh average Vp/Vs ratios (≥ 1.8) at elevated regions (≥ 1 km) and relatively thin crusts (≤ 40\nkm), indicating upper mantle-driven non-isostatic support.”\n- Lines 403-406 Can you elaborate on the observed changes in mantle composition?. You\nmentioned this in the “Tectonic context and motivation” section, but this is an opportunity to\nlink the specific changes to your modelling results.\n“This interpretation is409 corroborated by relatively local crustal thinning, thermal heating, and magmatism\nwith mixed410 geochemical signatures documented in the Armenian volcanic highland (Zabelina et al.,\n2016;411 Neil et al., 2015; Bewick et al., 2022). Although region-specific geodynamic models (e.g., Curie412\net al., 2015) provide additional kinematic constraints, our drip hypothesis establishes a useful413 approach for\ninterpreting both plumelet and Pontide-cold anomaly and possible related414 magmatic arc cycle with\nobservational (geological-magmatic-tectonic) signatures changing415 from arc to an intraplate deformation\nabove former subduction zones within continental416 collision systems (e.g., Burg, 2011; Allen et al., 2013;\nDucea et al., 2015; Maierová et al., 2025).”\nResponse: We agree with Reviewer 2, but since we have exceeded the 5000-word limit and\nhave already made the relevant explanation in the Introduction, expanding on this part is just\nnot feasible.\n- Lines 416-420 Shortening this sentence would improve clarity\nLines 416-420 in text with shifted lines: “Although region-specific geodynamic models (e.g., Curie et al., 2015)\nprovide additional kinematic constraints, our drip hypothesis establishes a useful approach for interpreting both\nplumelet and Pontide-cold anomaly and possible related magmatic arc cycle with observational (geological-\nmagmatic-tectonic) signatures changing from arc to an intraplate deformation above former subduction zones\nwithin continental collision systems”\nResponse: As suggested by the Reviewer, we have shortened this sentence to make the\npassage more succinct:\nDiscussion: “Although region-specific geodynamic models (e.g., Curie et al., 2015) provide\nadditional kinematic constraints, our drip hypothesis establishes a useful approach for\ninterpreting both a plumelet and Pontide-cold anomaly. We therefore claim that there is a\nPontide-cold-related magmatic arc cycle here supported by observational (geological-\nmagmatic-tectonic) signatures changing from arc to an intraplate deformation above former\nsubduction zones within continental collision systems….”\n- Line 469 - advection / diffusion equation is scrambled. It would be helpful to just write\nthese out with equation numbers as it is really difficult to read these equations in-line.\nResponse: This section has been revised with the expressions on separate lines (with\nequation numbers) to make the governing equations clearer.\n## Recommendation\nThis paper covers a novel and interesting topic and is potentially publishable in Nature\nCommunications. The writing is currently somewhat breathless and that seems to get in the\nway of readability. Similarly, the figures are elegant but they are cluttered in a way that gets\nin the way of understanding. We are happy to recommend publication upon addressing\ncomments\nResponse: We appreciate Reviewer 2 for the providing a careful review and constructive\ncomments.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This\nis part of the Nature Communications initiative to facilitate training in peer review and to\nprovide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We commend Nature Communications' for its initiative to train Early-Career\nResearchers.\nReviewer #4 (Remarks to the Author):\nReview of “Anatomy of a post-subduction collision” by Uluocak et al.\nThe manuscript by Uluocak and co-authors presents an interesting study on the relationship\nbetween mantle thermal structure and various observables from the Arabian–Eurasian\ncollision zone. Overall, the paper is well written and addresses a topic of broad relevance.\nThe potential implications are significant. However, I believe the current modeling approach\nrequires substantial strengthening before the conclusions can be fully supported.\nMajor comments\n1) Modeling strategy\n-To my understanding, only a single simulation is presented. This limits the robustness of the\nconclusions. Sensitivity tests exploring variations in the initial thermal structure (inferred\nfrom geophysical constraints) are essential. Similarly, the role of mantle rheology, crustal\narchitecture, inherited sutures, and other weak zones should either be explicitly investigated\nor at least discussed in detail.\n-A direct comparison between model outputs and observations is missing. For example,\nmodeled surface velocities should be plotted against GPS data in a single figure. Such\ncomparisons are crucial to convincingly demonstrate the model’s validity.\n-L328-L330. This is rather difficult to see both modelled velocities and GPS velocities are\nnot super-imposed on a figure. Having a figure showing how well the model and the data fit\nwould really help to convince the reader.\nResponse: We thank the reviewer for these important comments. We acknowledge that the\nmanuscript presents a single numerical model to support its conclusions. This approach was\nnecessary due to the space constraints of Nature Communications. However, the chosen model\nis the culmination of an extensive testing process, including comparisons with re-designed\nmodels from literature, the details of which we have already discussed in our previous studies\n(i.e., sensitivity of the model). Besides, as detailed in the Results and Discussion, our modeling\nresults are strongly corroborated by independent datasets, including a new seismic tomography\nmodel and various geophysical and geodetic data, each sensitive to different physical properties\n(e.g., Discussion: The calculated mantle flow vectors align closely with observations, such as modern kinematics\n(Fig. 1a), seismic anisotropy orientations and a low seismic velocity zone (Fig. 1c) that extends to the East African\nplume-magmatic belt at large-scales (Ershov & Nikishin, 2004; Civiero et al., 2023). These results reconcile\npreviously disparate observations - from shallow (< 300 km), fast seismic velocities (at ~250 km depth; Koulakov\net al., 2012) to overlying slow seismic anomalies (e.g., Koulakov et al., 2012; Civiero et al., 2023; the plumelet\nUluocak et al., 2024)- into a unified geodynamic framework validated by different numerical models (Uluocak et\nal., 2021) based on recent seismic tomography data (Kounoudis et al., 2020).”).\nOur modeling strategy, as highlighted by Reviewer 2, further, involves correlating the results\nof our 3D thermomechanical model with key geophysical observables (i.e., direct comparison\nbetween model outputs and observations). We present a systematic analysis, for the first time\nin the region, showing positive correlations between high Vp/Vs ratios and magmatic arcs with\nhighlands; between cold mantle anomalies and regions of thick lithospheric roots/deep basins,\nand between estimations and observed seismicity patterns. It is the strength of these consistent\nrelationships that validates our model predictions and forms the core evidence for our\ninterpretations. In response to Reviewer-4, we have improved the Method section to clarify\nvalidation of the numerical model:\nMethod: “Since different numerical models with various temperature configurations and\nrheological features were tested (Uluocak et al., 2021; 2024), in this study we followed a\nconvection modeling approach with inputs that have already been confirmed/validated by\nconsidering large sets of geophysical data and analyses documented in the Arabian-Eurasian\ncollisional region.”\nWe evaluated the Reviewer’s suggestion to include a direct comparison of surface velocities\nwith GPS data, which would require further discussion on spatially resolved geodetic\nmeasurements in the region. However, given the stringent space and reference limitations, we\nhave chosen to focus the manuscript on the core geodynamic analyses. We believe the existing\ncorrelations presented between our model results, major kinematics of plates, seismicity, and\nother geophysical observables provide robust and sufficient evidence for our arguments.\nTherefore, we feel that adding such a figure with references and discussion, while potentially\ninteresting, is not essential to the central conclusions of this study.\n2) Model setup and methodology\n-The modeling approach is not sufficiently detailed in its current form. Model dimensions and\nresolution should be specified.\n-Even in instantaneous flow models, runs should be evolved long enough to achieve isostatic\nequilibration. What timescales were used?\n-Initial conditions must be clearly introduced: what is the starting thermal and density\nstructure?\n-A table summarizing all material properties should be provided, ideally in the appendix.\n-What is your model size and resolution?\n- You also need to provide a table with all used material parameters.\nResponse: We thank the reviewer for their comments on expanding the explanation of the\nmodel. Following the suggestions, we have added a new Table 1 to the revised manuscript that\nlists the material parameters. To accommodate this, we have removed some references from\nthis section (Glerum et al., 2018; Uluocak et al., 2016, 2019) and added the relevant ones to\nthe table (Ranalli, 1995; Naliboff & Buiter, 2015). Furthermore, we have improved the\nMethods section with additional explanations regarding the model's spatial resolution, time\nscale, and the initial condition of the model.\nMethod: “…We used seismic data (Piromallo & Morelli, 2003) to constrain the geometry of\nlithospheric structures and mantle features for the model (Uluocak et al., 2024: Uluocak,\n2024a). Since estimations of instantaneous thermomechanical models are data-dependent,\nsensitivity analyses for the modeling approach used here were performed and documented in\nprevious work on the collisional zone …The temperature model (e.g., Figs. 3b, 4d & 5e) was\nobtained by using seismic data with a density scaling approach based on the thermal-expansion\nexpression: ρ(T)= ρ (1- α ΔT) (reference densities, ρ , are listed in Table 1; a coefficient of\n0 0\nthermal expansivity, α = 2x10-5 1/K was used) with the sampling intervals of 20 km in each\ndirection for the investigated region (Uluocak et al., 2024). The time-stepping in the model was\ncarried out using a Courant–Friedrichs–Lewy (CFL) number of 0.1, yielding a timestep of up\nto ~1x105 years and relatively small velocity errors (i.e., RMS ≤ 2 cm/year) with the effective\ncomputation time (Uluocak et al., 2021).”\n3) Figures and presentation\nFigure 2: The contours are unclear. What variable is being contoured? This needs\nclarification.\nResponse: Based on Reviewers’ suggestions, color scales, and missing location names have\nbeen renewed with an additional map to make more visible of deformation patterns in revised\nFigure 3 (previous Fig. 2). The caption has also been improved to more clearly indicate\nwhich variables are shown:\nFigure 3. “Iso-surfaces are imaged only for ≤ 1000 K (see color scale in Fig. 4d).”\n-Figure 3A: The method used to calculate dynamic topography must be described.\n-Fig. 3a. How do you compute the dynamic topography? This is quite critical to explain how\nthis is computed in your model as this is not trivial especially for non-modelers.\nResponse: We have improved the revised text with additional explanation of how we\ncalculated dynamic topography to make it more accessible to non-modelers and the broad\nreadership of Nature Communications. Due to journal space limitations, we also refer readers\nto specialized studies for a more detailed discussion of the modeling techniques. That is, we\nused the traditional approach to obtain the surface deflection of the top of the model solution\nspace as reported in previous studies:\nMethod: “Dynamic support in the topography- dynamic topography (e.g., Faccenna et al.,\n2014)- is defined along the free surface based on the vertical stress exerted on the lithosphere\nas a response to normal fluid stresses in the underlying mantle.”\n4) Data availability.\nIf I recall correctly, Nature Communications requires that the data necessary to reproduce the\nmodels and generate the figures be made publicly available, for example via Zenodo.\nResponse: Agreed. ASPECT is an open-source code available to the general public. As is\nstandard in publications with models, our input files and ASPECT (with open-source functions\nin ASPECT) can be compiled from Uluocak et al., (2024) and Bangerth et al. (2022) as cited\nin the manuscript. Seismic data used in modeling is also available in Piromallo & Morelli\n(2003), as stated in the text (e.g., Introduction, Method). Besides, following the reviewer’s\nsuggestion, we have enhanced the Method and Data Availability sections with an additional\nreference (Uluocak 2024, below) to clarify accessibility and complicitly of our modeling input\ncode via Zenoda, ensuring the same level of openness as the ASPECT code itself.\nUluocak, E. Ş. (2024). Dataset for “The role of upper mantle forces in post‐subduction\ntectonics: Plumelet and active rifting in the East Anatolian Plateau” (G‐Cubed). [Dataset].\nZenodo. https://doi.org/10.5281/zenodo.12568287\nIn accordance with our commitment to open science, the corresponding author (ESU) will\nprovide assistance to researchers seeking to compile the data upon reasonable request. The\nData Availability Statement has been now updated to the following standard explanation from\nNature Communications to formally confirm our compliance:\nhttps://www.nature.com/documents/nr-data-availability-statements-data-citations.pdf\nData availability: “Open-sourced functions/codes and data used in the numerical model can\nbe compiled from Uluocak, (2024) that are publicly available. The data that support the findings\nof this study are available from the corresponding author upon reasonable request.”\n-The manuscript presents a substantial body of work and demonstrates the authors’ strong\nunderstanding of the Arabian–Eurasian collision zone. However, the modeling component\nrequires further work and sensitivity testing to ensure the robustness of the conclusions.\nAddressing these points would significantly strengthen the study and enhance its impact.\nResponse: We would like to clarify that a sensitivity analysis for a similar modeling approach\nwas performed in previous work (Uluocak et al., 2021, 2024). The current manuscript expands\non that foundation to provide a deeper insight into deformation mechanisms in collisional\nsettings. The Method section has been revised to explain this:\nMethod: “We used seismic data (Piromallo & Morelli, 2003) to constrain the geometry of\nlithospheric structures and mantle features for the models. Since estimations of instantaneous\nthermomechanical models are data-dependent, sensitivity analyses for the modeling approach\nused here were performed and documented in previous work on the collisional zone (Uluocak\net al., 2021, 2024).”\nLine per line comments:\nL22. I am not quite sure to know/understand what is a “plumelet” at that stage (I am not on\ntop of this litterature). As this is likely a key concept of the paper you may want to clarify this\nfor a general audience, right from the start. Is that a small plume? Is it super-adiabatic? Is its\nbuoyancy chemically controlled? (in Abstract: collision—an archetype of post-subduction tectonics. Our\nkey findings reveal that plumelet23plate interactions drive deformations both within and at the margins of\nconvergent plates..)\nResponse: To avoid repeating our earlier work where the plumelet concept was first\nestablished in this region, we provide a succinct explanation -not in the Abstract but -in the\ntext:\nTectonic context and motivation: “…including a large-scale low-velocity zone extending\nfrom Afar to the Greater Caucasus (≤ ~300 km depth, Fig. 1c, e.g., Civiero et al., 2023),\nthe\nlatest numerical model by Uluocak et al. (2024) implies a lithospheric channel with SW-NE\nmantle flows (Ershov & Nikishin, 2004) –– a plumelet (i.e., upper-mantle plume migration\nwithout significant tail and a mushroom head; Uluocak et al., 2024)––here.”\nL54. “and the geodynamic model of the upper mantle and lithosphere structures”. Is that\nbased on previous work and reconstructions? You may need a reference here. (“the geodynamic\nmodel of the upper mantle and lithosphere structures, we resolve the 55 three-dimensional (3D) mantle\ncirculation patterns, including large- and small-scale convective 56 cells, and their dynamic interactions with\noverlying plates. Different from previous studies, 57 including numerical experiments in this part of the\ncollision (Fig. 1b, e.g., Faccenna et al., 2014; 58 Uluocak et al., 2021; 2024 and references therein),\nResponse: Thank you for the reviewer, this part has been rephrased to avoid\nmisunderstanding as follows:\nIntroduction: “Based on seismic imaging (Piromallo & Morelli 2003) and our geodynamic\nmodel of the upper mantle structures, we resolve the three-dimensional (3D) mantle circulation\npatterns with large- and small-scale convective cells, and their dynamic interactions with\noverlying plates.”\n§ L51-L65. It is not clear how the model is built from geophysical and geodynamic data at\nthat stage.\n(51 Here, we focus on the Arabian-Eurasian collision region (Fig. 1) as a plate collision\n52 archetype where a multitude of post-subduction mantle dynamics are actively shaping the\n53 surface geology and tectonic configuration. Based on seismic imaging (Piromallo & Morelli\n54 2003) and the geodynamic model of the upper mantle and lithosphere structures, we resolve the55 three-\ndimensional (3D) mantle circulation patterns, including large- and small-scale convective56 cells, and their\ndynamic interactions with overlying plates. Different from previous studies,57 including numerical experiments\nin this part of the collision (Fig. 1b, e.g., Faccenna et al., 2014;58 Uluocak et al., 2021; 2024 and references\ntherein), our 3D thermomechanical model, extending59 east to the Caspian Sea, is sensitive to both large (e.g.,\nplumelet, Uluocak et al., 2024) and60 regional-scale changes, which may correspond to specific zones with\nthinned/missing61 lithospheric mantle (e.g., the proposed delamination-modified regions of eastern Anatolia\nand62 western Greater Caucasus), high plateaus (e.g., Turkish-Georgian-Armenian-TGA, the East63 Anatolian\nand the NW Iranian plateaus), deep basins (e.g., Mesopotamian-Zagros, Kura,64 eastern Black and western\nCaspian seas basins) and plate boundaries (e.g., Bitlis-Zagros suture65 and fold-thrust belt) (Fig. 1).)\nResponse: According to the journal’s format we introduce our aim, approach, and results\nbriefly in the Introduction and expand more in the Method. And as cited from the various\nsources, we use geophysical data such as seismic to constrain the geometry of lithospheric\nstructures and mantle features for the model. The citations from previous studies reference the\nbasis for various configurations and assumptions we make for the model. In the revised\nmanuscript we have made a number of changes in response to specific points from all three\nReviewers to clarify aspects of the setup and design of our model. We believe that these\ncitations and expanded explanation together provide a clearer explanation of how the model\nwas built.\nL68-L69. “By coupling mantle flow models with thermal 68 and isostatic constraints” this\npartially answer my previous comment. Maybe this should be moved in previous section.\n(67 Our comprehensive geologic-geophysics-geodynamic analyses illuminate the tectonic\n68 anatomy of this complex post-subduction system. By coupling mantle flow models with thermal69 and\nisostatic constraints, we reveal a SW-NE oriented plumelet from the Arabian foreland70 injecting into the\nsubduction system, 3D thermomechanical erosion of the upper plates and the71 active Neotethys slab\ndeformation along the southern collisional front.)\nResponse: According to the journal’s format we introduce our aim, approach, and results\nbriefly in the Introduction, respectfully we choose to retain this flow.\nL73-L74 . “Our results for the first time establish new constraints”, remove “for the first time\n“, as you have “new constraints” in the same sentence, and also because you use “for the first\ntime” L78. “73 Our results for the first time establish new constraints on drip-74 like removal beneath the\npost-orogenic magmatic arc under75 convective forces. Compared to the well-studied seismotectonics beneath\nthe Greater Caucasus 76 Range, the poorly defined tectonic implications of crustal and subcrustal earthquakes\nalong the77 Arabian collision front between the Arabian-Anatolian-Iranian plate boundary are clarified for 78\nthe first time in this study.”\nResponse: Thank you for the suggestion, it has been done.\nFig. 1c. Shouldn’t the caption be delta_Vs (or delta S-wave and not just S-wave) as you\ndisplay % variation? I am not a geophysicist, but S-wave should be in km/s. If you have % it\nlikely means that this is a variation from a reference model.\nResponse: Both versions can be accurate, we utilize the same units as in the original study-\nCiviero et al. (2023).\nL147-148. “The northern branch of the slab subducted beneath Eurasia by creating arc and\nback-arc volcanism”. Something is wrong in this sentence. I suppose you meant “by creating”\n-> “created”?\n“147 The northern148 branch of the slab subducted beneath Eurasia by creating arc and back-arc volcanism,\nmostly149 dating from the Early Cretaceous to Eocene”\nResponse: The Reviewer is correct. We have fixed this grammatical error to “The northern\nbranch of the slab subducted beneath Eurasia created back-arc volcanism, mostly dating from\nthe Early Cretaceous to Eocene.”\nL165. I don’t see the Erzurum-Kars plateau “T” on figure 1D.\n(165 the Erzurum-Kars Plateau (“T”; Fig. 1d))\nResponse: The letter 'T' in Figure 1d indicates the location of the Erzurum-Kars Plateau to\nsimplify and make clear the image. The letter T in the revised Figure 1d was bolded, and the\ntext was clarified as follows:\n“the Erzurum-Kars Plateau (locations shown by “T” in Fig. 1d)”\nL192-L195. This should go in method and be expanded (and partly moved to appendix too).\nHow did you integrate the data? Did you invert solely for temperature and you that to define\nthe thermal structure of the model?\n(192 We determine the regional mantle circulation and lithosphere deformation with a 3D\n193 instantaneous thermomechanical model based on seismically defined upper mantle structures194 (Fig. 2).\nThe results define the anatomy of the structures and dynamics that interact in a complex195 way in this active\npost-subduction system of the Arabian-Eurasian collision zone.)\nResponse: We employ a modeling method previously developed and discussed (Uluocak et\nal., 2021, 2024). We thank Reviewer 4 for their comment and confirm that the inversion was\nindeed performed only for the temperature field to define the model's thermal structure. It is\nalso important to note that the core objective of this study is not to present a novel\ncode/modeling approach, but to utilize proven and robust methodology to investigate complex\ndeformation mechanisms. In response, we have revised the manuscript to make this modeling\nstrategy more explicit for a broader audience in the Method section.\nWhat are the densities of your material? How do you choose them? Are they computed from\nphase equilibrium modelling? How sensitive is you modeling results when propagating error\nof the inversion? Are you only imposing the initial temperature structure? What are the\nassumptions in terms of rheologies and material properties for the different crustal and\nlithospheric units?\nL298. The temperature model should be presented way earlier in my opinion as this is the\nbasis of your modeling!\nResponse: As suggested by Reviewer 4, Table 1 and additional references have been added to\nthe revised manuscript, and this includes material properties, including density. Model\nsensitivity was already discussed by our previous papers cited clearly in the text (and ASPECT\nis benchmarked, as described in these references). As we do not introduce new code in this\nwork, and also because of space limitations, there is not an opportunity to repeat such\ndiscussion here. We have expanded description of the temperature model in the revised Figure\n3 (former Fig. 2); at a stage well before Figs. 5 and 6 (modelling results). The Method section\nhas been revised as follows for further description:\n-Variation in the temperature model should also be explored to test how sensitive the\nmodelling outcome is.\nResponse: Because of space and reference constraints, we could not include all previously\ntested models for this collisional zone. Our aim is to focus specifically on deformation\nmechanisms, having already explored various parameters and thermal configurations in prior\nwork. In response to the Reviewer's comment, we have expanded the revised text with\nadditional explanation concerning the modeling methodology (see response above, Method and\nResults in the tracked change file). The full suite of input parameters is detailed in the following\nfigures from our previous studies: Figure 5 in Uluocak et al. (2021), Figure 3 in Uluocak et al.\n(2024), and Figure 6 in Uluocak et al. (2021). We have cited relevant references clearly in the\ntext for detailed input parameters.\n-Do you define fault, shear zones and weak zones? How are they distributed? Although it does\nnot look like you impose such structures, how important can they be in controlling the stress-\nstrain state of the region? I would suspect that they are rather important.\nMoreover, I understand you perform instantaneous flow models. Do you account for elasticity?\nHow is isostatic equilibrium achieved? How many timesteps are performed to reach elastic and\nisostatic equilibrium?\nResponse: We did not impose any pre-existing weaknesses in the model that could cause\ndeformation, such as shear zones. Any inputs from a priori knowledge should be important in\nevolutionary models, but this is not the case here. Besides, we don’t consider elasticity,\nmass\ntransfer, and radioactive heating. The Method section has been improved to more clearly define\nour approach and modeling criteria, and to clarify possible misunderstandings raised by the\nreviewer, as follows:\nMethod: “The time stepping in model was carried out using a small Courant–Friedrichs–Lewy\nnumber of 0.1, yielding a timestep of up to ~1x105 years and the relatively small‐velocity errors\n(i.e., RMS ≤ 2 cm/year) with the effective computation time (Uluocak et al., 2021)......Any\npossible responses of elasticity, mass transfer, and radioactive heating were not considered in\nthe incompressible solution domain. Since different numerical models with various\ntemperature configurations and rheological features were tested (Uluocak et al., 2021; 2024),\nin this study we followed a convection modeling approach with inputs that have already been\nconfirmed/validated by considering large sets of geophysical data and analyses documented in\nthe Arabian-Eurasian collisional region. For a more detailed discussion of the modeling\nmethodology and parameters, the reader should refer to Uluocak et al., (2024) and Bangerth et\nal. (2022).”\nL197-L199. At that stage I am still not quite sure to understand what is a plumelet. It seems\nto me that this is simply a buoyancy contrasts driving upward motion. An upwell?\n“197 Mantle dynamics in the region are dominated by a focused rapid SW-NE directed inflow198 from the\nArabian foreland to the eastern Greater Caucasus—a plumelet— penetrating through199 the complex\nmorphology of the post-collision environment (Fig. 2).”\nResponse: A plumelet can be explained as a regional upper-mantle plume migration with\nneither a significant tail extending to the lower mantle nor a mushroom head reaching the hot\nspots on the surface as stated in the Tectonic context and motivation before Results:\nTectonic context and motivation: ... “ a plumelet (i.e., upper-mantle plume migration\nwithout significant tail and a mushroom head; Uluocak et al., 2024)”\nFig. 2. What are the iso-contours you are showing on the figure?\nResponse: The caption of revised Figure 3 has been corrected to clarify the variables.\nFigure 3.: “...Iso-surfaces are imaged only for ≤ 1000 K (see color scale in Fig. 4d).”\nFigs 3b,c Shear stress and Strain rate. Are those invariants? Second invariant the shear stress\nand strain rates? Otherwise they should be tensors…\nResponse: Agreed, we thank for the correction; the revised text has been improved by\ncorrected definition of parameters:\nRevised Fig. 4: “(b) Amplitudes of the shear stress in a slice at 20 km depth”\nMethod: “…the square root of the second invariant of the strain rate tensor…\n...and the shear stress tensor...”\nTable 1.: `… the square root of the second invariant of the strain rate tensor is shown as strain\nrates (1/s) for simplicity in the study.”\nFig. 3d. 3D temperature model with a slice. Are you plotting an iso-surface of the\ntemperature? + a slice? Moreover, there are no yellow dashed lines in Fig. 3d. It needs be to\nmoved in caption 3c.\nResponse: The caption has been corrected, thank you for the comment. In Fig. 4d, we imaged\nthe 3D temperature field for a specific depth range ([40-200] km with a temperature slice at a\nspecific depth of 200 km.\nFig. 5. Is that realistic to have temperature reaching 2000°C at 300 km depth? Also it seems\nthat large regions of very low temperature are predicted at great depth ~1000°C. Without\ncontouring I have trouble to exactly see the temperature structure, but it seems rather extreme.\nResponse: For unity, we use scientific coloring for the same features in the text. It is suggested\nthat average temperature can reach 1628 K at 100 km depths in the Anatolian Plateau (e.g.,\nArtemieva, & Shulgin, 2019, https://doi.org/10.1029/2019TC005594). These extreme\ntemperatures can be linked to lithospheric deformation mechanisms and plumelet-induced\nthermomechanical erosion, as we claim in this study. We thank the Reviewer for this comment;\nthe method has been improved with additional explanations:\nL279-L284. Where is the cold anomaly on the figure 3d. I don’t see it. Actually, the figure\nshows really well the trench parallel hot anomaly due to slab break off.\nResponse: We thank the reviewer for the interpretation although we are not sure we can\nunderstand how they reached this conclusion. The text has been corrected as follows:\nResults: “arc-parallel relatively cold anomalies (< 1500 K) beneath the Transcaucasian\nzone...”\nL288-L290. The “Sinking Britis Slab” of figure 5b?\n(288 In the active Arabian collisional front, the results discover an interaction between the\n289 plumelet and cold upper mantle structures, namely a mantle anomaly beneath the thick northern290\nArabian lithosphere (~250 km depths; Fig. 5b).)\nResponse: We don’t think we fully understood the question here; the very next sentence\nexplains the Bitlis slab anomaly. This part has been improved by the previous Reviewers’\nsuggestions (see above) and these are documented in the tracked change file.\nL301-L305. First, this slab sinking and related upwell and potential erosion are very close to\nwhat seems the boundary of the model. How is this constraining the results? Or is it because\nof how the figure is clipped?\nResponse: We aim to make visible all contrasts that might indicate deformations in the study\narea by choosing a consistent color scheme and view of figures, including revised Figure 6\n(former Fig. 5). Interpretations in the text regarding the deformation mechanisms (sinking-\ntearing-breaking of the slab, erosion etc.) do not change based on the full-sized cross-section\nshown below. We should note that this interpretation: “First, this slab sinking and related\nupwell and potential erosion are very close to what seems the boundary of the model” made by\nthe Reviewer can be accurate based on their models, but we are not suggesting such a type of\ndeformation here based on our model.\nFull-sized cross-section along P2.\nOriginal figure in the text.\nL337-L340. This is also difficult to visualize this based on figure 6. Would be better to have a\nfigure super-imposing seismicity and modeled second invariant of the strain rate for instance.\nResponse: Revised Figure 7 (previous Fig. 6) has been improved in the revised text to clarify\nlocations and seismicity by following the Reviewers’ constructive comments.\nFig. 6c Given that you are using non-linear T dependent rheology, the discretized viscosity\ncolorbar is an odd choice\nResponse: We show this image in revised Figure 7c to emphasize the deformation patterns\nsuch as sinking Bitlis Slab, active plate tearing and mantle intrusion along the cross-section.\nL349-L355. How can the southeastern propagation of the deformation be discussed in the light\nof an instantaneous flow model? I may be misunderstanding here, but from my point of view\nthis can only be explored when performing forward models.\nResponse: We appreciate the reviewer's perspective (i.e., interpretations provided in our work\ncan only be done by forward models). We respectfully note, however, that the prevailing\nevidence from instantaneous modeling studies does not support this conclusion.\nL390-L392. It seems like a bit of a circular reasoning here. As the small scale anomaly are\nthe outcome of using the tomography model as input for the model.\n(We identify small-scale cold anomalies that are also visible in independent seismic 391 tomography models\nbeneath the eastern Pontides-Lesser Caucasus arcs and the East Anatolian392 Plateau)\nResponse: This is not circular reasoning, but the result of the flow modelling (which is, in\nturn, based on the seismically-derived mantle structures) based on the data imaging upper\nmantle structures.\nL435- Where is the plumelet coming from? What drives its ascent?\nResponse: This is an important question that has been/is being considered by many researchers;\nplease see Fig. 1c and relevant explanation in the text, e.g., “In the tectonic anatomy of the\ncomplex plate boundary, during its migration from deep mantle source(s) in the south…”\nFurther discussion is beyond our scope in this study and cannot be conducted without additional\ndata, investigations, and analyses.\nRESPONSE TO REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nDear authors and editors,\nI have read carefully the new version of the manuscript and the corresponding answers and it seems\nthe authors have carefully addressed all comments and suggestions given by the reviewers, so in my\nopinion this version is acceptable for publication in Nat. Comm.\nwith kind regards\nResponse: Thank you for your constructive and insightful comments that improve our manuscript.\nBest regards.\nReviewer #2 (Remarks to the Author):\nWe are pleased to see the extensive responses of the authors to the comments / recommendations that\nwe provided for this manuscript (and to the other reviews). We are satisfied that the authors have\naddressed the points we raised in a thoughtful manner and, given the constraints of space and the need\nto respond to other reviews, we are satisfied with their improvements. We endorse the paper for\npublication if the editors see fit.\nResponse: Thank you for your contributions and comments on this manuscript. Sincerely yours.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of\nthe Nature Communications initiative to facilitate training in peer review and to provide appropriate\nrecognition for Early Career Researchers who co-review manuscripts.\nResponse: We wish you success in your career. Best regards.\nReviewer #4 (Remarks to the Author):\nReview of “Anatomy of a post-subduction collision” by Uluocak et al.\nIn the revised version of the manuscript, the authors have addressed all of my comments. The study\nrepresents a very interesting contribution, and I thank the authors for their careful and thorough\nrevision.I have only one minor comment remaining.\nRegarding Figures 6b and 6d, I would like to clarify my comment from the previous review, which may not have been\nsufficiently clear. When I wrote, “First, this slab sinking and related upwelling and potential erosion are very close to what\nseems to be the boundary of the model,” I was questioning the extent to which the lateral boundary\nconditions may influence the model behavior near the margins.\nFor instance, in Figure 6e, the flow in the top left and top right corners appears to rotate rapidly and become parallel to the\nleft and right model boundaries. My question is whether this behavior reflects a boundary effect. If so, how does this\ninfluence or constrain the interpretation of the modeled flow patterns and associated processes near the boundaries?\nNicolas Riel\nResponse: Thank you for your comment regarding the patterns in the corners of the cross-sections in\nFig. 6. We agree that modeling results are dependent on chosen inputs including boundary conditions.\nHence, to ensure our interpretations are robust and not unduly affected by boundaries or specific\nparameters used in the model, we have systematically validated/confirmed our estimates against\nmultiple independent lines of evidence at different scales. These include various seismic tomography\ndatasets, patterns of magmatism, geochemical data, and results from large-scale numerical models that\nencompass our study region. We believe this multi-disciplinary approach strengthens the reliability of\nour findings. Thank you again for your insightful comment."}]

Decision letter