An Analysis of the Eurasian Tectonic Plate Motion Parameters Based on GNSS Stations Positions in ITRF2014

The article is the fourth part of our research program concerning an analysis of tectonic plates’ motion parameters that is based on an observation campaign of an array of satellite techniques: SLR, DORIS, VLBI, and now GNSS. In this paper, based on the International Terrestrial Reference Frame 2014 (ITRF2014) for observations and using the GNSS technique, the Eurasian tectonic plate motion was analyzed and the plate motion parameters Φ, Λ (the position of the rotation pole), and ω (the angular rotation speed) were adjusted. Approximately 1000 station positions and velocities globally were obtained from the GNSS campaign over a 21-year time interval and used in ITRF2014. Due to the large number of data generated using this technique, the analyses were conducted separately for each tectonic plate. These baseline data were divided into a number of parts related to the Eurasian plate, and are shown in this paper. The tectonic plate model was analyzed on the basis of approximately 130 GNSS station positions. A large number of estimated station positions allowed a detailed study to be undertaken. Stations that agree with the plate motion were selected and plate parameters were estimated with high accuracy. In addition, stations which did not agree with the tectonic plate motion were identified and removed. In the current paper, the influence of the number and location of stations on the computed values and accuracy of the tectonic plate motion parameters is discussed. Four calculation scenarios are examined. Each scenario contains 30 stations for the common solution of the European and Asiatic part of the Eurasian plate. The maximum difference between the four calculation scenarios is 0.31° for the Φ parameter and 0.24° for the Λ parameter, indicating that it is at the level of the value of the formal error. The ω parameter has the same value for all the scenarios. The final stage of the analysis is the estimation of parameters Φ, Λ, and ω based on all of the 120 stations used in the four calculation scenarios (i.e., scenario 1 + scenario 2 + scenario 3 + scenario 4). The following results are obtained: Φ = 54.81° ± 0.37°, Λ = 261.04° ± 0.48°, and ω = 0.2585°/Ma ± 0.0025°/Ma. The results of the analysis are compared with the APKIM2005 model and another solution based on the GNSS technique, and a good agreement is found.


Introduction
The earth, due to its dynamics, requires constant observation. Continuous changes occurring inside the earth caused by various geophysical and geological phenomena are reflected on the surface of the globe [1], such as in the form of tectonic plate movements. The theory of the movement of tectonic plates (the so-called continental drift hypothesis) was first published in 1915 by A. Wegener in his fundamental work "Die Entstehung der Kontinente und Ozeane" [2]. In the following years of the 20th century a number of scholars addressed these issues, e.g., geologists O.C. Hilgenberg [3], F.A. Vening-Meinesz [4], S.W. Carey [5], and X. Le Pichon [6], until the adoption by the scientific community in 1968 of the so-called Theory of Plate Tectonics, whose assumptions were described in [7].
During the 1980's and 1990's of the last century a dynamic development of satellite observation techniques (GPS, SLR, DORIS, VLBI) occurred, providing an unprecedented opportunity to conduct research on a global scale with very high accuracy. These techniques were also applied in the study of tectonic plate motion. The accuracy of determining the position of points on the earth's surface using satellite methods has continuously increased, and this increase in accuracy has allowed the numerical values describing the movement of tectonic plates to be determined. Continuity of research in this field is a noticeable trend, e.g., GNSS technique: [8][9][10][11][12]; SLR technique: [13][14][15][16][17][18]; DORIS technique: [10,19,20]; and VLBI technique: [21][22][23][24][25][26].
Initially, geological methods based on the description of the mechanism of tectonic plate motion were used to study the movement of tectonic plates; only the satellite techniques mentioned above made it possible to precisely estimate the numerical values of this movement, as discussed in [27]. Thorough knowledge and the ability to mathematically describe the phenomena causing the movement of the tectonic plates, and the possibility of taking precise measurements using satellite techniques, allows points on the earth's surface and their annual movements to be determined with high accuracy. This is the basis for determining the parameters that describe the movement of tectonic plates. Tectonic plates move in relation to each other on the astenospheric surface at a speed ranging from a few to several centimeters per year. The lithosphere is constantly destroyed in the subduction zones and is renewed on mid-ocean ridges. Continents usually form a section of a plate and are moved with it. The following types of boundaries between the plates exist: divergent, convergent, and horizontally sliding. Issues related to these boundaries have been previously addressed, e.g., [28,29]. Large stresses can be present on the edges of the plates. These are discharged in the form of earthquakes. In the case of tectonic plates, it is difficult to speak of rigid boundaries between plates and to indicate where one plate ends and another begins; at the plate boundaries, there is often an area in which displacements are incompatible with the movement of the whole plate.
The movement of the tectonic plate is described by the rotation vector Ω, whereas the parameters of this movement are described by the geographical position of the rotation pole-Φ and Λ-and the angular velocity of rotation-ω (or by individual components of the rotation pole) [30,31]. According to [30], the displacement of the observational station as a function of the parameters of plate motion is expressed by Equation (1): where: ∆φ, ∆λ-displacement of the observational station position in latitude and longitude; φ, λ-observational station position (GNSS, SLR, DORIS, VLBI); Φ, Λ, ω-plate motion parameters (the position of the rotation pole in latitude and longitude, and the angular rotation speed, respectively).
The aim of this study is to estimate and analyze the parameters describing the motion of the Eurasian plate based on GNSS station positions, and to identify and eliminate from the calculation those stations whose motion is not consistent with the motion of the plate, hence "polluting" the results of the motion parameters and increasing the value of the formal error.

Materials and Methods
There are approximately 1000 GNSS stations located around the world whose positions are provided in the International Terrestrial Reference Frame 2014 (ITRF2014) [32]. From this baseline data, about 130 stations located on the Eurasian continent were analyzed. These stations are not uniformly distributed. Approximately 70% of the stations are located in the European part of the Eurasian plate.
detail. It also allows stations that are not compatible with the Eurasian tectonic plate motion to be identified and removed, and only stations that are consistent with the solution to be selected.
The baseline data were divided into four calculation scenarios. Each scenario contained 30 randomly distributed stations for a common solution of the Eurasian plate. The stations used in the four scenarios of the calculations are tabulated in Appendixes A-D, separately for each scenario, and are shown in Figure  1. Detailed information on individual stations, including their positions and velocities, is presented on the ITRF website: http://itrf.ensg.ign.fr/ITRF_solutions/2014/. The final solution of the calculation was determined from all of the 120 stations (scenario 1 + scenario 2 + scenario 3 + scenario 4). The locations of the GNSS stations that are not consistent with the Eurasian plate motion are shown in Figures 2 and 3. All calculations were carried out using the authors' own software. The description of the method applied in this paper to determine the plate motion parameters, and the theoretical basis for this motion, has been extensively presented by the authors in their earlier papers on the analysis of the plate motion based on movements of observation stations: SLR [33,34]; DORIS [34,35]; and VLBI [34,36]. Thus these aspects were omitted herein.

Tectonic Plate Theory for the Eurasian Plate
The Eurasian plate, which is the subject of this analysis, is the most complex part of the earth in terms of geological structure [37]. On the western, northern, and north-eastern sides, the plate borders the North American plate. In the west and north, this boundary is a divergent one, on which the Mid-Atlantic Ridge was formed. In the east it is a convergent boundary [38]. In the north, the area of Sweden and southern Norway has been subject to continuous uplift movements for about 1.5 billion years to date [39]. The north-eastern part of the Asian continent belongs to the North American plate, whereas the area comprising the Koryak Mountains and the Kamchatka area belongs to the Eurasian plate; however, the course of the boundary itself has not yet been sufficiently studied. The Eurasian plate also borders on the Bering and Pacific plates, which are being pulled under the Eurasian plate, The description of the method applied in this paper to determine the plate motion parameters, and the theoretical basis for this motion, has been extensively presented by the authors in their earlier papers on the analysis of the plate motion based on movements of observation stations: SLR [33,34]; DORIS [34,35]; and VLBI [34,36]. Thus these aspects were omitted herein.

Tectonic Plate Theory for the Eurasian Plate
The Eurasian plate, which is the subject of this analysis, is the most complex part of the earth in terms of geological structure [37]. On the western, northern, and north-eastern sides, the plate borders the North American plate. In the west and north, this boundary is a divergent one, on which the Mid-Atlantic Ridge was formed. In the east it is a convergent boundary [38]. In the north, the area of Sweden and southern Norway has been subject to continuous uplift movements for about 1.5 billion years to date [39]. The north-eastern part of the Asian continent belongs to the North American plate, whereas the area comprising the Koryak Mountains and the Kamchatka area belongs to the Eurasian plate; however, the course of the boundary itself has not yet been sufficiently studied. The Eurasian plate also borders on the Bering and Pacific plates, which are being pulled under the Eurasian plate, constituting is a convergent boundary. Two small plates are positioned between the Pacific and Eurasian plates: the Okhotsk plate, including the island of Sakhalin on which tectonic movements are ongoing, and the Amur plate. These plates are often treated as part of the Eurasian plate because their motion is similar to that of the latter plate [40]. A zone of islands with heterogeneous crust structure extends from Hokkaido in the north to Taiwan in the south. This string lies at the meeting point of the Pacific, Philippine, Okhotsk, and Amur plates. The boundary between the Philippine and Pacific plates is a convergent one, running through the island of Honshu. In Indonesia, eastern Asia joins the Pacific range. The southern boundary of the Eurasian plate with the African, Anatolian, Arabic, and Indo-Australian plates is defined by approximately 12,000 km of the Alpine-Himalayan range, the Alpine fold zone connecting through the Strait of Gibraltar with the Alpine structures of the Atlas mountains, and the Alpine ranges in Asia. This boundary includes the Atlas, the Pyrenees, the Alps, the Apennines, the Carpathians, the Rhodopes, the Caucasus, the Tien Shan, Tibet, the Himalayas as far as Burma, and the Andaman Sea. The Apennines were created in the process of subduing the Tyrrhenian and Adriatic plates. The Caucasus was created as a result of a collision between the Anatolian and European plates. The area is characterized by high seismicity and is still subject to tectonic processes. In this region, the plates overlap and are squeezed and folded, the continents collide, and the phenomenon of orogenesis takes place. Now, this phenomenon occurs at the place where the Indian and Eurasian plates meet. This collision began about 40 million years ago. Since that contact, the Indian Peninsula has moved 2000 km into Eurasia. The consequences of this collision include the creation of the Himalayas, Tibet, and Tien Shan, in addition to earthquakes in China and north-western Iran [41,42]. The processes associated with the formation of the continental structure continue to this day. The Indian Peninsula moves northwards at a speed of 3.3-4.8 cm per year, resulting in a continuous uplift of the Himalayas. The Arabian Peninsula moves to the northeast at a speed of 1.4-1.8 cm per year. The Pacific plate is subdued by the Asian continental lithosphere at a speed of 6.7-7 cm per year. These processes result in strong earthquakes, tsunamis, and volcanic eruptions [43].
For the Eurasian, African, and Arabic plates, the whole area of southern Europe and the Mediterranean area is considered a boundary zone. The lithosphere here is very cracked and forms a large number of micro plates [44]. The African plate is rubbing against the Eurasian and Arabic plates. The area is highly seismically active and the lithosphere is strongly cracked and divided into the Anatolian, Black Sea, Aegean, Macedonian, Adriatic, Corsican, Atlas, and Iberian micro plates. The north-eastern part of Corsica belongs to the Alpine structures and came into existence as a result of collisions with the micro plate of the Tyrrhenian Sea (the Apulian plate) and the African plate.

Results and Discussion
This section presents the results of the estimation of the Eurasian tectonic plate motion parameters based on GNSS station positions in ITRF2014 [32]. Initially, the solution was derived for two stations located on a plate. Then, single stations were added individually (sequential solution) for a maximum of 29 steps for 30 stations (in each of the four scenarios of calculations). The selection of the stations in the scenarios was based on the following criterion: first stations located in the stable region of the plate were included; then, stations located near the plate boundary but which did not disturb the solution accuracy were included. Results for the four calculation scenarios are given in       The stations with velocity directions and values that are not consistent with the Eurasian plate motion are not included in the solution (these are shown in Figures 2 and 3). Inclusion of these stations results in a description of the Eurasian plate motion that is incompatible with the real movement of the whole plate rotation around poles Φ and Λ.
The final stage of the analysis was the estimation of parameters Φ, Λ, and ω based on all of the 120 stations used in the four calculation scenarios (scenario 1 + scenario 2 + scenario 3 + scenario 4). The following results were obtained: Φ = 54.81 • ± 0.37 • , Λ = 261.04 • ± 0.48 • , and ω = 0.2585 • /Ma ± 0.0025 • /Ma. These results were compared with the APKIM2005 model by H. Drewes [45] and with another solution based on the GNSS technique given in Larson et al. [8]. The comparison is shown in Table 2. The APKIM2005 model [45] is based on weekly solutions for SLR, DORIS, VLBI, and GNSS techniques observed for the time interval from 1993 to 2004. The database for these techniques allows estimation of the velocity vectors caused by plate motion and plate motion parameters for the following 17 plates: African, Amur, Antarctic, Arabian, Anatolian, Australian, Caribbean, Eurasian, Indian, Nazca, North American, South American, Okhotsk, Pacific, Somalia, Sunda, and Yangtze. In [8] authors analyzed the GNSS data from January 1991 to March 1996. All of the presented data were analyzed using the GIPSY/OASIS II software. On the basis of the computed network, velocities for 38 sites located on the African, Antarctic, Australian, North American, South American, Pacific, and Eurasian plates were estimated. For the Eurasian plate, the parameters were estimated on the basis of eight stations; six of these are located on the European part of the plate, and two are located on the Asian plate. of eight stations; six of these are located on the European part of the plate, and two are located on the Asian plate.

Conclusions
The computations and results of the analysis given in this paper allow the following conclusions to be drawn:

•
The convergence and stability of the solutions for the four calculation scenarios based on 30 randomly distributed stations for each scenario are obtained for about 18-20 stations, as presented in Figure 4 for latitude Φ, Figure 5 for longitude Λ, and Figure 6 for rotation speed ω.

•
The selection of suitable stations for determining the parameters allows the determination on the basis of approximately 20 stations to be made, as shown in the four calculation scenarios. Adding more stations to the calculations results in a change in the value of the determined parameters by a value that does not exceed the formal error. • The geological model of the Eurasian plate is the most complicated on earth. The Eurasian plate includes a subduction zone in the convergent eastern boundary. Part of the North American plate, on which Bilibino (BILI) station is situated, is located in the eastern section of the Eurasian continent. The annual shift of this station is not consistent with the Eurasian plate motion. This station cannot be included in the solution because of a change in the pole rotation of the plate parameters by approx. 2 degrees in latitude, which corresponds to 200 km, and 1.5 degrees in longitude, which corresponds to 150 km. Magadan (MAGO) station is located on the boundary of the American and Okhotsk plates. The annual shift of the station is not consistent with the Eurasian plate motion by approx. 1 degree in latitude, which equals 100 km, and 1 degree in longitude, which also equals 100 km, therefore this station cannot be used in the solution. Petropavlovsk (PETP) station is located on the Okhotsk plate on Kamchatka Peninsula. The shift of this station is also not consistent with the Eurasian plate motion, and amounts to 4.5 degrees in latitude, which corresponds to 450 km, and 6.8 degrees in longitude, which equals 680 km, hence the station cannot be used in the solution. Teheran (TEHN) station is also not used in the solution because it is located near the boundary of the Iran plate in a high seismic activity region, despite a relatively small change in the pole rotation of plate parameters, by approx. 0.2 degrees in latitude, which corresponds to 20 km, and 0.4 degrees in longitude, which equals 40 km.

•
The Eurasian plate contains areas of high seismic activity and cracked boundaries, both convergent and divergent, with other tectonic plates, such as the Pacific, Deccan, Iran, Arabian, and Anatolian plates. Often, shifts of the analyzed stations located on the Eurasian plate are not compatible with the tectonic plate motion; for example, Lhasa (LHAS) station, which is located on the boundary of the Eurasian plate and the Deccan plate in the South Tibet and Himalaya mountains, is in a very active and cracked region. This station increases the error of the estimated plate parameters significantly, by approx. 3 degrees in latitude, equaling 300 km, and 4 degrees in longitude, equaling 400 km. The station is not included in the solution.

•
The Japanese stations of Mizusawa (MIZU), Koganei (KGNI), Usuda (USUD), Tsukuba (TSKB), and Abashiri (P202) are located in a very active and cracked area on the boundary of the Philippine plate. In this region, three cracked plates come into contact: the Pacific plate with the Okhotsk and Amur plates. Shifts of these stations are not consistent with the Eurasian plate motion and are not included in the solution.

•
Wakkanai station (used in calculation scenario 4) is located on the boundary of the Amur and Okhotsk plates. The shift of this station is consistent with the Eurasian plate motion to a high degree, therefore it is included in the solution.

•
The application of the sequential calculation method allows stations whose movement is not consistent with that of the entire plate to be identified and eliminated from the solution.

•
Analysis and elimination of the selected stations, as shown in Table 1, is indispensable because it allows results to be obtained only on the basis of the stations with shifts that are consistent with the Eurasian plate.

•
The Mediterranean Sea and surrounding areas, i.e., the Anatolian and Arabic plates, are excluded from the analysis in this work. These regions will be analyzed separately because they include a significant number of micro plates, which each has a characteristic motion that differs from that of the Eurasian plate. In the Asian part, the Deccan plate is omitted because it does not belong to the Eurasian plate.