Unique Pre-Earthquake Deformation Patterns in the Spatial Domains from GPS in Taiwan

: Most earthquakes are considered to be caused by stress accumulating in and subsequently releasing from the crust. To extract non-linear and non-stationary earthquake-induced signals associated with stress accumulation, the Hilbert–Huang transform was utilized to ﬁlter long-term movements, short-term noise, and frequency-dependent (annual and semi-annual) variations from surface displacements measured by the global positioning system (GPS) in Taiwan. Earthquake-related surface displacements were expressed as horizontal directions (i.e., GPS azimuths) using the north–south and east–west components of residual GPS data to bypass inﬂuences resulted from the inhomogeneous nature of the crust. Analytical results showed that the relationships between earthquake occurrence and the aligned GPS azimuth passed the statistical test of the Molchan’s error diagram. Aligned GPS azimuths were in agreement with direction of earthquake-related P axes for 81% (26 / 32) studied events. Areas with the highest paralleling orientations of GPS azimuths appeared around epicenters several days to weeks before earthquake occurrence. Durations from aligned GPS azimuths to earthquake occurrence are roughly proportional to earthquake magnitude. Similar variations of the GPS azimuths were observed in GPS data containing or excluding co-seismic dislocation (i.e., one day before) in the temporal and spatial domain. These suggest that the aligned GPS azimuth could be a promising anomalous phenomenon for studying crustal deformation before earthquakes.


Introduction
Reid [1] proposed the elastic rebound theory based on the observations of surface deformation during the 1906 San Francisco earthquake. The theory held that stress accumulation would be stored in the crust via deformation due to the elasticity of rocks. When stress accumulation exceeds thresholds of fault rupture, subsequent release of partial stored energy results in earthquakes, stratum dislocation and the thrust fault. We make no distinction between earthquakes dominated by different fault types due to insufficient numbers for classification. Once agreements between the aligned horizontal azimuths and directions of the P axes were obtained, surface displacements were considered in relation to particular earthquakes. Locations and timings of the aligned horizontal azimuths were correlated with epicenters and earthquake magnitudes, respectively, to further advance our seismic knowledge during earthquake preparation periods. The solid (open) stars denote earthquakes for which aligned GPS azimuths agreed (disagreed) with the P axes. Also labeled is the HUAL (Hualien) station, which was illustrated as an example for the Hilbert-Huang transform (HHT) analysis.  [2006][2007][2008][2009] and were used in this study. NaN denotes that a phase of stress loading was not detected. ID with * reveals earthquakes utilized in Figure 8. The LT is the leading time that reveals the time of aligned horizontal azimuths before the earthquakes. The strikes, dips, and slips are dislocation parameters of the faults.

Theory and Methodology
Earthquake-related strains are generally considered to be non-stationary changes based on long-term variations of GPS data. Clarifying the earthquake-related strains in the GPS data is very difficult due to the varying velocities of plate movements and unknown timing of the earthquake-related stress that begins to obviously dominate faults before dislocation. Scientists have tried to retrieve the earthquake-related strains by investigating unusual variations of strain change rates. These investigations are generally based on an assumption of fixed velocities of plate movements and homogeneous materials in the Earth's crust (as illustrated in the upper plot of Figure 2).

Figure 2.
Diagrams of displacements resulting from stress accumulation in homogeneous and inhomogeneous materials. When the same stress accumulates in homogeneous materials in (a), the quantity of displacements is associated with the Young's modulus. In reality, because the surface crust is comprised of materials with different values of Young's modulus, characteristics of stress (as shown by blue arrows) are misunderstood when only the quantity of displacements is taken into account. On the other hand, displacements resulting from stress accumulation exhibit a similar direction regardless of the Young's modulus of materials (as shown by blue arrows) in inhomogeneous materials in (b). The direction of stress-inducing displacements is a good substitute for the quantity of the displacements and further approaches the intrinsic qualities of stress-inducing phenomena involving inhomogeneous structures.
The simplification of homogeneous materials in the crust can affect scientific understanding of earthquake-related stress accumulation. Different quantities of displacements can be generated by the same stress accumulating in rock layers due to discrepancies in the Young's modulus of materials underground (as illustrated in the bottom plot of Figure 2). In reality, surface displacements occur in response to inhomogeneous structures underground, such as lateral variations of elastic properties in rock layers [28,29]. Orientations of deformation and/or displacements most likely adapt the major axis of loading stress despite the quantities of displacements related to the Young's moduli of materials underground. When orientations of the major axes associated with earthquakes and large-scale tectonic stress loading in the crust owing to plate movements are simultaneously taken into consideration, a significant discrepancy between their orientations can often be observed. Thus, if orientations of displacements are utilized to study earthquake-related stress accumulation, both physical issues (i.e., underlying inhomogeneous structures and the source of strain contributions) can be simultaneously at play, obscuring the desired clarification. In short, orientations should be normally random as long as no significant earthquake-related stress accumulates in the crust, owing to the removal of the effects of the long-term plate movements. Once orientations were aligned in a similar direction, the major axes of loading stress associated with earthquakes were taken into comparison to clarify factors of the alignments.
Relative positions of 100 stations were calculated with respect to the Kinmen site (KMNM, located at 118.39 • E, 24.46 • N) by assuming a static relative distance for the daily epoch using Bernese software to process the double-differenced GPS phase observations. We corrected the tropospheric delay using the Saastamoinen model [30] with an elevation-dependent weighting of cos 2 (z), where z is the zenith distance. To estimate the wet delay of GPS signals in the zenith direction, an hourly wet tropospheric parameter was applied [31]. The ionospheric delay was removed using a combination of GPS data at L1 and L2 frequencies, and satellite clock errors are removed using the International GNSS Service (IGS) precise ephemerides [32]. The coordinates of the KMNM station were determined precisely using the international terrestrial reference frame for geocentric coordinates.
Because displacements associated with earthquakes are generally expressed as non-linear and non-stationary signals, the HHT [23][24][25] was used to remove long-term movements, short-term noise, and frequency-dependent variations from GPS data of the 100 stations. The HHT process consisted of two steps: the empirical mode decomposition (EMD) and the Hilbert transform [23]. Initial data were processed with EMD to generate several intrinsic mode functions (IMFs) depending on energy characteristics; these IMFs then underwent the Hilbert transform. The EMD process involved subtracting the average envelopes of the data computed based on the cubic spline constructed from the local maxima and minima of analyzed data (for more details, see Reference [23]). Here, IMF was determined as long as the mean of the difference between the analyzed data and the envelopes was less than 10 −6 m. If the mean of the difference was greater than 10 −6 m, the data were replaced by the difference, and the process was repeated accordingly. IMFs obtained from the EMD process were subsequently removed from the initial data one by one, and the residuals were used in repeating the EMD process. Once there were no longer sufficient local maxima and minima to generate the envelopes, the EMD process is terminated. Following the EMD process, the instantaneous frequency and amplitude of all points in the derived IMFs were calculated via the Hilbert transform. For each IMF, the quantitative time-domain displacements, instantaneous frequencies and frequency-domain amplitudes can be obtained. Characteristics of derived IMFs were examined and a frequency band is further chosen to filter out unwanted influences (i.e., noise, long-term movements, co-seismic dislocations, etc.). Orientations of GPS-azimuths were computed from the north-south and east-west components of the residual displacements (i.e., filtered displacements) at each station.
To examine whether GPS-azimuths were related to earthquakes or not, θ w values, which were determined by inverses of daily average difference from GPS-azimuths distributed within a moving spatial window in a particular day (t), were computed via where k is a total number of stations located within a moving spatial window, and A (i, t) and A (j, t) are GPS azimuths (A) at the i and j stations on day t, respectively. The spatially varying θ w values were utilized to construct θ w maps in Taiwan in a particular day. Notably, the small θ w values of about 0.011 (~= 1/90 • ) implied that GPS azimuths at each station were generally independent of each other (i.e., no preference in orientation among stations) in the moving spatial window. The small values indicated that either no significant stress accumulated in the crust, or the accumulation stress approached the threshold of the fault rupture (see also References [22,33]). In contrast, the relatively large θ w values suggested that the disordered GPS azimuths were aligned in a similar direction primarily due to adaptation of loading stress.
On the other hand, we determined θ(r, t) by using the inverses of average differences from GPS azimuths of l stations located within epicentral distances, r, from 20 to 200 km, and the temporal period, t, from 70 to 1 days prior to the earthquake, to examine relationships between GPS azimuths and epicentral distances, if any.
The duration of the 70 day period was determined to roughly avoid interplay between two contiguous events, because M5 earthquakes often occur with an average interval of about 50 days (estimated from 68 events in 9 years, 2004 to 2012) in high-seismicity Taiwan. To examine the relationships between the relatively large θ(r, t) and fault parameters of these studied earthquakes, the azimuths from θ(r, t) and directions of the P axes were taken into consideration simultaneously.

Examples and Results
We  Figure 3m shows the residual that represents the long-term movements. The slope of the residual (about −25 mm/year) indicated the relatively moving speed of east-west motion between Kinmen and HUAL during the study period and was consistent with previous studies [34].    We took 5.787 × 10 −7 Hz (period = 20 days) as an upper bound of the filter because high-frequency IMFs intrinsically represented the noise in this work. Figure 4a shows variations of the displacements derived from IMFs with instantaneous frequencies higher than 5.787 × 10 −7 Hz (period < 20 days). Large fluctuations appeared in every summer period due to disturbances in the atmosphere and ionosphere [35]. We subtracted the high-frequency noise from displacements to generate noise-free data (Figure 4b). The long-term effects (including long-term movements, interseismic deformation, and semi-annual and annual variations in Figure 4c), estimated from the sum of the residual and points in IMFs with instantaneous frequencies lower than 7.716 × 10 −8 Hz (period > 150 days), were further subtracted from the noise-free data again (shown as Figure 4d). The entire process therefore involved a critical step of using a band-pass filter with frequencies between 7.716 × 10 −8 Hz and 5.787 × 10 −7 Hz to retrieve the targeted GPS displacements (i.e., the residual displacements).
Here, a land earthquake, EQ28 (5 November 2009 M L = 6.2; 23.789 • N, 120.719 • E, depth of 6.9 km, Figure 5 and listed in Table 1), was taken as an example to examine the pre-earthquake changes of θ w in the spatial domain, due to its large magnitude and good coverage of GPS measurements. In this study, the θ w maps were constructed using a moving spatial window of 80 × 80 km 2 , which was able to cover a sufficient number of GPS stations with a moving step distance of 5 km. The daily θ w maps for the period from 30 (EQD-30) to 1 (EQD-1) days prior to the earthquake day (EQD) of EQ28 are shown in Figure 5. It was clearly that the entirety of Taiwan was characterized by: (1) small θ w (~0.011; pink color in Figure 5) between EQD-30 and EQD-22, (2) large θ w (>0.02; yellow color in Figure 5) between EQD-21 and EQD-11, and (3) small θ w between EQD-10 and EQD-1, a short transient stage with low-signals right before the occurrence of EQ28.
Before EQD-22, the small θ w values suggest that surface crust did not experience perceptible stress associated with EQ28. After that period, southeastward movements (~145 • ) were obviously observed in the subsurface crust and large θ w (~= 0.031; also see Figure 5) values appeared between EQD-21 and EQD-12. The direction of the southeastward movements was approximately orthogonal to the strike of the fault (228 • in Table 1). The persistent stress accumulation gradually approached the threshold of rock rupture, leading affected areas to drop to small θ w , which was initially observed on the southeastern side of Taiwan. The small θ w values gradually expanded to all directions between EQD-10 and EQD-1, eventually encompassing the entirety Taiwan. We also found that an area (120.8 • N, 23.5 • E) in central Taiwan displayed relatively large θ w between EQD-8 and EQD-1. Orientations of GPS azimuths in the area with relatively large θ w of about 110 • were consistent with the azimuth of the P axis (100 • ) retrieved using the method proposed in Robinson and McGinty [27]. The existence of this spot was caused by slight displacements and/or cracks [36][37][38][39] that would have developed before the main rupture of EQ28. In short, stress accumulation in the crust led to parallel alignments of those GPS azimuths to a large and detectable extent. When loading stress continuously approached to the threshold of crustal rupture, parallel features yielded by aligned GPS azimuths tentatively resumed disordered and disarranged patterns in the stressed area over a relatively short period prior to the earthquake. Areas characterized by disordered GPS azimuths (low θ w values) gradually extended due to the approach of the threshold of earthquake-related loading stress. Similar characteristics of the parallel GPS azimuths could also be observed before the Chi-Chi earthquake (M = 7.6 on 20 September 1999; Chen et al., 2013a) and the Tohoku-Oki earthquake (M = 9.0 on 11 March 2011 [33]). For the aforementioned thrust-type events, the alignments were orthogonal to the strikes of earthquake faults and/or were similar to the direction of the P axis, which prompted to further examination of the relationships between GPS azimuths from the residual surface displacement data and fault parameters of earthquakes.
On the other hand, the P axis of EQ32 was 111 • , which agreed with the 110 • result obtained from the aligned GPS azimuths of θ(50, 58) (i.e., in an area with the epicentral distance ≤ 50 km and appeared about 58 days before EQ32; also see Figure 6). In contrast, the P axis of EQ1 was either 141 • or 321 • . The azimuths of about 310 • , which was retrieved from functions via the relatively large θ(50, 34), was consistent with the P axis of EQ1. Azimuths determined from relatively large θ(r, t) agrees with the P axis of earthquakes could be observed in the 26 events of the 32 studied earthquakes. Disagreements included five M < 5.5 offshore earthquakes (i.e., EQ4, EQ8, EQ21, EQ25, and EQ30) which would have been caused by GPS stations located on land far away from the epicenters. In contrast, the other M = 5.2 land event (i.e., EQ 13) showed different characteristic changes in the GPS azimuths before the earthquake for unknown reasons (see also Reference [40]). This suggests that the analytical methods used in this study are most suitable for large and land earthquakes. The aligned azimuths from the residual displacements could be linked with the P axis, which was neither noise nor artifact, but was dominated by earthquake-related stress.   The relationship between the P axis and θ(r, t) for the 32 studied earthquakes. Colors ranging from pink to yellow denote quantities of θ(r, t). If θ(r, t) > 0.0211, the black arrows are shown on the diagrams and their directions denote orientations of the GPS azimuths derived from the residual displacements. Notably, the value of 0.0211 is the so-call anomalous band (i.e., mean+2σ) computed from the inverse of the difference of the GPS azimuths between every pair of stations during the entire study period. The fault plane solution of each earthquake is shown in the upper right diagram. The P axis agreement (disagreement) with the aligned GPS azimuths is shown by using blue (black) color and marked above on associated days (below the fault plane solutions).
To further examine the aforementioned relationships, the 12 large earthquakes with magnitudes between 5.4 and 6.9 were qualified in this regard (shown as ID with * in Table 1). Figure 7 shows that changes in θ(r, t) along with varying r from 20 km to 200 km were associated with the 12 events on the t day, with the azimuths of residual surface displacements being most representative of the P axis of earthquakes. The θ(r, t) values were inversely proportional to r only within a relatively short range, but became relatively constant and smooth farther away. This observation suggested that surface crust near epicenters undergoes substantial displacements due to the high loading stress associated with earthquakes. Epicentral distances constrained by θ(r, t) values could be considered as the affected range of accumulating stress before earthquakes. In general, the affected epicentral distances before earthquakes with magnitudes between 6.9 and 5.4 were about 100 and 50 km, respectively (i.e., relatively high θ(r, t) values in Figure 7a,l). The affected distances roughly agreed with the radius of the earthquake preparation zones estimated in Dobrovolsky et al. [41]. Meanwhile, the likely locations of epicenters were determined by using the top 0.5% of daily θ w values from the θ w maps associated with EQ28 during stress accumulation periods to compare with the real (as published by CWB) one of EQ28 via the inversely proportional relationships between θ(r, t) and r. The likely locations of epicenters associated with EQ28 determined from high θ w values were continuously found on days between EQD-21 and EQD-12 except on EQD-13 ( Figure 5). On days of EQD-20, EQD-19, EQD-17 and EQD-12, the likely locations of epicenters were shown at two possible locations, but on other days, only a single spot (120.35 • N, 23.75 • E) appeared. The distance difference between the CWB-published and the likely locations of epicenters of EQ28 on the EQD-12 was about 38 km.

Discussion
The aforementioned observations showed definitive relationships between the residual displacements and seismic parameters (i.e., strikes of faults, directions of the P axes, and areas of earthquake preparation zones) in physical nature. The values of θ(r, t) and leading times, which were the duration from the initial appearance of the aligned GPS azimuths that agreed with the P axis of the subsequent earthquake occurrence, were further related to earthquake magnitudes using those 26 direction-matched events. Clearly, no apparent relationship could be found between the values of θ(r, t) and earthquake magnitudes (Figure 8a) because some studied earthquakes occurred in the ocean, while the GPS measurements were primarily taken on land. On the other hand, a proportional relationship between the leading time and earthquake magnitude could be roughly obtained (Figure 8b). Note that the Chi-Chi (M = 7.6) and Tohoku-Oki earthquakes (M = 9.0) on 20 September 1999 and 11 March 2011, with the leading times of 54 and 65 days, respectively [33,42], were included to extend the studied scale.
The constructed θ w maps for Taiwan showed that loading stress has led to developments of aligned GPS azimuths prior to many earthquakes, and the leading time was related to earthquake magnitude. Likely locations of epicenters could be roughly determined in these maps based on the aligned residual displacements that occurred several days before the earthquakes. However, the aforementioned results, which were obtained from the time-series data covering earthquake occurrence, were dominated by co-seismic changes. To clarify the effects of co-seismic changes on our analytical results, EQ 28 was taken as an example again. The daily θ w maps for the period from 30 to 1 days prior to EQ28 were constructed using the GPS data from 1 January 2006 to the day before EQ28 to avoid the co-seismic effects (Figure 9), and were also compared with Figure 5. It was obvious that θ w patterns of GPS azimuths between Figures 5 and 9 were very similar. Slight discrepancies existed due to co-seismic (or step) changes, but did not affect the earthquake-related patterns of GPS azimuths. This comparison suggested that HHT analyses can effectively mitigate the effect of co-seismic (or step) changes from other minor to medium earthquakes with insignificant surface ruptures. However, this observation will not be valid if co-seismic (or step) changes are very large, as was the case in the Chi-Chi and the Tohoku-Oki earthquakes. Relationships between the maximum of θ(r, t) and earthquake magnitude are shown in (a). In contrast, relationships between the leading time and earthquake magnitude are shown in (b). The black circles denote the 26 earthquakes detected during the study period. Note that the maximum of θ(r, t) and the leading time associated with the Tohoku-Oki and Chi-Chi earthquakes were computed using different r values (i.e., r 50 km) because the Tohoku-Oki earthquake occurred far away from the land and due to the relatively-low density GPS array surrounding the Chi-Chi earthquake in 1999.
Although aligned residual displacements in agreement with the P axis direction of earthquakes were obtained in this study, we also computed the θ w maps in Taiwan every day in the entire studied period. We found no major earthquake occurrence following significant aligned GPS azimuths during the periods of May-June 2007 and December 2007-January 2008. To examine whether the aligned azimuths could be useful precursory information in Taiwan, we performed a statistical test by using Molchan's error diagram. Molchan's error diagram is a conventional method by which to evaluate the prediction performance of given observations for earthquake forecasting [43,44]. The diagram plots the rate of earthquakes detected against the rate of alarmed cells. In practice, for given series of precursory signals and related earthquake events, each prediction strategy was characterized by the lead time of alarms and the length of the alarm window. The lead time was the time length between a detected anomaly and its alarm, and the alarm window was the duration that an alarm lasted. A modified area skill score measuring the area between actual prediction curve and random prediction line in the Molchan's error diagram [43] was used to assess the efficiency of different prediction strategies, characterized by the lead time of alarms and the length of alarm windows [44]. A score of 0.5 indicated perfect skill and a score of −0.5 indicated complete non-skill. The expected score for a random prediction was 0. Scores above 0 indicated that the prediction was better than random and the observation contained precursory information [44].  Figure 10 shows the area skill score of prediction in color gradation from blue (low score) to yellow (high score), based on the observations using different lead times and alarm windows. There were two main high score regions (higher than 0.03): one with a lead time around few days and a long alarm window, and the other with lead time of 14-16 days and an alarm window of less than one week. The maximum S, which suggests the optimal prediction strategy, has a value of 0.068, and was located at lead time 1 day and alarm window 3 days. Figure 11 demonstrates the Molchan's error diagram of the optimal prediction strategy based on the observations (lead time = 1 day and alarm window =3 days). The prediction curve shown by the bold red line was clearly above the random prediction line given by the black diagonal line. Most predictive curves here clearly exceeded the 90% confidence interval, and quite a few even exceeded the 95% confidence interval. These results highly suggested that the observations contained precursory information for sizeable earthquakes and could be potentially used in short-term earthquake forecasting. Figure 10. Area skill scores for predictions with different lead times and alarm windows based on observations. We performed a forecasting experiment and predicted earthquakes with depths of less than 100 km and magnitude ≥ 5.0 during the studied period in and around Taiwan.

Conclusions
Time-series GPS data contain an abundance of factors. Once the effects of unrelated signals and artificial noise have been mitigated in long-term continuous GPS data using HHT, earthquakes occur following short-term crustal displacements with the aligned orientations in agreement with the P axes of earthquakes, as seen in in Taiwan. The aligned orientations before earthquakes could be consistently obtained from GPS data before earthquake occurrence and passed the statistical test of Molchan's error diagram. The alignments of the short-term crustal displacements could provide potential insight for discriminating and/or determining the timing of earthquake-related stress which significantly accumulates in the crust before the dislocation of faults. The highest paralleling alignments revealed the epicenter locations of forthcoming earthquakes. The durations of the disordered orientations after the alignments were proportional to the magnitudes of the forthcoming earthquakes. This suggested that short-term crustal displacements could be utilized as a promising indicator of seismo-crustal displacement anomalies.
Meanwhile, the alignments of seismo-crustal displacements anomaly can be utilized as references to examine other promising seismo-anomalies [40,45]. Negative (positive) strain changes can relate to uplift (depression) of seismo-ground water levels. Earthquake-related stress accumulates in the crust and can increase or decrease electrical conductivity underground, driving changes of seismo-electromagnetic anomaly. Once GPS data recorded by a dense array over a wide area can be utilized, regions of earthquake preparation could be determined by stations with aligned orientations. The large-scale stress accumulation would expose the causal mechanism of the GPS TEC (total electron content) and thermal anomaly that generally distributes in an area of more than ten thousand square kilometers (100 × 100 km) before earthquakes. In short, these examinations shed light on the integration of multiple physical parameters to construct a potential model that can encompass most pre-earthquake anomalous phenomena.