# Trans-lithospheric diapirism explains the presence of ultra-high pressure rocks in the European Variscides

Petra Maierová, Karel Schulmann, Pavla Štípská, Taras Gerya, Ondrej Lexa

► To cite this version:

Petra Maierová, Karel Schulmann, Pavla Štípská, Taras Gerya, Ondrej Lexa. Trans-lithospheric diapirism explains the presence of ultra-high pressure rocks in the European Variscides. Communications Earth & Environment, 2021, 2, 10.1038/s43247-021-00122-w. insu-03708186

HAL Id: insu-03708186

https://insu.hal.science/insu-03708186

Submitted on 29 Jun 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Trans-lithospheric diapirism explains the presence of ultra-high pressure rocks in the European Variscides

Petra Maierová<fnref n="1" />✉, Karel Schulmann<fnref n="1" />,², Pavla Štípská<fnref n="1" />,², Taras Gerya ID ³ & Ondrej Lexa ID ⁴

The classical concept of collisional orogens suggests that mountain belts form as a crustal wedge between the downgoing and overriding plates. However, this orogenic style is not compatible with the presence of (ultra-)high pressure crustal and mantle rocks far from the plate interface in the Bohemian Massif of Central Europe. Here we use a comparison between geological observations and thermo-mechanical numerical models to explain their formation. We suggest that continental crust was first deeply subducted, then flowed laterally underneath the lithosphere and eventually rose in the form of large partially molten trans-lithospheric diapirs. We further show that trans-lithospheric diapirism produces a specific rock association of (ultra-)high pressure crustal and mantle rocks and ultra-potassic magmas that alternates with the less metamorphosed rocks of the upper plate. Similar rock associations have been described in other convergent zones, both modern and ancient. We speculate that trans-lithospheric diapirism could be a common process.

<footnote n="1">¹Center for Lithospheric Research, Czech Geological Survey, Prague 1, Czech Republic. ²EOST, Institute de Physique de Globe, Université de Strasbourg, Strasbourg, France. ³Institute of Geophysics, Department of Earth Science, ETH-Zurich, Zurich, Switzerland. ⁴Institute of Petrology and Structural Geology, Faculty of Science, Charles University in Prague, Prague 2, Czech Republic. ✉email: petra.maierova@geology.cz</footnote>
Since 1970s, it has been accepted that oceanic subduction beneath a continent can be followed by arrival of continental mass into the subduction zone¹ and result into collisional deformation and formation of a wide mountain belt². It was later shown that subduction of continental mantle beneath progressively accreted crust results in an asymmetric wedge-shaped architecture dominated by high-pressure (HP) metamorphism typical for Alpine type orogens³. The large Himalaya-Tibet type orogens also develop due to subduction of continental lithosphere beneath a continent, but lead to formation of a wide elevated plateau flanked by oppositely dipping thrust wedges. These large hot orogens are characterized by strong upper crust underlain by hot and mobile crust which flows sub-horizontally due to a pressure gradient⁴. This channel flow transports ductile partially molten material towards a zone of efficient erosion where the high grade rocks are finally exhumed.

The European Variscan belt is a large fossil and deeply eroded bivergent orogen characterized by high-temperature (HT) metamorphism and widespread crustal melting⁵. This led some authors to compare it to the Himalayan-Tibetan orogen⁶,⁷. However, here the channel flow exhumation scenario fails to explain the occurrence of (ultra-)high pressure (U)HP continent-derived rocks and mantle relics in the axial zone of the Variscan orogen several hundred kilometers away from the plate interface⁸. This rock association was interpreted as klippen of “root-less” nappes of unknown origin overlying the accreted crust⁹.

However, an alternative scenario was proposed¹⁰ to explain their position in the orogen based on the general architecture of the Variscan Bohemian Massif¹¹ (Figs. 1a and 2). This scenario begins with Devonian continental subduction of the Saxothuringian plate beneath Late Proterozoic-Cambrian accretionary complexes unconformably overlain by an Ordovician-Devonian basin. A subsequent Carboniferous collisional event led to the thickening of the upper plate and formation of the Moldanubian orogenic root together with the Teplá-Barrandian upper crustal lid. In the current configuration, the UHP rocks occur along the plate boundary and within the Moldanubian orogenic root. The scenario further assumes that subducted continental crust was relaminated beneath the Moho of the upper plate and finally emplaced into the crust of the orogenic root in the form of giant elongated diapiric bodies composed of UHP felsic rocks and mantle fragments. This is supported by a remarkable similarity of the protolith geochemistry and geochronology of the UHP felsic rocks with Saxothuringian Ordovician felsic volcanics and granitoids pointing to their lower-plate origin¹⁰,¹²⁻¹⁴. In the present work we further develop this scenario and suggest that the lower-plate rocks were transferred through the asthenospheric mantle and upper-plate lithosphere by a process called trans-lithospheric diapirism. Such a mechanism of exhumation was described and studied numerically¹⁵, but it was not compared with the real pressure-temperature-time (PT-time) datasets and the complexity of the exhumed mantle-crust rock associations.

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

**Fig. 1** Simplified geological maps. **a** A tectonic map of the Bohemian Massif modified after¹⁰. (U)HP felsic rocks, peridotites, mafic complexes, and ultra-potassic plutons are highlighted. Color rectangles show positions of inset maps (**b**, **c**), color circles show positions of rocks in Fig. 3. The thin dashed line shows position of the profile in Fig. 2. From the NW to the SE, the main domains are: the lower plate (Saxothuringian), the upper plate composed of the upper-crustal lid (Teplá-Barrandian) and lower-to-mid-crustal orogenic root (Moldanubian), and the continental back-stop (Brunovistulian). Positions of suture and magmatic arc are shown. **b** A part of the metamorphic complex along the plate interface⁶⁰,⁶¹,¹⁶. Spatial reference was corrected according to the General Geological Map of the Federal Republic of Germany 1:200 000 (GÜK200) available online (https://www.bgr.bund.de). Felsic granulites crop out within felsic orthogneiss and migmatite next to an accretionary complex including bodies of eclogites. Several nappes with contrasting metamorphic conditions can be distinguished. The nappes are rooted beneath the upper plate and lie on top of weakly metamorphosed rocks of the lower plate. **c** A representative granulite massif in the orogenic root domain (St. Leonhard Complex)⁶² showing a granulitic core with bodies of garnet peridotite, a mafic rim of garnet amphibolite and serpentinized peridotite surrounded mostly by medium-pressure gneiss.
<figure><img src="image_4.png" /></figure>

**Fig. 2** Interpretative NW-SE cross-section through the Bohemian Massif. The profile is representative for the 200 km-wide central part of the Bohemian Massif and the strike of the section generally coincides with the direction of Variscan subduction. The overall coaxial convergence was later disturbed by underthrusting of the easterly continental backstop towards SW at ~330 Ma, and post-collisional strike-slip motions along crustal-scale shear zones (e.g.⁶³). Structures highlight polyphase deformations related to subduction-driven processes in the western part, followed by regional cross-folding and exhumation during continental backstop underthrusting along the eastern margin. Numbers in brackets are ages (in Ma) of dominant metamorphic/deformational structures or pluton emplacement (geochronology data are summarized by ref. ¹¹ and regional studies²⁸,²⁹,³¹,⁶⁴⁻⁶⁶). Extrapolation of surface geology is constrained by seismic reflection and refraction studies (solid and dashed red lines)⁶⁷,⁶⁸. The extents of the lower plate (laminated lower crust, inclined reflectors tracking the plate interface) and the continental backstop (seismic velocity gradient in the lower crust) are visible. The subsurface extent of the continental backstop was also determined from the gravity data.

In the next sections, we first review the geology of the (U)HP felsic rocks and their rock associations. Then we introduce our numerical models of trans-lithospheric diapirism and describe their typical features. Finally, we discuss the similarities and discrepancies between the natural observations and the models, and outline possible implications.

## Geological evidence for trans-lithospheric diapirism

In the Bohemian Massif, we distinguish two contrasting domains containing UHP felsic rocks:¹⁰ (i) the plate interface and (ii) the region between the magmatic arc and the continental backstop (the orogenic root).

**Plate interface association.** The interface between the lower and upper plate is characterized by schist and orthogneiss nappes underlying HP granulite and migmatitic orthogneiss units all derived from the subducted lower plate (Fig. 1b; PT paths and ages in Fig. 3, blue color). The mineral assemblage and mineral chemistry of the granulites commonly record HP metamorphism at ~800 °C and ~20 kbar¹⁶, but locally found UHP minerals like coesite and micro-diamond point to pressures of 40 kbar or more¹⁷. Lower plate schist and orthogneiss record conditions of ~650 °C at ~20 kbar¹⁸ and ~700 °C at 13 kbar¹⁶. All these rocks are characterized by cooling during exhumation.

**Orogenic root association.** (U)HP granulite bodies in the orogenic root (Fig. 1c; PT paths and ages in Fig. 3, green, orange, and red colors) are commonly rimmed by mafic complexes containing meta-gabbros, retrogressed eclogites, and peridotites surrounded by medium-pressure schists and gneisses (7–10 kbar¹⁹,²⁰). These bodies are typically accompanied by ultra-potassic plutons that formed by mixing of metasomatized mantle- and crust-derived melts¹². The (U)HP granulites contain numerous xenoliths of garnet and spinel peridotites with pressures estimated to 20–40 kbar⁸. The metamorphic conditions of the hosting granulites are ~800–1000 °C at ~20 kbar²¹⁻²⁴ but the presence of coesite and micro-diamond indicates early UHP conditions²⁵,²⁶. Granulites in the vicinity of the magmatic arc, i.e., closer to the plate interface, show lower and/or short lived peak temperature conditions compared to those near the backstop²⁷. The HP-HT conditions are typically followed by nearly isothermal decompression²²,²⁸,²³.

In both regions, modern U-Pb petro-chronological studies on zircon and monazite dated the prograde PT path to ~360–350 Ma and peak temperature conditions followed by retrogression in the middle crust to 345–338 Ma²⁹⁻³².

Altogether, the plate interface association reflects metamorphic evolution of continental rocks that return from different levels of the subduction interface. In contrast, the orogenic root association records mixing of the lower plate-derived rocks with the mantle at (U)HP conditions. This rock association and mantle-derived magmas were emplaced into the upper-plate crust. This view is further supported by the current structure of the crust as seen in deformation fabrics and seismic studies (Fig. 2).

## Magmatic-thermal-mechanical model of trans-lithospheric diapirs

We ran a series of 2D numerical models of oceanic and continental subduction and compared them with the data from the Bohemian Massif. In most of the presented models, trans-lithospheric diapirs develop. This is because we explore a limited part of the parameter space where the diapric exhumation is the preferred evolutionary style¹⁵. In the following text we describe the typical model evolution exemplified by model 8 (Fig. 4a, PT paths in Fig. 4b). In addition, we use other model cases to illustrate relationships between different model characteristics (input and output parameters). More details about the model setup and its limitations are in Methods. Full descriptions of the individual model cases are in Supplementary Discussion.

In the first stage, the models show steeply dipping oceanic subduction accompanied by the formation of a sedimentary wedge, hydration and melting of the mantle, and melting of the oceanic crust and sediments. After the closure of the oceanic domain, the continental margin is subducted. At ~100–180 km depth (after ~7 Myr of continental subduction), the felsic continental crust weakens sufficiently to separate from the rest of the plate. The felsic material returns back along the subduction zone, and accumulates at a depth of ~60–90 km, i.e., below the upper-plate lithosphere (model 8, 15 Myr, Fig. 4a). It is buoyant with respect to the surrounding mantle, but its motion is restricted by the overlying high-viscosity lithosphere. It flows laterally towards the retro-lithosphere forming a sheet that is up to several hundred kilometers long.

During this flow, the crustal material can mix with the sub-lithospheric mantle, water is released from the subducted
<figure><img src="image_5.png" /></figure>

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

**Fig. 3** Pressure-temperature-time data. **a** Selected pressure-temperature paths of (U)HP crustal rocks from the Bohemian Massif. Paths from the plate interface association (blue)¹⁶⁻¹⁸,⁶⁹,⁷⁰ from the orogenic root association near the magmatic arc (green)¹⁴,²⁵,²⁷,³⁰ from the central part of the root domain (orange)²²,²³ and near the continental backstop (red)¹⁹⁻²¹,²⁴,²⁶,²⁸,⁷¹,⁷². Locations of the studied rocks are shown in Fig. 1. For the plate interface association (blue) a large spread of peak temperatures and pressures is typical. For the orogenic root association (red), high temperatures at ~15–20 kbar are typical showing a significant deviation from the cold subduction environment. **b** Metamorphic ages show a similar pattern for all domains: ages of ~360–355 Ma for the prograde path or early UHP metamorphism²⁸,³⁰⁻³² ~340 Ma for the main HT metamorphic event¹⁴,²⁹,³² and ~340–335 Ma for retrogression in the middle crust²⁸,²⁹,³¹. An exception is the age of ~340 Ma reported for all the three phases in granulites near the plate interface⁷³. More details on the PT-time data are available in Supplementary Data 1.

sediments leading to hydration of the nearby mantle and its melting. Mixing is particularly efficient in models with a high viscosity of the felsic crust (model 3, subducted felsic crust has the viscosity of felsic granulite). The subducted felsic material partially melts as a result of progressive heating in contact with the mantle. Eventually, after ~10 Myr of continental subduction, it is abruptly exhumed as a diapir through the mantle lithosphere weakened by rising melt, and through the lower crust into the middle crust (model 8, 20 Myr, Fig. 4a). The resulting exhumed body typically contains sediments, pieces of oceanic crust and mantle, and is surrounded by mantle- and lower-crustal envelopes. In some models, the subducted felsic material is also exhumed along the plate interface forming a subduction channel (model 8, 26 Myr, Fig. 4a). The flow in the channel can be active

in several pulses and contain sediments and rare pieces of oceanic crust and mantle.

In the diapir, the PT conditions are rather homogeneous. The pressure during the subduction stage reaches 30–50 kbar corresponding to the maximum depth of subduction. The maximum pressure correlates with the viscosity of the subducting crust (compare PT paths in models 1, 2, and 3, Fig. 4c). The subsequent lateral flow in the mantle takes place at pressures around 20–25 kbar depending on the thickness of the lithosphere. The temperature increases during this process due to gradual thermal equilibration of the subducted crust with the surrounding mantle. The maximum temperature therefore depends mainly on duration of this flow. If the diapir exhumes soon after continental subduction, the peak temperatures are relatively low (~800 °C) (Fig. 4c, model 15). In the models where the sub-lithospheric flow is active for more than ~15 Myr, the material usually exhumes far from the trench and the peak temperature can reach 1000 °C (Fig. 4c, model 2). The decompression path is nearly isothermal (adiabatic) due to quick exhumation.

In the subduction channel, the peak temperature reaches ~600–900 °C at maximum, but can be as low as 500 °C in a similar final position in the same model; the peak pressures are heterogeneous as well (Fig. 4b–c, blue PT paths). In contrast to the diapirs, cooling is more pronounced during exhumation along the plate interface. The maximum temperature above 800 °C in the channel is rare in our models. Such a high temperature is attained only by the material that was in contact with the hot mantle for a prolonged time (models 9 and 21; Supplementary Movies 9 and 21).

## Discussion

**Comparison with the Bohemian Massif.** The stages of the model evolution correlate well with the metamorphic conditions of the Bohemian (U)HP rocks. The modeled deep continental subduction explains the occurrence of relict UHP minerals (Fig. 5a) and the conditions during the subsequent sub-lithospheric flow match the ubiquitous HP–HT metamorphism (Fig. 5b). Most of the published PT paths start at 16–20 kbar where the high temperature caused major chemical re-equilibration and erased the early history. The PT path along which the rocks got to this point is extremely uncertain and rarely studied. Chemical zoning in garnet with diamond and coesite inclusions suggests relatively low temperatures for UHP peak and early exhumation, followed by heating²⁵,²⁷. Such a PT path correlates with heating during the lateral flow stage in models 1, 2, and 22 (Fig. 4c). Another detailed study of coesite- and diamond-bearing rocks suggests that UHT was already attained at UHP, followed by almost isothermal decompression¹⁷, comparable to model 3 (Fig. 4c).

The trans-lithospheric diapiric exhumation explains both PT-time characteristics and rock associations in the root domain: high peak temperatures (800–1000 °C), isothermal decompression from HP conditions, and eclogite and garnet peridotite along the margins and inside the granulite massifs (Fig. 5c, right part). Conversely, the modeled subduction channel mimics formation of the metamorphic complex along the plate interface showing heterogeneous PT conditions, lower peak temperatures and cooling during decompression (Fig. 5c, left part). Temporal evolution of the model mimics the time lag of 10–20 Myr between subduction (~360–350 Ma) and exhumation (~345–335 Ma) of the granulites. In the model, the shorter time lag corresponds to lower peak temperatures (model 15, shorter sub-lithospheric flow, exhumation relatively near the trench). Longer residence of the lower plate rocks in the mantle (model 2, exhumation far from the trench) leads to higher peak temperatures. Duration of the flow is influenced by the solidus temperature of the material. The
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Results of numerical modeling. a Time evolution of model 8: continental subduction and sub-lithospheric flow (snapshot at 15 Myr from the beginning of continental subduction), exhumation in a diapir far from the trench (20 Myr), and final architecture including exhumation along the plate interface (26 Myr). b PT paths in model 8 for representative samples (stars and circles in the panels with the time evolution). c Final architecture in other model cases (for discussion see text). Small inset panels show typical PT paths in the subducted felsic crust (red for diapirism, blue for exhumation along the plate interface) compared to the PT paths in model 1 (gray). Color coding of materials: yellow - felsic crust of lower plate, orange - sediments, brown - felsic crust of upper plate, gray - mafic crust, green - oceanic crust, blue - mantle, light blue - hydrated mantle, pink - molten felsic crust, violet - quenched felsic crust, red - molten mantle. Tics on the xy-axes are every 100 km. Times are given from the beginning of continental subduction. The time evolutions of the individual model cases are captured in Supplementary Movies 1-23.</figcaption></figure>

lower solidus (model 15) leads to earlier melting, melt-induced weakening of the upper plate, and exhumation. Mixing with the mantle and formation of small diapirs can further increase the peak temperature (model 3, model 10).

A more detailed look on the PT paths reveals several discrepancies between the models and the data. The models predict somewhat lower temperatures than observed, especially in the subduction channel. A similar problem was previously reported in other numerical studies of subduction³³. Another discrepancy is in the depth of diapir emplacement. Observed PT

paths show that the level of emplacement was in the middle-lower crust, although at least some of the HP rocks reached the surface and were eroded and deposited already at ~330 Ma³⁴. In contrast, the diapirs in the model exhume very quickly along vertical weak zones and in most cases reach shallow crustal depths. There are several factors that may influence the level of emplacement. According to previous studies of magma intrusion, the depth of emplacement crucially depends on the viscosity of the crust: intrusion into the lower crust occurs only if it has a relatively low-viscosity and shows a ductile behavior³⁵.
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Proposed mechanism of formation of (U)HP felsic rocks.</figcaption></figure>

a Continental subduction, b sub-lithospheric flow and formation of mixed crust- and mantle-derived magma, c growth of a trans-lithospheric diapir and exhumation along the plate interface. Stars show positions of samples whose PT-time evolution is depicted in the inset panel. Red path shows long residence at 20–25 kbar and quick exhumation to the crust of the upper plate. Blue paths show subduction to various depths with possible residence beneath the lithosphere (dark blue path), followed by exhumation along the plate interface.

This implies that the middle-lower crust of the upper plate is probably too strong in our models. If we modify the stratification of the upper plate and prescribe a thicker felsic crust at the expense of the strong mafic lower crust, a significant part of the diapir remains at the mid-crustal level (Fig. 4c, model 22). In line with that, the PT data along the eastern margin of the orogenic root show thickening (6–9 kbar) and warming (600–700 °C) of the upper-plate middle crust prior to exhumation of the granulites²⁰. Another factor is the viscosity of the diapirs, where the diapirs with the strong rheology (Fig. 4c, model 3) remain trapped at deep crustal levels. Finally, the fast vertical exhumation in the model is mainly the result of the prescribed melt-induced weakening. In nature, we can expect a more complicated path, more efficient cooling, and a deeper level of emplacement.

The Bohemian Massif architecture is three dimensional and shows more complexities than our models. Properties of the plates that vary along the plate interface cause deviations from two-dimensionality. Besides that, the buoyant material separated from the downgoing plate may form columnar rather than sheet-like upwellings in 3D³⁶. As a result, the sub-lithospheric flow might be irregular or laterally discontinuous, and the diapirs could have complex shapes³⁷ varying between ellipsoidal massifs and linear zones parallel to the plate interface as observed in nature.

The implementation of (de)hydration, melting, melt percolation, and melt-induced weakening in the model (Methods), although highly simplified, allows us to resolve the inflow of crustal and sedimentary material underneath the mantle lithosphere, its mixing with the mantle, and hydration of the

surrounding mantle (Fig. 5b). These processes mimic formation of mélange of HP granulites, eclogites and peridotites, and associated mantle metasomatism. Subsequent melting of such a heterogeneous source, as seen in the model, would be a straightforward way to form the ultra-potassic granitoids and explain their association with the granulite massifs.

Our results indicate that the classical concepts of the supramantle collisional wedge and exhumation along the plate interface are insufficient to describe formation of the (U)HP rock assemblages in the Bohemian Massif. While the data from its western part are compatible with the classical concept, more stages are needed to explain the key characteristics of the orogenic root. There, reworking by addition of material from the lower plate and from the mantle is a first-order feature: the HP granulites including peridotites form ~6 percent of surface exposure and adjacent HP gneiss and ultra-potassic plutons another ~25 percent. In addition, the rapid vertical transfer of hot diapirs affected the thermal state of the upper-plate middle crust significantly, as evidenced by massive migmatization of metasediments close to the granulite massifs. Besides efficient advection of heat from the mantle, the diapirs and ultra-potassic plutons added a large amount of radiogenic elements into the crust thereby increasing its long-term heat production, as witnessed by voluminous post-orogenic magmatism³⁸.

**Ancient and modern mountain belts perspective.** Our work suggests that the continental subduction-to-exhumation process can lead to formation of two specific rock associations characterized by distinct PT histories. The first type is commonly developed along the plate interface and is well known to mark suture zones in a range of orogens. The second, orogenic root type, is characterized by a unique crust–mantle interaction and a specific metamorphic and magmatic history. Besides the Bohemian Massif, it was reported from the Schwarzwald in Germany³⁹, the Vosges in Eastern France⁴⁰, the Monts du Lyonnais in the French Massif Central⁴¹ and in the Variscan Ulten zone in the Alps⁴². In most of these European Variscan regions, the ages of the metamorphic peak, potassic magmatism, and exhumation correspond to those obtained in the Bohemian Massif⁴³ and modeled in this work. If we assume that the orogenic root association described above is indicative of trans-lithospheric diapirism then a significant part of the European Variscan orogen was shaped by this newly recognized crustal recycling process.

We can further speculate that a similar process operated also in other orogens. In the archetypal Tibetan-Himalayan system, the rock association characteristic for the trans-lithospheric diapirism is incomplete. In the western Himalaya, occurrences of (U)HP crustal xenoliths in ultra-potassic volcanics point to deep subduction of the continental crust⁴⁴, but an incorporation of the mantle into the UHP crustal rocks was not described there. In southern Tibet, interaction of crustal and mantle rocks is recorded by felsic HP granulite and ultramafic xenoliths entrapped in (ultra-)potassic magmas⁴⁵. In southern and eastern Tibet, the ultra-potassic magmatism was interpreted to result from hydration and metasomatism of the mantle wedge above deeply subducted continental crust⁴⁶,⁴⁷ similar as in our models. Further south the UHP rocks of Tso Morari and Kaghan units⁴⁸⁻⁵⁰ are developed entirely within strata becoming to the lower Indian plate and may therefore correspond to the plate interface association described in the Bohemian Massif. Altogether, these indices suggest that the common “supra-mantle” view of the Tibetan-Himalayan orogen should be reconsidered.

In addition, previous studies showed that trans-lithospheric diapirism can operate during oceanic subduction. In the Early
Paleozoic Qilian orogen, the rock association of UHP gneiss, eclogites, and garnet peridotites was interpreted to result from diapiric uplift of a mélange complex above an oceanic plate⁵¹. In the Neogene D'Entrecasteaux islands, gneiss domes containing eclogites and surrounded by amphibolites and ultramafics formed as a result of diapirism above a subducted continental margin through the (ultra)mafic crust⁵². These examples show that trans-lithospheric diapirs may form in various tectonic settings where buoyant material is efficiently subducted.

## Methods

We use a two-dimensional thermo-mechanical model based on the software I2VIS⁵³. It solves the equations of continuity, motion, and heat transport using a finite volumes scheme on a fully staggered computational grid. It adopts visco-plastic temperature-dependent rheology, percolation of fluids with constant fluid-to-matrix velocity, melting, melt extraction, and melt emplacement. Heterogeneous material composition is treated via the particle-in-cell method. The topmost part of the domain represents the air and water and it is treated as a "sticky air"⁵⁴ that mimics the free-surface condition at the top of the rocky earth surface and deformation of this surface. The surface is subject to altitude-dependent erosion and sedimentation. The model setup is detailed in the Supplementary Methods, the initial setup is depicted in Supplementary Fig. 1, input parameters of the model, material properties, and flow-law parameters are in Supplementary Tables 1-3.

Our model builds on the work by Maierová et al.¹⁵, who studied the fate of the crust during subduction of continental lithosphere. They showed that trans-lithospheric diapirs ("exhumation near the arc or in the back-arc" in¹⁵), which are the main focus of the present paper, form in models with a thin upper-plate lithosphere, quick subduction of the lower plate, and young oceanic plate. In the case of the Bohemian Massif, the thin upper plate is justified by its extensional setting prior to and during the Devonian¹¹. The rate of subduction and the ocean age are only weakly constrained. In our models, the oceanic subduction is only a preparatory stage for the continental subduction. We therefore choose the age of the ocean and velocity of subduction in such a manner that the bending of the plate due to slab pull is approximately balanced by the convergence of the plates and therefore no significant trench retreat occurs. In most models, we prescribe the initial age of the oceanic plate of 20 Myr (age at the beginning of continental subduction ~40 Myr) and the velocity of 4 cmyr⁻¹. The choice of the high convergence velocity results into low-temperature gradient with depth along the subducting plate.

Diapirs can exhume only through a sufficiently weak material. In the model, the most important hypothesized source of weakening is the melt rising from the subducted felsic crust. This weakening is prescribed in a simple manner: during melt extraction, the column of material above the site of extraction is weakened by a constant factor⁵⁵. The melt weakening hypothesis is justified by both analytical derivations and numerical models of lithospheric deformation assisted by frequent dyke propagation⁵⁵. It is further supported by geodynamic models that show a link between partial melting of the mantle and weakening of the lithosphere⁵⁶. It is also an essential physical mechanism for reproducing broad variability of Venus coronae shapes formed by tectono-magmatic plume-lithosphere interactions⁵⁷. In most models, the weakening factor is prescribed to $10^{-3}$, which may in some cases overestimate the effect of melt-induced weakening but is within the range of $10^{-4}$–$10^{-2}$ estimated numerically for deforming melt-bearing rocks (cf. discussions in ref. ⁵⁵, Methods). In nature, the melt-induced weakening may have been supplemented by other sources of weakening neglected here: addition of fluids into the mantle wedge during prolonged oceanic subduction, thermal weakening due to back-arc extension, and inherited zones of weakness in the upper plate.

We performed a parametric study, where we focused mainly on the properties of the converging continental plates. We tested the effect of the solidus temperature and density of the subducting felsic crust, thickness of the felsic crust, density of sediments and felsic melt, and parameters influencing the strength of the upper plate. Importantly, we varied the flow law for the subducting felsic crust (wet quartz, dry quartz, felsic granulite) to demonstrate the effect of its weak or strong rheology on the subduction–exhumation process. The natural (micro)structure of the HP felsic rocks, however, shows that their behavior was more complicated than pure dislocation creep. These rocks record a switch to deformation by grain-boundary sliding that was associated with significant weakening during cooling at ~800 °C, and subsequent hardening at ~700 °C⁵⁸,⁵⁹. Such a complexity is beyond the scope of this study, which imposes certain limitations on the comparison between the models and the Bohemian Massif.

It should be noted that the present model is not tuned to mimic formation of the Bohemian Massif. It is rather a model of trans-lithospheric diapirism that shows several robust features that can be compared to natural examples. The model is non-linear, has multiple stages and evolution during each of these stages influences the subsequent stages and the final properties of the model. For that reason, correlation of input and output characteristics is difficult, and we mostly focus on a comparison of different output characteristics (distance of diapirs from the trench, duration of sub-lithospheric flow, PT conditions). More details about different model cases are in Supplementary Table 4 and Supplementary Discussion.

## Data availability

No datasets were generated during the current study. The pressure–temperature–time data supporting the findings of this study are available in Supplementary Data 1. The time evolutions of the presented models are in Supplementary Movies 1–23.

## Code availability

The numerical modeling codes associated with this paper are available from the authors upon reasonable request.

Received: 6 May 2020; Accepted: 5 February 2021;
Published online: 11 March 2021

## References

1. Wilson, J. T. Did the Atlantic close and then re-open? *Nature* **211**, 676–681 (1966).

2. Dewey, J. F. & Bird, J. M. Mountain belts and the new global tectonics. *J. Geophys. Res.* **75**, 2625–2647 (1970).

3. Platt, J. P. Dynamics of orogenic wedges and the uplift of high-pressure metamorphic rocks. *Geol. Soc. Am. Bull.* **97**, 1037–1053 (1986).

4. Beaumont, C., Jamieson, R. A., Nguyen, M. H. & Lee, B. Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation. *Nature* **414**, 738–742 (2001).

5. Vanderhaeghe, O. & Teysier, C. Partial melting and flow of orogens. *Tectonophysics* **342**, 451–472 (2001).

6. Dewey, J. F. & Burke, K. C. Tibetan, Variscan, and Precambrian basement reactivation: products of continental collision. *J. Geol.* **81**, 683–692 (1973).

7. Maierová, P. et al. European Variscan orogenic evolution as an analogue of Tibetan–Himalayan orogen: insights from petrology and numerical modeling. *Tectonics* **35**, 1760–1780 (2016).

8. Medaris, L. G.Jr., Beard, B. L. & Jelínek, E. Mantle-derived, UHP garnet pyroxenite and eclogite in the Moldanubian Gföhl Nappe, Bohemian Massif: a geochemical review, new PT determinations, and tectonic interpretation. *Int. Geol. Rev.* **48**, 765–777 (2006).

9. Matte, P. The Variscan collage and orogeny (480–290 Ma) and the tectonic definition of the Armorica microplate: a review. *Terra Nova* **13**, 122–128 (2001).

10. Schulmann, K., Lexa, O., Janoušek, V., Lardeaux, J. M. & Edel, J. B. Anatomy of a diffuse cryptic suture zone: an example from the Bohemian Massif, European Variscides. *Geology* **42**, 275–278 (2014).

11. Schulmann, K. et al. An Andean type Palaeozoic convergence in the Bohemian massif. *Comptes Rendus Geosci.* **341**, 266–286 (2009).

12. Janoušek, V. & Holub, F. V. The causal link between HP-HT metamorphism and ultrapotassic magmatism in collisional orogens: case study from the Moldanubian Zone of the Bohemian Massif. *Proc. Geol. Assoc.* **118**, 75–86 (2007).

13. Lexa, O. et al. Heat sources and trigger mechanisms of exhumation of HP granulites in Variscan orogenic root. *J. Metamorph. Geol.* **29**, 79–102 (2011).

14. Nahodilová, R. et al. Exhumation of subducted continental crust along the arc region. *Gondwana Res.* **80**, 157–187 (2020).

15. Maierová, P., Schulmann, K. & Gerya, T. Relamination styles in collisional orogens. *Tectonics* **37**, 224–250 (2018).

16. Závada, P. et al. Role of strain localization and melt flow on exhumation of deeply subducted continental crust. *Lithosphere* **10**, 217–238 (2018).

17. Haifler, J. & Kotková, J. UHP–UHT peak conditions and near-adiabatic exhumation path of diamond-bearing garnet–clinopyroxene rocks from the Eger Crystalline Complex, North Bohemian Massif. *Lithos* **248**, 366–381 (2016).

18. Konopásek, J. Eclogitic micaschists in the central part of the Krušné hory Mountains (Bohemian Massif). *Eur. J. Miner.* **13**, 87–100 (2001).

19. Štípská, P., Schulmann, K. & Powell, R. Contrasting metamorphic histories of lenses of high-pressure rocks and host migmatites with a flat orogenic fabric (Bohemian Massif, Czech Republic): a result of tectonic mixing within horizontal crustal flow? *J. Metamorph. Geol.* **26**, 623–646 (2008).

20. Racek, M., Štípská, P., Pitra, P., Schulmann, K. & Lexa, O. Metamorphic record of burial and exhumation of orogenic lower and middle crust: a new tectonothermal model for the Drosendorf window (Bohemian Massif, Austria). *Miner. Petrol.* **86**, 221–251 (2006).

21. Carswell, D. A. & O'Brien, P. J. Thermobarometry and geotectonic significance of high-pressure granulites: examples from the Moldanubian Zone of the Bohemian Massif in Lower Austria. *J. Petrol.* **34**, 427–459 (1993).

22. Franěk, J., Schulmann, K., Lexa, O., Tomek, Č. & Edel, J. B. Model of syn-convergent extrusion of orogenic lower crust in the core of the Variscan belt: implications for exhumation of high-pressure rocks in large hot orogens. *J. Metamorph. Geol.* **29**, 53–78 (2011).
23. Štípská, P., Powell, R., Racek, M. & Lexa, O. Intermediate granulite produced by transformation of eclogite at a felsic granulite contact, in Blanský les, Bohemian Massif. *J. Metamorph. Geol.* **32**, 347–370 (2014).

24. Tajčmanová, L., Konopásek, J. & Schulmann, K. Thermal evolution of the orogenic lower crust during exhumation within a thinned Moldanubian root of the Variscan belt of Central Europe. *J. Metamorph. Geol.* **24**, 119–134 (2006).

25. Perraki, M. & Faryad, S. W. First finding of microdiamond, coesite and other UHP phases in felsic granulites in the Moldanubian Zone: implications for deep subduction and a revised geodynamic model for Variscan orogeny in the Bohemian Massif. *Lithos* **202**, 157–166 (2014).

26. Kobayashi, T., Hirajima, T., Hiroi, Y. & Svojtka, M. Determination of SiO₂ Raman spectrum indicating the transformation from coesite to quartz in Gföhl migmatitic gneisses in the Moldanubian Zone, Czech Republic. *J. Miner. Petrol. Sci.* **103**, 105–111 (2008).

27. Jedlička, R., Faryad, S. W. & Hauzenberger, C. Prograde metamorphic history of UHP granulites from the Moldanubian Zone (Bohemian Massif) revealed by major element and Y+ REE zoning in garnets. *J. Petrol.* **56**, 2069–2088 (2015).

28. Tajčmanová, L., Soejono, I., Konopásek, J., Košler, J. & Klötzli, U. Structural position of high-pressure felsic to intermediate granulites from NE Moldanubian domain (Bohemian Massif). *J. Geol. Soc.* **167**, 329–345 (2010).

29. Friedl, G., Cooke, R. A., Finger, F., McNaughton, N. J. & Fletcher, I. R. Timing of Variscan HP-HT metamorphism in the Moldanubian Zone of the Bohemian Massif: U-Pb SHRIMP dating on multiply zoned zircons from a granulite from the Dunkelsteiner Wald Massif, Lower Austria. *Miner. Petrol.* **102**, 63 (2011).

30. Nahodilová, R., Štípská, P., Powell, R., Košler, J. & Racek, M. High-Ti muscovite as a prograde relic in high pressure granulites with metamorphic Devonian zircon ages (Běstvina granulite body, Bohemian Massif): consequences for the reclamation model of subducted crust. *Gondwana Res.* **25**, 630–648 (2014).

31. Štípská, P., Powell, R., Hacker, B. R., Holder, R. & Kylander-Clark, A. R. C. Uncoupled U/Pb and REE response in zircon during the transformation of eclogite to mafic and intermediate granulite (Blanský les, Bohemian Massif). *J. Metamorph. Geol.* **34**, 551–572 (2016).

32. Závada, P. et al. Monazite geochronology in melt-percolated UHP meta-granitoids: an example from the Erzgebirge continental subduction wedge, Bohemian Massif. *Chem. Geol.* **559**, 119919 (2021).

33. Penniston-Dorland, S. C., Kohn, M. J. & Manning, C. E. The global range of subduction zone thermal structures from exhumed blueschists and eclogites: rocks are hotter than models. *Earth Planet. Sci. Lett.* **428**, 243–254 (2015).

34. Hartley, A. J. & Otava, J. Sediment provenance and dispersal in a deep marine foreland basin: the Lower Carboniferous Culm Basin, Czech Republic. *J. Geol. Soc.* **158**, 137–150 (2001).

35. Gerya, T. V. & Burg, J. P. Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation. *Phys. Earth Planet. Inter.* **160**, 124–142 (2007).

36. Hasenclever, J., Morgan, J. P., Hort, M. & Rüpke, L. H. 2D and 3D numerical models on compositionally buoyant diapirs in the mantle wedge. *Earth Planet. Sci. Lett.* **311**, 53–68 (2011).

37. Kaus, B. J. & Podladchikov, Y. Y. Forward and reverse modeling of the three-dimensional viscous Rayleigh-Taylor instability. *Geophys. Res. Lett.* **28**, 1095–1098 (2001).

38. Finger, F., Roberts, M. P., Haunschmid, B., Schermaier, A. & Steyrer, H. P. Variscan granitoids of central Europe: their typology, potential sources and tectonothermal relations. *Miner. Petrol.* **61**, 67–96 (1997).

39. Marschall, H. R., Kalt, A. & Hanel, M. P-T evolution of a Variscan lower-crustal segment: a study of granulites from the Schwarzwald, Germany. *J. Petrol.* **44**, 227–253 (2003).

40. Skrzypek, E. et al. The significance of Late Devonian ophiolites in the Variscan orogen: a record from the Vosges Klippen Belt. *Int. J. Earth Sci.* **101**, 951–972 (2012).

41. Gardien, V., Tegyey, M., Lardeaux, J. M., Misseri, M. & Dufour, E. Crust-mantle relationships in the French Variscan chain: the example of the Southern Monts du Lyonnais unit (eastern French Massif Central). *J. Metamorph. Geol.* **8**, 477–492 (1990).

42. Langone, A., Braga, R., Massonne, H. J. & Tiepolo, M. Preservation of old (prograde metamorphic) U-Th-Pb ages in unshielded monazite from the high-pressure paragneisses of the Variscan Ulten Zone (Italy). *Lithos* **127**, 68–85 (2011).

43. Lardeaux, J. M. et al. The Moldanubian zone in the French Massif Central, Vosges/Schwarzwald and Bohemian Massif revisited: differences and similarities. *Geol. Soc. Lond. Spec. Publ.* **405**, 7–44 (2014).

44. Hacker, B. et al. Near-ultrahigh pressure processing of continental crust: Miocene crustal xenoliths from the Pamir. *J. Petrol.* **46**, 1661–1687 (2005).

45. Chan, G. N. et al. Probing the basement of southern Tibet: evidence from crustal xenoliths entrained in a Miocene ultrapotassic dyke. *J. Geol. Soc.* **166**, 45–52 (2009).

46. Tian, S. H. et al. Subduction of the Indian lower crust beneath southern Tibet revealed by the post-collisional potassic and ultrapotassic rocks in SW Tibet. *Gondwana Res.* **41**, 29–50 (2017).

47. Campbell, I. H. et al. The origin of shoshonites: new insights from the Tertiary high-potassium intrusions of eastern Tibet. *Contrib. Miner. Petrol.* **167**, 983 (2014).

48. de Sigoyer, J. et al. Dating the Indian continental subduction and collisional thickening in the northwest Himalaya: multichronology of the Tso Morari eclogites. *Geology* **28**, 487–490 (2000).

49. Parrish, R. R., Gough, S. J., Searle, M. P. & Waters, D. J. Plate velocity exhumation of ultrahigh-pressure eclogites in the Pakistan Himalaya. *Geology* **34**, 989–992 (2006).

50. Wilke, F. D. et al. The multistage exhumation history of the Khanag Valley UHP series, NW Himalaya, Pakistan from U-Pb and 40Ar/39Ar ages. *Eur. J. Miner.* **22**, 703–719 (2010).

51. Yin, A. et al. Early Paleozoic tectonic and thermomechanical evolution of ultrahigh-pressure (UHP) metamorphic rocks in the northern Tibetan Plateau, northwest China. *Int. Geol. Rev.* **49**, 681–716 (2007).

52. Little, T. A. et al. Diapiric exhumation of Earth's youngest (UHP) eclogites in the gneiss domes of the D'Entrecasteaux Islands, Papua New Guinea. *Tectonophysics* **510**, 39–68 (2011).

53. Gerya, T. V. & Yuen, D. A. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. *Phys. Earth Planet. Inter.* **140**, 293–318 (2003).

54. Crameri, F. et al. A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the 'sticky air' method. *Geophys. J. Int.* **189**, 38–54 (2012).

55. Gerya, T. V., Stern, R. J., Baes, M., Sobolev, S. V. & Whattam, S. A. Plate tectonics on the Earth triggered by plume-induced subduction initiation. *Nature* **527**, 221–225 (2015).

56. Bahadori, A. & Holt, W. E. Geodynamic evolution of southwestern North America since the Late Eocene. *Nat. Commun.* **10**, 1–18 (2019).

57. Gürcher, A. J., Gerya, T. V., Montési, L. G. & Munch, J. Corona structures driven by plume-lithosphere interactions and evidence for ongoing plume activity on Venus. *Nat. Geosci.* **13**, 547–554 (2020).

58. Franěk, J. et al. Origin of felsic granulite microstructure by heterogeneous decomposition of alkali feldspar and extreme weakening of orogenic lower crust during the Variscan orogeny. *J. Metamorph. Geol.* **29**, 103–130 (2011).

59. Maierová, P., Lexa, O., Jeřábek, P., Schulmann, K. & Franěk, J. Computational study of deformation mechanisms and grain size evolution in granulites-implications for the rheology of the lower crust. *Earth Planet. Sci. Lett.* **466**, 91–102 (2017).

60. Mlčoch, B. & Konopásek, J. Pre-Late Carboniferous geology along the contact of the Saxothuringian and Teplá-Barrandian zones in the area covered by younger sediments and volcanics (western Bohemian Massif, Czech Republic). *J. Geosci.* **55**, 81–94 (2010).

61. Sebastian, U. *Die Geologie des Erzgebirges*. (Springer Spektrum, 2013).

62. Racek, M., Štípská, P. & Powell, R. Garnet-clinopyroxene intermediate granulites in the St. Leonhard massif of the Bohemian Massif: ultrahigh-temperature metamorphism at high pressure or not? *J. Metamorph. Geol.* **26**, 253–271 (2008).

63. Edel, J. B., Schulmann, K., Lexa, O. & Lardeaux, J. M. Late Palaeozoic palaeomagnetic and tectonic constraints for amalgamation of Pangea supercontinent in the European Variscan belt. *Earth-Sci. Rev.* **177**, 589–612 (2018).

64. Collett, S. et al. Combined Lu-Hf and Sm-Nd geochronology of the Maríánské Lázně Complex: new constraints on the timing of eclogite-and granulite-facies metamorphism. *Lithos* **304**, 74–94 (2018).

65. Peřestý, V., Lexa, O. & Jeřábek, P. Restoration of early-Variscan structures exposed along the Teplá shear zone in the Bohemian Massif: constraints from kinematic modelling. *Int. J. Earth Sci.* **109**, 1189–1211 (2020).

66. Racek, M. et al. Re-evaluation of polyphase kinematic and 40 Ar/39 Ar cooling history of Moldanubian hot nappe at the eastern margin of the Bohemian Massif. *Int. J. Earth Sci.* **106**, 397–420 (2017).

67. Tomek, Č., Dvořáková, V. & Vrána, S. Geological interpretation of the 9HR and 503M seismic profiles in western Bohemia. *Sborník Geologických Věd* **47**, 43–50 (1997).

68. Hrubcová, P. et al. Crustal and uppermost mantle structure of the Bohemian Massif based on CELEBRATION 2000 data. *J. Geophys. Res.* **110**, B11305 (2005).

69. Schmädicke, E., Okrusch, M. & Schmidt, W. Eclogite-facies rocks in the Saxonian Erzgebirge, Germany: high pressure metamorphism under contrasting PT conditions. *Contrib. Miner. Petrol.* **110**, 226–241 (1992).

70. Klápová, H., Konopásek, J. & Schulmann, K. Eclogites from the Czech part of the Erzgebirge: multi-stage metamorphic and structural evolution. *J. Geol. Soc.* **155**, 567–583 (1998).

71. Cooke, R. A. High-pressure/temperature metamorphism in the St. Leonhard Granulite Massif, Austria: evidence from intermediate pyroxene-bearing granulites. *Int. J. Earth Sci.* **89**, 631–651 (2000).
72. Schantl, P., Hauzenberger, C., Finger, F., Müller, T. & Linner, M. New evidence for the prograde and retrograde PT-path of high-pressure granulites, Moldanubian Zone, Lower Austria, by Zr-in-rutile thermometry and garnet diffusion modelling. *Lithos* **342**, 420–439 (2019).

73. Kotková, J., Whitehouse, M., Schaltegger, U. & D'Abzac, F. X. The fate of zircon during UHT-UHP metamorphism: isotopic (U/Pb, δ18O, Hf) and trace element constraints. *J. Metamorph. Geol.* **34**, 719–739 (2016).

## Acknowledgements

This study was funded by the GAČR (Grantová Agentura České Republiky) project 19–25035 S (P.M., K.S., and P.S.). The Center for Geosphere Dynamics (UNCE/SCI/006) is acknowledged for the support of O.L. Software GMT (https://wwwridic-mapping-tools.org/) was used for visualization of the modeling results. The maps were processed using QGIS 3 and Inkscape softwares.

## Author contributions

P.M. designed the study, conducted the numerical experiments, and interpreted their results. K.S. designed the study and interpreted the natural data in a broader perspective. P.S. compiled and interpreted the natural PT-time data for the orogenic root association. T.G. designed the study, designed the numerical code, and the model setup. O.L. analyzed and interpreted the natural data for the plate interface association. All authors discussed the results and problems, and participated in writing the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-021-00122-w.

**Correspondence** and requests for materials should be addressed to P.M.

**Peer review information** Primary handling editor: Joe Aslin

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021