An optimal geodetic dynamic reference frame realization for Greece: Methodology and application

In this study the theoretical and practical approach of an optimal geodetic reference frame realization is described in detail. The Hellenic area due to its strong geodynamic behavior creates a rather inhomogeneous velocity field. This issue plays a critical role in stable and long-term geodetic datum. The methodology of the Minimization of the Kinetic Energy Criterion (MKEC) is defined and analyzed on a permanent GNSS network consisting of 151 stations. The site velocity values were estimated from seven years of continuous data. The application formula shows that the kinematic energy is minimized more than 60 percent in two different approaches for the International Terrestrial Reference Frame and the European Terrestrial Reference System 1998 (ETRS89). In addition, the realization of the ETRS89 approach in Greece is not offering any significant advantage. However, we proposed a new strategy, which minimizes the total kinetic energy and it is found effective in local areas with strong geodynamic activity like the Greece case.


Introduction
South-eastern Mediterranean and especially the Hellenic area show an intense geodynamic behavior.The principal geomorphological features of the Hellenic area [Papazachos et al. 1998] listed from South to North as: (1) Mediterranean ridge (East Mediterranean chain) which extends from the Ionian Sea to Cyprus and Levantine Arabic areas.
(2) Hellenic trench begins from Ionian Islands, continuing at south Crete and east to Rhodes where splits at the Ptolemy, Pliny and Strabo trenches.There detected the deepest point of the Mediterranean Sea (~5 km).
(3) Hellenic arc which includes the outer sedimentary arc (Hellenic island arc) and parallel the inner volcanic arc (southern Aegean volcanic arc), which consists of Santorini, Nisyros, Milos and Kos islands.
Between of these two approximately parallel arcs is the Cretan back-arc basin with depth less than 2 km.Some of the remaining tectonically active features are the Kefallinia Transform Fault (KTF) in the Ionian Sea, the Corinthian Rift (CR) in Central Greece and finally, the North Anatolian fault (NAT) at the northeast part of Hellenic area.Figure 1 shows the major seismogenic sources of the border area of Greece as described in detail at the Greek Database of Seismogenic Sources (GreDaSS) [Caputo et al. 2012].
Since the last decades, various studies focused on tectonics of the Hellenic subduction zone and the geodetic velocity field [Davies et al. 1997, Rossikopoulos et al. 1998, Cocard et al. 1999, Hollenstein et al. 2008, Vernant et al. 2014].The horizontal velocity field characterized by inhomogeneous behavior which reflecting the geophysical complexity of the Hellenic area.Since the early 1970s several studies [McKenzie 1972, Le Pichon 1995, McClusky et al. 2000, Nyst and Thatcher 2004, Reilinger et al. 2006, Floyd et al. 2010] attempted to determine the number and the block boundaries of the microplates in the Aegean Sea.The number of the estimated blocks increases at the most recent studies and varies from two [McKenzie 1972] to fifteen [Floyd et al. 2010] microplates.The block boundaries cannot be defined with high accuracy for geodetic applications and consequently the realization of the local reference frame (based on microplates) it is not feasible.Corinth gulf is one of the most seismically active areas in Europe and located at the collision boundaries between two microplates [Avallone et al. 2004], thus the accurate definition of the borders is cumbersome.
For that reason, alternative methodologies have been already proposed, like the implementation of a se-mi-kinematic datum for Greece, based on seven stable crustal blocks [Chatzinikos et al. 2015].The block definition does not reflect the geophysical structure of the crust, but it is only based on geometrical interpretation by applying the Euler Pole model and statistical criteria.An improvement of the Hellenic velocity field, study using seven years of continuous GNSS data, was estimated and confirmed the inhomogeneity of this specified area [Bitharis et al. 2015].
The main aim of the present study is to propose a new methodology of an Optimal Reference Frame (ORF) realization in a local sense.This new methodology is most beneficial on areas with strong and inhomogeneous velocity field as Greece.This need led us to suggest an alternative technique, which minimizes the kinetic energy of the Reference frame and providing the sustainability, and also the reliably cartographic and geodetic products.

The ITRF approach:
In modern International Terrestrial Reference Frame (ITRF) estimation scheme (from ITRF96 and for the upcoming released versions), the derived coordinates and velocities for the contributed stations are determined from a cumulative solution, combining data from four different space techniques (VLBI, SLR, GPS and DORIS).The latest released versions (ITRF2005, ITRF2008 and recently ITRF2014) are based on long period time series analysis [Altamimi et al. 2007, Altamimi et al. 2011].Velocity determination relies on the estimation of the linear trend of each station's time series analysis, taking into account some rigorous statistical criteria in order to remove the outliers; [Altamimi et al. 2002].Nevertheless, time series analysis also reveals the existence of some other geophysical phenomena [Nikolaidis 2002, Williams 2008], like the ocean loading and atmospheric loading [van Dam et al. 2010] and episodic events (earthquakes).It must also be noticed that the selected stations for the ITRF realization must fulfill some special criteria such as at least 3 years of data and stability in their solution [Altamimi et al. 2002].
Stations coordinates and velocities estimated from each technique are optimally combined by the so called time-dependent Helmert-type transformation; see [Altamimi et al. 2002, Altamimi et al. 2011].From a pure mathematical point of view, time-dependent transformation implemented in an ITRF cumulative solution, has a total rank deficiency of 14.For more details on reference system rank deficiency, see [Davies andBlewitt 2000, Altamimi andDermanis 2012].The rank deficiency refers to the problem defining ITRF's origin, scale and orientation and their rates respectively.To overcome this problem, proper constraints must be imposed.Space techniques do not sense the orientation (and its rate) and so the associated to them deficiency should be solved using another strategy.
Particularly, the rank deficiency for the orientation rate is treated introducing the Tisserand's condition [Tisserand 1889].Tisserand's condition is realized on a network of points fulfilling the criterion of the minimal relative angular momentum and is mathematically expressed through the following relation (Lambeck 1988): (1a) in integral form and in discrete form (which is used in the ITRF practice), (1a) yields: (1b) where h r the relative angular momentum of the network, x i =[x y z] T i the 3D position vector for a point i, the infinitesimal mass of a point i and n the number of discrete points covering Earth's surface.The criterion of minimizing the relative angular momentum implies the realization of a reference frame whose total kinetic energy is minimal (Dermanis 2001).Hence, the implementation of Tisserand's condition does not only stand as an algebraic constraint but at the same time provides a solution with a physical meaning (minimization of the kinetic energy).
The practical implementation of Equation (1b) requires a-priori information for a set of stations velocities, distributed all over the earth's surface.This could be done either using velocities from a global tectonic motion model e.g.NNR-NUVEL1A [DeMets et al. 1994] or AM-02 [Minster and Jordan 1978] or to constrain the orientation rate of the new ITRF to be the same as the exact previous one.The latter treatment is implemented from ITRF96 and for the upcoming ITRF realizations; see [Boucher et al. 1998, Altamimi et al. 2002].It is also worth to mention that the first ITRFs (e.g.ITRF88 and ITRF89) stations velocities were not be estimated through optimal combination of space techniques, but they were derived directly from AM-02 [Boucher and Altamimi 1989].
ITRF's velocity field serves a special role: It can be considered as a conventional choice realizing a reference frame with the minimal kinetic energy trying in the same time to compromise the global geodynamic behavior variation.As a matter of fact, ITRF velocities do not offer an explicit geometrical interpretation like the intraplate velocities (with respect to a major tectonic plate), which are commonly used in geodynamics.

ETRS89 approach:
The European Terrestrial Reference System of 1989 (ETRS89) is created to follow the motion of the stable part of the Eurasian plate [Gubler et al. 1992].The ETRS89 is realized by the so-called European Terrestrial Reference Frame (ETRF YY ).ETRF-based velocities are obtained by reducing the initial ITRF velocities, following [Boucher and Altamimi 2011], point wise: (2) T the 3D position vector with respect to ITRF YY , the 3D vector of velocity in ETRF YY , the 3D vector of velocity in ETRF YY and is the 3x3 anti-symmetric matrix that contains the rotation rates (ṙ x ,ṙ y ,ṙ z ) about each axis of a (clockwise) Cartesian system respectively.Rotation rates are either derived from a global tectonic motion model or by Least Square estimation from properly selected European sites in Central and northern Europe; see [Boucher and Altamimi 2011].The initial reference frame for the ETRS89 realization is the most recent ITRF [Boucher and Altamimi 2011].The ETRS89 establishes a rather stable Terrestrial Reference Frame (TRF) for any geodetic and cartographic use all over Europe.This methodology leads to horizontal velocities determination for the most part of the European continent, at the level of few mm/yr [Gurtner et al. 2008].To be more specific, the ETRS89 explicitly uses intraplate velocities, with respect to the stable part of the Eurasian plate as quantities which describe its own dynamic behavior.The rest regional TRFs, North American Reference Frame (NAREF) and Sistema de Referencia Geocéntrico para Las Américas (SIRGAS) do not follow any particular methodology on their velocity field determination.The velocities (in both cases) referred to an existing ITRF.For more details about NAREF and SIRGAS see Craymer et al. [2007] and Seemüller et al. [2008], respectively.By its fundamental concept, the ETRS89 relies on the assumption that the entire European continent moves uniformly with respect to the Eurasian plate.Actually, this is not true; southern part of Europe and especially the Hellenic area does not seem to be consistent with this assumption [Clarke 1996].This is mainly caused due to area's geodynamic behavior as it was briefly described on section 1.The velocity field is strongly inhomogeneous and there are cases that the horizontal velocity can reach 35 mm/yr with respect to ETRS89 [Bitharis 2015].
Many countries in the European continent adopt an ETRS89 realization for their official reference frame.E.g the CHTRF2010 in Switzerland [Brockmann 2010], is also providing both horizontal and vertical velocities, in Czech republic the S-JTSK/05 [Simek 2011], the ITALREF 2008 in Italy [Baroni et al. 2011], the SWEREF in Sweden [Lidberg et al. 2011], EUREF IE/ UK 2009 in United Kingdom [Greaves et al. 2011] are also aligned to an ETRS89 realization.For Central and northern Europe ETRS89 serves as a stable reference system, providing minimal horizontal velocities.
As far as Greece is concerned, is maintaining (from 2007) a permanent network of GPS called HE-POS (Hellenic Positioning System).HEPOS consists of 98 stations all over the country.Its reference system is HTRS07 [Hellenic Terrestrial Reference System 2007]which is practically a local densification of ETRF 2005 [Katsampalos et al. 2010].HEPOS coordinates are fixed to 2007.5 epoch.HTRS07 is a static reference frame taking no consideration area's strong and inhomogeneous velocity field.This fact inextricably leads to severe network distortions, taking into account area's strong velocity field.It should also be reminded that ETRF 2005 is a withdrawn ETRF due to some problems related to coordinates' shifts, [Boucher and Altamimi 2011].

General concept
The main idea for the ORF realization relies on the Helmert type velocity similarity transformation between two different reference frames.The general mathematical expression for the velocity transformation is given by [Altamimi et al. 2002], pointwise: where v A i is the 3-D velocity vector in reference frame A, v B i is the 3-D velocity vector in reference frame B, Ṫ=[ṫ x ṫ y ṫ z ] T the vector of the shift rates with respect to each axis, Ḋ is the differential scale rate Ṙ as defined in Equation (3).Should be noted that the rotations and scale factor terms are negligible.In a more compact form, Equation (4) could alternatively be expressed, pointwise: where the design matrix.The Cartesian coordinates for each point could be referred to any arbitrary epoch t and the vector that contains all the transformation rate parameters.In the case of the transformation between a specified TRF to the ORF, Equation ( 4) can be modified accordingly to the following relation, pointwise: where v i ORF is the 3-D velocity vector in ORF, v i TRF is the 3-D velocity vector in the TRF.In order to realize the ORF it is necessary to introduce a specified optimization criterion which should be satisfied.In our case, this particular criterion is the minimization of frame's total kinetic energy, according to the Least Square principal [Ampatzidis 2011, Ampatzidis et al. 2015]: where P is a properly chosen weight matrix, the vector that contains all points velocities.A straightforward choice for P is the inverse covariance matrix of velocity error P=C e -1 .The transformation rate parameters are estimated through the following relation: and the velocities of the ORF will be: where Q is the projection matrix (Q 2 =Q) The covariance matrix of the LRF velocities is: The optimization criterion described in Equation ( 8), leads to ORF realization which its total kinetic energy is minimal.It could be interpreted as a physical equivalent of Tisserand's criterion implementation, in terms of minimizes frame's kinetic energy.Nevertheless, the Tisserand criterion is imposed to ITRF realization in order to treat the problem of orientation rate deficiency, having simultaneously a physical meaning (minimizes the total angular momentum and kinetic energy of the underlying reference frame).This is not the case for ORF realization, which is actually a full rank problem (the coordinates and the velocities are already referred to a global or regional TRF). (5) (8)

The geometrical meaning of the ORF
Another critical issue for the practical implementation of this procedure is the choice of the optimal rate transformation parameters θ .
. In fact, a scale-rate parameter should not be used; otherwise an artificial scaling distortion is introduced into the ORF.Moreover, in small geographical areas, the use of shift-rate parameters should be avoided due to their high correlation with the three orientation rate parameters.
In case of estimating only the three rotation rates we conduce to the same physical meaning like the estimation of a single set of Euler Pole Parameters (EPPs); see [Drewes 1982], for the underlying area.EPPs are related to the rotational rate parameters according to the following mathematical conversion formulas (spherical approach,): where ω EP , φ EP , λ EP , are the three EPPs (angular velocity, latitude and longitude of the Euler Pole reference point respectively).Hence, the ORF is directly associated with an Euler pole estimation which describe the motion of the underlying area, with respect to the initial reference frame.Finally, it must be noticed that the estimated rotation rates are high correlated, especially when they referred to a relatively small area.This practically means that a slight change in the number of stations and/or the distribution could modify significantly the estimated parameters.

The spatial connection between existing TRFs and the ORF
Until now our analysis scheme was concentrated to the dynamic part of local reference frame (LRF) and its derivation from a modern TRF.Assuming that in an arbitrary epoch t 0 the two reference frames (regional/ ITRF and ORF respectively) coincide, we have the following condition, pointwise: The 3-D position vector to any epoch t in the ORF is computed through the following expression: Take into account that (see Equation 8): Combining ( 14) and ( 15) the final mathematical relation becomes: Equation ( 16) allows the direct spatial connection between the initial TRF and the ORF (and vise-versa) at any epoch t.Thus the ORF approach does not cancel out the global or the regional TRF realization.On the contrary, the velocities of ORF can be transformed in any of the existing TRFs.

The ORF realization using exclusively 2-D geodetic velocities
In the most cases, the most dominant deformation is found to the horizontal plane.Thus, our strategy should also take into account for only the horizontal velocities.The 2D case of the ORF where the design matrix referring to the horizontal topocentric velocities, the horizontal curvilinear velocities (E for East and N for North components, respectively) in the ORF, the horizontal curvilinear velocities in the TRF and the vector contains the rotation rates with respect to each axis.For n points we have the following mathematical expression: (18) where where are the curvilinear velocities of an arbitrary point i and R the mean earth radius, respectively.
One might raise doubts when are used horizontal curvilinear velocities exclusively, claiming that omitting the vertical velocity part could not be provided the full 3-D position information.This possible disadvantage can be overcome introducing conventional values for the vertical velocity component e.g. from estimated global or regional TRF application in Hellenic area.

GPS Data and Analysis
The GPS data was analyzed on a previous work [Bitharis 2015], based only on permanent GPS network and carried out between 2008 and 2014.In this extended Hellenic Permanent GNSS Network (HPNet) we included 155 continuously GPS sites which are operated under the responsibility of different agencies (see Acknowledgements), as illustrated in Figure 2. HPNet included six stations which also contribute in the EUREF Permanent GNSS Network (EPN) [Bruyninx et al. 2012, Fotiou andPikridas 2012].The daily GPS data were processed using GAMIT/ GLOBK software [Herring et al. 2010], and following the methodology which is described in [Dong et al. 1998].The reference frame realization has been accomplished in ITRF2008, which was conducted by means of the 7-parameter similarity transformation using 24 IGS stations, at Eurasian and African tectonic plate.The IGS reference stations [Dow et al. 2005] selected carefully in order to obtain some crucial criteria as: i) The good geometrical distribution of GPS network.ii) Constant and continuous operation without discontinuities in their time-series for a period exceeding at least 7 years.iii) if it is possible, to collocated with more than one geodetic technique and included at the IGS fiducial (Core) network.
In the first step, daily raw GPS data were processed, with observation rate of 30 sec in order to estimate station coordinates, zenith atmospheric delay for each station (2hr interval), orbital and Earth orientation parameters.We applied some necessary corrections and models for the ocean [Lyard et al. 2006] and atmospheric loading tides [Tregoning and van Dam 2005] as recommended.For the orbital a-priori information IGS precise final orbits were used.
Before estimating site coordinates and velocities, we free our data from outliers or discontinuities, equipment changes (antenna replacement), jumps depend on strong EQ events.Secondly, we use the loosely constrained daily solutions of site positions, orbit and EOP, to estimate station coordinates and velocities estimated by a Kalman filtering sequential approach.In most sites the data span exceeds the 4 years, and the standard deviation of velocities determination in the horizontal plane is better than 0.8 mm/yr.In the present study, we choose to remove four sites from our analysis, which are located in one of the most geophysically active areas of Greece (Santorini Island), due to its strong local deformation activity.
In this case, we provide the ORF using only the horizontal velocities for the ITRF2008 and ETRF2000, as presented in the following chapters, respectively.

Results in ITRF
Firstly, the implementation of the ORF based on the horizontal velocities in the ITRF and the graphical representation of our test is illustrated at Figure 3.It should be noted that the magnitude of the velocities in the northern part of the Hellenic area is the largest one, however, they are close to the observed velocities of the Eurasian stable part.Unfortunately, the major discrimination occurs in the Aegean Sea as extensively described in many previous publications, [e.g.Floyd et al. 2010].
The reduction of kinetic energy is clear from the graphical representation of the velocities vectors, with more than 60 percent of total kinetic energy from the initial RF for the whole GPS network.In particular, the statistics of the 2D optimization scenario in the ITRF are given in Table 1.Bias indicated as the mean average of the sample.The bias of the ORF velocity magnitudes is less than 10 mm/yr.In a more detailed view, we show the statistics per topocentric component (East, North), as presented at Table 2.We found that the most significant reduction is at the East component, which is the dominant (with largest velocities than the North component) in ITRF at European region.
The estimated transformation parameters θ between the initial (ITRF2008) and the new ORF are given in the Table 3.The three rotation rates expressed in the appropriate Eulerian Pole parameters as described in detail at previous subchapter (see section 3.2.1).
We should also underline that the covariance matrices are not multiplied by the variance factor.The present procedure does not refer to the classical Gauss-Markov model; here we minimize the velocity norm and not the errors.Thus, the variance has not any physical meaning.

Results in ETRF2000
Secondly, we obtain the results in the ETRF2000, which is the most recent and recommended realization of the ETRS89.The main idea of the ETRS89 it has been to design a TRF with criterion the minimum velocities at whole the European region.Nevertheless, the statistics reveal that the ETRF2000 does not present any significant advantages with respect to the ITRF2008.Greece is an exception on this idea having velocities larger than everywhere in Europe.It should be mentioned that the median velocities are amplified in ETRF2000 from 15.5 (in ITRF) to 26.0 mm/yr (see Table 1, Table 4).In Figure 4 the 2D velocities in the ETRF2000 and the new ORF are illustrated.According to the statistics showed in Table 4, AN OPTIMAL GEODETIC REFERENCE SYSTEM FOR GREECE  the improvement is significant in ETRF, because the total kinetic energy was reduced more than 70 percent for the whole test area.The mean average of the velocity norm in ORF is less than 10 mm/yr.In Table 5 we present the statistic values for the velocities in ETRF2000 and the ORF per topocentric component, respectively.
Figure 5 shows a cumulative diagram of the total kinetic energy per initial RF and their optimal implementation and it is clear that ORF decreases the kinetic energy in both RFs.The results of the assessment between the two different approaches pointing out that ETRF and their implementations like HTRS07 could not be recommended as fundamental -official RF, in order to generate the cartographic and topographic   between the initial and the optimal ETRF2000 (units: mas/yr for the orientation rates, degrees for EPPs latitude and longitude).

Strategies Advantages Disadvantages
ITRF YY Represents the nature of the crust displacements, globally.
Not recommended in areas with inhomogeneous velocities.
The most models and products (orbits) are generated in ITRF.
A dynamic TRF which is not easily understand in order to generate the cartographic materials and intended for expert users.
Is a multi-technique combination of high accurate geodetic techniques (GNSS, SLR, VLBI, DORIS).

ETRF2000
Directly transformation with ITRF YY .
Strong and inhomogeneous velocities field in south Greece.
Homogenous velocity field in central Europe.
Recommended datum for regional mapping and surveying applications in Europe.

Local Euler Pole
In small areas gives good results and represents the local characteristics of the inhomogeneous velocity field.
The results are high correlated with microplates boundaries.
The accuracy of the estimated EPPs is low in small areas.

Minimization of Kinetic Energy Criterion
The velocity is minimal and led to a more stable TRF.
The velocity field in the most cases remain inhomogeneous.
Directly transformation between global/regional TRFs and ORF.
Do not reflect the geodynamic behavior of the area and not recommended for geophysical studies.
Could be applied both in 2D and 3D.
Recommended for national cartographic and geodetic purposes.materials, due to their velocities values.The estimated three-rotation rate parameters and the related EPPs are given in Table 6.The comparison between of two ORFs shows that, in both cases i) The total kinetic energy is severely reduced ii) The amount of kinetic energy is almost equal (see Table 1, Table 4) iii) The inhomogeneous velocity field still remains but with the minimum velocity norm.This could be confirmed due to their EPP which are different in the ITRF2008 and ETRF2000 (see Table 3; Table 6).Fundamentally, each different Terrestrial Reference Frame realization strategy has some advantages and disadvantages, in Table 7 summarized the main characteristics of each approach.

Conclusions
In this study the specific geodynamic activity of Greek area was mentioned.Also, the various methods for reference frame realization and minimization of kinetic energy are described in detail.In order to evaluate the effectiveness of the optimal reference frame, an application using the methodology of kinetic energy minimization is performed to a network of 151 permanent GNSS stations which covers all the Greek area.The well distributed network provides significantly results of the key role of ORF strategy.
As it was found the criterion generally reduces the velocity values by means of 60 percent.Applying the minimization algorithm with ETRF no significant advantages were came up for the study area.Instead of the application of optimal reference frame provides the advantage of using its coordinate system at various activities as the National Cadastral applications, where small velocities are of importance.
In addition, our proposed methodology gives flexibility, because it is possible to provide a direct transformation between any modern RF to an optimal RF.Thus, besides of some especially global and regional geodetic applications, which appropriately need ITRF and ETRS, our proposal would have some important advantages.
Finally, the ORF realization is independent of microplates boundaries definition and should be mentioned that do not have any obvious physical interpretation This strategy could be adopted from countries with active geodynamic as Turkey, Italy and Japan.

Figure 1 .
Figure 1.Principle tectonic setting of the Greek area and seismogenetic sources (Mw >5.5), from: Caputo et al. (2012), KFZ, Kefallinia Transform Fault, GC, Gulf of Corinth, GS, Gulf of Saronic NAT, North Anatolian Fault velocities are initially expressed in rad/yr.It is easy to convert them to a more comprehensive form (m/yr) by multiplying them by earth radius, as the following relation yields (pointwise):

Figure 3 .
Figure 3. Horizontal velocities at Greek area in the ITRF2008 (black color) and the optimal RF (red color)

Figure 4 .
Figure 4. Horizontal velocities at Greek area in the ETRF2000 (black color) and the optimal RF (pink color).

Table 2 .
Statistics of velocities per topocentric component in the Greek study area for ITRF approach (units: mm).

Table 3 .
Estimated rotation parameters and Euler Pole parameters between the initial and the optimal ITRF (units: mas/yr for the orientation rates, degrees for EPPs latitude and longitude).

Table 4 .
Statistics of horizontal velocities in the Greek study area for ETRF approach (units: mm, (mm/yr) 2 for kinetic energy).

Table 5 .
Statistics of velocities per topocentric component in the Greek study area for ETRF approach (units: mm).

Table 7 .
Characteristics of each Terrestrial Reference Frame realization strategy.