The Caucasus Territory Hot-Cold Spots Determination and Description Using 2D Surface

Waves Tomography

Seyed Hossein Abrehdari ^{1, 2, ✉}, Jon K. Karapetyan³, Habib Rahimi ⁴, Eduard Geodakyan³

¹ PhD Student, Institute of Geophysics and Engineering Seismology, National Academy of Sciences; Gyumri, Armenia

² Seismology Research Division, Institute of Geophysics, University of Tehran; Tehran, Iran 

³ Candidate of Geol. Sci., Institute of Geophysics and Engineering Seismology, National Academy of Sciences; Gyumri, Armenia

⁴ Associate Professor, Department of Earth Physics, Institute of Geophysics, University of Tehran; Tehran, Iran

^✉ Corresponding author e-mail: abrehdari@ut.ac.ir (S. H. Abrehdari)

Abstract

Many questions have been raised about the thermal-mechanical development of plate tectonics boundary interactions, lithospheric processes, mantle activity, movement of faults, continental thinning, and generally the heat beneath our feet. The earthquake waves are originating in the Earth’s crust or upper mantle, which ricochet around the earth's interior and traveling most rapidly through cold, dense regions, and more slowly through hotter rocks. In this paper, in order to identify and describe the Caucasus territory Hot-Cold spots and better understand the regional tectonic activities based on the fast and slow wave velocity anomalies, the 2D tomographic maps of Rayleigh wave dispersion curves were imaged. To obtain these maps in the ever-evolving collision zone of the Eurasian-Arabic plates, we performed a 2D-linear inversion procedure on the Rayleigh wave in a period ranging from 5 to 70 s (depth= ~200 km). To conduct this, ~1500 local-regional earthquakes (M≥3.7) recorded by the 48 broadband-short period stations from 1999 to 2018 were used. In this study, we assumed that the low-velocity tomography images or dark red-orange shades indicate hot spots (slow-regions) and high-velocity or dark blue-green-yellow shades imply cold spots (fast-regions). Therefore, by using the technique of increasing-decreasing the velocity anomaly in a wide area with complicated tectonic units the hot-zones and extensive cold-aseismic areas were described and investigated. Hence, for short-periods (5≤T≤25 s; 6.6≤depth≤30.8 km) 15 hot spots were determined. The result for medium-periods (30≤T≤45 s) show two hot spots with a depth of ~108 km. In long-periods (depth= ~200 km), most part of the study area has covered by ultra-low-velocity anomaly as a permanent hot spots.

Keywords: Caucasus territory, Hot-Cold spots determination, 2D surface wave tomography, Geothermal resources, 2D linear inversion

1. Introduction

The geologists to map earth’s interior, use the seismic waves chart movement (traces) generated by earthquakes that has been recorded by a worldwide network of seismometers. Each recorded wave trace reveals many useful informations from the extent and density of Earth’s deepest regions.

According to the seismic tomography images of [Suzan van der Lee et al., 2019] studies  from within the earth, the colors show anomalies in rigidity, which correlate with temperature anomalies (also visit this website: [https://www.iris.edu/hq/inclass/fact-sheet/seismic_tomography]). The dark blue-green-yellow shades mean colder and stiffer rock (Cold spots) that are the remnants of an old tectonic plate that has been subducted underneath the Earth plates (large cold and aseismic area during million years) and dark red-orange shades mean warmer and weaker regions (Hot spots).

Geologists believe that a hot spot is a location on the Earth's surface that has experienced active volcanism for a long period of time and to be fed by underlying mantle that is anomalously hot compared with the surrounding mantle. The origins of the concept of hot spots lie in the work of Wilson [1963], who postulated that the formation of the hot spot is from the slow movement of a tectonic plate across a hot region beneath the surface.

In effect of the chemical interactions, volcanic zones, tectonic activities, high temperatures, and pressures in the earth's interior some of the rocks melted, the solid behavior becomes plastically and geothermal zones (hot spots) arise. Geologists have identified some 40-50 such hot spots around the globe (Fig. 1).

Figure 1. Mantle plume locations (yellow circles) are the world's famous (Hot spots). {Figure retrieved from: [https://www.sciencelearn.org.nz/images/350-tectonic-plate-boundaries], University of Waikato}.

Estimates for the number of hotspots postulated to be fed by mantle plumes have ranged from about 20 to several thousands, over the years, with most geologists considering a few tens to exist. Hawaii, Réunion, Yellowstone, Galápagos, and Iceland are some of the most active volcanic regions to which the hypothesis is applied [Foulger, 2010; Wright, 2000]. However, keep in mind these are just theories and nobody really knows the answer. The honest answer is that lots of folks are working on it but haven't come up with the answer yet.

Hot spots (close to the earth's surface) are the same geothermal resources and is useful for identifying sources of geothermal energy. For example, some countries such as Iceland (Nesjavellir Geothermal Power Plant), New Zealand (Wairakei Geothermal Power Plant), Turkey (Efeler geothermal power plant) and NW Iran (using water from hot springs or hydrothermal) are using this energy.

This paper attempts to generate the 2D tomographic inversion of Rayleigh wave velocity dispersion maps in order to identify and describe Hot-Cold spots inside the earth based on fast and slow wave velocity anomalies. To do this, the Rayleigh wave group velocity dispersion curves for each source-station path (single-station method) are estimated using the Hermann's do_ mft package [Herrmann, 2013].

Then, using a 2D-linear inversion method developed by [Ditmar and Yanovskaya, 1987] and [Yanovskaya and Ditmar, 1990], the 2D group velocity maps were generated. Some studies of the crustal structure have been conducted related to the estimation of the crustal thickness beneath this region [e.g., NW Iran: Gheitanchi, 1996; Caucasus: Randolph et al., 2007; EAAC: Skobeltsyn, et al., 2014]. These maps show excellent agreement with results of the previous studies and many of the geological features of the Caucasus territory.

The main purpose of this study is to investigate the Hot-Cold spots in the Caucasus based on scattering characteristics of the wave velocity increasing and decreasing using 2D tomographic maps. In the Caucasus dominion, no study has been conducted entitled hot-cold spots using surface wave tomography in large-small scale yet. This study benefits from a rich earthquake database (1999-2018) and new permanent seismic stations installed in NW Iran, Russia, Armenia, Turkey, Azerbaijan, and Georgia, which provides much better ray path coverage in the Caucasus for the resolution of tomographic images. Also, to provide more useful information to researchers about the hot-cold zones, the results of the performed previous tomographic studies in the Caucasus were used. Also, it is for the first time, that we use terms such as temporary, unstable, and permanent hot-cold spots in geology.

2. Caucasus general information

The Caucasus is a very famous and enigmatic landmark from various aspects of ongoing tectonic evolutions. The structure of the Caucasus mountain region and surrounding areas is primarily controlled by the collision and continuing convergence of the Arabian-Eurasian plates. Over time, this motion led to subsequent collision stages between Arabia and smaller continental blocks resulted from the break-up of Gondwana until the final closure of Neo-Tethys Ocean.

Because of this geodynamic evolution, a complex geological structure has formed that characterized by important lateral variations in age, composition and tectonic style [Hatzfeld and Molnar, 2010]. Greater Caucasus (GC), Lesser Caucasus (LC), East Anatolian Accretionary Complex (EAAC), Bitlis Massif (BM), Pontide (PN), NW Iran, Tlesh (TAL), South Caspian Basin (SCB), Kura Basin (KB), Rioni Basin (RB) and Eastern Blake Sea Basin (EBSB) are some part of the Arabia-Eurasia collision.

The Greater Caucasus mountains are an orogenic belt that was raised as a result of the collision, with altitudes of more than 5 km above sea level and 1,300 km in the NW-SE direction between the Black Sea and the southern basin of the Caspian Sea. The convergence between Arabia and Eurasia began in Late Cretaceous [Golonka, 2004].

Also, it is for the first time, that we use terms such as temporary, unstable, and permanent hot-cold spots in geology.

According to some studies [e.g., Copley and Jackson, 2006; Jackson, 1992] this continent-continent collisional tectonics processes has begun at about 12 Ma. Caucasus region is compressed between Arabian-Eurasian plates and due to N-S compression expanded the main seismo-active structures in NW Iran, Greater Caucasus (GC), Lesser Caucasus (LC), Eastern Anatolian Accretionary Complex (EAAC). The fault zones include reverse strike-slip, strike-slip sinistral, strike-slip dextral, wrench and major thrust faults with WNW-ESE direction are developed (Fig. 7). Further, several large Neogene-Quaternary strato-volcanoes [< 0.4 Ma- 6.5 Ma; Bavali et al., 2016] such as Elbrus is located in this region.

According to some studies [Ismail-Zadeh et al., 2020], the Caucasus ranges incload the Greater Caucasus, is consisting mostly of Paleozoic metasedimentary rocks and granitoids, Jurassic sediments, Mesozoic and Cenozoic volcanism and the Lesser Caucasus consist Paleozoic granitoid metamorphic basement overlain unconformably by shelf carbonates of Paleozoic Triassic age, respectively. In this enigmatic area, there are complicated geological structures and large volcanic complexes and basins.

3. Data and study area

The study area is the territory of the Caucasus (Longitude: 38°- 53°; Latitude: 37°- 44°). The dataset used for tomography consists of ~1500 local-regional events with magnitude M ≥ 3.7, which have been collected by 48 broadbands and short-period Seismic Network Incorporated Research Institutions for Seismology (IRIS) including Armenia (GNI), Georgia (GO), Turkey) (TK, TU), Azerbaijan (AB), as well as earthquake data recorded by the Iranian Seismological Center (IRSC), International Institute of Earthquake Engineering and Seismology (IIEES), and temporary network of the Institute for Advanced Studies in Basic Sciences (IASBS) from 1999 to 2018 (Table 1). The study area, distributions of earthquakes and seismic stations used in the analysis are shown in Figs. 2 and 7. The Geological units, faults, volcanoes, rivers, and basins that used for interpretation the 2D tomography maps in the study area was shown in Fig. 7.

Figure 2. a) The seismic events (red circles) and b) locations of the stations used in this study.

Table 1. Station and network codes, coordinates, and data center website considered in this study.

3.1 Ray paths distribution and resolution parameters

The resolution of group velocity maps depends more on the density of the paths and their balance distribution (crossing paths). In tomographic studies related to seismic events data, these two parameters depend on the distribution of the earthquakes and geometry of the seismic array that can limit the number of available paths for some directions which is almost beyond our control.

According to Yanovskaya [1997 and 1998], the stretching of the averaging area and the mean size parameters are used to estimate the lateral resolution. For 2D tomography problems, a function S (x, y) for different orientations of the coordinate system is used in order to determine the sizes of the averaging area along different directions. The averaging area which gives us an idea of the obtained resolution can be approximated by an ellipse centered at a point, with axes equal to the largest Smax (x, y) and to the smallest Smin (x, y) values of S (x, y). The smallest Smin (x, y) and largest Smax (x, y) axes of the ellipse are calculated, and the resolution in each point is given the mean size of the averaging area:

L=Sminx,y+Smaxx,y2                                                              (1)

As the resolution is closely correlated to the density of the crossing ray paths in each cell, it is clear that small values of the mean size of the averaging area L  (corresponding to high resolution) should appear in the areas that are crossed by a large number of ray paths and vice versa.

The second parameter is the stretching (ex or ε)  of the averaging area, which provides information on the azimuthual distribution of the ray paths and is given by the ratio:

ex or ε=(Smaxx,y-Sminx,y(Smaxx,y+Sminx,y         (2)

Where Smax  and Smin  are the large and small elliptical focal lengths. According to study [Yanovskaya, 1997], large values of the stretching parameter (usually > 1) imply that the paths have a preferred orientation and along this preferential direction is likely to be quite small. On the contrary, small values of the stretching parameter imply that the paths are more or less uniformly distributed along all directions; hence the resolution at each point can be represented by the mean size of the averaging area. By calculating the Smax ، Smin parameters and L in this study, the value of ɛ  = 0.6985 ≅ 0.7 was obtained.

Figure 3. a) The number of ray paths used in the tomography inversion respective to the period in this study. b) Distribution of seismic stations and interstation path coverage. c) The appropriate value of α (0.2) selected for performing tomography at period of 10 s. d) shows the sample of events-stations rays coverage (yellow lines) at period of 40 s.

Figure 3 a-c shows the number of paths in each period, distribution of stations-earthquakes (T=10 s, a= 0.2), and ray coverage between seismic events epicenters-stations in the study area, respectively. Distribution of stations and earthquakes controls the amount of stretching parameter and data density. The mean size of the averaging area of our tomographic results is of the order of 183.57 (Table 2) in most of the study region. The values of the stretching of the averaging area are between 0.5 and 0.95 in most part of the study area at the fourteen periods. This indicates that the azimuthal distribution of the paths is sufficiently uniform and the resolution is almost the same along any direction. The averaging area value is larger than ~150, with its maximum equal to ~2400. The dense rays path distribution (Figs. 4 and 6, a= 0.2) controls the reliability and the high resolution of tomography results (red shads in fig. 4- averaging area (L)). For more details in this case, see 4.2. sub-section.

Therefore, stretching ε  and averaging area (L) values are two parameters that indicate the orientation and resolution of the different areas within the study area for each period and at any latitude (Y) and longitude (X) direction. Figure 4 illustrates the variation of these two parameters for some periods. Figure 4a-f shows the cell size of 0.2° × 0.5° (20 × 50 km²) for the stretching, data density, and averaging area parameters at period 5-70 s (generated by MATLAB software). The Yanovskaya’s methodology [1997] is used to calculate the spatial resolution, which varies from 20 to 50 km in our study region. We constructed the 2D tomography velocity model in Caucasus by inverting the pure path Rayleigh wave dispersion curves at (35×30 =1050) nodes using a grid with cell size of 0.2° × 0.5° (20 × 50 km²).

Figure 4. Resolution of parameters map: stretching value (ɛ), data density and average value (L) for short-periods of 5, 10, 15, and 20 s; medium-period of 30, 40, and 45 s; and long-periods of 50, 65, and 70 s in this study. Figures a, c-f show the size of the cells 0.2° × 0.5° (20 × 50 km²) for parameters of stretching, averaging area, and data density for the shortest and the longet periods (generated by MATLAB software). Fig. b shows the zoomed cells.

In this study the various α values (0.1, 0.2, and 0.3) were tested and the number of ray paths passing through each cell were observed and selected. So, the checkerboard test was ignored and it is not required to perform. This method is much more useful than checkerboard test (Fig 6).

Table 2. Fluctuations in the parameters of resolution: stretching value (ɛ ), averaging area (L ) value (km), and data density. Also the table shows the velocity (km/s), depth (km), and possible Moho-LAB-LVZ depth (km) discontinuities in each period of 5 to 70 s.

Also, in order to have a good quality data bank in the time domain, we used SNR = 3 for initial processing of data in different periods. Then, by choosing the proper regularization parameter (a)  in the computer code, we improved the resolution of the tomography images by observing the most data and the beam passing through each cell.

4. Methodology

4.1. Dispersion estimation

After preparing the earthquakes waveform and preliminary corrections on it, for each station-earthquake pair (single-station method), the group velocity dispersion curve of Rayleigh waves by applying the Herrmann's do_mft package [Herrmann, 2013] to the vertical component (Z) of motion on each event is estimated. Modified sacmft96 work around problems with improper station and component specifications in Sac files. Sacmft96 is called do_mft for interactive analysis of group velocities and spectral amplitudes. SAC (Seismic Analysis Code) is a general purpose interactive program designed for the study of sequential signals, especially time series data.

In fact, the frequency-time analysis of surface waves is used to estimate the dispersion curves. This method is used for estimating phase and group velocity of surface waves. It passed the preprocessed signal through a system of narrow-band filters in which the central frequency is varying and the amplitude of filter outputs is visualized in time and frequency domains. Then, on the do_mft diagram, the group velocity dispersion curve for each path is obtained. Figure 5 shows an example of determining group velocity dispersion curve for the vertical component (Z) of the Oni (ONI) station using do_mft processing.

To conduct this, we applied Herrmann's do_mft package on waveforms of ~1500 earthquakes recorded by the 48 stations in the Caucasus region. We then processed more than ~34000 vertical component (Z) of dispersion curves (Fig. 5 e). For this purpose, first, in Ubuntu system the earthquake data (miniSEED format) was converted to a SAC file format and then the fundamental mode of the Rayleigh wave for each vertical component (Z) using the do_mft package was determined.

Figure 5a shows the raw, radial component, and the processed waveform. Figure 5b shows the dispersion curve measurement by do_mft to separate the fundamental mode of earthquake data waveform. Figure 5c shows the picked dispersion curve related to the energetic part of the signal (red area) in b for an earthquake recorded in the Oni seismic station (ONI, Georgia).

In the single-station method, since estimating the dispersion curves depends on the basic parameters of earthquakes such as magnitude, epicentral distance, depth, and etc.; different period ranges by applying do_mft to every epicenter-station pair are attained; hence, for different periods we have various path numbers. Finally, a set of dispersion curves for the fundamental mode Rayleigh wave in the period range from 5 to 70s to create 2D tomography maps was estimated (Fig. 5 e). The periods greater than 70 s in this study (see Fig. 5e), due to poor coverage of dispersion curves, which leads to poor coverage of the rays to conduct the tomography images were ignored.

Figure 5. a) Raw, radial component and cleaned seismogram waveform (up to down). b) An example of determining group velocity dispersion curve using Herrmann's do_mft package for the vertical component (Z) of the Oni (ONI- Georgia) seismic station. Earthquake source parameters and the recording station name are mentioned on the dispersion curve plot. The raw waveform is given in the right-hand rectangular block. c) The picked dispersion curve related to the energetic part of the signal (red area) in b. The vertical red lines (with the alphabet ‘O’ above them) show the onset of chosen pickfile by SAC software (automatic default of start reading of arrival time). e) Dispersion curves (~34000 curves).

4.2. Two-dimensional surface wave tomography 

When the ground shakes at the onset of an earthquake, seismic waves race outwards from it and travels most rapidly through cold, dense regions and more slowly through hotter rocks. Waves travel faster through cold-rigid material (like a subduction plate inside the mantle) and passes through warmer materials more slowly (like hot rocks rising to the surface). Seismic tomography is like taking a Computed Tomography or CAT scan of the Earth. In a method similar to CT scans, scientists instead use seismic waves to make images of Earth's interior.

With reference to the interpretation of tomography images [Suzan van der Lee, 2019] from within the earth, the colors show anomalies in rigidity, which correlate with temperature anomalies. Hence, the dark blue-green-yellow shades mean colder and stiffer rock (cold spots areas) and dark red-orange shades mean warmer and weaker regions (hot spots areas). Thus, in this study, the tomography images with dark red-orange shades show low-velocity (slow-hot) zones and dark blue-green-yellow shades indicate high-velocity (fast-cold) areas.

In order to construct group velocity distribution maps in period ranges from 5 to 70 s, a 2D-linear tomographic inversion technique developed by Ditmar-Yanovskaya [1987] and Yanovskaya-Ditmar [1990] was used. This methodology is a generalization of the classical 1-D method of Backus-Gilbert [1968].

The study region was parametrized into grids with cell size of 0.2° × 0.5° and with proper regularization parameters to provide relatively smooth maps with small data misfits. The same regularizations parameters were used for producing the maps, before and after the impoundment. The resolution of the data set is controlled by the average path length, density of paths, and the azimuthal coverage.

In addition to group velocity maps, the corresponding resolution information will be provided at each period. The Yanovskaya’s methodology [1997] is used to calculate the spatial resolution at each point and in different directions. Therefore, it is clear that checkerboard testing is not required and the averaging area maps (red shades) show the highest resolution in each period by choosing the proper regularization parameter(a). The main advantage of this method is that in the cases of uneven distribution of surface wave paths, it works well. The dataset in this method are the travel times along different paths at each period that were calculated by the Herrmann's do_mft package. The method estimates the lateral variation of group velocity Vx  at each period.

The lateral group velocity distribution could be estimated by minimizing the following function:

d-GmTd-GM+aÑm(X)2dX=min              (3)

Where

mX=(V-1(X)-V0-1V0                                                (4)  

Where

di=Ti- Ti0                                                                         (5)

(Gm)i =GiXm(X)dX= Loim(X)dsV0                        (6)

GiXdX= LoidsV0= ti0                                                   (7)

In the relations (3–7), X = (q, f) is the position vector, V0  is the velocity corresponding to a starting model, ti  is the observed travel time along the ith  path, tio is the travel time calculated for the starting model, a is a regularization parameter, tio  is the length of the ith  path and s is the segment along which the inversion is performed.

The regularization parameter (a) that depends on the accuracy of the data, is the trade-off between the fit to the data and the smoothness of resulting velocity distribution. Decrease in a gives a sharper solution region with an increase in solution error, whereas increase in a leads to smoothing of the solution region with decrease in solution error [Yanovskaya et al., 1998].

The parameter a controls the trade-off between the fit to the data and the smoothness of the resulting group velocity maps. Therefore, to improve the resolution and for having a real model, we tested various values of the regularization parameter (a). We chose α = 0.2 that gives relatively smooth maps with small solution errors which was conducted by testing different α values and observing the number of rays passing through each cell size of 0.2° × 0.5° by running the specialized computer codes in MATLAB software which used in this study (Fig. 6); (Also, see 3.1. sub-section).

Fig. 6. Calculations of group velocity maps are imaged for several regularization parameter (a ). Decrease in a gives a sharper solution region with an increase in solution error, whereas increasing α results in smoothing of the solution region with decreasing solution error. The small solution errors by testing different α values (a = 0.1, 0.2, and 0.3) and observing the number of rays path passing through each cell size of 0.2° × 0.5° by running the specialized computer codes in MATLAB software, which was determined in this study. The σ is an estimate of the standard error of the data.

As shown in Figure 6, choosing the proper regularization parameter (α = 0.2), reflects the uniformly distributed ray paths for conducting the tomography (for averaging area and stretching red shades). Finally, by selecting the compatible regularization parameter and using the GSAC, GMT software, and computer specialized codes in the Ubuntu operating system and MATLAB software, the 2D tomography group velocity, stretching, data density, and averaging area maps were plotted and estimated in period ranges from 5 to 70 s (Fig. 4).

Yanovskaya [1997 and 1998] proposed to use two parameters to estimate the lateral resolution: the mean size and the stretching of the averaging area. The resolution is directly controlled by the coverage of the ray paths and the distribution of stations and earthquakes. Therefore, determining the compatible regularization parameters (a) plays a very important role in the resolution of tomography images (Figs. 4 and 6).

Moreover, by determining the compatible regularization parameter (a) and thus correctly to estimate the mean size and the stretching of the averaging area, there is no need to a checkerboard test for resolution and gives relatively smooth maps with small solution errors.

Another criterion that controls the quality of the solution is the comparison between the initial mean square travel time residual and the remaining (unaccounted) residual σ. It is assumed that the unaccounted residuals are random, so σ can be accepted as an estimate of the standard deviation (error) of the data, which allows a standard error of the solution σ_m to be computed. Therefore, in this study, the value of σ is used for the selection of the appropriate data: if the travel time residual for one path is larger than 3σ, this path is eliminated from the dataset and the solution is recalculated [Yanovskaya et al., 1998]. The standard deviation (σ) with selecting the regularization parameter α = 0.2 is reasonably low which showing the stability of the method. (Fig. 6).

5. Hot-Cold spots determination and description using 2D tomography velocity maps

As mentioned earlier, the wave travels more rapidly through cold, dense regions, and more slowly through hotter rocks and zones. Thus, in this study, we assumed that each low-velocity (slow) region with a dark red-orange shade is a hot spot and each high-velocity (fast) region with dark blue-green-yellow shade is a cold spot. Hence, in order to identify and describe Hot-Cold spots inside the earth of the study area based on the increase and decrease of the wave velocity anomalies, the 2D tomographic velocity maps obtained in Fig. 8 were generated.

Schematic diagram of Fig. A1 has depicted to better understand the hot-cold spots procedure inside the earths. Figure A1 shows the physical processes within the Earth’s upper mantle that lead to the generation of magma in steps A to D for different plate tectonic settings. Tomographic maps with distinct velocities over short-periods of 5 to 25 s (equivalent to a depth of 6 to 53 km), are more sensitive to the structure of the upper to lower crust, and Moho. These short-periods represent sediments in the basins, chemical interactions of hydrocarbon resources, molten material and magma chambers beneath volcanoes and Moho discontinuity which these areas can be considered temporary and unstable hot and cold spots.

These periods, which is also known as the crust, include soil, vegetation growth, construction, surface-groundwater, oil-gas resources, magma chambers, metallic, non-metallic mines, and chemical interactions. Although, some researchers believe that slow velocities in the crust or upper mantle under regions of active volcanism do not require "Hot Spot" (i.e, plume-related) magmatism.  These could simply reflect decompression melting and/or crustal melting following slab-rollback, delamination, or breakoff (e.g., Keskin et al., 2008).

The lithosphere-aesthenosphere and upper mantle is the source of volcanic lava and the origin of some earthquakes in the mantle and remnants of the old tectonic plate. So, short-periods (T= 5-25 s) can be considered temporary and unstable hot and cold spots.

Dark red spots that are seen below the chain of volcanoes (e.g., Elbrus, Ararat, Aragats, Kazbek), upon some segments of faults (e.g., PSSF, NTF) and basins (e.g., SCB, KB, EBSB, RB) in the study area indicate the hot spots. And these hot spots are located in an appropriate depth of shallow area of the earth's crust, which can be considered as geothermal energy resources (e.g., Iceland-Nesjavellir Geothermal Power Plant). The rest of the areas on the tomography maps with various periods that has shown with dark blue-green-yellow shades, include colder and more rigid rocks and stones and remnants of an old tectonic plate and is calm (aseismic), which represent cold spots. Cold spots usually cover a wide area and can cover tectonic plates and continents and may even include the mantle core, which is the source of cold (old) lava volcanic and some deep earthquakes. The core controls the Earth planet balance in the solar system and its magnetic property.

Figure 8 illustrates the labelled major geological units (e.g., GC, EAAC, LC, BM); Volcanoes (e.g., El., Ar., Arg., Sab.); Basins (e.g., EBSB, SCB, RB, KB) to describe in this study.  Figure 9 shows all the hot spots identified in this study for different periods.

Tomography velocity map at period 5 s (Fig. 9) shows 8 low-velocity hot spots (dark red shade). Hot spot number 1, 4, and 5 in our study are located in the East Anatolian Accretionary Complex (EAAC) which is known for its thin lithosphere [Sengör et al., 2003] in the east of Turkey. The depth of these hot spots is ~6 to 13.33 km (Fig. 9). 

Hot spots number 1 is a small area near Kars mountain in NE Turkey known as the Erzurum-Kars Volcanic Plateau (EKVP). This part of the plateau has been formed by the eruptions during the Zanclean (~4.5 Ma) period, related an earlier continental collision event between Eurasian and Arabian continents ~15 Ma ago. The EKVP is composed mainly of andesitic and dacitic lavas and their trachytic equivalents intercalated with acidic ignimbrites and tuffs. In the northwest of Kars, an eroded stratovolcano is present which is possibly coeval with the plateau. It consists of a thick sequence of rhyolitic lavas, tuffs and perlitic-obsidian [e.g., Duru Olgun et al., 2020].

Hot spot number 2 and 3 are located in the centeral Armenian block (CAB) and north Armenia, respectively. Hot spot number 2 is approximately located on the northern slope of Aragats volcano (depth of ~7 to 13.66 km), which may be the reason for the existence of magma {(but at 10 s period this volcano is covered with a high-velocity anomaly [as a cold spot], which this feature might be explained by the old age of this volcano which does not express any activity for more than half million years)}. Based on some studies [e.g., Chernyshev et al., 2002; Milukov et al., 2018] the Aragats center, one of the largest Quaternary volcanic centers in the Caucasus, is confined to the Aragats neovolcanic area located in the western part of Armenia, at the intersection of tectonic zones of a general Caucasian extension and the sublongitudinal Transcaucasus uplift. The development of the Pliocene-Quaternary volcanism of the Aragats area is defined by complex late collisional Geodynamics, which is related to global processes of the convergence of Eurasian and Arabian continental plates. Also, perhaps the emergence of a hot spot under the northern slope of the Aragats volcano in tomography image of our study indicates new volcanic developments under this strato-volcano.

Hot spot number 3 in our study covers a wide area such as Garni, Shoraghbyur, Yerevan, Avan salt dome, and Harazdan which introduced from oil and gas resources [e.g., Jrbashyan et al., 2001]. The Paleocene and Lower Eocene of the subthrust section yielded oil-saturated cuttings and oil-cut mud and oil-stained cuttings were reported as features of the upper Eocene section in these area (thermogenic chemical interactions). By [Milanovski, 1962] about this hot spot, detailed explanation is given as Sevan and Central Troughs. Due to the chemical interactions of in-earth materials in oil-rich areas (contains hydrocarbons) and gas resources, the temperature inside the earth is high and spots relevant with gas plumes is anomalously hot compared to the surrounding. The depth of this hot zone varies from 6.6 to 13.66 km. On the other hand, according to [Kagramanov et al., 2001], there are only a few oil and gas indications in the Sevan Trough. Gas shows and cuts of oil encountered during drilling on the southwestern shore of Lake Sevan (Yeranos-1) suggest that the molasse is prospective for hydrocarbons and country energy in the future.

As well as, in central Armenia, two Paleocene-Miocene flysch-molasse troughs, Central and Sevan, have been identified and several oil and gas shows have been reported from these troughs in Paleogene and Neogene strata. The Shoraghbyur High in the south-eastern part of the trough was the structural feature where shows were recorded from several intervals in the Shoraghbyur-1, where included an oil flow from a Paleocene horizon (3,474 to 3,589 m) and oil-cut mud from 3,640 to 3,634 m. The Paleocene and Lower Eocene of the subthrust section yielded oil-saturated cuttings. Oil-cut mud and oil-stained cuttings were reported as features of the upper Eocene section in Garni-1. Gas shows are reported in a number of wells over the Hrazdan structure. The Avan salt-dome area, between Hrazdan and the Vokchaberd Plateau, had oil shows in the Middle Miocene section [e.g., Klett, 2013]. Therefore, the low-velocity anomaly in the Central Armenia Block in our study is located around these oil-gas areas (Yerevan, Sevan, Gavar, Martuni cities) and should be taken into consideration for the energy requirements of the country in the future.

Figure 7. Geological units, faults, volcanoes, rivers, and basins that used for interpretation the 2D tomography maps in the study area. Abbreviations: F= Fault, Sab.= Sabalan, Sah.= Sahand, Sup.= Suphan, Ara.= Ararat, PT= Pontide, BM= Bitlis Massif, EAAC= East Anatolian Accretionary Complex, Arg.= Aragats, El.= Elbrus, Kaz.= Kazbek, Baz.= Bazarduzu, Ya.= Yanardag (natural gas fire on a hillside), Na.= Nakhchivan, GO.= Georgia, AR.= Armenia, AZ= Azerbaijan, LV= Lake Van, LU= Lake Urmia, LS= Lake Sevan, MD= Mingachevir Dam, TAL= Talesh, KB= Kura Basin, RB = Rioni Basin, SCB= South Caspian Basin, EBSB= Eastern Blake Sea Basin, PSSF= Pambak-Sevan-Syunik Fault, and NTF= North Tabriz Fault. The seismic sources (faults) of the Caucasus are retrieved from [Adamia et al., 2011].

Figure 8. Rayleigh wave group velocity tomography maps for short-periods (5, 10, 15, 20, 25 s); medium-periods (30, 35, 40, 45, 50 s), and long-periods (55, 60, 65, 70 s) used in this study. The white lines are faults and the gray triangles are volcanoes and mounts.

As mentioned earlier, the hot spot number 4 is situated in the EAAC, approximately beneath Ararat strato-volcanic structures (depth= ~6.6-13.66 km, velocity= 2 km/s) in the Julfa region, which could be due to the presence of a magma chamber beneath this volcanic complex. East-north foothills of this volcanic complex are affected by sediments of Aras river, which is limited by the uplifted basement of the Ararat volcanos to the south and by the Hrazdan Trans verse Fault Zone to the west.

Also, some seismic tomography study indicates a magma reservoir at great depths (20-30 km) below the Ararat volcano [e.g., Karaoğlu, 2017]. Geochemical constraints on some of the later-formed rocks suggest an interaction between a shallow chamber (8-10 km) and the deep reservoir approximately 0.5 Ma. This depth is consistent with the result of our study in period of 5 s (depth= 7-13 km; Fig. 8). Although, slow velocities under regions of active volcanism do not require hot spot magmatism.

Low-velocity (dark red area) that we assume as the hot spot number 5, is approximately located in the northeast of the Lake Van includes Tendurek, Suphan and Nemrut mountains. Study of [Vural Oyan et al., 2018] shows collision related to Quaternary Mafic Volcanism to the north of Lake Van has been occurred by eruptions from both volcanic centers and extensional fissures trending approximately north-south. Low-velocity anomaly and the existence of the hot springs around these mountains signify high temperature rocks. Also, the volcanic products in this area consist of mildly alkaline lavas and calculations based on crustal temperatures and Curie point depths indicate that the magma chamber might have been located at a depth of around 6-8 km, within the upper crust.

For the period of 5 s in our study, this property has been shown at a depth of 6 to 13.66 km. We infer, that perhaps, the molten material beneath the Ararat volcano complex and the mountains around Lake Van are quite interconnected. As well as, the pattern of concepts of hot spots 1, 4, and 5 is almost the same, as these spots are situated in the EAAC which is known for its thin and hot (shallow) lithosphere structure [Angus et al., 2006]. So, there is a possibility of the hot rocks and the high temperature inside the earth of this region as a hot spot.

The hot spot number 6 and 7 are located in NW Iran near the north part of Sahand volcano and southeastern segment of the north Tabriz fault. This fault is responsible destructive earthquakes in Tabriz (e.g., M= 7.7, 1721). The epicenter of this earthquake, is located right into the hot zone of number 6. It is clear that these hot spots are perhaps due to the interactions of the rocks of this famous active fault or magma chamber beneath the Sahand volcano. In a study has reported limited volcanic eruptions evidences in the South of the Tabriz fault (Sahand block) that are characterized by ages ranging from 11 Ma to present (era 4). The 11 Ma lavas have an alkaline potassic to ultrapotassic composition [Aghazadeh et al., 2010]. Also, existence the low-velocity zones beneath the NW Iran volcanoes in our results at period of 5 s could be due to the high temperature of the volcanic rocks or shallow magma chamber beneath this region. The existence of the hot springs around these volcanoes signify this high temperature rocks. In contrast, beneath the Sahand volcano a high-velocity zone is observed that could be due to the low temperature volcanic rocks or a deeper magma chamber at a depth of ~30.8 km.

Hot spot number 8 is located near a segment of Salvard fault in the north east of Nakhichevan (Fig. 9). According to several studies [e.g., Danelian et al., 2014; Sokolov, 1977], exposures of Jurassic sequences are located in Nakhichevan and in Iran, where a 500 m -thick Lower and Middle Jurassic sedimentary sequence overlies Upper Triassic strata. Lower Cretaceous deposits are absent on the south Armenian block and the Triassic-Jurassic deposits are unconformably overlain by Cenomanian reefal limestones that are covered by marls. Upper Devonian (the fourth period of the Paleozoic era) and Permian (the fifth period of the Paleozoic era) rocks could be petroleum source rocks [Sosson et al., 2010]. Silurian and Lower and Middle Devonian marine clastic and carbonate rocks crop out in Nakhichevan and are presumed to be present in Armenia. Our study has shown this property in a period of 5 s at a depth of 6 to 9.5 km (hot spot number 8). Presence of petroleum and carbonate rocks resources and chemical-thermal activity related to it can be a sufficient reason for the hot temperature inside the earth in this area as a hot spot.

Tomography velocity map of period of 10 s shows the hot spot number 9, 10, 11, and 12, which are located in the Greater Caucasus. These hot spots are including the eastern Black Sea basin (EBSB) and a segment of the Odishi fault in the Rioni basin (number 9); the Terek basin in Russia and NSW of Kazbek neovolcanic center (number 10); the Chatma region in east of Georgia (number 11), and South Caspian Basin (number 12). These hot spots follow the same pattern as described for the hot spots 1 to 8.

Figure 9. Labeled areas number 1 to 15 are the major low-velocity with dark red shade (slow) that we assumed as hot spots and the rest of the areas with dark blue-green-yellow shades are high-velocity (fast) or cold spots, which represents the remnants of an old tectonic plate that has sunk beneath the Earth's plates in our study. Greater Caucasus hot spots are shown with the number 13. The white lines are faults and the gray triangles are volcanoes and mount.

As mentioned, the low-velocity anomaly beneath the volcanoes in the depth associated with this period (14 to 28.66 km), reveals the presence of magma and the magmatic reservoir and sediments. Although, slow velocities under regions of active volcanism complexes do not require hot Spot magmatism. A broad low-velocity zones are observed in the SCB-Kura basins, and Baku region (hot spot number 12), which according to a study by [Bochud, 2011], the presence of abundant major oil and gas fields (hydrocarbons) in the Baku-Kura region could be the reason for the low surface wave velocity in this area. Sediments reaches 5-7 km in more in part of the -Rioni-Kura foredeep basin [Jrbashyan et al., 2001]. The Black Sea is generally thought to have a basement of oceanic crust that is overlain with 10-20 km of sediment. Similarly, the basement of the Southern Caspian Sea basin has geophysical attributes like that of thick oceanic crust and is overlain by ~20 km of sedimentary cover [Mangino and Priestley, 1998]. The presence of abundant major oil and gas fields (hydrocarbons) in the Baku-Kura-Terek region could be another reason for the low surface wave velocity [Martin Bochud, 2011].

Dark red Low-velocity area which we assume as the hot spots number 13, are include the Elbrus volcanic complexes, Dzirula Massif, Kazbek neovolcanic center, and Yanardag (natural gas fire on a hillside) in the Greater Caucasus. The low-velocity anomaly beneath the volcanoes in the depth associated with ~30.8 km and 158.6 km in these periods reveals the presence of magma and the magmatic reservoirs and mantle plume. In our tomography map with long-period of 70 s (equivalent to a depth of 158.6 km); Azerbaijan, Kura-South Caspian basins, and Talesh heights are covered with high-velocity, which we suggest cold lithosphere roots for deep areas. On the contrary, the low-velocity in the Greater Caucasus, eastern Black Sea basin and EAAC are resulting in very thin (shallow) lithosphere and hot asthenosphere (in great depth). The hot spots in period of 15 s follow the pattern described in periods of 5 and 10 s.

The results at period of 40 s (depth= ~100 km) show a low-group velocity for most parts of NW Iran, Talesh heights, South Caspian Basin (SCB), and Astara region, that could be due to the warm upper mantle in these regions. We determined the hot spot number 14 in NW Iran, where, according to some studies [e.g., Sugden et al., 2018], the mid-lithosphere magma source has a distinct composition compared to the base of the lithosphere, that is argued to be the result of the increased retention of metasomatic components in phases such as apatite and amphibole, that are stabilized by lower temperatures prior to magma generation. Also, partial melts of the deep lithosphere ~120 km (in our study 111 km) and mid-lithosphere sources to give a composition intermediate between magmas from the northern Lesser Caucasus and NW Iran could be the reason for this extensive hot spot and in addition, leads our mind to onset the low velocity zone (LVZ, depth of ~104, Fig. A2) discontinuity.

A wide very low-velocity anomaly at period of 45 s is observed in Sahand, Sabalan, Bitlis, Nakhchivan, Astara, South Caspian Basin, lake Urmia, and especially in the South Armenia Block (SAB). We infer that the hot spot number 15 is situated in the South Armenia Block (SAB) and as the hot spot number 15 was determined.  About the low-velocity in Nakhchivan, Syunik-Gegham-Vardenis highlands or SAB, we interpret that due to the upper Devonian (the fourth period of the Paleozoic era) and Permian (the fifth period of the Paleozoic era) rocks, these rocks could be petroleum source rocks. Silurian and Lower and Middle Devonian marine clastic and carbonate rocks crop out in Nakhichevan and are presumed to be present in Armenia [e.g., Sosson et al., 2010] (also, see the explanation of hot spot number 8).  

The Gegham volcanic group in Armenia also match with the location of the low-velocity anomaly. Besides this, in a study of [Sugden et al., 2018], diagram of depth vs. temperature of melting shows that after the depth of ~80 (1100 ℃) and 150 km (1400 ℃), the temperature has a significant increase in Gegham, Syunik, and Vardenis. So, we propose the interactions and intrusion of very hot molten material from the upper mantle into the asthenosphere-lithosphere discontinuities for creating hot spot number 15 is not unexpected.

The rest of the areas with dark blue-green-yellow shades are cold spots, which represents the remnants of an old tectonic plate that has sunk beneath the Earth's plates.

The results for the long-periods are different, and in this case (deep areas with depth= ~180 km), the Rayleigh waves are more influenced by the velocity structure of the uppermost mantle. We argue at these periods the low-velocity anomalies are mainly due to the absence (thinning) of lithospheric mantle or a thin mantle lid, while high-velocities can be related to the presence of a stable continental mantle lid or of an oceanic-like lithosphere. Our interpretation is in a good agreement with some studies [e.g., Zabelina et al., 2016; Koulakov et al., 2012].

Seismic waves pass through the lithosphere-asthenosphere very slowly and wave velocity reduction from lithosphere to asthenosphere, could be caused by the presence of a very small percentage of melt in the asthenosphere. The upper mantle low velocity zone (LVZ) is a depth interval with slightly reduced seismic velocity compared to the surrounding depth intervals. The zone is present below a relatively constant depth of 100 km (in our study depth of 104 km, Fig. A2) in most continental parts of the world [Thybo 2006]. And the LVZ, extends from about 65 to 220 km depth in the ocean basins [Presnall et al., 2011].

Therefore, the hot wide area is not unexpected. Figure A2 has depicted the approximate depth of Moho, LAB, and LVZ for this study. In the Lesser Caucasus, there is the link between the volcanic manifestations and low-velocity patterns, but it is not as clear as in the Great Caucasus. The result of some studies also confirm our interpretation [e.g., Kearey et al., 2009; Condie, 1997].

At long-period (55-70 s and approximate depth = 200 km), velocity structures of our tomographic maps indicate ultrahigh-velocity anomalies (5.04 km/s) and ultralow-velocity (1.4 km/s) zones (Fig. 8). We infer the deep ultrahigh-velocity anomalies may be the broken off cold lithosphere generated slabs were sinking into the mantle transition zone and very-hot upper mantle with low mantle lid (cap). In contrast, for ultralow-velocity regions, is thought that the upper mantle has been rejuvenated by a phase of the upwelling hot mantle, and this metasomatic refertilization of the upper Cratonic mantle has increased its density and reduces seismic velocity and rocks experience temperatures above 1300-1600 ℃ at these depths.

Also, accordance with depth-temperature diagram [e.g., Sugden et al. 2018], at these depths, some interactions such as onset of dry melting in the convecting mantle cause an increase in temperature and density conflicts tension. So, we propose that at the depths joint between the lithosphere-asthenosphere-upper mantle anomalies accumulation, inhomogeneities, and antagonistic behaviors are prevalent in surface wave velocity variations. In these regions, due to continuous changes in temperature caused by the plate tectonic activity, the effect of active liquids penetrated by the asthenosphere, subsidence, uplifts, hot asthenospheric diapirs intrusion, the velocity of waves experiences many fluctuations [e.g., Sugden et al. 2018].

Also, as mentioned, seismic waves traverse slowly from the lithosphere-asthenosphere boundary (LAB), which known as the low-velocity zone (LVZ), and then enter the upper mantle (Fig. A2). Poor coverage of ray paths in this part of the study area (T= 60, 65 and 70 s) leads to stretching and smearing (butterfly-shaped areas) by this feature toward the northwest and southeast of the study area. In these periods, due to some reasons such as plates tectonic activities, hot asthenospheric diapirs intrusion, the effect of active liquids penetrated by the asthenosphere, subsidence, uplifts and mantle plumes the temperature changes constantly and therefore, the surface waves have variable behavior.

6. Discusion

Hot spots and related structures in the mantle to understand the dynamics of Earth and the modes of heat transfer inside the planet is very important and there are evidences for subduction or underplating crust in these regions. In here, we’re looking for cold and hot spots inside the Earth. Therefore, using the 2D tomography technique and increasing-decreasing the wave velocity anomalies in different geological units of the study area, hot and cold spots were determined.

In other words, using the commentary seismic tomography images results from within the earth, we are looking for the following results in the period of 5 to 70 seconds: 1. Dark blue-green-yellow shades mean colder and stiffer rock (Cold Spots- areas with fast wave velocity), which are the remnants of an old tectonic plate that has been subducted underneath the Earth’s plates (large cold and aseismic area during million years). 2. Dark red-orange shades mean warmer and weaker regions (Hot Spots- areas with slow wave velocity).

So, our tomography maps show the hot zones with dark red shade, where, there is diapirs intrusion, mantle plumes, chain of volcanoes, and fault activity. And the maps with dark blue-green-yellow shades show the cold zones, where, oceanic plates have sunk into the earth’s interior in the past, and supporting the idea that dense slabs of oceanic crust may penetrate to the lower mantle (underlie areas).

In the study area, the highest number of hot and cold spots are corresponding to the short-periods (T= 5 to 20 s) and approximately 13 major hot spots (6.6 < depth < 30.8 km; 1.4 < velocity < 2 km/s) for different geological units of the Caucasus were identified. The location of these hot and cold spots are in good agreement with the results of mentioned tomographic studies in this region.

Since the geothermal resources are shallow-crustal phenomena (a few km above the crust), so these short-period hot spots, are located in an appropriate depth of shallow area of the earth's crust as geothermal energy for humans. In fact, this is the first study to determine and interpret the hot-cold spots of the Caucasus region using Rayleigh surface wave velocity, and so far, no study has been conducted that directly examines the hot-cold spots in the Caucasus by using decreasing and increasing the surface wave velocity.

As well as, determined hot spot and cold spot (15 regions) follows physical processes within the earth and magma generation in Fig. A1 steps. Based on tomographic maps in periods of 5 to 70 s and geological evidences, 15 hot spots with dark red shades in the study area were determined and analyzed and the rest of the areas with dark blue-green-yellow shades are cold spots, which represents the remnants of an old tectonic plate that has sunk beneath the Earth's plates.

As an overall result for three groups of short, medium and long periods in the Caucasus, tomography maps with distinct velocities for the short-periods (T= 5-25 s; equivalent depth of 6.6-53.33 km) are more sensitive to the structure of the upper-middle crust. Based on the reduction and increasing the surface wave velocity, these periods represent sediments in the basins, hydrocarbon resources, molten material and magma chambers beneath volcanoes and Moho discontinuity which these areas can be considered temporary and unstable hot and cold spots inside the earth.

Of course, as mentioned, slow velocities under regions of active volcanism complexes do not require hot spot (i.e., plume-related) magmatism. These could simply reflect decompression melting or crustal melting following slab-rollback, delamination, or breakoff.

In tomographic maps with medium-periods (T= 30-45 s; equivalent depth of 68 to 108 km) the Azerbaijan, Kura-South Caspian basins, and Talesh heights are covered with high-velocity, which according to the results of the previous studies this condition is due to the cold lithosphere roots. On the contrary, due to the presence of a very thin lithosphere and hot asthenosphere, the low-velocity anomaly in the Eastern Greater Caucasus, Eastern Black Sea Basin, Eastern Anatolia, and NW Iran are observed, that this can be considered temporary and semi-stable hot and cold spots inside the earth.

In long-period tomographic maps (50 ≤ T ≤ 70 s; depth of ~180 km), there are the ultrahigh-velocity anomalies (dark blue) under the South Caspian Sea-Kura Basins, Baku, and even has been spread to the Talesh heights, NW Iran, and the Bitlis Massif, while, a wide area is covered with ultralow-velocity (dark red). These deep ultrahigh-velocity anomalies may be due to the broken off cold lithosphere generated slabs were sinking into the mantle transition zone and very-hot upper mantle with thin mantle lid (cap). Whereas the subduction system (e.g., Eastern Greater Caucasus and South Caspian Sea basin) and asthenosphere with significant amounts of melt (e.g., Armenia, Georgia, EAAC, and Western Greater Caucasus) is the major factor in creating the ultralow-velocity zone, which these areas can be considered permanent and stable hot and cold spots inside the earth.

7. Conclusions

In this study, we have performed 2D tomography maps of Rayleigh wave for the entire Caucasus using developed method by Yanovskaya-Ditmar. The derived 2D tomography velocity anomaly maps of Rayleigh wave dispersion curves, were carefully verified using fundamental mod of vertical component (Z) of earthquake waveform energy in order to identify hot-cold spots, better understand the regional tectonic activities, faults activities, and lithospheric blocks interactions as geothermal resources and surface waves velocity variations in the ongoing collision-compressed edge zone of the Eurasian-Arabic plates. These maps show excellent agreement with many of the geological features of the Caucasus territory, such as Volcanoes Complex, Troughs, Uplifts and Basins.

Also, 15 low-velocity area (hot spots) were identified that the Hot Spots Number 10, 11 and 13 (in Greater Caucasus); Number 1, 4 and 5 (EAAC); Number 2, 3 and 15 (Lesser Caucasus-Armenia); Number 6, 7 and 14 (NW Iran); Number 12 in the South Caspian Sea Basin; Number 9 (Rioni-Eastern Black Sea Basin) and Number 8 are located NW of Nakhchivan and is in good agreement with previous studies and geological evidences. The rest of the areas with dark blue-green-yellow shades are cold spots, which represents the remnants of an old tectonic plate that has sunk beneath the Earth's plates. The hot spots close to the earth's surface (beneath our feet) during the short periods can be considered as geothermal resources to provide the heat energy of cities and power plants (e.g. hydrothermal of Sabalan volcano in Iran-Ardabil and Iceland-Nesjavellir Geothermal Power Plant).

Acknowledgments    

Thanks to the Iranian Seismological Center, Incorporated Research Institutions for Seismology (IRIS) Data Management Centre, Iranian National Seismological Network, National Strong-Motion Network of Turkey (TR-NSMN), National Seismic Network of Turkey (DDA), National Seismic Network of Azerbaijan and National Seismic Network of Georgia who provided some of the seismic data used in this study. Many of the figures in this paper were prepared using GMT (Wessel & Smith 1995), which thank them. We also thank the Institute of Geophysics of the University of Tehran, which managed and supported this study in the framework of the educational mission (2016-2018). Special thanks to Michael Smith (USA) & Natasha Lewis (USA) for editing the grammar of the article.

Appendix

Figure A1. (a) Schematic diagram showing the physical processes inside the Earth that lead to the generation of magma that partial melting begins above the fusion point b) Shows a cross section through the Earth's lithosphere (yellow) with magma rising from the mantle (red). Figure retrieved from Christian Schluchter [2010]. The graphs show the geothermal curves (the temperature curve inside the Earth, red) and the solidus (temperature where rock starts to melt, green). When the two curves cross each other, magma is generated by partial melting. A) The curves do not cross- no magma is generated B) at mid-ocean ridges magma generation occurs at quite shallow depths due to high temperatures and very thin lithosphere C) over mantle plumes magma generation occurs at larger depths due to even higher temperatures, but thicker lithosphere D) over subducting slabs magma generation occurs at larger depths due to lowering of melting temperature of the rock by fluids released from the slab.

Figure A2. Velocity variations vs. depth in different periods and approximate location of Moho, Lithosphere- Asthenosphere Boundary (LAB) and Low Velocity Zone (LVZ) respect to the group velocity in this study.

References

Abrehdari S. H. et al. (2022), The ray paths distribution significant role in the resolution of tomographic images, Biomed J Sci & Tech Res 43(4)-2022. BJSTR. MS.ID.006938, DOI: 10.26717/BJSTR.2022.43.006938.

Aghazadeh Mehraj et al. (2010), The gabbro (shoshonitic)-monzonite-granodiorite association of Khankandi pluton, Alborz Mountains, NW Iran, DOI: 10.1016/j.jseaes.2010.01.002.

Angus D. A. et al. (2006), Lithospheric structure of the Arabian and Eurasian collision zone in eastern Turkey from S-wave receiver functions, Geophys. J. Int., doi: 10.1111/j.1365-246X. 2006.03070.x

Backus G.E., Gilbert F. (1968), The resolving power of gross earth data, Geophys. J. 16, 169-205.

Bavali K. et al. (2016), Lithospheric structure beneath NW Iran using regional and teleseismic travel-time tomography, Physics of the Earth and Planetary Interiors, 253, 97-107.

Bochud M. (2011), Tectonics of the Eastern Greater Caucasus in Azerbaijan, Corpus ID: 129649434, Geology.

Chernyshev I.V. et al. (2002), Isotopic geochronology of Quaternary volcanic eruptions in the Greater Caucasus.

Copley A., Jackson J. (2006), Active tectonics of the Turkish-Iranian Plateau, Tectonics, 25, TC6006.

Danelian T. et al. (2014), Reconstructing Upper Cretaceous paleoenvironments in Armenia based on Radiolaria and benthic Foraminifera, http://dx.doi.org/.10.1016/j.palaeo.2014.03.011.

Ditmar P. G., Yanovskaya T. B. (1987), Generalization of Backus-Gilbert method for estimation of lateral variations of surface wave velocities: Phys. Solid Earth, Izvestia Acad. Sci. USSR, 23 (6), 470–477.

Duru Olgun et al. (2020), Magmatic Evolution of the ALADAĞ Volcanic System and Southern Edge of the Erzurum-Kars Volcanic Plateau (SARIKAMIŞ, City of Kars, NE Turkey), Bibcode: 2020 EGUGA. 2221626D.

Foulger G.R. (2010), Plates vs. Plumes: A Geological Controversy; Wiley-Blackwell.

Gheitanchi M.R. (1996), Crustal structure in Nw in Iran, revealed from the 1990 Rudbar aftershock sequence, J. Earth Space Phys., 23, 7-14.

Golonka J. (2004), Plate tectonic evolution of the southern margin of Eurasia in the Mesozoic and Cenozoic, Tectonophysics, 381, 235-273.

Hatzfeld D., Molnar P. (2010), Comparisons of the kinematics and deep structures of the Zagros and Himalaya and of the Iranian and Tibetan Plateaus and geodynamics implications, Rev. Geophysics., 48, RG2005, doi:10.1029/2009RG000304.

Herrmann R. B. (2013), Computer Programs in Seismology: An Evolving Tool for Instruction and Research, Seis. Res. Lettr., 84, https://doi.org/10.1785/0220110096.

Ismail-Zadeh Alik et al. (2020), Geodynamics, seismicity, and seismic hazards of the Caucasus, https://doi.org/10.1016/j.earscirev.2020.103222.

Jackson J. (1992), Partitioning of strike-slip and convergent mountain between Eurasia and Arabian in Eastern Turkey and the Caucasus, J. Geophys. Res., 97, 12471-12479.

Jrbashyan R.  et al. (2001), Geology of Meso-Cenozoic Basins in Central Armenia, with Comment on Indications of Hydrocarbons, Search and Discovery Article #30007.

Kagramanov Yu.R. et al. (2001), www.searchanddiscovery.com/documents/armenia/.

Karaoğlu Ozgür et al. (2017), Volume lava flows fed by a deep magmatic reservoir at Ağrı Dağı (Ararat) volcano, Eastern Turkey, DOI: 10.1007/s00445-016-1098-0.

Klett T.R. (2013), Geology and Assessment of the Undiscovered, Technically Recoverable Petroleum Resources of Armenia, U.S. Department of the Interior, U.S. Geological Survey.

Koulakov I. et al. (2012), Nature of orogenesis and volcanism in the Caucasus region based on results of regional tomography, https://doi.org/10.5194/se-3-327-2012.

Mangino S., Priestley K., (1998), The crustal structure of the southern Caspian region, https://doi.org/10.1046/j.1365-246X.1998.00520.xCitations.

Milanovski E.E. (1962), Sevan Basin. In: Geology of Armenian SSR, v. 1, Geomorphology, Ar. Ac. Sci. Publ., Yerevan, p. 115-136 (in Russian).

Milukov Vadim et al. (2018), Contemporary State of the Elbrus Volcanic Center (The Northern Caucasus), DOI: 10.1007/s00024-017-1595-x.

Presnall Dean C. et all. (2011), The Source of Oceanic Volcanism-Plumes vs the Low-Velocity Zone, Journal of Petrology, doi:10.1093/petrology/egq093.

Randolph J. Martin et al. (2007), Caucasus Seismic Information Network: Data and Analysis, New England Research Institutes.

Ritzwoller M.H., Levshin A.L. (1998), Eurasian surface wave tomography: group velocity, J. geophys. Res., 103, 4839–4878.

Ryan C. Porter, Suzan van der Lee, Steven J. (2019), Whitmeyer Synthesizing EarthScope data to constrain the thermal evolution of the continental U.S. lithosphere, Geosphere) 15 (6): 1722–1737, https://doi.org/10.1130/GES02000.1.

Schluchter Christian (2010), University of Bern, Institute of Geological Sciences.

Sengör A. M. C. et al. (2003), East Anatolian high plateau as a mantle supported, north-south shortened domal structure, Geophys. Res. Lett., 30, 8045, doi:10.1029/2003GL017858.