Nature of orogenesis and volcanism in the Caucasus region based on results of regional tomography

In the paper, we discuss the problem of continental collision and related volcanism in the Caucasus and surrounding areas based on the analysis of the upper mantle seismic structure in a recently derived model by Koulakov (2011). This model, which includes P and Svelocity anomalies down to 1000 km depth, was obtained from tomographic inversion of worldwide travel time data from the catalogue of the International Seismological Center. It can be seen that the Caucasus region is squeezed between two continental plates, Arabian to the south and European to the north, which are displayed in the tomographic model as high-velocity bodies down to about 200–250 km depth. On the contrary, a very bright low-velocity anomaly beneath the collision area implies that the lithosphere in this zone is very thin, which is also supported by strong horizontal deformations and crustal thickening indicating weak properties of the lithosphere. In the contact between stable continental and collision zones, we observe a rather complex alternation of seismic anomalies having the shapes of sinking drops. We propose that the convergence process causes crustal thickening and transformation of the lower crust material into the dense eclogite. When achieving a critical mass, the dense eclogitic drops trigger detachment of the mantle lithosphere and its delamination. The observed high-velocity bodies in the upper mantle may indicate the parts of the descending mantle lithosphere which were detached from the edges of the continental lithosphere plates. Very thin, or even absent, mantle parts of the lithosphere leads to the presence of hot asthenosphere just below the crust. The crustal shortening and eclogitisation of the lower crustal layer leads to the dominantly felsic composition of the crust which is favourable for the upward heat transport from the mantle. This, and also the factors of frictional heating and the radioactivity of felsic rocks, may be the origin of volcanic centres in the Caucasus and surrounding collisional areas.


Introduction
Caucasus is a part of the Alpine-Himalayan orogenic belt which is the largest continental collision zone in the world.In the Caucasis segment of the belt, the collision occurs due to convergence of the Arabian and European continental plates in a zone located between two basins of presumably oceanic nature, Black Sea and South Caucasian Plate (Fig. 1).This collision determines active recent tectonic processes manifested in intensive mountain building, seismicity and Cenozoic volcanism.High level of seismic hazard in this densely populated region is one of the main reasons for vital interest to the tectonics of this region.
Mechanism of continental collision is presently not completely understood and it is actively discussed by specialists in different domains of geosciences.Considering most examples of continental collision (e.g., Dewey et al., 1986), one can see that convergence of continental blocks causes considerable crustal thickening which is roughly proportional to the value of shortening.At the same time, the fate of the lithosphere beneath the continental belts remains a disputable topic.While doubling the crust, the collision hardly results at thickening of the lithosphere: strong thick lithosphere would make impossible the observed active deformations in orogenic belts.Active mountain building and strong deformations imply that the lithosphere in the collision zones is weak, and this means that a part of the mantle lithosphere disappears.However, the details of the mantle lithosphere recycling are still not well understood.Is it subducted to the mantle similarly as in cases of oceanic subduction or sink in another way?These questions are actively discussed in the scientific community (e.g., Dewey and Bird, 1970;England and Houseman, 1989;Ershov and Nikishin, 2004).Seismic tomography, which allows imaging the structures at great depths, is one of the most powerful tools to clarify these and other geodynamic questions.
Deep seismic structure beneath the Caucasus region has been investigated in many geophysical studies, mainly using seismic tomography.Mantle structure beneath Caucasus and the surrounding areas has been studied using global and regional seismic modelling based on travel times of body waves (e.g., Neprochnov et al., 1970;Hearn and Ni, 1994;Al-Lazki et al., 2004;Gök et al., 2003), surface waves (e.g., Maggi and Priestley, 2005;Sandvol et al., 2001;Villasenor et al., 2001) and seismic attenuation (Sarker and Abers, 1998).
Most of the studies display generally consistent features in the lithosphere depth intervals (50-250 km) with higher seismic velocities in areas corresponding to the continental blocks and low velocities beneath the folded belt.Here we base our discussion on a recent seismic model of Asia by Koulakov (2011) which enables higher frequency features compared to most of the previously published regional and global tomographic studies.In order to verify several issues and to check the reliability of the proposed geodynamic scenario, in this study we provide additional synthetic tests oriented specially to the target region.Based on considering the tomographic model by Koulakov (2011) we will provide our answers to the following two questions: 1. What happens to the mantle part of the continental lithosphere during the continent-continent collision in Caucasus?
2. What is the nature of the active Cenozoic volcanism in Caucasus and surrounding areas?

Geodynamics and volcanism of Caucasus and surrounding areas
The geological evolution of the Caucasian region (Fig. 1) mostly controlled by convergence of Eurasian and Africa-Arabian continental lithosphere plates.According to geodetic data, the total rate of the convergence is ∼20-30 mm yr −1 (e.g., DeMets et al., 1990).More detailed analysis of regional deformations shows that about 60 % of this rate is taken by the Lesser Caucasian suture, and the rest is accommodated in crustal shortening in the Southern Caucasus (Allen et al., 2004;Forte et al., 2010).Pre-Cenozoic evolution of this area was connected with the closing of the Tethys Ocean between Eurasian and Gondwana continental parts (e.g., Khain 1975;Adamia, 1975;Adamia et al., 2008;Zakariadze et al., 2007).Convergence of the ocean with continents during the Late Proterozoic-Early Cenozoic pre-collisional stage resulted at accretion of island arcs, intra-arc rifts and back-arc basins etc.Thus, very large variety of arc volcanism age can be found in the collision zone around Caucasus (Adamia et al., 2011).Foldthrust belts in the Great and Lesser Caucasus and, in between, the Transcaucasian intermontane depression were formed after definitive closure of Tethys in this segment of the collision belt during syn-collisional (Oligocene-Middle Miocene) and post-collisional (Late Miocene-Quaternary) stages of the Late Alpine tectonic cycle (Adamia et al., 2011).High intensity of the collision processes in Caucasus are possibly due to its location in a gap between two presumably oceanic type basins, Black Sea and South Caspian, which are thought to be the remnant parts of Tethys (see for example, Fig. 21 in Adamia et al., 2011).
Late Cenozoic intrusive and extrusive volcanism is widespread throughout the collision zone from Anatolia to Iran.In Turkey the outcrops of calc-alkali volcanic rocks are observed along two suture zones in different sides of the Anatolian peninsula.The southern volcanic branch extends to Iran along the Urmieh-Dokhtar magmatic arc which is oriented parallel to the Zagros belt.The northern branch appears to be adjoining with intensive post-collisional magmatic manifestations in Caucasus.
Quaternary volcanism in Caucasus and surrounding areas is mostly represented by andesites-to-dacite series.The dacitic lavas were actively erupted in a time period from 760 000 a to 30 000 a in the Javakheti highlands (Lebedev et al., 2004).Products of the most recent and mostly uplifted segments of Caucasus including post-collisional volcanoes of the Elbrus, Chegem and Keli-Kazbegi are represented by lavas of calc-alkaline subalkaline andesite-basalt, andesitedacite rhyolite composition (Tutberidze 2004;Koronovsky and Demina, 2007).There are many evidences of Pliocene-Quaternary ages of eruptions for some of these volcanoes, for example, ∼6000 years near the Kazbegi volcano (Djanelidze et al., 1982).It is interesting that some authors (e.g., Lebedev et al., 2008) observe a "dominoes effect" when the magmatic activity is migrated northward from one volcanic center to another.
Several geodynamic models have been proposed to explain the Late Miocene-Quaternary calc-alkaline volcanism of Caucasus, such as the detachment model of the last piece of subducted oceanic lithosphere (e.g., Innocenti et al., 1982) or the lithosphere delamination (Pearce et al., 1990;Keskin et al., 1998).In this paper, we will provide additional arguments for the second concept based on recent tomographic images of the upper mantle.

Tomographic model
This study is based on analysis of P and S velocity models beneath Asia down to 1000 km depth computed by Koulakov (2011).This model was constructed using arrival times of seismic body waves reported in the worldwide ISC catalogue (ISC, 2001) in the time period from 1964 to 2007 based on the tomographic approach developed in Koulakov and Sobolev, (2006).All the data from the catalogue were initially reprocessed which resulted at relocation of the events and rejection of large amount of outliers.
The inversion was performed separately in 32 overlapped circular areas which covered most part of Asia.All data with ray paths traveling, at least partly, through the study volume, were considered in this study.This included, the data from events located in the study area recorded by the worldwide station network and picks from long-distant events recorded by stations in the study region.It was important that free inversion parameters were determined separately for each circular window based on the results of synthetic modelling.This allowed minimising a problem existing in global studies when the same damping parameter caused loss of information in densely covered areas and artificial instabilities in parts with insufficient data amount.In the presented model, for each window, the value of damping was separately tuned to enable the optimal reconstruction of a synthetic model.
For the selected area, we used four circular windows of 8 • radius with approximately 300 000 rays in total.S-data took about 10 % of the total amount.The distribution of stations and events used in this study is shown in Fig. 2. Number of the parameterisation nodes in each window varied from 6000 to 9000 for the P-data and from 4000 to 5500 for the S-data, depending on the data coverage.
Here we provide several horizontal and vertical sections of the considered model shown in Figs. 3, 4 and 5 which correspond to the Caucasus region.Note that one should be careful with the interpretation of absolute values of anomalies given in this model.The amplitudes of anomalies in seismic tomography studies are strongly affected by damping values which are used for the inversion.This problem is especially important in the case of using the ISC data which contain strong noise.To extract a coherent signal, one should apply strong damping which reduces the amplitudes of retrieved heterogeneities.It can be seen, for example, that in the presented results, the amplitudes of P-anomalies are stronger in shallow layers than those of the S-anomalies, and this can be explained by stronger noise level in S-data and, correspondingly, higher damping used in inversion.In most cases, true amplitudes cannot be achieved, as the reducing of damping causes the inversion instability.This is a fundamental problem which is actual not only for the model considered here, but for any tomographic studies.This should be taken into account when interpreting these results on a quantitative level, and especially, when converting P and S velocities into petrophysical parameters (temperature, composition, density, etc.).One of the approaches to estimate the realistic amplitudes of heterogeneities is synthetic modelling simulating realistic patterns, noise level and the main workflow used for real data processing.Comparison of the initial and recovered structures gives an idea about reduction of anomalies due to damping.
The model by Koulakov (2011) has been verified using many different tests.For example, the contribution of the random noise in the data was estimated based on the "odd/even" test with independent inversions of two data subsets with odd and even numbers of events.The spatial resolution was evaluated using several checkerboard tests with different sizes of patterns.In order to ground the approach with inversions in overlapping windows, the synthetic modelling was performed using a model with realistic shapes of structures.The travel times for this test were computed in the entire area, whereas the inversion was performed in separate windows.
Here we provide some additional tests.In Fig. 6, we present a series of checkerboard models: two models Fig. 4. Same as Fig. 3, but for the S-velocity anomalies for the P-data with the size of 2 • × 1.5 • × 300 km and 1.5 • × 1 • × 200 km and two models for the S-data (3 × 2.5 × 300 km and 2 × 1.5 × 200 km).In all cases, these models were finer than considered in tests in Koulakov (2011).When computing the synthetic data, we added the random noise with the rms of 0.3 s.Furthermore, to simulate the existence of blunders in the ISC catalogue, we also added 5 % of "outliers" for which the noise was multiplied by ten.As a result, the variance reduction after the inversion for the synthetic data was about 25-35 % which is significantly lower than in the case of real data (45-50 %).Despite these "pessimistic" simulations, it can be seen that in areas with sufficient amount of data, the checkerboard patterns can be correctly resolved.For the coarser model, the resolved area covers the entire central part of the study area.For the finer board, having a size of about 100 km, the satisfactory recon-struction is achieved only in areas of Turkey, Caucasus and Zagros; in most other areas the anomalies are not visible.Fortunately for us, the areas with the highest resolution are the most attractive from the geological point-of-view and mostly discussed in the next section.
Another test shown in Fig. 7a consists of the recovering of a model with realistic shapes of anomalies defined in vertical Sect. 2. The conditions of synthetic modelling were the same as in the cases of the checkerboard tests.Because of larger size of synthetic patterns, close to the real anomalies, the value of variance reduction in this case was about 50 %, which is similar to one observed for the real data inversion.The reconstruction results show that the shapes and locations of most features are generally correct both for P and S data.However, at greater depths, where the amount of rays is much lower, the resolution of recovering is much poorer, the amplitudes of anomalies are much weaker than in the "true" model and they are strongly smeared, especially for the S model.For example, for the high-velocity "drops", in the reconstruction results we cannot separate them and say exact the number of anomalies.At the same time, it is important for the interpretation that we can detect the existence of these drops in the mantle, though without resolving their details.This should be taken into account while constructing a geodynamic interpretation.
To check the possibility of vertical smearing and leakage of the crustal anomalies to the mantle, we made another test shown in Fig. 7b.The configurations of synthetic anomalies are the same as in the previous case, except for the low-velocity anomaly in the middle part of the profile defined down to 50-60 km which represents the thick crust.It can be seen that this anomaly is correctly resolved in both P and S velocities; no vertical leakage is observed and the lower boundary of the "crust" is reconstructed at the correct depth.From this test, we can conclude that the low-velocity anomaly, which is observed beneath the Caucasus mountains, really represents the mantle structure.

Discussion
In shallower depth sections down to 220 km depth, we can clearly observe higher P and S velocities associated to the south with the Arabian plate and to the north with the European plate which consists of several microplates in the contact zone, such as Scythian and Turan plates.Lower seismic velocities are observed beneath the collision zone in the areas of the major mountain belts.P and S models display generally consistent structures; however, the amplitude of P anomalies in the shallower sections is considerably higher.However, as was discussed in the previous section, this reflects rather the damping issues that the real relationships of amplitudes.It can be seen that all recent volcanic centres exactly fit to the low-velocity patterns of P and S anomalies in shallower sections.In vertical sections, we can see that the thickness of highvelocity layers related to the Arabian and European lithospheric plates is about 200-250 km which is a little bit higher than estimated by other authors based on different methods (e.g., Artemieva, 2003).However, it should be kept in mind that the lithosphere related anomalies might be smeared downward due to the limited vertical resolution.In the transition zones between high-velocities in continental blocks and low-velocities in the collisional belt, the structure of anomalies is rather complex with alternating high-and low-velocity anomalies.It can be seen that high-velocity anomalies form drop-shaped bodies which seem to sink to a greater depth mantle.Taking into account the results of synthetic test with vertical anomalies in Fig. 7, we can propose that the real data inversion smears and reduces the amplitudes of anomalies at greater depths.Thus, the true amplitudes of anomalies in these drops might be stronger than observed after inversion of real data.
The obtained results allow us to propose mechanisms of the lithosphere recycling due to collision and origin of volcanism illustrated in Fig. 8.It can be seen that the Arabian and Eurasian parts are represented by approximately similar lithosphere type of about 200-250 km thick.It can be proposed that the lithosphere of these plates have the standard continental type structure which includes upper felsic (gran-ite) and lower mafic (basaltic) crustal layers and a rigid mantle layer which dominates in total strength of the lithosphere (e.g., Burov and Diament, 1995).When the plates collide, the crust thickens in the shortening areas between these plates.In this case, the lower mafic crustal layer appears at greater depth.Temperature and press increase lead to phase transformation of the mafic layer into denser eclogite (e.g., Sobolev et al., 2006).The drops of eclogite are united into larger bodies which, after reaching a critical mass, descend to the mantle.These changes in the lower crust may lead to detachment of the gravitationally unstable mantle lithosphere.The presence of dense eclogite drops may trigger the lithosphere detachment and accelerate sinking of separate pieces of the lithosphere as shown in Fig. 8.As proposed by Burov and Watt (2006), this may lead to the "crème-brûlé" behaviour of the mantle lithosphere in the collision zone instead of "jelly sandwich" rheology which is characteristic for nondeformed continental lithosphere.This process abruptly decreases the total strength of the lithosphere and leads to its fast degradation through active delamination (e.g., Kay and Kay, 1993).According to this hypothesis, the mantle lithosphere at the edges of the collided continental plates should be gradually destroyed and delaminate together with dense eclogite produced in the lower crust.These sinking drops are probably visible as high-velocity bodies in vertical sections of our seismic model beneath the edges of the Arabian and European continental parts and beneath the collision zone.
It is important to mention that a similar processes (tectonic shortening -crustal thickening -eclogitisation of the lower crust that triggers delamination of the mantle lithosphere) was previously suggested for Central Andes based on petrological arguments (Kay and Kay, 1993) and numerical modelling (Sobolev and Babeyko, 2005).Rate of this process was analysed by Jull and Kelemen (2001).
After detachment of the mantle lithosphere, it is replaced by hot asthenosphere which may appear directly beneath the crust.Note that the crust in the collision zone is mostly composed of thick felsic rocks, whereas the lower mafic crust was transformed to eclogite and sank together with the mantle lithosphere.The existence of thick felsic crust is supported by the observed very low velocities in the collision zone at 50 km depth, which can be considered as an integral layer for the crustal properties.It is known that the felsic layer is composed of mechanically weaker rocks which facilitate circulation of hot materials in the thick crust (e.g., Babeyko et al., 2002;Babeyko and Sobolev, 2005).This favours active heat transport from the asthenosphere to the surface that explains the existence of active volcanic fields in Caucasus and surrounding collisional areas.At the same time, additional heating can come from the frictional effects due to the strong compressional deformations.Some authors sug-gest that this factor cannot be ignored in the process of delamination (e.g., Schott et al., 2000) and in the origin of volcanism (e.g., Stüwe, 1998).Some contribution in heating can also be related to radioactivity of felsic rocks.

Conclusions
Based on our seismic model, we can conclude that mantle lithosphere beneath the collision zone of Caucasus and surrounding areas is very thin or absent.That is why this segment squeezed between rigid Arabian and European blocks behaves as weak lithosphere and is strongly affected by tectonic shortening and orogenesis.Tomography results show that the continental lithosphere in the contact with the collision area is actively destroyed by the process of delamination and sink as separate drops.We propose that a big role in triggering the detachment of the mantle lithosphere and its more active descending is played by eclogitisation in the lower crust affected by strong shortening.We believe that delamination mechanism of such type is the major candidate for the lithosphere recycling in all continent-to-continent collision zones of the world.
Important conclusion of this research consists in the definition of a collisional type of volcanism which is principally different of the intraplate and subduction types of volcanism.We propose that the volcanic activity in Caucasus and surrounding collisional areas is presumably due to direct heating of the crust from the asthenosphere which is possible due to lack of the mantle lithosphere and thinning of the mafic lower crust layer.This type of volcanism might be expected in most areas of continental collision, however, in practice it is observed in a limited number of regions.This might be explained by a hypothesis that the collisional type of volcanism occurs only under the condition of complete detachment of the mantle lithosphere which is not realised in all cases of continental collision.However, this topic needs additional investigations based on thermo-mechanical modelling and analysis of geological data in different collisional belts.

Fig. 1 .
Fig. 1.Main tectonic units in Caucasus and surrounding areas overlaid on a shaded relief map.Yellow stars depict the recent volcanoes in Caucasus (Adamia et al., 2011).Major volcanoes are named with yellow characters.Areas of Cenozoic volcanism in Iran compiled from Nezafati (2006), Verdel et al. (2007) are marked with yellow.TP is Turan Plate; KD is Kopeth-Dagh; SCB is South Caspian Basin.White arrow marks the direction of the Arabian Plate displacement.

Fig. 2 .
Fig. 2. Distribution of data from the ISC catalogue: triangles depict stations, red dots are the events.

Fig. 3 .
Fig. 3. P-velocity anomalies in six horizontal sections, yellow stars and polygons mark the locations of recent volcanoes and folcanic fields in Caucasus and surrounding areas.

Fig. 5 .
Fig. 5. P-and S-velocity anomalies in three vertical sections.The locations of the profiles are shown in maps in Figs. 3 and 4. Relief along the profile is shown above each plot.

Fig. 6 .
Fig. 6.Checkerboard tests with different sizes of anomalies for P-and S-models in three horizontal sections.Depths of the sections correspond to the middle level of the checkerboard patterns.The sizes of synthetic anomalies are indicated above each column.

Fig. 7 .
Fig. 7. Two synthetic tests with realistic patterns defined in a vertical Sect.2, same as indicated in Fig. 5. Upper plots show the configurations of the synthetic models; middle and lower plots are the reconstruction results for the P and S anomalies.

Fig. 8 .
Fig. 8. Schematic representation of the delamination mechanism in the Caucasus region.The crust is composed of the upper felsic (orange) and lower mafic (green) parts.Blue areas indicate the mantle parts of the lithosphere.Background is the distribution of P-velocity anomalies in vertical Sect.2, same as in Fig. 5. Relief along the profile is shown above the plot.Green ellipses schematically mark the possible locations of the eclogite drops which were transformed from the lower mafic crust in the shortening zone.Red arrows mark the asthenosphere upwelling.