Model design and material properties

A broad suite of drip tectonics models that explore a large range of modeling space (e.g., density, activation energy, viscosity, initial size and the width of the drip) are conducted in this work. By showing a representative suite of experiments in the context of orogenic evolution of Central Anatolia (and possibly to other areas where dripping lithosphere has been proposed), we have focused on cases tailored to the region and highlighting the influence of the important parameters. The choice of major model parameter approximations for Central Anatolia are given below (please see Supplementary File ¹ for other model parameters).

The wet quartzite rheology used in these models is an approximation for the general lithological properties of the Central Anatolian basins (e.g., Tuz Gölü basin) that are made up of thick layers of Oligo-Miocene conglomerate and Pliocene sandstones²⁹. When weak lower crust is inserted in the models (e.g., felsic granulite rheology) the dense layer starts to peel away/delaminate depending on the lateral extent of the weak layer⁹. Based on the same flow law used for the mantle (see Model set-up section), the viscous deformation of the continental crust is controlled by the material parameter (A = 1.1 × 10²⁸ Pa⁻⁴ s⁻¹), power law exponent (n = 4), and the activation energy (Q = 223 kJ mol⁻¹) based on wet quartzite⁶⁴ (Supplementary Table ¹). In addition to the viscous response, the crust is able to deform by frictional plastic yielding in which the Drucker–Prager yeld criterion is used, equivalent to the Coulomb criterion in plane strain, σ _y = p sin ϕ + Cc (crust). For the crust, an empirical weakening is imposed with the internal angle of friction varying from ϕ = 15−2° dependent on the strain. Here, ϕ = 15° is an effective internal angle of friction that implicitly includes the effects of pore fluid pressure P _f in the crust. The weakening/softening process of the crust represents the increasing fluid pressure through partial melt and fluid infiltration into the Kırşehir arc by subduction processes. This is a regular approach in these types/scales of models (e.g., ref. ⁵⁰) Furthermore, the crustal weakening employed in these models implicitly takes into account the shear zone related deformations (e.g., cataclastic flow, fault gouges)⁶⁵. In these experiments, the isostatic thickening of the crust is also effective in increasing the elevation of the surface topography, even without asthenospheric mantle upwelling. Further, such elevation change is amplified compared to a perfect Airy isostasy most likely because of the strain weakening implemented in the crust. For instance, when comparing the results of EXP-3 against EXP-1 at later stages the crustal thickening of 7 km corresponds to change in elevation from 1 to 2 km in the center of the drip where asthenospheric mantle upwelling does not occur. For a perfect Airy type isostasy, 1 km of surface topography is associated with ~5.6 km of crustal thickening (assuming ρ _cont =2800 kg m⁻³ and ρ _m =3300 kg m⁻³).

A locally pre-existing perturbation (280 km wide and 160 km thick) represents a shortened and thickened Kırşehir arc root at the center of the model that instigates the removal process. Inferences made by geodynamic models suggest that such plate shortening may not only effect the crust and but also the underlying mantle part of the lithosphere (sub-arc lithosphere) which may lead to lithospheric instabilities^6,49,50. The initial width of the arc (~280 km) is based on the distribution of arc-related granitoids in Central Anatolia, and the 160 km thickness is chosen for 25% thickening, inferred by paleomagnetic restorations³¹. The geometry of the instability is chosen based on available observational constraints for the arc root under Central Anatolia. The duration of the models is kept in the 0–10 My interval because this was when the plateau uplift started to occur and possibly switched from shortening to extension at ~6 Ma.

In these experiments, the reference density of the sub-arc mantle lithosphere (arc root) is initially set to be higher than the underlying asthenospheric mantle. Petrological studies infer the opposite buoyancy conditions especially for older lithospheric plates (e.g., cratons)⁴⁷ due to chemical depletion. However, owing to several geologic factors a higher density of the arc root is suitable for the ~80 my old Central Anatolia (Kırşehir) continental magmatic arc. First, sub-arc mantle lithosphere (arc roots) in thickened/shortened mature arcs, (e.g., Late the Cretaceous arc in Sierra Nevada) may become denser because of refertilization/enrichment due to the entrapment of the melts induced by subduction fluids⁶⁶. The sub-batholitic arc root for various regions is considered to be of high relative density^47,51,67, as with the Kırşehir arc root. In support of this interpretation, it has been suggested that the Central Anatolian basalts (Cappadocia area) carry components of enriched/refertilized mantle, interacted with slab derived fluids²³. A second geologic factor for the dense unstable arc root in Central Anatolia is potential eclogitization of the lower crust. Shortening of lower crust (as in the region) can transform gabbroic rocks into eclogitic facies conditions⁶⁸. Subsequently, eclogites in the lower crust can sink into the underlying arc root and promote the instability with other mafic residues⁶⁹. Though the mechanism is plausible for Central Anatolia, the evidence for the presence of eclogite bearing xenoliths has not been yet been substantiated.

The effects of varying viscosities of the dripping arc root are tested and an alternative model with 10 times higher viscosity (ranging from μ = 5.10²⁰‒5.10²²Pa s for the same strain rates as the preferred experiment EXP-1) is shown in Supplementary Fig. ². The EXP-1 is associated with lower arc root viscosity which controls the initiation of drip which started at ~10 Ma. This timing is and the evolution of surface tectonics is consistent with the geological evolution of Anatolia.

The numerical (width) and (depth) resolution is 201 × 101 Eulerian nodes and 601 × 301 Lagrangian nodes. Half of the Eulerian and Lagrangian elements are concentrated in the top 160 km in order to enhance resolution in the lithosphere. The model has a free top surface, allowing topography to develop as the model evolves. The mechanical boundary conditions at the other three sides are defined by zero tangential stress and normal velocity (e.g., “free slip”). We have extended the depth of the solution space into the lower mantle so that the sinking mantle lithosphere material moves away from the lithosphere. The initial geotherm for the experiments is laterally uniform and is defined by a surface temperature of 25 °C, an increase to 550 °C at the Moho, an increase to 1350 °C at the base of the mantle lithosphere, and an increase to 1525 °C at the bottom of the model. We tested the influence of Moho temperature (350–750 °C) in a series of models and the general model predictions discussed here (e.g., surface elevation, crustal thickness) are only minorly affected by this.

The surface and bottom temperatures are held constant throughout the experiments and the heat flux across the side boundaries is zero. The initial temperature profile is the same in all experiments. Thermal properties (thermal conductivity k = 2.25 W m⁻¹ K⁻¹, heat capacity c _p =1250 J kg⁻¹ K⁻¹) are the same for all materials and we ignore radioactive heat production and shear heating in the model. In this work we do not take into account explicitly petrological processes of decompression and/or hydrated mantle melting.