# Angular velocity of tectonic plates determined by non-rigorous combination of space geodetic techniques.

1. INTRODUCTION

Plate tectonics is a theory which describes the large scale motions of Earth's lithosphere. The theory builds on the older concepts of continental drift, developed during the first decades of the 20th century by Alfred Wegener, and seafloor spreading, developed in the 1960s. The lithosphere is broken up into what are called tectonic plates.

Earth scientists have identified eight major and several minor plates. The plates are in motion either with respect to the fixed mantle or to a chosen reference plate.

We may distinguish between three types of plate boundary zones: convergent or collisional boundaries, divergent boundaries, also called spreading centers, and transform boundaries. The latter movements of the plates are about 10-100 mm annually.

In addition to these individual plate movements there are rotations of the whole Earth's body with respect to the inertial and Earth's fixed system: Firstly, it rotates about its axis in the inertial celestial system (GCRS). This motion is called precession-nutation and is due to the action of the Sun, Moon and other celestial bodies. The International Astronomical Union recommends using the model of precession IAU2006 (Capitaine et al., 2003) and nutation IAU2000 (Mathews et al., 2002). Secondly, the orientation of the axis with respect to the Earth's body is also changing. This so-called polar motion is mainly affected at long time scales e.g. by changing distribution of land-masses and by displacements of water and air-masses on short time scale. Finally, the Earth rotates around its axis with irregular velocity that is often described as a difference between Universal time and Atomic time.

Thus, the orientation of the Earth body can be described by five parameters, Earth orientation parameters (EOP): two coordinates, [x.sub.p], [y.sub.p], of the pole with respect to the International Terrestrial Reference Frame (ITRF), the Earth Rotation Angle (ERA) that is given as a linear function of Universal Time (UT1) and two components of the celestial pole offset, dX, dY, which denote the observed corrections to the adopted precession-nutation model.

Monitoring of EOP and movements of tectonic plates are provided by modern geodesy techniques: Global Position System (GPS), Very Long-Baseline Interferometry (VLBI), Satellite and Lunar Laser Ranging (SLR, LLR) and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS). Only VLBI and LLR techniques can connect the terrestrial system with the celestial one, consequently only these techniques provide data associated with precession and nutation.

Each technique has several analytical centers, associated to the corresponding service, which generates an intra-technique combined product, primarily EOP and station coordinates. The best way how to get the most representative EOP is to combine all solutions together, respecting the advantage of each technique by appropriate weighting.

We used "non-rigorous approach" that derives an appropriate solution by combining the results of individual techniques omitting co-variances, i.e. omitting interrelations between the input parameters, which are treated as independent.

2. NON-RIGOROUS COMBINATION METHOD

The basic idea of this method that was proposed by Kostelecky and Pesek (2006) is based on combining station position vectors in the celestial reference frame. It has been improved a few times within the past years. Below is a brief description of the applied method. For a more detailed description the reader is referred to (Stefka et al., 2009). The transformation from ITRS to GCRS (i.e. [x.sub.T] [right arrow] [x.sub.C]) was modified to serve our purpose:

[x.sub.C] = Q(t)[R.sub.3](ERA)[R.sub.3](-s')[R.sub.1]([y.sub.P])[R.sub.2]([x.sub.P])R(p) [x.sub.t]0000, (1)

where Q(t) is the precession-nutation matrix, s' shifts from ITRF x axis to the terrestrial non-rotating origin, TIO, along the intermediate equator, and [x.sub.P], [y.sub.P] are coordinates of the CIP pole. [R.sub.i] is the matrix of rotation along the i-th axis. Finally, R(p) is the matrix of the seven-parametric transformation, which is considered be a linear function of time.

Since satellite techniques do not provide precession and nutation, the method treats them as being known functions in time and so only the coordinates of the pole and proper rotation are combined. Thus the input data for the combination consists of M sets of EOP [([x.sub.P], [y.sub.P], ERA, X, Y).sub.m] and corresponding sets of station coordinates [([x.sub.T]).sub.m], m = 1,..., M, as derived by analysis centers for individual techniques. Values of UT1-UTC were firstly derived by connection of integrated LOD-estimates provided by GPS and SLR with UT1-UTC measurements from VLBI. Secondly, they were converted into ERA using the standard formula (e.g. McCarthy, 2003).

Partial derivatives of the formula (1) with respect to any unknown, U, yields observation equations of the form:

[summation over (j)] [partial derivative][x.sub.C] / [partial derivative][U.sub.j] d[U.sub.j] = [x.sub.C][|.sub.obs][- x.sub.c][|.sub.0] + r, (2)

where the "observed" vectors [x.sub.C] are calculated from the respective input data and [x.sub.C][|.sub.0] are functions of adopted a priori values of the unknowns and r is residual.

Additionally, two constraints must be added to the system, which is composed of observation equations (2), to stabilize it. The first constraint comes from the smoothing method (Vondrak, 1969 and 1977) and ties the values of the respective EOP at four adjacent epochs. By weighting, a smoothness of the solution was controlled in order to retain as much as 99 % of the signal with periods greater than 5 days. Secondly a no net-rotation is applied which removes the singularity of the system and preserves the system as a whole.

The system of observation equations including all applied constraints is solved using modified Cholesky decomposition proposed by Cepek and Pytel (2005). It solves the sparse matrix of normal equations very effectively.

The unknowns are:

* daily values of [x.sub.P], [y.sub.P]

* daily values of (UT1-UTC)

* values of seven-parameter transformation, p([t.sub.0]), and their time derivatives, [??], for each technique, determined once for a whole period of processing. Station positions change very slowly so that the rates [??] are included only in the case of longer data span, when they can be derived reliably.

3. DATA AND NUMERICAL SOLUTION

The data covers an eight-years period (2000-2008) and was obtained by three types of space geodesy techniques, namely, GPS, VLBI and SLR. The GPS- and VLBI-data have been retrieved from the IERS Combination Pilot Project database (iersl.bkg.bund. de/projects/combination/intra-technique/). The last one, the constrained ilrsb solution was used, as published by ILRS analysis center (ftp://cddis.gsfc.nasa.gov/pub/slr/products/pos+eop/). The GPS and SLR data were weekly SINEX (Solution in Independent Exchange format) solutions, from which the EOP and station coordinates were extracted. The VLBI data consists of per observation singular normal equation matrices. To solve them, the constraints have to be added to tie station coordinates to the VTRF 2005 frame (Nothnagel, 2005) with the a priori precision of 5 mm. Since GPS and SLR do not provide celestial pole offset, only the [x.sub.p], [y.sub.p], and UT1-UTC are combined.

The latter data are successively processed in order to obtain monthly solution of EOP and p. The techniques entered the adjustment with the following weights: 1.44, 0.8 and 1.0 for GPS, SLR, and VLBI, respectively. The values of weights were taken over from the fathers of the present method to keep continuity. Moreover, the mutual ratio of the VLBI to the SLR weight is nearly same as that applied in the IERS Dynamo program (Richard et al., 2008).

From the monthly solutions, final coordinates for each collocation station were obtained as a weighted mean of transformed coordinates of the techniques contributing to the particular solution. Then, biases and linear trends (coordinates and velocities of stations) were computed relative to the epoch J2000. Some comparisons of our solutions with ITRF 2005 have been already published (Stefka et al., 2009 and 2010).

Further, the computed velocities were sorted by their assignation to tectonic plates they are located on. They were compared with the ones published by ITRF 2005 (Altamimi at al., 2007) and depicted as differences, [dv.sub.i], transformed from the terrestrial system to the local system (S--south direction, E--east direction, R--vertical direction) using following formula:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

where [phi], [lambda] are ellipsoidal coordinates of the involved stations. The transformation matrix is obtained by multiplication of two rotation matrices: [R.sub.Z]([lambda]) and [R.sub.Y](90 - [phi]).

The well known formula (4) gives us the mathematical relation between

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

the [??] is angular vector of tectonic plate and [??] and [??], which are vectors of collocation station position and velocity, respectively. Low indices, i and j, indicate the appropriate collocation station and tectonic plate, respectively.

[FIGURE 1 OMITTED]

The equation (3) above was used to formulate the system of observation equations to compute the angular velocity vector of tectonic plates, on which collocation stations are situated. The system had to be constrained by the additional condition:

0 = [v.sub.R] = [[??].sub.1]([x.sub.2]sin[phi] - [x.sub.3]cos[phi]sin[lambda]) + [[??].sub.2]([x.sub.3]cos[phi]cos[lambda] - [x.sub.1]sin[phi]) + [[??].sub.3]([x.sub.1]cos[phi]sin[lambda] - [x.sub.2]cos[phi]cos[lambda]), (5)

where [[??].sub.j] (j = 1,2,3) are components of the angular vector [??]. This constrain assured the primary condition of the tectonic model: total sum of vertical velocities must be equal to zero ([summation][v.sub.R] = 0).

Finally, the introduced system was solved by the least square method. The results as well as comparisons with the original [[??].sub.j] (published by ITRF 2005) are displayed in subsequent tables. The Tables (2 through 8) show the components and their formal errors of the angular velocity vector of the individual plates as computed by the combination method. For comparison the corresponding ITRF 2005 rotation vectors are provided as well.

4. CONCLUSION

The combination method of the results of space geodesy techniques was used to compute the positions and their time derivatives for collocation stations, which are always equipped by at least two of the following techniques (GPS, SLR and VLBI). Station velocities were compared with the ones issued by ITRF 2005 (see Figure 1). The RMS of the comparison is 2.7 mm/y.

Moreover, the station velocities were used to form a system of observation equations with the constraint satisfying the main tectonic model condition. The system was solved by means of the least squared method to compute angular velocity vectors of six major and two minor tectonic plates. The results are shown in detail in Tables (1 trough 8), where the values of angular velocity vectors published by ITRF 2005 are also provided.

By comparing the results, we can see that computed angular velocity vectors of five plates (Australian, European, North American, Nubian, and Pacific) are in good agreement, taking into account the given standard deviations. Other angular velocity vectors of three plates (Okhotsk, South American, and Yangtze) are different. The differences are probably originating from errors in computed velocities of particular collocation stations, which are located on those plates. The errors are over the level of 5 mm/y caused mainly by using not very precise collocation ties that affect combination process. These results could be improved by using the data measured also at non-collocation stations and by including more precise ties in the future.

AKNOWLEDGEMENTS

The author greatly appreciates the support of grant LC506 awarded by Ministry of Education, Youth and Sports of the Czech Republic.

REFERENCES

Altamimi, Z., Collilieux X., Legrand, J., Garayt, B. and Boucher, C.: 2007, ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters, J. Geophys. Res., 112, B09401, doi:10.1029/2007JB004949.

Capitaine, N., Wallace, P., and Chapront, J.: 2003, Expressions for IAU 2000 precession quantities. Astron. Astrophys., 412, 567-586.

Cepek, A. and Pytel, J.: 2005, Progress Report on Numerical Solutions of Least Squares Adjustment in GNU Project Gama. Acta Polytechnica, Czech Technical University in Prague, 45(1), 12-18.

McCarthy, D.D. and Petit, G. (eds.): 2003, IERS Conventions (2003), IERS Technical Note, No. 32.

Mathews, P.M., Herring, T.A. and Buffet, B.A.: 2002, Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior. J. Geophys. Res. 107, doi: 10.1029/2001JB000390.

Nothnagel, A.: 2005, VTRF2005--A combined VLBI Terrestrial Reference Frame. http://miro.geod.uni bonn.de/vlbi/IVS-AC.

Pesek, I. and Kostelecky, J.: 2006, Simultaneous determination of Earth orientation parameters and station coordinates from combination of results of different observation techniques, Stud. Geophys. Geod., 50, 537-548.

Richard, J.Y, Gambis, D., Bizouard, C., Lemoine, J.M. and Cancale, R.: 2008, Global combination of stations coordinates & Earth rotation parameters. Presented at EGU General Assembly held in April 13-18, 2008 in Viena.

Stefka, V., Kostelecky, J. and Pesek, I.: 2009, Combination of different space geodesy techniques for EOP and terrestrial reference frame, Acta Geodyn. Geomater., 6, No. 3 (155), 239-246.

Stefka, V., Pesek, I, Vondrak J. and Kostelecky J.: 2010, Earth orientation parameters and station coordinates from space geodesy techniques, Acta Geodyn. Geomater., 7, No. 1 (157), 29-33.

Vondrak, J.: 1969, A contribution to the problem of smoothing observational data, Bull. Astron. Inst. Czechosl., 28, 84-89.

Vondrak, J.: 1977, Problem of smoothing observational data II, Bull. Astron. Inst. Czechosl., 28, 84-89.

Vojtech STEFKA

Center for Earth Dynamics Research (CEDR), Astronomical Institute of the Czech Academy of Sciences of the Czech Republic, Bocni II, 141 31 Prague 4, Czech Republic

Corresponding author's e-mail: stefka@ig.cas.cz

(Received January 2010, accepted August 2010)
```Table 1 The table shows components
of angular velocity vector of
Australian plate and their standard
deviations, all computed from the
combination results. There are also
in the right part of the table.
There are 4 collocation stations
on the plate.

Australian plate

[[??].sub.computed]   [[??].sub.ITRF2005]
[deg / My]            [deg / My]

0.42128               0.42138
[+ or -] 0.26894
0.33213               0.32179
[+ or -] 0.23046
0.34422               0.33656
[+ or -] 0.25965

Table 2 The table shows components
of angular velocity vector of
European plate and their standard
deviations, all computed from the
combination results. There are also
in the right part of the table.
There are 23 collocation stations
on the plate.

European plate

[[??].sub.computed]     [[??].sub.ITRF2005]
[deg / My]              [deg / My]

0.03002               -0.01507
[+ or -] 0.07065
-0.16651               -0.14391
[+ or -] 0.04218
0.27112                0.21722
[+ or -] 0.08505

Table 3 The table shows components
of angular velocity vector of North
American plate and their standard
deviations, all computed from the
combination results. There are also
in the right part of the table.
There are 9 collocation stations
on the plate.

North American plate

[[??].sub.computed]     [[??].sub.ITRF2005]
[deg / My]              [deg / My]

-0.03125                 0.00874
[+ or -] 0.06381
-0.17341               -0.19126
[+ or -] 0.15479
-0.06494               -0.01437
[+ or -] 0.15151

Table 4 The table shows components
of angular velocity vector of Pacific
plate and their standard deviations,
all computed from the combination
results. There are also components
part of the table. There are 6
collocation stations on the plate.

Pacific plate

[[??].sub.computed]     [[??].sub.ITRF2005]
[deg / My]              [deg / My]

-0.06388               -0.12212
[+ or -] 0.18941
0.31692                0.28948
[+ or -] 0.13131
-0.62717               -0.60532
[+ or -] 0.10503

Table 5 The table shows components
of angular velocity vector of South
American plate and their standard
deviations, all computed from the
combination results. There are also
in the right part of the table.
There are only 2 collocation stations
on the plate.

South American plate

[[??].sub.computed]     [[??].sub.ITRF2005]
[deg / My]              [deg / My]

-0.02640               -0.07388
[+ or -] 0.21893
-0.30772               -0.08921
[+ or -] 0.26798
0.04393               -0.03497
[+ or -] 0.17606

Table 6 The table shows components
of angular velocity vector of Nubian
plate and their standard deviations,
all computed from the combination
results. Although the estimated
standard deviations are quite big,
probably because of only 2 collocation
stations being on the plate, the
computed values fit well with ones
the right part of the table.

Nubian plate

[[??].sub.computed]   [[??].sub.ITRF2005]
[deg / My]            [deg / My]

0.02027           0.02259
[+ or -] 1.35490
-0.18167          -0.17159
[+ or -] 0.62431
0.20327           0.20593
[+ or -] 0.83142

Table 7 The table shows components
of angular velocity vector of
Okhotsk plate and their standard
deviations, all computed from the
combination results. The standard
deviations are too big so that
these results are unrealistic and
will be excluded. There are only 2
collocation stations on the plate.

Okhotsk plate

[[??].sub.computed]   [[??].sub.ITRF2005]
[deg / My]            [deg / My]

-0.81317          -0.04790
[+ or -] 14.82910
0.64348          -0.05153
[+ or -] 12.53821
0.74007          -0.04403
[+ or -] 13.96711

Table 8 The table shows components
of angular velocity vector of Yangtze
plate computed from the combination
2005. One can see that the estimated
standard deviations are quite big
especially one related to second
component of computed angular
velocity vector. It is mainly due
to one of participated stations
(21602M001) whose vertical velocity
differs from the result of ITRF 2005

Yangtze plate

[[??].sub.computed]   [[??].sub.ITRF2005]
[deg / My]            [deg / My]

0.01186           -0.05325
[+ or -] 0.80084
-0.24928           -0.14842
[+ or -] 1.32241
0.19951            0.26690
[+ or -] 0.92527
```