Local Deformation Precursors of Large Earthquakes Derived from GNSS Observation Data

Research on deformation precursors of earthquakes was of immediate interest from the middle to the end of the previous century. The repeated conventional geodetic measurements, such as precise levelling and linear-angular networks, were used for the study. Many examples of studies referenced to strong seismic events using conventional geodetic techniques are presented in [T. Rikitake, 1976]. One of the first case studies of geodetic earthquake precursors was done by Yu.A. Meshcheryakov [1968]. Rare repetitions, insufficient densities and locations of control geodetic networks made difficult predicting future places and times of earthquakes occurrences. Intensive development of Global Navigation Satellite Systems (GNSS) during the recent decades makes research more effective. The results of GNSS observations in areas of three large earthquakes (Napa M6.1, USA, 2014; El Mayor Cucapah M7.2, USA, 2010; and Parkfield M6.0, USA, 2004) are treated and presented in the paper. The characteristics of land surface deformation before, during, and after earthquakes have been obtained. The results prove the presence of anomalous deformations near their epicentres. The temporal character of dilatation and shear strain changes show existence of spatial heterogeneity of deformation of the Earth’s surface from months to years before the main shock close to it and at some distance from it. The revealed heterogeneities can be considered as deformation precursors of strong earthquakes. According to historical data and proper research values of critical deformations which are offered to be used for seismic danger scale creation based on continuous GNSS observations are received in a reference to the mentioned large earthquakes. It is shown that the approach has restrictions owing to uncertainty of the moment in the beginning of deformation accumulation and the place of expectation of another seismic event. Verification and clarification of the derived conclusions are proposed.


Introduction
Since the middle of the last century, the study of precursors of strong earthquakes has become an important task of modern geophysics. Repeated traditional geodetic measurements, such as highprecision levelling [1] and observations of linear-angular networks, were used to study this issue. Many examples of the detection of precursors of strong seismic events using classical geodetic methods are presented, for example, in [2]. According to the recording frequency, deformation precursors (horizontal and vertical deformations of the earth's surface, slopes and deformations in boreholes) in the Rikitake report are presented as the most numerous. Considering creep anomalies at the fault and groundwater levels, they accounted for 38% of the total number of geophysical anomalies recorded before the earthquakes. One of the first researchers of deformation precursors of earthquakes recorded by geodetic methods was Meshcheryakov [3]. The model proposed by Meshcheryakov [3], for example, is confirmed by the results of analysis of multiple levelling data on the Sheki-Kurdamir line, located on the southern slopes of the Greater Caucasus [1]. The normativity of the forerunners of this type is limited by the absence of the possibility of predicting the more or less accurate position of the future epicentre. This can be explained by the low density of precise levelling networks.
Today, sufficiently dense GPS observation networks are operating in seismically active regions. Analysis of the results of these observations is of great interest for studying the behaviour of the earth's surface near the epicentres of earthquakes during the process of their preparation.
Evidence of the existence of a deformation precursor of a strong earthquake based on GPS observations in a regional network is presented in [4,5]. Anomalous behaviour of the earth's surface before the strongest earthquake of Japan, remote from the regional observational network [6], was also discovered.
A review of the evidence for the existence of deformation precursors is presented in [7]. In particular, the article gives several cases of their registration using GPS. The author makes an encouraging conclusion regarding the possibility of earthquake predictions using deformation precursors.
This study is devoted to the analysis of continuous GPS observations before and after major earthquakes to identify local deformation precursors.

Territories under study, observation networks, data and coordinate time series collection
The object of the study was the western coast of North America, which is the boundary of the Pacific and North American global tectonic plates. The territory is covered by a rather dense network of GPS observations, within which three strong earthquakes occurred during the last decades: Parkfield (September 28, 2004, Mw = 6.0); El Mayor Cucapah (April 4, 2010, Mw = 7.2); Napa (August 24, 2014 Mw = 6.0). CORS geodetic networks, such as the Plate Boundary Observatory (PBO) GPS network, the Southern California Integrated Geodetic Network (SCIGN), the Bay Area Regional Deformation (BARD) network, and others [8][9][10], operate in the region. US state services, municipalities, and universities participate in a single complex of network development and functioning. The local fragments of the GPS network covering the areas of postseismic ruptures of the earth's surface, as well as the Napa epicentral zone of the latest seismic event (August 24, 2014 Mw = 6.0), were selected. The Delaunay triangulation was used to form triangles representing the finite elements of deformation analysis.
For the cases of the Parkfield and El Mayor Cucapah earthquakes, daily observation files were selected from the open data archive [http://sopac.ucsd.edu/], as well as navigation files with a frequency of once every five and ten days, respectively, to the first and second observational networks. The interval of registration of GPS radio signal was 30 s.
The preliminary processing of observational data was performed using the Topcon MAGNET Tools software. The processing yielded three-dimensional vectors of the baselines and their covariance matrices Q xyz . The components of the baseline vectors were increments of rectangular geocentric coordinates Δx, Δy, and Δz.
The values of horizontal displacements in the cases of the first two GPS networks were calculated from the results of free adjustment of time differences of geocentric coordinate increments ΔΔx i = Δx i - Δx 1 , ΔΔy i = Δy i -Δy 1 , and ΔΔz i = Δz i -Δz 1 , where the increments of the initial cycle 1 are subtracted from the increments of coordinates of the current measurement cycle i. The system of correction equations of the differences in the coordinate increments was solved using the least square method, taking into account the covariance matrices Q xyz and using the principal pseudo-inverse matrix of the coefficient matrix of the normal equations. This approach is called a free adjustment [11]. As a result, the displacement vectors of all the GPS network points without exception and the mean square errors of their determination were obtained. These results served as the basis for calculating the space-time sequences of horizontal deformations of the earth's surface. The GNSS point displacements were obtained in the local topocentric reference frame.
For the Napa earthquake, the time series of horizontal displacements obtained by the Nevada Geodetic Laboratory [http://geodesy.unr.edu/index.php] were used directly. In this case, the processing strategy was a "point precise positioning" (PPP) mode using the precise GNSS products of the Jet Propulsion Laboratory. The daily GNSS point data were processed using the GIPSY OASIS II package. The displacement values are obtained in the IGS08 global reference frame.

Earth crust strain calculation
The horizontal deformations were calculated for the Delaunay triangulation triangles produced.
To calculate deformations of finite elements, the strain tensor (1) the elements of which are , and , correspondently, wherein , was used.
To study the spatiotemporal distribution of horizontal deformations, the following invariant characteristics were obtained. 1) Principal strains ε 1 and ε 2 .

Graphical representation and visualization of the deformation process
The obtained values of the maximum shear and dilatation were interpolated to a regular grid with Hermitian splines using the standard Matlab software package procedure. In this way, the pictures (frames) of the horizontal distribution of dilatation deformations and the maximum shear were obtained for dates specified. Data coverage intervals and strain determination spacing are presented in table 1. To study the nature of the temporal changes, animation was produced from each sequence of frames, demonstrating the behaviour of deformation in time, both before and after the earthquakes. Corresponding sequences of frames with a rarer reproduction frequency than indicated in the table, as well as for time intervals from the beginning of observations to the moment of the main shock, are presented in figures 1-6.

Interpretation of the results and discussion
The obtained images of the space-time distribution of deformations show the following features.
Over some years before the studied earthquakes within the observational networks, zones of deformation inhomogeneities were formed and expressed in changes in dilatation and maximum shear. Values of deformations reached 10 -5 or more. Obviously, over time, these zones continue their development in the direction of deformation rising.
The outlines of zones of inhomogeneous deformation are similar to the corresponding patterns of coseismic deformations as evidenced by the last frames of the presented sequences. Thus, coseismic deformations seem to slowly appear in the forms of interseismic deformation, gradually reaching a critical level.
In his studies, Rikitake [2] showed that the strain level from 0.5 * 10 -5 to 1.7 * 10 -4 is characteristic of seismic ruptures for earthquakes of 5.9 < M < 8.4. This result is based on the analysis of repeated measurements in geodetic networks. As studies of rock samples show, the threshold of their destruction is reached in the range of values 5 * 10 -4 -10 -4 [12]. These levels of deformation can be considered as signals predicting the occurrence of seismic events near the localization of these deformations.
According to the hypothesis of elastic rebound during the earthquake, elastic deformations accumulated during the preparation of the event are realized in its focus [13,14]. On this basis, the level of coseismic deformation can be taken as a critical level. Achieving this level will indicate a high seismic danger. To assess this critical level, data were collected on the coseismic displacements near the epicenters of strong earthquakes in the regions under study (see table 2). From the values of the maximum coseismic horizontal displacements, the values of the dilatation and the maximum shear were calculated. In this case, an equilateral triangle with a side length equal to the average length of the triangle of each geodetic network was taken as a finite element. The maximum coseismic displacement was assigned to one of the vertices of this triangle and is oriented from the centre to periphery. The other two vertices were assumed to be stationary due to a rapid decrease in displacements as they moved away from the seismic rupture zone. The values of the seismic displacements and the corresponding deformations are presented in the table 2.
To obtain a homogeneous data set, the dilatation values were reduced to a uniform average area of the triangle P m with the side equal to the averaged side of the triangle of all geodetic networks used. The scale factor was m = P i / P m . The scaled dilatation was calculated as Δ m = ΔP m . A plot was produced for the dependence of magnitude of earthquakes and scaled dilatations. It is shown in figure  7.

Conclusions
The study of spatiotemporal changes in horizontal deformations due to strong earthquakes on the northwestern coast of North America (Parkfield M6.0, USA, 2004; El Mayor Cucapah M7.2, Mexico, 2010; and Napa M6.1, USA, 2014) revealed the peculiarities of the occurrence of anomalous deformation zones years before the future seismic event. The pattern of coseismic deformations repeats the preceding pattern of deformation accumulation as the deformation process develops, but with a higher intensity.
The revealed features may be local deformation precursors of strong crustal earthquakes occurring on seismically active strike-slip faults.
Continuous GPS observations are an important element of the seismic event prediction system. With the accumulation of empirical data on critical deformations and temporal features of the development of the deformation process, the success of developing a seismic hazard scale with the use of continuous GPS observations in networks covering well-known seismogenic zones is very likely.