GPS Seismology for a moderate magnitude earthquake : Lessons learned from the analysis of the 31 October 2013 ML 6 . 4 Ruisui ( Taiwan ) earthquake

The 31 October 2013 ML 6.4 Ruisui earthquake was well recorded by twelve 50-Hz, four 20-Hz and thirteen 1-Hz GPS receivers, and twenty-five strong motion stations located within the epicentral distance of 90 km in eastern Taiwan. Kinematic positioning solutions estimated by four GNSS software (TRACK, RTKLIB, GIPSY, VADASE) are used to derive the seismic waveforms and the co-seismic displacements for this event; strong motion accelerometers are used to verify the capability of high rate GPS to detect seismic waves generated by this earthquake. Results show that the coordinate repeatability of the GPS displacements time series are ~6 mm and ~20 mm standard deviation in the horizontal and vertical components respectively, after applying spatial filtering. The largest co-seismic displacement derived from high-rate GPS is nearly 15 centimeter at 5 km northeast of the epicenter. S waves and surface waves are successfully detected by motions of high-rate GPS and double-integrated accelerometers within the 15 km epicentral distance. For the first time twelve 50-Hz and four 20 Hz GPS observations for seismological study were used and analyzed in Taiwan; a clear benefit was evidenced with regard to the seismic waves features detection, with respect to the 1-Hz GPS data, so that ultra-high rate (> 1-Hz) observations can compensate the sparse coverage of seismic data, provided proper monuments for the GPS permanent stations are realized. Spectra analysis between co-located GPS and strong motion data further suggests that the optimal sampling rate for high-rate GPS Seismology study is 5 Hz. The 2013 Ruisui Taiwan earthquake recorded by the high-rate GPS permanent stations network in Taiwan demonstrates the benefits of GPS Seismology for a moderate size earthquake at a local scale.


Introduction
High-rate GPS has become an important sensor for the seismic wave detections and earthquake warning systems for moderate and large earthquakes.The true ground displacements retrieved from the processing of high-rate GPS observations provide precise seismic waveforms to identify the wave property, the arrival time of the body waves and the surface wave [Kouba 2003, Larson et al. 2003, Bock et al. 2004, Langbein and Bock 2004, Kouba 2005, Ohta et al. 2006, Larson et al. 2007, Bilich et al. 2008, Larson and Miyazaki 2008, Miyazaki and Larson 2008, Davis and Smalley 2009, Larson 2009, Shi et al. 2010, Avallone et al. 2011, Ohta et al. 2012, Branzanti et al. 2013, Li et al. 2013, Benedetti et al. 2014, Li et al. 2014a,b].For large disastrous earthquakes, such as the 2004 Sumatra and the 2011 Tohoku-Oki earthquakes, highrate GPS gives the possibility to provide information about fault rupture processes, focal mechanism determinations, and tsunami warning and monitoring in near real-time [Ohta et al. 2006, Blewitt et al. 2009, Crowell et al. 2009, Colosimo et al. 2011a,b, Ohta et al. 2012], also considering the remarkable advantage of the no-clipping and no-tilting features with respect to the strong motion and broadband seismometers [Bilich et al. 2008, Zheng et al. 2012].Significant surface wave propagation and attenuation are observed by regional dense Continuous GPS Permanent Stations (CGPSs) network [Davis and Smalley 2009, Grapenthin and Freymueller 2011, Hung and Rau 2013].
However, the traditional sampling interval of the high-rate GPS (1-Hz) cannot always describe the seismic waveforms completely due to waveforms of higher frequency in near-field earthquakes [Avallone et al. 2011].Efficient sampling rate of the position time series for the waveforms monitoring are controlled by the Nyquist-Shannon sampling theorem, which states that the recording sampling rate should be twice larger (in practice, it should be at least four times larger) than the main frequency for the interesting signals.Some moderate earthquakes recorded in near-field with 5-Hz and 10-Hz sampling rate demonstrate the benefit of the very-high-rate GPS for the waveforms recover [Avallone et al. 2011, Zheng et al. 2012, Lou et al. 2014].On the other hand, differently from the reconstruction of the entire waveforms, it was showed that co-seismic displacements can be accurately detected by the combination of the 1-Hz GPS solutions and co-located accelerometer [Bock et al. 2011, Melgar et al. 2013].
As regard the present study, since 2010 a dense CGPSs network of nearly four hundred sites has been established at 10-20 km spacing in Taiwan from various institutions.About 250 stations established mainly by Central Weather Bureau (CWB), Central Geological Survey (CGS) and National Land Surveying and Mapping Center (NLSC) are available for high-rate GPS recording for both geodetic and geophysical applications.The receivers installed more recently are characterized by a larger storage and a higher sampling rate capacity [Yeh et al. 2012, Hung andRau 2013].
These features were exploited to collect GPS observations at high-rate and very high sampling rate for this study, where we used 1-Hz, 20-Hz and 50-Hz CGPS observations for the investigations on the moderate Ruisui earthquake (M L 6.4) occurred on 31th October, 2013 in eastern Taiwan, which was already object of other seismological researches [e.g., Lee at al. 2014, Wen et al. 2016].In details, thanks to the availability of both CGPS network, broadband seismometers and strong motion accelerometers, the goals of the work was to develop a deep comparison about the solutions supplied by these different instruments (also considering different GPS software) as regards seismic waveforms and the co-seismic displacements, and to verify the capability of high-rate GPS to detect body waves and surface waves.
In section 2 the used data and software are introduced; in section 3 the main results are presented and discussed; in section 4 some conclusions are outlined.

Data and method
CGPS dual-frequency observations at 1-Hz, 20-Hz and 50-Hz sampling rate were collected to estimate the epoch-by-epoch position time series during the 2013 Ruisui earthquake in eastern Taiwan (Figure 1) for the seismic waveforms and co-seismic displacements monitoring.Detailed information for the four 20-Hz GPS stations from National Cheng Kung University and twelve 50-Hz GPS stations from Central Weather Bureau, including the types of the GPS receivers and antennas are listed in Table 1.Observations of broadband seismometers (BB) were collected from BATS (Broadband Array in Taiwan for Seismology, http:// bats.earth.sinica.edu.tw/) and strong motion data were collected from the Central Weather Bureau, Taiwan.Note that some near-field seismic waveforms obtained from broadband seismometer stations HGSD and YULB are clipped (Figure 2).GPS epoch-by-epoch displacements were estimated from phase measurements in kinematic positioning [Hofmann-Wellenhof et al. 2006], following three different approaches and using four GNSS software: -Differential Positioning (DP), using TRACK software [Herring 2009].The reference station used, PNHU, as shown in Figure 1 is located 200 km west of the epicenter, outside the deforming region with stable coordinate time series.-Precise Point Positioning (PPP), using GIPSY-OA-SIS II [Zumberge et al. 1997] and RTKLIB software [Takasu 2011]; precise orbits and clock error corrections were generated by Jet Propulsion Laboratory for GIPSY-OASIS II, and from final International GNSS Service (IGS) solutions for RTKLIB.-Variometric approach, using VADASE software [Colosimo et al. 2011a, Branzanti et al. 2013, Benedetti et al. 2014]; here solutions were obtained in post-processing mode, but in order to simulate its real-time capabilities standard broadcast orbits and clocks were used.
To reduce the position errors due to the multipath effect in the cycles of the sidereal day, modified sidereal filtering (MSF) was constructed and applied to the original positions [Choi et al. 2004, Hung andRau 2013].Subsequently, a spatial filter (SP) were applied to each position time series to remove the common-mode error resulted from unmodeled effects during data processing [Wdowinski et al. 1997, Bilich et al. 2008, Yin and Wdowinski 2013], where the stations MOV1 and WUST (Figure 1) located about 150 km away from the epicenter are used to construct the spatial filter.These two stations were checked for their data quality and if there is any disturbance induced by the earthquake.The displacements after applying spatial filtering are analyzed for the investigations of their precisions and waveform propagations.The main features and parameter settings for the processing with all GNSS software are listed in Table 2.As regard the strong motion data (SMT), the 200-Hz sampling rate observations were transformed into displacements by double integration using software TSPP established by USGS [Boore 2010].Baseline corrections were applied to SMT using a zeroth-order-correction by removing the mean of 8-seconds segment before the earthquake origin time.Instrument responses were removed by deconvolution of the low-cut (0.02 Hz) and high-cut (100 Hz) filters in the acceleration space.For the broadband observations, the displacements were generated by single integration after removing the instrument response given the poles and zeros of the instrument [Hung and Rau 2013].
Therefore, a comparison among the different GNSS software, the strong motion and broadband integrated solutions was possible at the level of epoch-by-epoch displacements.

Results and discussion
Several aspects were investigated considering the results obtained using different software, with different observations sampling rate and at different CGPSs.Hereafter the main results are presented and discussed.

Seismic waveforms vs GPS observations with different sampling rate
To investigate the reconstruction of the seismic waveforms from GPS at different sampling rate, we firstly considered the effectiveness of the waveforms sampling interval at 0.2-, 1-, 5-, 10-, and 20-Hz.Figure 3a shows that the clear seismic shaking occurs after 5 seconds with respect to the earthquake origin time for all the ground motion time series.The sampling rate at 0.2-Hz and 1-Hz truncate the seismic signals and cannot retrieve the complete waveforms; also 1-Hz sampling rate apparently is not high enough to retrieve the complete patterns of the ground vibrations.As a matter of fact, spectral analysis show that the major energy in this earthquake is located at 0-1.5Hz interval (Figure 3b).This is also why waveforms retrieved from 10-Hz and 20-Hz sampling rates (Figure 3a) indicate similar ground motions, therefore in the following we refer to the 10-Hz GPS solutions only instead of the original sampling rate.Notice that the frequency spectra of sampling rate of 5 Hz, 10 Hz and 20 Hz become flat at the frequency lower than about 2.5 Hz.

Software comparison
The seismic displacements derived from 10-Hz observations at CGPS SIFU for software TRACK, GIPSY-OASIS II, RTKLIB and VADASE are displayed in Figure 4, as an example.The coordinate repeatability of these 10-Hz solutions for the 20 seconds interval before the earthquake was estimated for the four software (TRACK, GIPSY-OASIS II, RTKLIB and VADASE (Table 3)).It is clear that similar waveforms were obtained with the different software at CGPS SIFU in all the three components during 5-30 seconds after the earthquake origin time (Figure 4).VADASE software, whose solutions were obtained using broadcast ephemeris and clocks only, was able to reconstruct waveform features similar to those of other software using precise orbits and clocks information.Albeit with its 2 times larger uncertainties (Table 3), solutions with VADASE software confirmed the approach reliability.Spectral analysis also indicates similar results in the frequency domain among the solutions obtained with the four different software (Figure 5).Similar results confirming the good agreement among the software were obtained for other CGPSs.It has to be recalled that here all the solutions were obtained in post-processing mode, but all the approaches are able to supply real-time solutions, albeit with different technological  6).Larger motion drifts occur for VADASE solutions, which can be reduced after applying SP.The noise reduction effect are similar for the four software, except for the higher drifts for TRACK solutions, which may be resulted from longer baseline estimations (around 200 km).

GPS vs strong motion
Figure 6 shows that the ground motions have excellent agreement between GPS and co-located double-integrated SMT data (Table 4).Correlation coefficients in horizontal components between both measurements are 0.6-0.75 after applying 0.2-Hz high-pass filter (Figure 7).The need to high-pass the SMT data is to remove the effects of low frequency bias, introducing trends after integration.Spectral analysis shows that both observations have similar patterns in the frequency span of 0.3-2.0Hz (Figure 8).The ground motion sensitivity limit is about -20dB for GPS observations for a frequency higher than about 2.5 Hz, and above this frequency the noises of GPS are always larger than those of the corresponding SMT data (Figure 8).This suggests that a minimum sampling rate of 5 Hz is required for optimum high-rate GPS observations.Energy in SMT solutions at the frequency larger than 2.5-Hz decreases rapidly and remains at the flat noise level for the GPS results.The noise level in the Up component is always larger than those in the horizontal component, especially for the station far away from the epicenter, such as station CHNT (Figures

Waves propagation and attenuation
Displacements for all the high-rate GPS of 10-Hz and 1-Hz sampling interval in the transverse and radial components show clear directivity effects (Figure 10).Apparent group velocities of about 4.0 km/s are shown in both components for waves propagating northward.The wave amplitudes in the northern direction are apparently larger than those toward South.Seismic wave are hardly detected at the epicenter distance larger than 90 kilometer.Displacements of 1-Hz GPS (green line) improve the spatial coverage of the network for the estimation of the apparent group velocity (Figure 10).

Co-seismic displacements
Co-seismic displacements for each CGPS in all the three components were determined by the difference between the coordinate averages of 100-150 seconds after the earthquake origin time and  50 seconds before the earthquake.Ground motions were estimated from the 1-Hz high-rate GPS measurements just using software TRACK software, considering the proven good agreement among the four software; after applying the spatial filtering, around 2-12 centimeters co-seismic displacements at the stations within 5-40 km to the epicenter were detected (Figure 11a and Figure 11b).The precision of the co-seismic displacements are 3 mm and 5 mm in the horizontal and vertical components, respectively.
Figure 12a and 12b show that the largest co-seismic displacements are located about 10 km NE of the epicenter.For the four GPS stations with noticeable co-seismic displacement (horizontal component > 16 mm), we compare their co-seismic displacements estimated from 1-Hz and daily solutions (the coordinate average of 3 days after the earthquake subtract the 3 days average before the earthquake), respectively (Figure 12a and 12b).We found that the two stations (DSIN and GUFU) closer to the epicenter have larger co-seismic displacements estimated from the daily solutions than those derived from 1-Hz solutions, while the other two stations (SIFU and FENP) showing similar results.This indicates that the nearer stations may have early post-seismic deformation.On the other hand, peak ground displacement (PGD) of the stations show the largest value of about 200 mm for both horizontal and vertical components (Figure 12a and 12b), and PGD distributions indicate clear attenuations of the amplitude with respect to the epicentral distances.

Ground vs. roof monuments
We selected three pairs of the CGPS to investigate how the waveforms computation is affected by the different types (ground, roof ) of monuments (Table 5, Figure 13).Ground-type monuments are established by four steel piles on the ground, and the roof-type antenna are mounted on the roof using the single-steel piles (Figure 13).Results show that the waveforms (1-Hz sampling rate) have no significant differences for pair height difference lower than 10 meters.Pair DSIN-GUFU, around 40 Figure 11a.Time series highlighting the co-seismic displacements of 1-Hz CGPSs using TRACK software, after applying the spatial filtering in all the three components.meter height difference, has around 3 cm motion difference at 7-10 second after the earthquake origin time (Figure 14).Since the inter-station distance of DSIN-GUFU is around 3 kilometers, so that the waveform differences can be due also to this different location.

Trapped waves in the Longitudinal Valley
Seismic waveforms obtained from CGPSs and SMTs display regional inconsistence in eastern Taiwan.Seismic waveforms triggered by the Ruisui earthquake are significantly different among stations in the Central Range Figure 11b.Time series highlighting the co-seismic displacements of 1-Hz CGPSs using TRACK software, after applying the spatial filtering in all the three components.
(CER), Longitudinal Valley (LV), and Coastal range (COR) (Figures 15,16,17).The later phases, ~15 seconds after the S waves, have clearly larger amplitudes in LV than in COR and CER (where basically cannot be detected), especially in the radial component.This is also confirmed by the spectral analysis of the later phases, highlighting a significantly stronger energy in the frequency range 0.5-1 Hz within LV than in the other two regions, with ~1-2.5 cm peak-to-peak amplitudes (Figure 18).Therefore, the joint analysis of the CGPSs and SMTs derived waveforms suggest that later phases are trapped within the Longitudinal Valley, being this likely due to the fact that this is a suture zone composed of Holocene thick sediment deposits.

Conclusions
In this study, the first twelve 50-Hz, four 20-Hz and thirteen 1-Hz CGPS epoch-by-epoch displacements in Taiwan are derived for seismological applications with respect to the 31 October 2013, M L 6.4 Ruisui earthquake, Taiwan.The importance of GNSS seismology, better if in real-time, can be well understood considering that, even for a moderate magnitude earthquake like the Ruisui one, waveforms generated from the broadband seismometer are clipped for some of the near-field stations within 80 km epicentral distances.These truncated waveforms may cause severe mistakes for the estimation of the focal mechanism and the earthquake magnitude.GPS derived waveforms can provide quite useful information to compensate   the weakness of the seismic measurements close to the epicenter; in particular, observations at sampling rate as high as 50-and 20-Hz demonstrated their capability in retrieving reliable seismic waveforms for the near-fault observations.

GPS
We used several types of positioning approaches and related software to evaluate their reliability for seismological applications.At first, results showed that the displacements estimated by PPP have a similar precision as those derived from the DP.Furthermore, the basic feature of PPP is the displacements/waveforms estimation without reference station(s); that can avoid spurious effects due to the motion of the reference station(s) required in DD approach, in case these stations are affected by the earthquake too.On the other hand, real-time PPP solution is achievable using dual-frequency (geodetic class) receivers only.Secondly, seismic displacements estimated using VADASE supplied similar waveforms in comparison with results derived using other software implementing DP and PPP approaches.The three approach solutions are consistent at the level of 1-2 centimeters in the frequency range 0.08-3 Hz.Moreover, VADASE approach supplies a more efficient data processing with respect to the DD and PPP, since no convergence time for initial phase ambiguities is required, so that the observations just collected during the earthquake are sufficient to provide the solution (in our case, 2 minutes instead of 2 hours observations were used).In addition, since VADASE is suitable to process single frequency observations, it actually represents a good alternative solution for the low-cost and real-time GPS seismology.On the other hand, VADASE solutions can suffer for unmodeled effects common to close CGPS, which have to be removed by a spatial filter; anyway, this requirement is not a severe drawback and can be acknowledged if (as standard) all the solutions of a CGPSs network are managed within a centralized data center.
The high-rate GPS waveforms were compared with those derived through double integration from approximately co-located SMT, within an inter-distance within 2 kilometers.A good agreement, with correlation coefficients 0.7 in the frequency span of 0.3-2 Hz, confirmed that 50-and 20-Hz GPS can be an alternative with respect to accelerometers in the moderate or large earthquake.However, we found that, as expected, the accuracy of the GPS waveforms relies on the type of the monument (ground or roof ); anyway, if the height above the ground is small (let say, within 10 meters), the effect is still limited, within 1-2 centimeters.This fact must be accounted for since in the urban areas of Taiwan, most of the CGPS stations are established on the roof of the buildings due to the poor GPS satellite geometry on the ground.Results of spectra analysis between co-located GPS and strong motion data indicate that the ground motion sensitivity limit for GPS is about -20dB for a frequency above 2.5 Hz, and thus we suggest that the optimal sampling rate for high-rate GPS is 5 Hz for seismology research, which is also shown by Galetzka et al. [2015].
The performed analyses proved that waveforms derived from 50-, 20-, and 1-Hz CGPSs and strong motions measurements are substantially coherent at very few centimeters level.Their combination was proven to be effective: when jointly analyzed for seismological interpretation, the waveforms generated from both high-rate GPS and SMT could give information about

Figure 1 .
Figure 1.Map of CGPS and seismic stations around the epicenter of the 2013 Ruisui earthquake, eastern Taiwan.Focal mechanism of the Ruisui earthquake is shown in lower-hemisphere projection.Yellow, green and red dots denote CGPS at 50-Hz, 20-Hz and 1-Hz frequency sampling, respectively, cyan dots are strong motion stations; blue diamonds are broadband seismic stations.Stations MOV1 and WUST located at the western Taiwan are selected as the base stations for the composition of the common mode filters.Station PNHU was the fixed station for the GPS processing using software TRACK.

Figure 2 .
Figure 2. Seismic waveforms generated from broadband seismometer station HGSD, YULB, and NACB in all the three components in Ruisui earthquake.Waveforms are clipped for stations HGSD (all the three components) and YULB (East and North components).

Figure 3 .Figure 4 .
Figure 3. (a) Displacements time series at CGPS SIFU in sampling interval from 0.2 to 20 Hz in the East component using software VADASE (t=0 sec is the earthquake origin time).(b) Power spectral density (PSD) of the displacements time series at the station SIFU with the 20-Hz (brown), 10-Hz (green), 5-Hz (purple) and 1-Hz (blue) sampling interval in the East component; the major energy content of the seismic waveforms is located at 0-1.5 Hz.

Figure 5 .
Figure 5. Power spectral density (PSD) of the displacements time series at CGPS SIFU generated from software TRACK (brown), GIPSY (green), RTKLIB in the PPP solution (blue), and VADASE (grey) in all the three components.Figure 6. GPS displacement time series retrieved after applying spatial filters using four GPS software.Stations MOV1 and WUST are used to compose spatial filters.Brown, green, blue, and purple lines show the motions of MOV1 and WUST calculated using VADASE, RTKLIB, TRACK, and GIPSY, respectively.Time series of station SIFU after applying spatial filtering are shown at the bottom.

Figure 6 .
Figure 5. Power spectral density (PSD) of the displacements time series at CGPS SIFU generated from software TRACK (brown), GIPSY (green), RTKLIB in the PPP solution (blue), and VADASE (grey) in all the three components.Figure 6. GPS displacement time series retrieved after applying spatial filters using four GPS software.Stations MOV1 and WUST are used to compose spatial filters.Brown, green, blue, and purple lines show the motions of MOV1 and WUST calculated using VADASE, RTKLIB, TRACK, and GIPSY, respectively.Time series of station SIFU after applying spatial filtering are shown at the bottom.

7
and 8).Spectrum coherence analysis between 60 seconds time series of GPS and SMT station pairs indicate large coherencies of near 1 in the frequency span of 0.3-2 Hz, especially for the station pairs SLIN-HWA020 and SHUL-HWA001 in the horizontal components (Figure9).

Figure 8 .
Figure 8. Power spectral density (PSD) of the displacements time series from CGPS and co-located strong motion data after double integration and application of 0.2-Hz high-pass filter.Blue is GPS, and red is SMT.

Figure 7 .
Figure 7. Displacements time series in all the three components from resampled 10-Hz CGPS and co-located strong motion data after double integration and application of 0.2-Hz high-pass filter.CC is the correlation coefficient.

Figure 9 .
Figure 9. Spectrum coherence in three components for the 60 seconds times series of station pairs (Figure 7) of 10-Hz GPS and co-located strong motion data after double integration and 0.2-Hz high-pass filter.

Figure 10 .
Figure 10.Resampled 10-Hz displacements (red line) and 1-Hz displacements (green line) time series in the transverse and radial components with respect to the epicentral distance (an epicentral distance is positive northward with respect to the epicenter).

Figure 12 .
Figure 12. (a).Horizontal co-seismic displacements derived from the 1-Hz and daily CGPS measurements.The 1-Hz displacements were determined by the difference between the coordinate averages of 100-150 seconds after the earthquake origin time and 50 seconds before the earthquake.The largest horizontal co-seismic displacement was detected at SIFU, approximately 6 cm toward SE.Peak Ground Displacements (PGD) of the stations are shown in horizontal and vertical components, respectively.(b).Vertical co-seismic displacements derived from the 1-Hz and daily CGPS measurements.See the figure caption of Figure 12a.

Figure 13 .
Figure 13.Examples of CGPS antennas in ground type and roof type monuments.

Figure 14 .
Figure 14.Comparison between approximately co-located groundtype and roof-type 1-Hz CGPSs (detailed information about these pairs are shown in Table5and Figure13).

Figure 16 .
Figure 16.Seismic waves in radial component generated from resampled 10-Hz CGPS (brown), 1-Hz CGPS (green) and strong motion observations (blue) at zone CER, LV, and COS (an epicentral distance is positive northward with respect to the epicenter).

Figure 17 .
Figure 17.Seismic waves in transverse component generated from resampled 10-Hz CGPS (brown) and strong motion observations (blue) at zone CER, LV, and COS (an epicentral distance is positive northward with respect to the epicenter).
significant differences among various tectonic settings in eastern Taiwan (Central Range, Longitudinal Valley, and Coastal Range), where the derived waveforms suggest that later phases are trapped within the Longitudinal Valley.So, dense seismic and GPS network can contribute to more precise results for the interpretation of changing of the waveform features with different tectonic regime.

Figure 18 .
Figure 18.Power spectral density (PSD) of the seismic waveforms at the locations of the CER (brown), LV (green), and COS (blue).Note the stronger energy in the frequency range 0.5-1 Hz within LV than in the other two regions.

Table 2 .
Main feature of different adopted GPS processing strategies.

Table 4 .
Global information of the approximately co-located CGPS and strong motion stations (average).