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

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. Deeply subducted continental crust can form partially molten buoyant diapirs and explain the presence of ultra-high pressure rocks far from plate boundaries, according to a comparison of numerical models and observations from the Bohemian Massif

S ince 1970s, it has been accepted that oceanic subduction beneath a continent can be followed by arrival of continental mass into the subduction zone 1 and result into collisional deformation and formation of a wide mountain belt 2 . 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 3 . 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 4 . 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 5 . This led some authors to compare it to the Himalayan-Tibetan orogen 6,7 . However, here the channel flow exhumation scenario fails to explain the occurrence of (ultra-)high pressure (U)HP continentderived rocks and mantle relics in the axial zone of the Variscan orogen several hundred kilometers away from the plate interface 8 . This rock association was interpreted as klippen of "root-less" nappes of unknown origin overlying the accreted crust 9 .
However, an alternative scenario was proposed 10 to explain their position in the orogen based on the general architecture of the Variscan Bohemian Massif 11 (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 10,[12][13][14] . In the present work we further develop this scenario and suggest that the lowerplate 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 15 , but it was not compared with the real pressure-temperature-time (PT-time) datasets and the complexity of the exhumed mantle-crust rock associations.  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 60,61,16 . 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) 62 showing a granulitic core with bodies of garnet peridotite, a mafic rim of garnet amphibolite and serpentinized peridotite surrounded mostly by medium-pressure gneiss.
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: 10 (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 16 , but locally found UHP minerals like coesite and micro-diamond point to pressures of 40 kbar or more 17 . Lower plate schist and orthogneiss record conditions of 650°C at~20 kbar 18 and~700°C at 13 kbar 16 . 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 19,20 ). These bodies are typically accompanied by ultra-potassic plutons that formed by mixing of metasomatized mantle-and crustderived melts 12 . The (U)HP granulites contain numerous xenoliths of garnet and spinel peridotites with pressures estimated to 20-40 kbar 8 . The metamorphic conditions of the hosting granulites are~800-1000°C at~20 kbar [21][22][23][24] but the presence of coesite and micro-diamond indicates early UHP conditions 25,26 . 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 27 . The HP-HT conditions are typically followed by nearly isothermal decompression 22,28,23 .
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 [29][30][31][32] .
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 mantlederived 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, translithospheric diapirs develop. This is because we explore a limited part of the parameter space where the diapiric exhumation is the preferred evolutionary style 15 . 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 upperplate 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 sublithospheric mantle, water is released from the subducted  63 ). 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. 11 and regional studies 28,29,31,[64][65][66] ). Extrapolation of surface geology is constrained by seismic reflection and refraction studies (solid and dashed red lines) 67,68 . 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.
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 25,27 . 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 17 , 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  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 33 . 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 34 . 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 35 . 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 20 . 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 sheetlike upwellings in 3D 36 . As a result, the sub-lithospheric flow might be irregular or laterally discontinuous, and the diapirs could have complex shapes 37 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 38 .
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 39 , the Vosges in Eastern France 40 , the Monts du Lyonnais in the French Massif Central 41 and in the Variscan Ulten zone in the Alps 42 . 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 43 and modeled in this work. If we assume that the orogenic root association described above is indicative of translithospheric 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 44 , 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 45 . 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 46,47 similar as in our models. Further south the UHP rocks of Tso Morari and Kaghan units [48][49][50] 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 51 . 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 52 . These examples show that translithospheric 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 53 . It solves the equations of continuity, motion, and heat transport using a finite volumes scheme on a fully staggered computational grid. It adopts viscoplastic temperature-dependent rheology, percolation of fluids with constant fluidto-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" 54 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. 15 , who studied the fate of the crust during subduction of continental lithosphere. They showed that translithospheric diapirs ("exhumation near the arc or in the back-arc" in 15 ), 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 11 . 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 −1 . 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 55 . The melt weakening hypothesis is justified by both analytical derivations and numerical models of lithospheric deformation assisted by frequent dyke propagation 55 . It is further supported by geodynamic models that show a link between partial melting of the mantle and weakening of the lithosphere 56 . It is also an essential physical mechanism for reproducing broad variability of Venus coronae shapes formed by tectono-magmatic plume-lithosphere interactions 57 . 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. 55 , 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 grainboundary sliding that was associated with significant weakening during cooling at 800°C, and subsequent hardening at~700°C 58,59 . 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.