UNCERTAINTY IN THE DETERMINATION OF FAULT LOCKING DEPTH AND STRIKE SLIP RATES BY GNSS MEASUREMENTS

Original scientific paper A GNSS network of 36 points was established in 2006 in the eastern part of the North Anatolian Fault Zone, which is located between Tokat and Erzincan. GNSS measurements were made in 2006 ÷ 2008 in the Kelkit valley and these points’ velocities were calculated. This study compares GNSS velocities estimated with different datum definitions and the effects on fault locking depth and strike-slip rates, which are parameters of elastic rebound theory. These two parameters contain preliminary information on possible earthquakes in the region. An average slip rate of 21 mm/year and a locking depth of 12,72 km are estimated.

Results from only one method are inadequate for tracking crustal deformations.Therefore, it is crucial to use up-to-date data as well as to carry out concurrent analyses with different methods.Geodetic data are important for defining earthquake areas and source parameters [2,16].
GNSS studies on crustal deformations in Turkey were first done in the early 1990s.In [20] were reported movements of 16.÷.18mm/year in the western part of the North Anatolian Fault Zone (NAFZ).In [11] was estimated a velocity of 22.÷.24±1 mm/year along the North Anatolian Fault.In [17] were carried out studies on the movement of the earth's crust around Turkey, concluding that it moves about 25mm/year along the fault.The tectonic setting of Turkey can be seen in Fig. 1.
The study area is located between Erzincan and Erbaa (Tokat) in the eastern part of the NAFZ.The area has suffered intense and devastating earthquakes in 1939, 1942 and 1992 (Fig. 1).The 1939 Erzincan earthquake is the strongest earthquake recorded in Turkish history (M s = 7,9), which created a 360-km surface breakage from Erzincan to Ezinepazari [10].The 1939 Erzincan earthquake was rapidly followed by a series of disastrous earthquakes in Niksar-Erbaa in 1942 (M s = 7,1), Ladik in 1943 (M s = 7,3) and Bolu, Gerede and Çerkeş in 1944 (M s = 7,3).The centre of the destructive earthquake of 13 March 1992 (M s = 6,8) was in eastern Turkey, in the densely populated Erzincan region, 700 km east of the capital city of Ankara.This study examines the effects of GNSS velocities derived from different datum definitions on fault locking depth and strike-slip rates, the parameters of elastic rebound theory.These two parameters offer preliminary information on possible earthquakes in the region.Due to the simplicity of elastic rebound theory, it is frequently used by geologists in tectonic studies [5,12,19,23,29].
In these studies, one-dimensional elastic rebound theory parameters (fault locking depth and strike-slip rates) were obtained from GNSS and radar interferometry.Fairly similar results can be obtained from a wide variety of geodetic data for GNSS-derived velocity values.For instance, in this study, the Eurasian Plate was kept fixed, represented by data from several International GNSS Service (IGS) points on the Eurasian plate.In this way, more than one combination of data points can be used to estimate GNSS velocities.If velocities are similar using different combinations and their differences are stochastic, it suggests that each combination is equally valid, statistically speaking.

Data collection and processing 2.1 The GNSS measurements
A number of geodetic studies focusing on tectonics have been conducted in the central and western parts of the NAFZ [2,11,12,14,16,17,20,27].However, there are few studies from the eastern part, in the regions of Tokat-Erbaa and Erzincan-Cayirli [23,29].In order to determine the movements in this region with the required level of accuracy, localized errors have to be minimized.Therefore, GNSS measurements must be made on a stable point and the number of satellites must be considered.Periodic GNSS measurements were made at 36 points from 2006 ÷ 2008 in the eastern part of the NAFZ, in an area that includes the interface between the Eurasian and Anatolian plates, along the Kelkit Valley between Tokat and Erzincan.As shown in Figs. 2 and 4, the fixed points were TUTGA (Turkish National Fundamental GNSS Network), Sivas (SIVA), Çördük (CRDK), Gürgentepe (GURE) and Kelkit (KLKT).At these fixed points observations were carried out for 15 days ÷ 24 hours in each period.The duration of observation in each day was about 10 hours with an interval of 30 seconds and the elevation mask angle for the campaign observations was taken as 15 degrees (Fig. 2).GNSS observation campaigns were only conducted between July and August in order to minimize systematic seasonal errors.The same receivers and antennas were used at all sites (Trimble 4000 and 5700 with choke-ring and geodetic Zephyr antennas) [23] in order to reduce antenna phase-pattern problems and errors in the computation of the vertical components.These measurements made data recurrence possible and data quality was enhanced.With these criteria, 36 points were used to characterize the tectonics of the region.

GNSS data processing
GNSS measurements were evaluated with the GAMIT/GLOBK software developed by the Massachusetts Institute of Technology (MIT) [7,8].GAMIT can estimate 3D coordinates, satellite orbits, atmospheric zenith delays and earth rotation parameters using carrier-phase measurements and pseudo-range observations.GLOBK provides the velocity of each point based on a KALMAN filter algorithm [7,8].Data were processed according to the GAMIT manual.
All data were obtained between the years of 2006 and 2008, which were then processed using the same processing strategies described in the GAMIT manual.GAMIT uses file extension such as SP3 files, USNO_bull_b as outside rinex files, broadcast ephemeris.The precise orbit information by IGS (International GNSS Service was obtained in SP3 (Standard Product 3) format from SOPAC (Scripps Orbit and Permanent Array Center).Earth Rotation Parameters (ERP) came from USNO_bull_b (United States Naval Observatory_ bulletin).The 9-parameter Berne model was used for the effect of solar radiation pressure on the GNSS satellites [11,14].
GNSS observations were processed in three stages, following [4].In the first stage, GAMIT was used to determine the coordinates and atmospheric zenith delays of each point and Earth Orientation Parameters (EOPs), based on daily observations of doubly-differenced GNSS phases.Next, this local network was tied to a global network in order to improve calculations, using data from five IGS stations (BUCU, GRAZ, KIT3, MATE and SOFI), which track the earth's orbit and rotation to the nearest millimetre.
In the second stage, a Kalman Filter was performed with EOP values, daily orbit coordinates and covariance to estimate site coordinates and velocities from the combined solutions.The regional h-files (adjustments and full variance-covariance matrix for input to GLOBK) from each GAMIT session were integrated with global IGS analyses that SOPAC publishes daily.This ensured that the regional orbit module and the global orbit module were in agreement, by using sensitive orbit and rotation parameters from IGS.The present study used IGS's general networks generated and igs1, igs2 and eura.
In the third stage, the Eurasian reference frame was defined and checked daily, based on a reliable set of global IGS stations, the ITRF2000 and ITRF2005 no-netrotation frames, in order to estimate velocity.Finally, GLORG provided local network velocity.At this step, the local network was defined relative to a general reference system, so translation, orientation and scale parameters were calculated [7,8,11].Origin can be determined with any one fixed site but any error in this site will be translated to other stations.Points defined as datum are chosen from the SOPAC archives by doing a time series analysis in order to select the points that are the least affected by crustal movements.Orientation is defined by EOP values but the errors are again going to be propagated into all sites.Scale may or may not be a problem in case of usage.If it is going to be used, GNSS should have well defined scale [7,8].The Eurasian plate was defined by minimizing horizontal velocities of sites (Table1).The velocities in the region according to the analysis of the third field trip were computed with the methods described in [11] (Fig. 3).At this point, reference frames of different point combinations were used (Tab.1).
Points were selected that represented Eurasia in the ITRF 2000 reference framework and also the ITRF 2005 reference framework.A regional velocity area was created by keeping the GURE and KLKT points at the upper parts of the fault fixed.Points representing the Eurasian plate were selected for stabilization.In this way, velocity values relative to the Eurasian Plate were estimated.
The first work on datum determination kept the Eurasia Plate fixed [11].Some points were excluded: ZWEN (Russia), BOR1 (Poland), ONSA (Switzerland) and NYAL (Norway) GNSS points.These points could have adversely affected the GNSS velocity area because field measurements were of poor quality or not made at the same time as those used in this study (Tab.1).On the other hand, points from other sources were suitable: BUCU, GRAZ, KIT3, MATE and SOFI.Velocity values were determined from different datum definitions (Tabs. 2 and 3).The six solutions in Tab. 2 refer to: [11] for solution 1, five IGS points for solution 2, KLKT and GURE for solution 3, which all use the ITRF00 reference frame.Solutions 4 ÷ 6 use the same points but in the ITRF05 reference frame (Tab.3).All solutions used Eurasia-fixed reference frames but have different results.Solutions 1 and 4 suggest velocities were approximately 2 mm/year faster than solutions 3 and 5, which are approximately 4 mm/year faster than those calculated with local points in solutions 3 and 6.
At-test was used to evaluate whether these differences were statistically significant (see Appendix A).There is no statistically significant difference between the different solutions' velocities and their standard deviations.This shows that all of the solutions and data sources are statistically acceptable for estimating velocity.

Locking Depth and Slip Rate Estimated from Geodesy
Studies on earthquake cycles and surface displacements have used elastic rebound theory [15].Relative movements of blocks separated by fault boundaries create linear (time-based) accumulations of tension.The areas close to the boundary cannot shift as much as the areas that are farther from the boundary as inter-seismic tension accumulates, from blocks that are locked above a certain depth under the fault boundary.Tension is released when earthquakes occur, depending on the geometry and geology of the crust, and the slip deficit near the fault boundary is cleared.Tension then begins to accumulate again and the cycle continues [24].
Considering Tabs. 2 and 3, the velocity vectors of the points are differenced slightly.Also in Fig. 3, similarity between the rotations directions of the vectors when reference system of ITRF00 is used to define 3 different data is indicated.Fig. 4a shows the initial situation before tension begins to accumulate, when line A-A' is straight and perpendicular to the fault, the vertical line.Line A-A' deforms as tension accumulates (Fig. 4b).Line B-B' represents the situation shortly before an earthquake.
After the earthquake, line A-A' is broken at the fault (Fig. 4c).A and A' move away from each other but remain straight.Line B-B', on the other hand, bends away from the fault line, based on distance from the fault line.In its simplest form, elastic rebound theory [18] can be expressed as V p is a constant rate equal to the long-term surface (i.e.geologic) slip rate.Locking depth is shown in Fig. 4. For this one-dimensional elastic rebound theory model, the inverse model finds the best V p and D (slip rate and locking depth) for a set of observed velocities (v is GNSS points' velocity and x is perpendicular distance from fault).Elastic rebound theory is discussed in more detail in other studies [12,19].Like other similar research, this study uses Monte-Carlo simulations [5,29].
Monte-Carlo simulations are used in earth sciences when samples or observations are few and no theoretical model can be proposed for the population of data [1].The best-fit values presented in Tab. 4 were calculated.Fault locking depth and strike-slip rates were obtained through 500 iterations of 3 ÷ 20 km locking depth and 10 ÷ 30 mm/year slip rate.These ranges were based on the analysis of earthquakes and previous research in the region [11,17,23,29].
The slip rate reported by [11] is 22 ÷ 24±1 mm/year for the NAFZ, whereas [17] estimated a slightly higher rate, 25 mm/year, based on block modelling.The average rate for this study varies between 19±3 mm/year and 22±3 mm/year, in agreement with [11,17,23].The slip rate of the eastern NAFZ increases westwards.Over a horizontal distance of about 400 km, the rate increases from 16,3±2,3 to 24,0±2,9 mm/year and locking depth increases from 8,1±3,3 to 12,8±3,9 km.Slip rates based on geological observations are 20,5±5,5 ÷ 27±7 mm/year, consistent with [6,9].Fault locking depth for the study area was obtained 15 km depth from the earthquake catalogue published by Kandilli Observatory and the Earthquake Research Institute.As shown in Tab. 4, six different solutions for locking depths and slip rates produced fairly similar results.Eq. (A2) demonstrates that the differences between the solutions are statistically not significant.When Tab. 4 is examined, it has been determined by applying t-distribution that statistically all obtained values can be used and that the differences between them can be random.According to a t-test, these differences are not statistically significant.The values of Solution 1 and Solution 4 have been obtained using the velocities obtained from the points that take the Eurasian Plate to be constant as used by McClusky et.al. in (2000).On the other hand, the fault locking and slip rates in Solution 2 and Solution 5 have been determined using the velocity area created by the points used to determine the rates relative to the Eurasian Plate.Solution 3 and 6 have been derived by using the rate values obtained according to the Eurasian plate defined by using the regional points.In general, the selection of the points representing Eurasia have been decreased from 12 to 2 and the same values with random errors have been obtained for the velocity areas obtained for each solution and the parameters of the elastic rebound.The comparison of results obtained from six different solutions in fault locking depths and strike slip rates in eastern part of NAFZ can be seen in Tab. 4.  Using elastic rebound theory, locking depth and slip rates, point velocities were recalculated and plotted relative to distance from the fault.In Fig. 5, the horizontal axis indicates the distance from the fault.The vertical axis indicates the velocities of the individual points.Velocities of the points south of the fault are 22 mm/year and velocities decrease as they approach the fault.Moving northwards, these velocities continue to decrease until reaching zero.Tension accumulation in the study area has a locking depth of approximately 12 km.There is an inter-seismic period of slow accumulation of elastic strain that coincides with frictional locking of a fault between earthquakes throughout the region, that is, pre-earthquake tension development.
According to the model, the high rate of GNSS velocity at Ataköy (ATKY) is because it is located between NAFZ's Kelkit segment and the Almus Fault Zone (Fig. 2) [3].As seen in Fig. 5, KMAH, ILIC and IMRN do not follow the expectations of the elastic rebound model.It is possible that the change of velocity at these points occurs because of their location between the NAFZ, Malatya-Ovacik Fault Zone (Fig. 2) [28] and the Central Anatolian Fault Zone (Fig. 2) [3] and because these fault zones are closing in on each other.

Conclusions
We have reported GNSS measurements from a geodetic network established for tectonic research.We analysed GNSS velocity data in the eastern NAFZ region and estimated locking depth and slip rates.Several solutions were tested and it was determined that the upper 12 km of the eastern NAF is locked and the nearby crust behaves elastically.
GNSS velocities across the fault segments between Niksar and Erzincan were calculated with GNSS field measurements between 2006 and 2008.In 2012, [23] examined tectonic movements in this segment by constructing three fault-perpendicular profiles.Slip rates of the NAFZ along these east-west profiles were 16,3±2,3, 18,5±2,2 and 24,0±2,9 mm/year, respectively, and increase westwards.Fault locking depth also increases westwards, from 8,1±3,3 in Erzincan to 12,8±3,9 km in Niksar.The slip rate in the central part of the NAFZ, immediately to the west of the study area, was estimated by [27] to be in the range of 18,7±1,6 ÷ 21,5±2,1 mm/year.To the east of the study area, slip rates were deduced by [13] to be in the range of 16 ÷ 24 mm/year.The slip rates obtained by geological observations are 20,5±5,5 ÷ 27±7 mm/year [6,9] and agree with the present study.It estimated an average slip rate of 21 mm/year and a locking depth of 12,72 km.
Stabilization of the base-vector components was carried out in order to determine the velocity area for different datum definitions.These were three groups of points: those used by [11], the IGS stations used for GNSS evaluation (BUCU, GRAZ, KIT3, MATE and SOFI) and two local stations.These three groups of points produced differences of 1 ÷ 4 mm in estimated velocity.According to a t-test, these differences are not statistically significant.
Velocity increases from north to south on the Eurasian and Anatolian plates and from east to west on the Anatolian plate.A westwards increase in velocity was determined for the Anatolian plate, relative to the Eurasian plate.In [11] this was calculated to be 24 mm/year along the NAFZ, while [17] reported the same value to be 25 mm/year on average.There is a counterclockwise angular rotation on the Anatolian Plate, from east to west, with an average velocity of 15 ÷ 22 mm/year (Fig. 3).This seems to be in agreement with the overall counter-clockwise rotation of Anatolia [11,17].
Eurasian-plate velocities decrease north of the fault until they reach zero (Fig. 5).The approximate velocity at the GURE point, located at the northernmost point of the fault, was calculated to be 4 mm/year.This indicates the width of the deformation zone on the Eurasian plate and shows that the elastic deformation in the region spreads over a wide area of approximately 140 km in length covering both plates.These findings indicate the presence of tension accumulating in the study area, which has a locking depth of approximately 12 km.It might also be determined that there is an inter-seismic period of preearthquake tension development, which is eminent throughout the region.
In the article there are differences between solutions of few millimetres.Few millimetres of point movements in seismic sense might be significant according to different geosciences.But, in the text there is also a statistical analysis that defines these movements as insignificant movements considered as random errors by means of Geodetic data usage.Still these insignificant movements and their effects should be inspected by other geosciences.
As can be seen in Fig. 5 the best-fit curve agrees well with the observations, indicating that deformation is due to fault locking and that the block behaves elastically.The horizontal region with major velocity variation is located at twice the locked depth (12 km).Whereas near-fault observations provide information about fault locking depth, observations farther from the fault define general block behaviour.Hence, both types of information are necessary in a fault-zone study.Data acquisition strategies should be developed along profiles perpendicular to the fault axis and cover an area about twice the locking depth in order to define the strain accumulation as a function of depth.

Appendix A: t-distribution
The t-distribution was developed by William Gosset.It is often applied to small samples with a humped distribution that do not closely follow the normal distribution [1,25,26].Two different stabilization points were used to obtain y i and y j forecast values along with s i and s j standard deviations for the magnitude of a point's velocity relative to the Eurasian plate.The statistical equality of y i and y j values is tested with the null hypothesis given below.The velocities estimated with

Figure 1 (
Figure 1 (A) To define Eurasian Plate are selected IGS points and local points; (B) Tectonic setting of Turkey and surrounding regions.The blue rectangle shows the location of the study area in the fault zone.Thin red lines are active faults in Turkey.The focal-mechanism solutions of the earthquakes are from [21].

Figure 2
Figure 2 The GNSS network established along the Kelkit valley.Black points were monitored for 24-hour periods.Red points were monitored for 10 hours per day over three consecutive days.Thin black lines show active faults in the Kelkit Valley.

Figure 3
Figure 3 Eurasia-plate velocity vectors according to points from the Kelkit GNSS Network in the ITRF00 reference frame.The red dashed line shows the approximate location of the study area.

Figure 4
Figure 4 Conceptual drawing of displacement around the strike-slip fault: (a) at initial stage, (b) during deformation and (c) after the fault slip.

Figure 5
Figure 5 Elastic rebound model of the Kelkit GNSS network.Vertical bars represent error ranges for the point velocities.

Table 1
GNSS Sites used to define the Eurasia-fixed reference frame

Table 2
Lateral velocities and standard deviations (mm/year) of the GNSS deviations, using the ITRF00 frame and a fixed Eurasian Plate

Table 3
Lateral velocities and standard deviations (mm/year) of the GNSS deviations, using the ITRF05 frame and a fixed Eurasian Plate

Table 4
Comparison of six solutions for fault locking depths and strike-slip rates in eastern NAFZ