Velocity Field Derived from Taiwan Continuous GPS Array ( 2007-2013 )

Data were collected from 281 Taiwan continuous Global Positioning System (cGPS) Array sites from 2007 2013 and processed with GAMIT/GLOBK software. Power spectral density stacking from cGPS position time series in Taiwan found the spectral index as -0.72, -0.77, and -0.57 for the E, N, U components, respectively. This indicates the cGPS data errors can be described as a combination of white noise and flicker noise. The common-mode errors are removed by stacking data from 50 cGPS sites with data periods greater than 5 years. By removing the common-mode errors the GPS data precision is further improved to 2.3, 1.9, and 6.9 mm in the E, N, U components, respectively. After strict data quality control, time series analysis and noise analysis, we derive a new Taiwan velocity field using cGPS data from 2007 2013. The general pattern of the newly derived 2007 2013 velocity field is quite similar to that from previous studies, but the station density is much larger and spatial coverage better. About 80 mm yr-1 plate convergence rate is observed, half of the rate is accommodated on the fold and thrust belt of western Taiwan and another half is taken up in the Longitudinal Valley and Coastal Range in eastern Taiwan. The velocities in western Taiwan generally show a fan-shaped pattern, consistent with the maximum compression tectonic stress direction. In northern Taiwan the velocity vectors reveal clockwise rotation, indicating the on-going extensional deformation related to the back-arc extension of the Okinawa Trough. In southern Taiwan, the horizontal velocity increases from about 40 mm yr-1 in the Chia-Nan area to 55 mm yr-1 in the Kao-Ping area with a counterclockwise rotation.


INTRODUCTION
Taiwan is situated in an active tectonic region with numerous thrust faults and folds due to on-going collision between the Luzon arc and Chinese continental margin.Numerous devastating earthquakes with magnitudes greater than 6 have occurred since 1900.The oblique convergence of the Eurasian and Philippine Sea Plates started at the beginning of Late Miocene period (Ho 1976; Barrier and Angelier 1986;Huang et al. 2006).
The Global Positioning System (GPS) has become an efficient tool for studying active tectonics and geodynamics (Dixon 1991;Hager et al. 1991;Feigl et al. 1993).Utilizing satellite positioning techniques each station can provide precise global coordinates for its antenna position which can be used to monitor the horizontal and vertical crustal movements at the site (Altamimi et al. 2002).We derived a velocity field using 2007 -2013 data from 281 sites in the Taiwan continuous GPS (cGPS) Array.The instrumentation and status of the cGPS array, data quality control and data processing procedures are well described.After strict and careful data quality control, times series analysis, noise analysis and common-mode error correction, we obtained a more reliable velocity field with more realistic uncertainty from the GPS position time series optimal model parameters.The results provide important information on the present-day crustal deformation in the Taiwan area.Further studies on the relation between seismic activity and strain accumulation or release can be carried out (Dixon 1991;Bock 1994;Hudnut 1995;Segall and Davis 1997).

TAIWAN CGPS ARRAY
In order to investigate the relationship between seismic activity and crustal deformation in the Taiwan area, the Central Weather Bureau (CWB) established a cGPS network composed of 16 sites in 1994 -1996.After the 1999 Chi-Chi earthquake, CWB increased the total number of cGPS stations to 150 sites from 2001 -2005.In the meantime, several agencies in Taiwan, including the Institute of Earth Sciences, Academia Sinica (IESAS), Central Geological Survey (CGS), and Ministry of the Interior (MOI), also installed many cGPS stations in the Taiwan area.Data collected from 2007 -2013 from a total of 281 cGPS stations operated by CWB, IESAS, and CGS are used in this study (Fig. 1).
All stations in the Taiwan cGPS array are equipped with dual-frequency geodetic Global Navigation Satellite Systems (GNSS) receivers and antennas.The sampling interval for data collection is either 30 or 1 sec.GPS data acquired at continuous stations are transferred to the central control system at main office by ADSL, wireless, or manual downloading.Some CWB cGPS stations are co-located with borehole broad-band seismic stations [e.g., I-Chu station (ICHU) and LONT] or borehole strainmeter sites (e.g., JSUI and NHSI).Furthermore, 77 cGPS stations are located within 1 km of a strong-motion seismic station.Thus, we can obtain coseismic displacements and ground motion from both GPS and seismic data at the same site.These data can therefore be used to compare high-rate GPS results with seismic wave observations.

DATA QUALITY CONTROL
Raw data from each station is converted to Receiver-INdependent-Exchange (RINEX) format.Sometimes RINEX files may give wrong information or format errors because of incorrect initial setup for receivers.The most common problems are illogical or garbled data format, incorrect receiver or antenna types.The Translation, Editing, and Quality Check (TEQC) software developed by UNAVCO (a non-profit university-governed consortium that facilitates geoscience research and education using Geodesy) is employed to fix these type of problems.TEQC can read original binary file format and output into RINEX common format.The program allows editing, modification, cut, and merge observational data and even reassign data sampling intervals.In addition, TEQC can also be used to check static and dynamic GPS observation data quality (Estey and Meertens 1999).
Most of the cGPS stations are located at relatively stable sites that exhibit little possibility of local movement and can be preserved for a long period of more than 10 years.Most sites have good sky visibility for elevation angles larger than 15 degrees.For example, the ICHU is a deeply anchored and braced monument station with very good data quality and sky visibility (Figs.2a and 3).In Fig. 2a the ticks outside the circle indicate the station's azimuth from 0 -360°, 0° is in the north direction.Blue dots describe the path of the received satellites on that day from elevation 0 -90°.The ICHU station sky-plots at the same day (May 31) in 2007 and 2013 reveal almost the same patterns and show good sky visibility (Fig. 2a).The proportion of ICHU receiving observations to predicted observations are 100 and 97% on May 31 in 2007 and 2013, respectively.Figure 2b shows the site photos taken toward the north, east, south, and west directions, respectively in 2013.The quality-control-time-series figure is another way to realize observation data qualities.Figure 3 presents the predicted and observed number of data (Fig. 3a), cycle slip numbers (Fig. 3b), and multipath errors for L1 and L2 (Fig. 3c) ICHU station.Basically, the ICHU site received more than 95% of the predicted number of observations.The multipath errors stay in the same level of less than 50 cm.Cycle slips occurred more frequently from 2010 -2013, which might be due to intensive solar activity.
Another example, the Shan-Da-Thi station (DASI), was installed in 2005 on a hill near the east coast of Taiwan.Figure 4a shows the DASI station sky-plots for 31 May in 2007 and 2013, respectively.We can see better sky visibility in 2007 and that fewer satellite signals were received in 2013 due to bad sky visibility, suggesting that the local environment at this site changed dramatically.The proportion of receiving observations to predicted observations for DASI also decreased from 82 -54% from 31 May in 2007 and on the same day in 2013.As is well known, data quality is usually relative to site sky visibility and good sky visibility comes from good site environment maintenance.Because the DASI site land owner grew trees near the antenna monument several years after the first installation, the surrounding tall trees explain the very bad sky visibility (Fig. 4b).Figure 5 presents observed data numbers (Fig. 5a), cycle slip numbers (Fig. 5b), and multipath errors of L1 and L2 (Fig. 5c) of DASI.The percentage of the total number of DASI received data is less than 60% (Fig. 5a) with a lot of cycle slips (Fig. 5b).Furthermore, multipath errors increased from 2007 -2013 and reached a maximum value in 2013.All of the QC results from Figs. 4a and 5 are consistent with the site environment changes shown in Fig. 4b.Therefore, the routine sky-plot and quality control are very important, convenient tools for the cGPS operation.

DATA PROCESSING
All cGPS data were processed by GAMIT/GLOBK software v.10.4 (Herring et al. 2010) with standard procedures.GAMIT uses the least-squares method to process data, and then obtain the site coordinates from daily solutions.GAMIT contains 5 modules, including ARC, MOD-EL, AUTCLN, CFMRG, and SOLVE modules.GLOBK uses Kaman filter to get a combined solution of several subnets or periods (Dong et al. 1998).The main GAMIT body functions and features are as follows: (1) ARC: ephemeris format conversion module, using partial differential to calculate orbit ephemeris.(2) MODEL: estimate theoretical value, the residual of observations and theoretical value will be considered as a reference parameter for the adjustment basis.(3) AUTCLN: fix cycle slip and correct model parameters.(4) CFMRG: generate the parameter configuration files.(5) SOLVE: solving equations.We first use the least-squares method to get the optimal model parameters and minimize residuals, then estimate the coordinates and orbital parameters of each station.These processes were generally executed 2 iterations (Fig. 6), and the normalized root-mean square value (NRMS) of the final results should be smaller than 0.25, otherwise we need to evaluate if any parameters are over-constrained in GAMIT processing.
Some GAMIT processing tables are very important and always need to be updated to maintain the correct information.For example: antenna mode table (antmod.dat),earth rotation parameters table (ut1., pole.), nutation table (nutabl.)…etc.The GAMIT and GLOBK processing control files need to be setup correctly.The main control files are as follows: (1) sites.defaults: the list of sites used for each different subnet.(2) Process.defaults:file for setup computation environment, including reference frame, data sampling rate, orbit parameters and other information.We use International Terrestrial Reference Frame 2008 (ITRF2008) as the reference frame.(3) Sittbl.: site control file, setup constrain of coordinate and velocity for each station.(4) Sestbl.: session control table, setup all kinds of parameters and modules.We set the orbit constrain as "RELAX."for loose orbit parameter constraint with an elevation cutoff angle of 15° to reduce any multipath effects and noises.We used the Vi-enna Mapping Function 1 (VMF1) for troposphere correction parameters, which is developed by the Victoria Vienna team.VMF1 is based on GMF (global mapping function) which developed from a numerical weather model (NWM) estimated using more than 20 years weather data (Boehm et al. 2006).The VMF1 can use different weights for different zeniths to conduct more accurate tropospheric error correction.(5) Station.info:site information table, includes site name, site code, receiver type, antenna type, data time span, and antenna height.This table needs to be updated when a station has an instrument change.(6) Apr file: a priori site coordinate files, set all initial station coordinates for GAMIT processing.The initial coordinates cannot have large bias from the site coordinates solution, or it will make the processing crash.All of the control files used for GAMIT processing in this paper can be found at the Geophysical Database Management System (GDMS) website (http://gdms.cwb.gov.tw/index.php) of CWB.
Usually we obtain loosely constrained parameter estimates and covariance solution files from the GAMIT processing results.These data are then passed to GLOBK for data combination to estimate station coordinates and velocities and orbital and Earth-rotation parameters (Feigl et al. 1993;Dong et al. 1998).GLOBK contains three major modules, globk, glred, and glorg.The functions of three modules are as follows: (1) globk: used for data combination in time domain or space domain.In other words, globk can combine decades of data into one solution, or merge the number of subnets into a daily solution, also estimate the site coordinates or velocity at the same time.It outputs four kinds of files, they are *.gcr (constrained, bias free solution), *.gcx (constrained, bias fixed solution),*.glr(loose, Fig. 6.GAMIT/GLOBK data processing flow chart.Including three parts, the pre-process and main-process for GAMIT part, and post-process for GLOBK part.Each section has its own main modules, input files and output files.bias free solution), and *.glx (loose, bias fixed solution).These files can be used for different purposes.(2) Glred: used for repeatability analysis.The major program is similar to globk, but glred takes the daily solution as independent and makes prominent of coordinate repeatability.(3) Glorg: merge GAMIT solution to get optimal coordinates by rotation, translation, and scaling.There are three steps in GLOBK processing (Fig. 6): (1) use glred first to obtain daily solution and time series of each site and remove outliers.(2) Without outliers, use globk/glorg to merge all solutions and estimate the optimal coordinates and velocities for each site.We use a loose constraint in this step to ensure that glorg has enough room to perform adjustment on rotation, translation and scaling.(3) Finally, glorg is used again to define a reference frame with tight constraints to achieve the optimal ITRF2008 coordinates.The International Terrestrial Reference System (ITRS) is realized and maintained by the International Earth Rotation and Reference Systems Service (IERS).The ITRS realization, through the International Terrestrial Reference Frame (ITRF) is updated regularly to take into account new accumulated data and also improved analysis strategies applied by the analysis centers of the contributed techniques.ITRF2008 is a refined version of the ITRF based on reprocessed solutions for the four space geodetic techniques: Very Long Baseline Interferometry (VLBI), Satellite Laser Ranging (SLR), GNSS, and Doppler Orbitography Radiopositioning Integrated by Satellite (DORIS), spanning 29, 26, 12.5, and 16 years of observations, respectively (Altamimi et al. 2011).
We processed a total of 281 Taiwan cGPS stations with time periods longer than 1.5 years and 25 International GNSS Service (IGS) stations.The precise IGS final orbits are used in the processing.Considering the computing capability and limitations of the number of sites for each subnet in GAMIT, we divided the 281 Taiwan cGPS stations into 12 sub-networks.In order to constrain the Taiwan subnet solution to the ITRF2008 reference frame, an additional 25 IGS station subnet was also used (Fig. 7).In order to connect all subnets we choose 5 local sites with longer data time span and good quality as common stations for all subnets.These common sites are S01R (Paisha, Penghu), S101 (Nankang, Taipei), S104 (Fushan, Taitung), LANY (Lanyu, Taitung), and HENC (Henchun, Pingtung) (Fig. 1).
The loosely constrained daily solutions from 13 GAMIT subnets are combined and constrained into the ITRF2008 reference frame through IGS station stabilization in GLOBK to obtain the precise position time series for each site.The coordinate constraints for 18 stabilization IGS stations are 5, 5, and 10 cm on X, Y, Z components, respectively, to avoid over-constraining in GLOBK processing.The correlation statistics between the solution precision and number of stabilization sites point out that the daily solution will usually be an outlier if less than 4 sites are used in GLOBK stabilization.This could be because GLORG needs at least 4 stabilization sites to adjust the optimal network geometry using translation, rotation and scaling.Figure 8 shows the statistics for 18 IGS sites used in GLOBK stabilization.The black, cyan and yellow color bars indicate RINEX archive data used in GAMIT, and the data used in GLOBK for stabilization, respectively.We can see that TSKB has long and complete raw-data and GAMIT solutions from 2007 -2013, but we do not use it for stabilization after the 11 March 2011 off-Tohoku earthquake.The FAIR, KOKB, and THTI sites are frequently excluded in stabilization due to their long baselines.Figure 9 is the statistics for the number of sites used in GLOBK stabilization from 2007 -2013.It reveals the minimal number of used sites is 8 in May 2012, and less than 10 days with only 9 stabilization sites.Twelve to seventeen sites are used for stabilization during the study period, suggesting the daily solutions are well constrained on the ITRF2008 reference frame.

Pre-Processing and Time Series Analysis
The velocity of each cGPS site is estimated from the position time series through a model that removes outliers and systematic errors.We transform the coordinates from ITRF2008 Cartesian coordinates into ENU (east, north, and up components) coordinates for easier site movement understanding in geophysical applications.The outliers need to be removed before time series analysis.All outliers were removed using the standard deviation (STD) and moving-average method (MA).The STD represents a mean data series values trend and most outliers can be easily removed using 3.5 times STD.For some special cases, e.g., the large offset caused by coseismic displacement or instrument change, the STD method may not work well because offsets will enlarge the STD value and reduce the capability of removing outliers.Therefore, combining STD with the MA method can make outlier removal complete.For the 7 years (2007 -2013) long GPS data the optimal window size is 30 days and shifts every 10 days.More than 98% of the 281 cGPS position time series had a removed outlier percentage less than 3% for whole data period.Nine percent of the 281 cGPS only removed less than 1% outliers, which means very good data quality, like the ICHU station (Fig. 10).Only a few sites presented large scattered position time series.For example, the DASI and FKDO stations show large scattered data distribution after the middle of 2011 caused by worse sky visibility (Figs.11 and 12).
The velocity of each cGPS site is estimated from the position time series through a model as follows (Nikolaidis 2002):  x is given different values for various stations and earthquake events to make the best time series fit.

Noise Analysis
In the past most studies took GPS time series observation errors as a random process with time-independent white noise.The time-independent errors can be reduced or eliminated through taking average of a large amount of data.Johnson and Agnew (1995) studied California cGPS data and found that if the time-dependent noises are not removed, we may take the noise as a strain variation with time signal and give rise to incorrect geophysical interpretation.Without considering the time-dependent color noise the error will be underestimated and the parameter precision overestimated (e.g., velocity field).The uncertainty of velocity components taking into account the time-dependent color noise can be four times the white noise only.Furthermore, the estimated velocity filed may also be slightly changed (Langbein and Johnson 1997;Zhang et al. 1997;Mao et al. 1999).In the low-latitude region the error is significantly enlarged, especially on the vertical component.The error could be underestimated as large as 10 times.The atmosphere in the low-latitude area contains abundant water vapor, which seri-ously affects GPS signal transmission.The water vapor is correlated with the seasons.The GPS position time series in Taiwan is usually more scattered during the summer.Therefore, we should consider the noise characteristics in the GPS position time series to be time-dependent.
GPS time series noise is a power-law process, its powerspectrum, P x (f), can be expressed as the following equation: Where f is temporal frequency, P 0 , f 0 are normalizing constant and k is spectral index.The spectral index k is used to define three kinds of color noise which is white noise (k = 0), flicker noise (k = -1), and random-walk noise (Brownian motion, k = -2).Noise can be a combination of different noise modes, i.e., white-noise-only model (WH), white noise + flicker noise model (FN), or white noise + random-walk noise model (RW).By stacking the power spectral densities from time series analysis residuals, we are able to estimate the spectral index to distinguish the noise characteristics (Tsai 2004(Tsai , 2013;;Tsai et al. 2007).Stacking the 281 Taiwan cGPS stations used in this paper, the spectra indexes for the east, north, and vertical components are -0.72,-0.77, and -0.57, respectively (Fig. 13).This result indicates the Taiwan cGPS data noise model should be white noise + flicker noise model (FN), consistent with the results from previous studies that k of GPS time series is usually between 0 and -1 (Zhang et al. 1997;Mao et al. 1999;Nikolaidis 2002;Yu et al. 2003).
A more realistic noise model can give a better estimate of the full covariance matrix and model parameters in GPS position time series.We used the Maximum-likelihood Method (MLE) to estimate the correlated noise amplitude (Williams 2003).The best-fit noise model (FN model) to the data is that which maximizes the probability function with the covariance matrix used on re-weighting the time series model parameters to obtain more realistic parameter errors.

Common-Mode Error
The GPS position time series can be used to better understand the spatial and temporal variations in crustal deformation.Therefore, it is important to identify the signal characteristics in the GPS time series.Removing the systematic errors can effectively improve the data precision and identify signals in the time series.Common-mode error is a kind of systematic error coming from poor constraint in data processing.Wdowinski et al. (1997) pointed out the common errors of all sites are the average from the summation of all station errors and it can be removed by subtracting the GPS position times series data deviation in three components.Tabei and Amin (2002) also indicated the GPS time series errors approach the same value when the total number of stations for common-error-estimation reaches a certain quantity.In other words, removing common error will not affect the velocity field result, but improve the data precision instead.For example, the seasonal GPS position time series signal errors can be effectually reduced by removing the common error.The common error can be obtained by stacking the GPS position time series model residuals for some specified stations with good data quality, uniformly disturbed, and long time period.The common-mode error i e can be represented as the following equations: Where ( ) and S i (t i ) are the time series analysis model residuals and number of stations at a certain time, t i .The common error is calculated only when the number of stations is greater than 3 at that day.Here i = 1, 2, 3…..N, indicates from day one to day N.
We chose 50 relatively stable stations with longer time period to estimate the common-mode error in the Taiwan area.Figure 14 present the data span for 50 sites for common-mode error estimation, with nearly all sites having data periods longer than 6.5 years.By stacking these 50 fiducial sites the common-mode errors are 1.3, 1.2, and 2.7 mm in the east, north, and vertical components, respectively.Figure 15 shows the common-mode error distribution in the E, N, U components.To evaluate if any common-mode error directionality occurs in the Taiwan area, we calculated the average distribution slope in Fig. 15 using linear regression.The results show the slopes are almost zero in the E, N, U components, indicating no significant common-mode error directionality.The root-mean-square (RMS) value time series for the common-mode errors show larger values during the summer time (Fig. 16).After common-mode error corrections we re-modeled all parameters for the 281 cGPS position time series.The results indicate the GPS data precision was been further improved by removing the commonmode errors.The standard errors for the 281 cGPS stations decreased from 2.7, 2.4, and 7.7 to 2.3, 1.9, and 6.9 mm in the east, north and up components, respectively.
Figure 17   for the other 281 stations used in this study can be found at the GDMS website (http://gdms.cwb.gov.tw/index.php) of CWB. Figure 18 presents the CHEN model residuals in the ENU, respectively.The random model residuals distribution indicates the CHEN time series position is well corrected.The model results linear rates represent a more accurate velocity field.Tables 1 and 2 give the 2007 -2013 ITRF2008 velocities and velocities relative to S01R (Paisha, Penghu) of cGPS sites.DASI and FKDO velocities are estimated using only data before 2011 to exclude the large scattering data due to bad sky visibility.

TAIWAN VELOCITY FIELD FROM 2007 -2013
We derived a 2007 -2013 velocity field based on the time series analysis results from 281 cGPS stations in the Taiwan area.Figures 19 -21 show the ITRF2008 velocities and velocities with respect to S01R in the horizontal and vertical components, respectively.The error ellipse on the tip of each velocity vector indicates 95% confidence interval of uncertainty for the horizontal velocity (Figs.19 and 21).On the other hand, the black cross represents various scales of rates and errors in vertical velocity.The red and blue color circles indicate uplift and subsidence, respectively (Fig. 20).It should be noted that stations with time periods less than 3 years and bad data quality are not shown.
The ITRF2008 horizontal velocity field in Fig. 19 is consistent with the plate motion direction which reveals the oblique collision between the Philippine Sea Plate and Eurasian Plate (Angelier 1986).The S01R station located at the Eurasia Plate are moving southeastward with a velocity of about 34 mm yr -1 .From the west coast to the Western Foothills (WF), the velocities show ESE direction motion at a rate of about 21 -33 mm yr -1 .The velocity decreases eastward to 0 -5 mm yr -1 in the Central Range.The LANY station, located at the Philippine Sea Plate, is moving northwestward with a rate of about 50 mm yr -1 .Velocities in the Coastal Range still reach 40 mm yr -1 in the NW direction.The significant velocity change is detected across the Longitude Valley Fault (LVF) which is the suture zone of two plates.Across the LVF there is about 20 -30 mm yr -1 velocity change.The velocities in northern and northeastern Taiwan are about 40 -60 mm yr -1 and reveal significant clockwise rotation.The velocity reaches a maximum value of about 60 mm yr -1 at SUAO (Suauo, Ilan).A dramatic change in velocities is found in the region between Hualien and Ilan, which is the place that the Philippine Sea Plate starts to subduct underneath the Eurasian Plate.In southern Taiwan the velocity field between the Chishan Fault (CSF) and Chaozhou Fault (CZF) shows obvious tectonic extrusion, moving in WSW-SSW directions with a rate of about 23 -44 mm yr -1 .
The ITRF2008 vertical velocity field in Fig. 20 indicates uplift sites are found mostly in the hill and mountain areas of the WF and Central Range.The Jiuchiunken-Muchiliao-Li-uchia Fault system (JMLF) seems to be a vertical velocities change boundary.There are about 3 -10 mm yr -1 uplift and 5 -14 mm yr -1 subsidence rates in the eastern and western

Site
Lon. (°) Lat.(°) V e (mm yr -1 ) V n (mm yr -1 ) V u (mm yr    sides of JMLF, respectively.Most of the large subsidence sites located at the western Coastal Plain could be caused by groundwater withdrawal (e.g., HUWE).Both the Longitudinal Valley and northern Coastal Range in eastern Taiwan show significant subsidence.In contrast, the southern Coastal Range is uplifting.The rates of uplift and subsidence are about 3 -11 and 2 -12 mm yr -1 , respectively.
Figure 21 shows horizontal velocity field with respect to S01R (Paisha, Penghu) which is located in the Chinese continental margin.The Luzon arc (Coastal Range, Lutao and Lanhsu) is moving northeastward at a rates of about 60 -84 mm yr -1 .The velocities in the Coastal Range are in the 307 -311° directions at a rate about 60 -69 mm yr -1 .Across the LVF, a major velocity discontinuity about 20 -30 mm yr -1 is observed.The movement direction is slightly changed from NW in the Coastal Range to WNW in the Central Range.Along the eastern side of the LVF the velocity decreases northward from about 68 mm yr -1 in Taitung to 38 mm yr -1 in Hualien.A clockwise velocity vector rotation is observed in the Ilan Plain and the velocity reaches its maximal value of 32 mm yr -1 at SUAO in the 141° direction.The Taipei Basin and surrounding area are relatively stable with minor crustal deformation and the velocities here are about 0 -5 mm yr -1 .
Approximately half of the plate convergence rate is accommodated on the fold and thrust belt of western Taiwan.The crustal motion shows nearly westward direction with rates of 35 -45 mm yr -1 in the Central Range and decreases to 0 -5 mm yr -1 in the Coastal Plain area.A significant velocity change of about 12 -20 mm yr -1 is detected across the fold and thrust belt where the JMLF and Chelungpu Fault (CLPF) zones are located.In southern Taiwan the horizontal velocity field becomes more complex.It is about 40 mm yr -1 in the Chiayi-Tainan area and then slowly increases to about 55 mm yr -1 in the Kaohsiung-Pingtung area with a counterclockwise velocity vector rotation.Since the ITRF2008 S01R vertical velocity is only 0.5 mm yr -1 , the relative vertical velocity field pattern with respect to S01R is quite similar to the ITRF2008 vertical velocities, as shown in Fig. 20.
Generally, the relative velocity field from 2007 -2013 exhibits the characteristics of different stages for the Taiwan collision process from south to north: pre-collisional with rapid and distributed convergence (southern Taiwan), collision and suturing (western-central Taiwan), and postcollisional with collapse and extension (northern-northeastern Taiwan) (Lallemand and Tsien 1997;Shyu et al. 2005).The fan-shaped velocity field (Fig. 21) indicates the ongoing collision is resisted by the Peikang High, a Pre-Miocene basement, which acts as a buttress for advancing thrusts and creates a frontal thrust fault system parallel to the basement margin (Tsan and Keng 1968;Mouthereau and Petit 2003).The velocity vector directions are perpendicular to the deformation front (the thrust and fold belt).

DISCUSSION AND CONCLUSION
We used GAMIT/GLOBK to process 2007 -2013 Taiwan cGPS data collected by the CWB, IESAS and CGS from 281 sites.After time series analysis, noise analysis and common-mode error correction, we obtained high precision GPS position time series and velocity fields with more realistic uncertainty.
The GPS position time series is used for velocity estimation and denotes the relationship between the site environment and data quality.Some geophysical signals related to the groundwater variation and seismic activity may be detected.The GPS position time series of the SUCH site (Shuangqi, Taichung) presents obvious periodic signals in the E, N, U components (Fig. 22a) and a scattering of data points changed with time.We found this periodic motion is a non-tectonic signal and is in fact due to the changes in site environment.The SUCH station was established in 2005 and a local government agency kept adding lots of decorations and grew trees nearby during the past few years.Figure 22b shows photos taken in 2007 showing the site environment was fine even with some decorative statues in NE GPS antenna monument direction.However, the site photos taken in 2013 (Fig. 22c) reveal big changes in site environment, with more buildings and trees surrounding the station.These obstructions reduce sky visibility and increase GPS position time series data scattering.Furthermore, the sakura trees next to the sites make sky visibility worse during the flower season, which is the obvious source that causes periodic signals in the SUCH position time series.The other two cGPS sites, CLON (Changlong) and NJOU (Nanzhou), are located in the Pingtung Plain.Both sites show remarkable periodic signals with a similar pattern in the GPS position time series U component (Figs. 23 and 24).These two stations are only about 8 km apart and both located at the footwall of the CZF.Comparing with the water level time series of the Qhishan groundwater station (operated by CWB; Fig. 25), the GPS position time series of CLON and NJOU in the U component show similar patterns with the groundwater data.The GPS site displacements and groundwater station levels present a decrease in the first half of the year and increase in the last half of the year.This result indicates a periodic variation in GPS time series in the U component for CLON and NJOU is a seasonal signal affected by groundwater level change.
The general pattern of the newly derived 2007 -2013 velocity field from the Taiwan cGPS Array is quite similar to that from previous studies (Yu et al. 1997(Yu et al. , 1999;;Hsu et al. 2009;Lin et al. 2010;Ching et al. 2011).The station density is much larger and spatial coverage is better.This indicates the complexity of tectonic structure in Taiwan.There are two major velocity discontinuities of about 20 -30 and 12 -20 mm yr -1 across LVF in eastern Taiwan and frontal thrust faults in western Taiwan, respectively.In northeastern Taiwan, the velocity vectors reveal a clockwise   rotation pattern.This indicates the area is undergoing extensional crustal deformation related to the back-arc extension of the Okinawa Trough (Kimura 1985).The velocity of western Taiwan generally shows a fan-shaped pattern, which is consistent with the maximum compression tectonic stress direction inferred from focal mechanisms (Yeh et al. 1991), borehole break-out data (Suppe et al. 1985), and geological strain field (Lee and Wang 1988;Lu et al. 1989;Angelier et al. 1990).
In southern Taiwan the velocities increase southward with the directions becoming southwesterly.It reaches 50 -55 mm yr -1 in the Kaoshiung-Pingtung area and deflects about 30° counterclockwise in the direction near the Kaoshiung-Pingtung coast.This movement could be related to tectonic extrusion.Widespread soft sediments and mudstone resulting in plastic deformation may be the reason that rapid crustal deformation but low seismic activity are found in this area.

Fig. 1 .
Fig. 1.Taiwan continuous Global Positioning System (cGPS) Array.The green triangle, yellow circle, and red square symbols indicate the cGPS sites operated by the Central Geological Survey (CGS), Central Weather Bureau (CWB), and Institute of Earth Sciences, Academia Sinica (IESAS), respectively.Pink lines are active faults.

Fig. 3 .
Fig. 2. (a) Sky plots of ICHU station on May 31, in 2007 and 2013, respectively.The ticks outside of circle indicate the station's azimuth from 0 -360°, 0° is in the north direction.Blue dots describe the received satellites path at that day from elevation 0 -90°.(b) Station photos taken toward the north, east, south, and west directions, respectively in 2013.

Fig. 5 .
Fig. 5.The quality-control time series of DASI station.(a) The green square and blue diamond symbols indicate predicted and observed number of observations, respectively.(b) Gray bars indicate number of cycle slips in logarithmic scale from 2007 -2013.(c) The multipath errors, MP1 and MP2 of L1 and L2, respectively.
Fig. 7. Distribution of IGS stations used in GMAIT/GLOBK processing.The triangular indicates the stabilization sites used in GLOBK, and the square represents the IGS sites only used in GAMIT processing.Black dash line denotes plate boundaries.

Fig. 9 .
Fig. 9.The statistics for the number of sites used in GLOBK stabilization from 2007 -2013.

Fig. 12 .
Fig. 12. GPS position time series of FKDO station in the ENU, respectively.The black solid circles and light gray bars indicate the observed daily solutions and their formal errors, respectively.
shows the position time series with FN model and common-mode error correction for the CHEN station (Chenkung, Taitung) in the ENU, respectively.The blue solid circles and light gray bars indicate the observed daily solutions and their formal error.The pink line indicates model predictions.The dashed line denotes the occurrence time for the 2008 Taitung earthquake (15 April 2008, M L = 5.1).The smoothing model predictions fit the observed GPS position time series very well.The GPS position time series figures

Fig. 13 .
Fig. 13.Stacking the power spectral density of 281 cGPS sites.The slopes of -0.72, -0.77, and -0.57indicate the spectral index for the east, north and up components, respectively.

Fig. 14 .
Fig. 14.The data time period of 50 fiducial stations used for stacking common-mode errors.The black dots indicate data are available for the cGPS sites at that time.

Fig. 16 .
Fig. 16.RMS distribution of common-mode errors in the east, north and up components, respectively.

Fig. 17 .
Fig. 17.CHEN station position time series analysis results in the east, north and up components, respectively.The blue solid circles and light gray bars indicate the observed daily solutions and their formal errors; the pink line indicates model predictions.The dashed line denotes the occurrence time of 2008 Taitung earthquake (M L = 5.1).

Fig. 19 .
Fig. 19.Horizontal ITRF2008 velocities of cGPS stations in the Taiwan area (2007 -2013).The green lines indicate the surface traces of active faults.JMLF and LVF denote the Jiuchiunken-Muchiliao-Liuchia Fault system and Longitude Valley Fault, respectively.CSF and CZF denote the Chishan Fault and Chaozhou Fault, respectively.WF denotes Western Foothills.The blue arrows show velocity vectors.A 95% confidence error ellipse is shown at the tip of each velocity vector.The purple dashed line describes the mirage of Peikang High.

Fig. 20 .
Fig. 20.Vertical ITRF2008 velocities of cGPS stations in the Taiwan area (2007 -2013).The circle size is proportional to the vertical motion magnitude.JMLF, LVF, CSF, CZF, and WF have the same definition as Fig. 19.Blue circles denote subsidence and red circles denote uplift.The black crosses are STDs and the purple dashed line describes the Peikang High mirage.

Fig. 21 .
Fig. 21.Horizontal velocities of cGPS stations in the Taiwan area w.r.t.Paisha, Penghu (2007 -2013).The green lines indicate the surface traces of faults.JMLF, LVF, CSF, CZF, and WF have the same definition as Fig. 19.CLPF denotes the Chelungpu Fault.The blue arrows show velocity vectors.A 95% confidence error ellipse is shown at the tip of each velocity vector.The purple dashed line describes the Peikang High mirage.

Fig. 23 .
Fig. 23.GPS position time series of CLON station in the east, north and up components, respectively.The black solid circles and light gray bars indicate the observed daily solutions and their formal errors, respectively.

Fig. 24 .
Fig. 24.GPS position time series of NJOU station in the east, north and up components, respectively.The black solid circles and light gray bars indicate the observed daily solutions and their formal errors, respectively.
Fig. 22.(a) GPS position time series of SUCH station from 2007 -2013.(b) and (c) station photos taken in 2007 and 2013, respectively.
Note: Lon. and Lat. are the longitude and latitude of GPS sites; V e , V n , and V u are the ENU of ITRF2008 station velocities with standard errors, respectively; period shows the start and end time of observation period.
Note: Lon. and Lat. are the longitude and latitude of GPS sites; V h and V u are the horizontal and vertical components of station velocities with standard errors, respectively; Azi is the azimuth of the horizontal velocity with standard errors; period shows the start and end time of observation period.