Tectonic effect for establishing a semi-dynamic datum in Southwest Taiwan

To accommodate the effects of crustal deformation in the current national static geodetic datum (Taiwan Geodetic Datum 1997 (TWD97)) in SW Taiwan, 221 campaign-mode global positioning system (GPS) stations from 2002 to 2010 were used in this study to generate a surface horizontal velocity model for establishing a semi-dynamic datum in SW Taiwan. An interpolation method, Kriging, and a tectonic block model, DEFNODE, were used to construct the surface horizontal velocity model. Forty-four continuous GPS stations were used to examine the performance of the semi-dynamic datum through exterior validation. The average values of the residual errors obtained using the Kriging method for the north and east components are ±1.9 and ±2.2 mm/year, respectively, whereas those obtained using the block model are ±2.0 and ±2.9 mm/year, respectively. The distribution of residuals greater than 5 mm/year for both models generally corresponds to a high strain rate area derived using the horizontal velocity field. In addition, these residuals may result from deep-seated landslide and active folding or mud diapir in a mudstone area. Similar exterior checking results obtained using the Kriging interpolation method and block model for SW Taiwan indicate a high station density and a relatively satisfactory station spatial coverage. However, the block model is superior to the Kriging method due to the consideration of characteristics of the geological structure in the block model. In addition, result from traditional coordinate transformation was used to compare with the semi-dynamic datum. The results indicate that a semi-dynamic datum is a feasible solution for maintaining the accuracy of TWD97 at an appropriate level over time in Taiwan.


Background
Taiwan is located at the present-day plate convergence boundary between the Eurasian plate and the Philippine Sea plate (Fig. 1); earthquake activities are abundant and the convergence rate is approximately 82 mm/year across the island, reflecting a high surface strain rate of up to 1 μstrain/year Bos et al. 2003;Chang et al. 2003;Byrne et al. 2011). However, the current Taiwan Geodetic Datum (Taiwan Geodetic Datum 1997 (TWD97)), announced in 2012, is considered a static geocentric datum, and it complied with the International Terrestrial Reference Frame 1994 (ITRF94) in the first epoch of 2010 (MOI 2012). The baseline accuracy in a horizontal component aimed at TWD97 is 2.0 mm + 1 ppm × k, where k is the length of the baseline, which corresponds to approximately 0.05 m horizontally. In other words, this datum will be distorted and loose its accuracy after 4-5 years because of the presence of a velocity gradient of up to 1 μstrain/year across the geodetic network. Therefore, the main concern in constructing a modern national geodetic datum in Taiwan is determining how to accommodate the effects of crustal deformation to maintain a spatial accuracy of its coordinates during a long period.
In general, establishing a semi-dynamic datum by assigning deformation models to a specific geodetic datum, such as the Japanese Geodetic Datum 2000 (Tanaka et al., 2007) and the New Zealand Geodetic Datum 2000 (Grant and Blick 1998;Beavan and Haines 2001;Beavan and Blick 2007), is an appropriate means for sustaining the coordinate spatial accuracy for a national-based geodetic datum at the plate boundary to account for the high surface strain (Tregoning and Jackson 1999). The coordinate accuracy requirements have been achieved over time by using this deformation model. The long-term surface velocity resulting from interseismic plate tectonic loading and permanent surface displacement caused by discrete events, such as an earthquake, are generally included in the deformation model to reflect the true deformation field with adequate accuracy and resolution (Beavan and Blick 2007). Therefore, constructing a semi-dynamic datum is a solution for maintaining the accuracy of TWD97 at an appropriate level over time in Taiwan.
More than 750,000 people reside in SW Taiwan, and a high velocity gradient has been observed in this area (Ching et al. 2007(Ching et al. , 2011b (Fig. 1). Previous studies have inferred a high contraction rate of approximately 1.0 μstrain/year with a clockwise rotation of 14.5°-27.1°/myr, and right-lateral shearing from sparse global positioning system (GPS) horizontal velocities in the fold-and-thrust belt of SW Taiwan (Bos et al. 2003;Chang et al. 2003;Ching et al. 2007Ching et al. , 2011bHsu et al. 2009). In addition, the earthquake record from the Central Weather Bureau of Taiwan indicates that no M L > 6.0 earthquakes have occurred in our study area since 1900, except for the 2003 M w 6.8 Chengkung earthquake in eastern Taiwan and the 2010 M w 6.4 Jiasian earthquake occurred close to the study area. Most of the coseismic displacements caused by those two earthquakes are smaller than 5 mm (e.g., Ching et al. 2011a). In other words, no major postseismic deformation has disturbed the secular motion in this region. Consequently, SW Taiwan is an excellent area for establishing a semi-dynamic datum by estimating a horizontal velocity field to maintain the coordinate accuracy of sites. The purpose of this paper is to provide a horizontal velocity model based on 265 GPS observations for estimating and predicting the coordinate changes associated with the horizontal crustal motion in SW Taiwan (Fig. 1). The interpolation method, the tectonic block model, and traditional coordinate transformation were used to analyze the accuracy in estimating and predicting the coordinate changes at arbitrary sites. This paper presents a comparison of efficiency between the interpolation method and block model.

Geological background
The major geological characteristic in SW Taiwan is the thick Plio-Pleistocene mudstone, Gutingkeng Formation. The Gutingkeng Formation dominates the study area, and its presence is proposed to be responsible for the nonoccurrence of large earthquakes in SW Taiwan since the last century. The Gutingkeng Formation consists of gray sandy siltstone and sandy mudstone intercalated with lenticular greywacke and subgraywacke with abundant Mollusca (Chou 1971).
Six major active faults, the Liuchia-Muchiliao fault system (LCMF), the Hsinhua fault (HHAF), the Houchiali fault (HCLF), the Hsiaokangshan fault (HKSF), the Chishan fault (CHNF), and the Fengshan transfer fault zone (FTFZ) from north to south ( Fig. 1), account for the surface motion observed in this study. For the LCMF (number 1 in Fig. 1), although this fault system has remarkable landscape feature based on studies of aerial photographs, the evidence of separation in strata is absence in the field (Sun 1971). The northern segment of the LCMF (the Muchiliao fault) is a blind thrust and might have reactivated during Late Quaternary (Shyu et al. 2005). The southern segment of the LCMF (the Liuchia fault) is also a thrust fault and cuts across Holocene sediments. For the HHAF (number 4 in Fig. 1), the length of unambiguous surface rupture associated with the 1946 M L 6.3 Hsinhua earthquake is about 6 km (Bonilla 1977). The high dip angle of about 70°to 90°is deduced by the surface trenching data (Lee et al. 2000) and the focal mechanism solution of the 1946 Hsinhua earthquake (Cheng and Yeh 1989). For the HCLF (number 3 in Fig. 1), a normal fault associated with the growth of diapiric fold has first been proposed using the seismic gravity, drilling well and shallow seismic survey data (Hsieh 1972;Kuo et al. 1998). In addition, a reverse fault has also been proposed using the geophotograph technique (Sun 1964) and the repeated geodetic survey and D-InSAR results (Fruneau et al. 2001;Huang et al. 2006Huang et al. , 2009Mouthereau et al. 2002). For the HKSF (number 6 in Fig. 1), an obvious fault scarp indicates that the length of the NNE-striking HKSF is approximately 8 km and it dips to the east (Hsu and Chang 1979;Sun 1964). For the CHNF (number 7 in Fig. 1), the NE-trending CHNF is a reverse fault with left-lateral components according to slickenside analysis along the fault plane (Chen 2005). However, the fault analysis suggests a normal displacement on the CHNF (Gourley et al. 2012). In addition, GPS measurements indicate that the CHNF is a reverse fault with a rightlateral strike-slip component (Ching et al. 2007(Ching et al. , 2011bHu et al., 2007;Lacombe et al., 2001). For the FTFZ (number 9 in Fig. 1), the N140°E-striking FTFZ has been detected by analyzing a geomorphic feature acquired from the Digital Elevation Model (DEM) (Deffontaines et al. 1994;1997;Lacombe et al. 1999) and by conducting GPS surface velocity analysis (Ching et al. 2007).

GPS data acquisition and processing GPS data
GPS observations from 221 campaign-mode GPS stations used in this study were installed by the Central Geological Survey and 44 continuously operating GPS stations established by the Central Weather Bureau, the Academia Sinica, and the Central Geological Survey in the study area in SW Taiwan from 2002 to 2010 (Fig. 1). GPS surveys are generally conducted annually. A benchmark is usually occupied by more than two sessions per year. Each session is 6-14 h, and all available satellites are tracked and rising higher than a 15°elevation angle. The sampling interval for data logging is 15 s. The 44 continuous GPS stations are used for examining the accuracy of the horizontal velocity model established by the 221 campaign-mode GPS stations in the study area.
The campaign surveying and continuously recorded GPS data were processed session by session and date by date, respectively, by using Bernese software v.5.0 (Dach et al. 2007) to obtain precise station coordinates. Precise ephemerides provided by the International GNSS Service (IGS) were used and fixed during processing. Five global IGS fiducial stations surrounding Taiwan (IRKT, TSKB, GUAM, PERT, and IISC) on the International Terrestrial Reference Frame (ITRF2008) (Altamimi et al. 2011) were used to determine the coordinates of all campaign-mode and continuous GPS stations in the study area. The horizontal uncertainties of coordinates were estimated to be 2-5 mm.

Horizontal velocity field in SW Taiwan
The GPS horizontal velocities were estimated on the basis of the coordinate time series in a time span of 9 years from 2002 to 2010. The horizontal velocities for the east and north components were estimated using least squares estimation (LSE). An empirical equation y = a + bt + Hc for the study area was used to determine the standard deviation of the horizontal velocity field, where y is the coordinate in each component, a is the offset from origin, b is the velocity in each component, c is the coseismic offset of the 2003 M w 6.8 Chengkung earthquake in eastern Taiwan and the 2010 M w 6.4 Jiasian earthquake, and H is a step function (Fig. 2). The formal uncertainties (σ ls ) obtained using the LSE of velocities were scaled by the amount k = (mis/2) 2 (Ching et al. 2011b), where k is the effect of daily coordinate variation presented in the coordinate time series and mis is the residual of calculations and observations for the time series for the east and north components, respectively. The velocity uncertainty σ for each component was re-estimated using the equation σ = (σ ls 2 × k) 1/2 (Ching et al. 2011b).
GPS horizontal velocities are relative to the station S01R at Penghu Island in the Chinese continental  A velocity gradient appears between the LCMF-HHAF-HCLF and the CHNF (Fig. 1).

Establishment of surface velocity model
In this study, the interpolation method and tectonic block model were used to construct a surface horizontal velocity model. In general, the spatial variation of surface horizontal velocity is caused by the interaction between fault coupling and plate motion (e.g., McCaffrey 2002). A tectonic block which is a physical model containing the conservation of momentum and geological constraints (e.g., McCaffrey 2002) is therefore adopted in this study. On the other hand, a continuous surface velocity field is shown in our study area because most of active faults are locked at the surface. Hence, an interpolation method which is mathematical method is also properly adopted to obtain a surface velocity model in this area.

Interpolation method
Various interpolation methodologies have been developed to construct new data points within the range of a discrete set of known data points using by curve fitting or regression analysis. Because the surface velocities containing a geographic feature are noisy data, a non-linear interpolation is needed not only for fitting an interpolant that passes through the given data points but also for regression. The Kriging interpolation which is mathematically closely related to regression analysis has been developed in geographic statistics for the estimation and prediction of spatial data (e.g., Roush et al. 2003;Samsonov and Tiampo 2006;van Dalfsen et al. 2006;Samsonov et al. 2007;Miura 2010;Majdański 2012). In the Kriging method, four frequently used semi-variance functions, exponential, Gaussian, linear, and spherical models, related to the distances among sites can be employed to express the spatial variation, and they minimize the error of predicted values which are estimated by spatial distribution of the predicted values (Wackernagel 2003). Therefore, the Kriging interpolation is not dependent on given data points but on the data configuration and variogram parameters (Goovaerts 1997;van Dalfsen et al. 2006). Analyzing the root mean square (RMS; Eq. 1) of the differences between modeled and observed velocities at the 44 continuous GPS stations revealed that the optimal function is the exponential model. The results of RMS analyses for the four aforementioned semi-variance functions are shown in Table 1. The optimal performance relative to other functions is obtained using the exponential model ( Eqs. 2 and 3).
where n is the number of sites and v model i and v observed i are the −modeled and observed velocities at the 44 continuous GPS stations, respectively.
where h is the lag distance between data point pairs, C 0 is a nugget effect, C is a scale, and a is a range. The horizontal velocities of the 221 campaign-mode GPS stations between 2002 and 2010 were employed with an assumption of equal weighting, and a horizontal velocity model with grids with 2.5-km spacing (approximately 90 s in both latitude and longitude) in the study area (70 km (E-W) and 80 km (N-S)) was generated using the Kriging method, and the exponential model with the range a = 15 km (Eq. 2) due to the average of distances among the 221 campaign-mode GPS sites is approximately 5 km (Fig. 3). The critical semi-variance was set as C 0 + C (Eq. 2 if h→∞) as the points were located at the exterior area of the observed sites, that is, 221 campaign-mode GPS stations. Moreover, to calculate the velocities at arbitrary points, the gridded velocities at the four points that form the cell in which the point lies were again interpolated using the bilinear method.

Tectonic block model
The interseismic velocity field at the plate boundary zone is generally dominated by tectonic block rotations and interseismic coupling on faults (Savage and Simpson 1997). Based on this concept, an elastic kinematic block model with the code DEFNODE, developed by McCaffrey (2002), was adopted in this study to construct a surface horizontal velocity model for SW Taiwan. The nonlinear The optimal performance is obtained using the exponential model in this study inversion in this program involved applying simulated annealing to downhill simplex minimization (Press et al. 1989) to invert GPS horizontal velocities simultaneously for Euler pole locations and the angular velocities of tectonic blocks and coupling coefficients on the blockbounding faults. The coupling coefficient is used to describe the velocity gradient, which is caused by the friction on the fault plane, across the fault. In this inversion, a constraint that decreases the coupling coefficient down-dip from being totally stuck at the surface to totally creeping at the bottom of the fault was imposed because most terrestrial faults in SW Taiwan exhibit no clear evidence of aseismic surface creeping. The contribution of fault coupling to the velocities was calculated using the formulations of Okada (1985) in an elastic half-space material. The best fit parameters were determined by minimizing data misfit, defined by the reduced chi-square statistic (χ 2 ) (McCaffrey 2002). A χ 2 value considerably greater than 1 indicates a relatively poor fit for the model, and a value close to 1 indicates an acceptable model fit.

Tectonic blocks
Eleven elastic tectonic blocks in SW Taiwan were identified according to the concept proposed by Ching et al. (2011b) (Fig. 4). Because the trends of the GPS velocity gradients are generally parallel to the strikes of active faults in this region (red lines in Fig. 1), the block boundaries (black lines in Figs. 4 and 5) are mainly selected by determining whether the locations of GPS velocity gradients and active faults match. The blocks used in this study are closed, spherical polygons on the Earth's surface and cover the entire model domain. All points within a block are assumed to rotate with the same angular velocity.

Fault configurations
The contribution of coupling on ten active faults to surface velocities in SW Taiwan was estimated according to the concept proposed by Ching et al. (2011b) in this study ( Figs. 1 and 5). Faults with a slip rate are represented by a series of node points with 3-D shapes within  Ching et al. (2011b), are fixed in this model (Fig. 5). Faults are assumed to extend to a depth of 5-10 km, below which they creep at the full relative plate velocity. The exact dip angles for most reverse faults remain undetermined, and the lower edges of east-dipping thrusts in SW Taiwan increase from west to east with a depth from 5 to 7 km and with dip angles from approximately 25°to approximately 45°. In addition, because the depths and dip angles of two nearly E-W-striking strike-slip faults are still unclear, the fault depths of the Hsinhua (number 4 in Fig. 1) and Fengshan transfer faults (number 9 in Fig. 1) were set as 10 km, and the nearly vertical fault planes (dip angles between 90°and 70°) were assigned on the basis of limited balanced cross sections (Huang et al. 2004).

Kriging method interpolating results
The RMS interior , denoting the residual errors of the horizontal velocities of the 221 campaign-mode GPS stations, was evaluated using Eq. 1. The average RMS interior is ±0.3 mm/ year for both north and east components (Fig. 6a). Model residuals are generally lower than ±2.0 mm/year (Fig. 6a). However, the residual for the CHNF area is greater than ±3.0-4.0 mm/year (number 7 in Figs. 1 and 6a).

Tectonic block modeling results
The preliminary locations of Euler poles for each tectonic block were calculated by inverting horizontal velocities without considering the effect of elastic strain because of the locking on the faults. The locations of Euler poles and the angular velocities of the aforementioned inversions serve as a priori information in inversions for motions of tectonic blocks and interseismic coupling on faults (Fig. 5). The reduced χ 2 statistic of the optimized model is 1.29. Model residuals are generally lower than 2.0 mm/ year (Fig. 6b). However, the residual around the CHNF region is greater than 4.0 mm/year (number 7 in Figs. 1  and 6b). The RMS interior residual errors of the 221 campaign-mode GPS stations are ±2.8 and ±2.0 mm/ year for the east and north components, respectively.

Coordinate transformation
Conventionally, the two-dimensional coordinate transformations are used to transform the coordinates of stations from the observed time t 0 to a specific epoch t 1 by applying a fourparameter transformation (Helmert transformation) (Ghilani 2010) or six-parameter transformation (affine transformation) (Chang 2004) to reduce the horizontal distortions caused by the crustal deformation in a local area. The six-parameter transformation is more commonly used than the fourparameter transformation because it uses more geometric parameters to absorb the horizontal deformation. Therefore, six-parameter transformation was adopted in this study.
The six-parameter transformation (Eqs. 4 and 5) uses 33 control points (Fig. 6a), which were selected from the 221 campaign-mode GPS stations in the study area ( Fig. 1) by using the LSE method.
where (x,y) and (x′,y′) are the coordinates of the control points in the epochs t 0 and t 1 , respectively, and a, b, c, d, f, and g are the six parameters of the coordinate transformation.
To avoid the contamination of coseismic displacements resulting from the 2003 M w 6.8 Chengkung earthquake and the 2010 M w 6.4 Jiasian earthquake, the period adopted for the coordinate transformation experiment was 2004-2009. Therefore, the six parameters were evaluated for five periods, 2004-2005, 2004-2006, 2004-2007, 2004-2008, and 2004-2009, by using the LSE method. In this study, the observed coordinates of the 44 continuous GPS stations in the study area (Fig. 1) were compared with their corresponding transformed coordinates to analyze the accuracy of the coordinate transformation. The RMS was evaluated using Eq. 1 for the residuals in different epochs (Table 2 and Fig. 7). The RMS of the residuals quickly increased from ±7.8 mm in 2004-2005 to ±19.1 mm in [2004][2005][2006][2007][2008][2009]. Based on these results, coordinate transformation is not an appropriate method for reducing the distortion caused by crustal deformation between different time epochs in a local area. The longer the time interval is, the greater the residuals of the coordinate transformation are. Moreover, the distribution of control points and surface deformation are the main limitations of coordinate transformation. In a local area, uniformly distributed control points and uniform surface deformation are essential for achieving a high-quality coordinate transformation. The 33 control points used in this study were not uniformly distributed (Fig. 7a). Therefore, coordinate transformation could not satisfactorily reduce the horizontal distortions in this study.

Precision of velocity models
In this study, the Kriging interpolation method and block model were used to establish the horizontal velocity model for SW Taiwan. The average residual RMS interior of the Kriging method is considerably lower than that of the block model. Therefore, performance evaluations of the two methods were conducted to determine which is more suitable for establishing a surface horizontal velocity model for SW Taiwan.
First, the signals in surface deformation are composed of the tectonic block motion, coupling coefficients due to friction on the fault plane, seasonal variation, deepseated landslide, volcanic activity, artificial groundwater withdrawal, and so on (e.g., Murray-Moraleda 2011). The temporal variations in the coordinate time series are stable from some of aforementioned mechanisms, such as the tectonic block motion and locking on the fault plane caused by tectonic loading. The first-order pattern of a horizontal velocity field generally represents tectonic movements. However, the temporal variation of the coordinate time series resulting from deep-seated landslide or groundwater pumping may not be stable (Fig. 8). Those mechanisms temporarily change the orientation of the velocity vector. Block model of velocity tend to consist of tectonic block motions and interseismic coupling on faults. In other words, modeling events caused by inelastic deformation, such as folding, deep-seated landslide, and artificial groundwater pumping, is difficult. By contrast, the Kriging interpolation method attempts to minimize prediction errors, which are themselves estimated (Oliver and Webster 1990). Therefore, the performance of the Kriging method in minimizing residuals is superior to that of the block model.
Next, to evaluate the precision performance of the two methods, 44 continuous GPS stations observed from 2002 to 2010 in the study area were adopted for exterior checking (blue arrow lines in Figs. 1 and 6). The RMS exterior values of the residual errors of horizontal velocities from 44 continuous GPS stations were calculated using Eq. 1. The RMS exterior values of the Kriging method for the north and east components are ±1.9 and ±2.2 mm/year, respectively, and those of the block model are ±2.0 and ±2.9 mm/year, respectively (Fig. 6). The RMS exterior of the two models is comparable in value and spatial distribution (Fig. 6). In addition, the RMS interior and RMS exterior of the block model are almost consistent. Compared with the Kriging method, the block model is a physical model and is constrained by geological parameters, such as the locations of faults, dips of faults, and sense of fault slips. Therefore, under a suitable geological constraint, the RMS interior and RMS exterior of the block model are consistent. However, for the Kriging method, favorable RMS exterior values of ±1.9 and ±2.2 mm/year for the north and east components, respectively, are based on dense spatial coverage of the horizontal velocity field.
Although the Kriging method can easily establish a velocity model, a high station density (e.g., the approximately 5 km spacing in this study area) is required for providing a satisfactory spatial coverage condition for interpolation in SW Taiwan. In addition, a long dataacquisition duration is also required to minimize contamination from inelastic deformation or artificial sources. Conversely, constructing a block model of velocity is relatively difficult because more geological information is required to build the block boundaries and fault parameters (McCaffrey 2002). However, fewer observations are required in the tectonic block model because of the physical constraints originating from geological investigations.
Finally, according to the patterns in the RMS exterior of the Kriging method and RMS interior and RMS exterior of the block model (Fig. 6), the distribution of residuals greater than 5 mm/year along blocks B, C, D, E, F, and H of the block model (Fig. 6b) generally corresponds to the high strain rate area derived from the horizontal velocity field, particularly for vectors in blocks D, E, F, and H (Fig. 6). Notably, orientations of some velocity vectors in blocks B and C, such as those of stations S336 and G126, are sub-perpendicular to the strikes of mountains (Fig. 4). Although there is no strong evidence proving that these vectors are contaminated by the deep-seated landslide, the poor data fitting of the block model in blocks B and C and the orientations of these vectors sub-normal to the strikes of mountains may imply that deep-seated landslide disturbs the horizontal velocities. In addition, folding or mud diapir in a mudstone area, which acts as a non-recoverable, inelastic strain, is a crucial structural feature in blocks D, E, F, and H (Hsieh 1972;Lacombe et al. 1997Lacombe et al. , 2004Pan 1968;Shih 1967) (Fig. 6). Therefore, modeling surface deformation caused by active folding or mud diapir in a mudstone area is difficult by using an elastic block model. In addition to the possibilities of deep-seated landslide and active fold or mud diapir in a mudstone area, these residuals may be attributable to poorly defined block boundaries or unclear fault parameters, because a detailed geological investigation is still required to clarify the present-day characteristics of active structures.

Conclusions
A national geodetic datum is crucial in studying Earth science, establishing basic infrastructure, developing technology, and conducting academic analyses in a country. Therefore, to consider the tectonic effect in the Taiwan national geodetic datum, 221 campaign-mode GPS observations from 2002 to 2010 were employed to establish a surface horizontal velocity model for SW Taiwan by using the Kriging method and block model. According to exterior checking of the 44 continuous GPS stations in the study area, the averages of the residual errors (RMS exterior ) obtained using the Kriging method for the north and east components are ±1.9 and ±2.2 mm/year, respectively, and those obtained using the block model are ±2.0 and ±2.9 mm/year, respectively. Similar exterior checking results regarding the Kriging interpolation method and block model in SW Taiwan indicated a high station density and a relatively satisfactory station spatial coverage in the study area. However, the block model is more favorable than the Kriging method while lack of observation sites due to the consideration of characteristics of the geological structure. Traditional coordinate transformation was used to compare with the coordinate precision of semi-dynamic datum. Comparing the Kriging method and block model revealed that the RMS of the residuals of the coordinate transformation rapidly increased from ±7.8 in 2005 to ±19.1 in 2009 (Fig. 8). Therefore, coordinate transformation is not an appropriate method for eliminating distortion over a long time in a deformed area.
In this study, only surface horizontal motions were considered to establish a horizontal velocity model for the semi-dynamic datum. The effects of coseismic and postseismic deformation were not considered. Therefore, in the future, crustal deformation models should include the surface vertical velocity field, coseismic displacements, and postseismic deformation caused by major earthquakes.