Realization approach of non-linear postseismic deformation model for Taiwan semi-kinematic reference frame

Surface displacements associated with earthquake-cycle deformation of active faults significantly influence the accuracy of geodetic datum and reference frames, especially for the Taiwan area with high plate convergence and deformation rates. Following the current architecture of the semi-kinematic reference frame in Taiwan, which does not particularly consider the non-linear behavior of postseismic deformation, we explored the methods to implement a non-linear postseismic deformation model using the 2003 Chengkung earthquake as an example. Together with linear interseismic function, we utilized a logarithmic function to approximate the non-linear postseismic decays. For the time series without preseismic observations, we removed the interseismic velocities by spatial interpolation and fitted the resultant time series with the logarithmic function. After estimating postseismic decays for all GPS stations, we conducted two grid models of accumulative displacements for only postseismic deformation and total deformation after the earthquake. The first grid model provides a useful prediction for tracking surface movements, and the latter model provides a straightforward view to access the timing and amount of displacements to correct the semi-kinematic reference frame. The implementation of the grid models can well approximate non-linear postseismic trends for the semi-kinematic reference frame.


Introduction
The semi-kinematic reference frames implemented with surface deformation models are one practical approach to maintain the accuracy of reference frames distorted by the earthquake-cycle deformation. (Blick et al. 2003;Tanaka et al. 2007;Pearson et al. 2010;Pearson and Snay 2012;Blick and Donnelly 2016;Crook et al. 2016;Klein et al. 2019). Grid models of surface deformation used for the semi-kinematic reference frame usually contain interseismic velocities (e.g., McCaffrey 2005; Meade and Hager 2005;Ching and Chen 2015;Li et al. 2019) and coseismic displacements (e.g., Pearson and Snay 2012;Tanaka et al. 2007;Li et al. 2019). However, different from the signals of coseismic and interseismic deformation, the postseismic deformation shows transient, non-linear decay patterns (e.g., Shen et al. 1994;Savage and Svarc 1997;Freed et al. 2006). Afterslip and postseismic relaxation are two major mechanisms of postseismic deformation, which can be represented by logarithmic (Marone et al. 1991;Scholz 1998) and exponential (Nur and Mavko 1974;Pollitz 1997) functions, respectively. The two mechanisms may exist simultaneously and generate more complex decay patterns after large earthquakes (Hetland and Hager 2006;Tobita 2016). For most crustal faults, it is suitable to use one decay function to describe postseismic deformation (e.g., Hearn 2003;Savage et al. 2005;Freed et al. 2006). In terms of practically maintaining the reference frames, however, the methods to describe the postseismic decay pattern are nonunique. Many studies used the piecewise linear functions to approach the nonlinear postseismic pattern (Altamimi et al. 2007(Altamimi et al. , 2011Blick and Donnelly 2016;Klein et al. 2019). Instead, some other studies used a ramp function (Crook et al. 2016) or proposed velocity varying with time after the earthquake (Li et al. 2019) to fit the non-linear trend. These linear functions can more or less correct the postseismic signals, but the fitting of these functions needs to wait until the non-linear decay is ceased or the complete postseismic deformation is finished. Therefore, a nonlinear (logarithmic or exponential) function can describe the postseismic pattern appropriately (Bevis and Brown 2014) and this concept was used in the latest version of the International Terrestrial Reference Frame, ITRF2014 (Altamimi et al. 2016).
The island of Taiwan is located at a plate boundary zone with rapid deformation rates, causing the urgent need to maintain the reference frame. Since 1979, successive static reference frames were announced such as TWD67 (based on the GRS67 reference ellipsoid), TWD97 (based on the GRS80 reference ellipsoid with respect to ITRF94 at the epoch of 1997.0) (Yang et al. 2001), andTWD97[2010] (refer to ITRF94 at the epoch of 2010.0). Because the current TWD97 [2010] datum was announced a decade ago with some major earthquakes occurred since then, a preliminary semi-kinematic reference frame was recently proposed and tested (Ching and Chen 2015;Li et al. 2019). Yet, these prior examinations did not particularly consider the non-linear postseismic decay. In addition, since the installation of continuously recorded GPS stations in Taiwan was conducted by different government agencies mainly including Academia Sinica, Ministry of Interior, Central Geological Surveys and Central Weather Bureau, the installation time of the GPS stations was diverse in a way that the stations were gradually installed after major earthquakes and the distribution of the most stations is not designed for specific faults (i.e., denser close to the faults and less dense away from the faults). Therefore, although current ~ 500 GPS stations in Taiwan provide a plausibly enough spatial density, the spatial and temporal resolutions of postseismic deformation around the seismogenic faults were generally insufficient.
Thus, to fully estimate the effects of earthquake-cycle deformation for accurately maintaining the reference frame in Taiwan, in the study we propose methods to characterize the non-linear postseismic deformation under the limitation of sparse spatiotemporal coverage. We used the 2003 Chengkung earthquake in eastern Taiwan as an example to examine the efficiency of nonlinear function fitting and to explore the ways to utilize GPS data for correcting postseismic deformation without enough spatial and temporal resolution. With the consideration of the non-linear postseismic decay, the deformation models can more accurately characterize earthquake-cycle deformation and therefore assist the maintenance of semi-kinematic reference frames. Following the work of Ching and Chen (2015) and Li et al. (2019), this study focuses on the horizontal component only, which has a strong influence in cadastral surveys, and neglects the vertical component at the current stage of examination.

Tectonic background and the 2003 Chengkung earthquake
As an active arc-continental collision between the Eurasian and Philippine Sea plates (Suppe 1981(Suppe , 1984Teng 1990;Byrne et al. 2011) with a rapid plate convergence rate of 82 mm/year (Yu et al. 1997), the island of Taiwan is consequently having high deformation rates across active structures (Shyu et al. 2016) and high seismicity (Wu et al. 2008). The longitudinal Valley in eastern Taiwan is considered as the main suture zone between the Eurasian and Philippine Sea plates (Shyu et al. 2005). Two major active faults are present in the Longitudinal Valley: the Longitudinal Valley fault and the Central Range fault (Fig. 1). The Central Range fault is the west-dipping fault at the eastern side of the Central Range and became buried in the northern valley (Shyu et al. 2006), probably served as a retro-wedge fault (Chuang et al. 2014). The Longitudinal Valley fault located at the boundary between the western edge of the Coastal Range and the eastern edge of the valley has the highest geologic and geodetic fault slip rate in Taiwan (Chang et al. 2016;Ching et al. 2011;Shyu et al. 2016) and has partially ruptured during the 1951 Taitung-Hualien earthquake sequence (Shyu et al. 2007). The Longitudinal Valley fault has several segments (Shyu et al. 2005) and is partitioned into two strands to the southernmost end (Shyu et al. 2008;Chuang et al. 2012). In addition, the Longitudinal Valley fault in the southern half of the valley has significant fault creep (Lee et al. 2001(Lee et al. , 2003Hsu and Bürgmann 2006;Champenois et al. 2012).
The dense GPS network in Taiwan provides the information on surface deformation associated with the earthquake cycle (Yu et al. 1997;Hsu et al. 2009b;Lin et al. 2010;Ching et al. 2011). The M w 6.8 Chengkung earthquake occurred in December, 2003 is the largest onshore earthquake after the 1999 M w 7.6 Chi-Chi earthquake with a discernible postseismic deformation Cheng et al. 2009). The earthquake occurred on the listric Chihshang reverse fault (Kuochen et al. 2007), which is one major segment of the Longitudinal Valley fault in the southern valley (Shyu et al. 2005). The coseismic fault slip is mainly concentrated at the depth between 5 and 20 km (Wu et al. 2006a;Ching et al. 2007;Hsu et al. 2009a). No shallow fault slip of the Chengkung earthquake and significant fault creep during the interseismic period suggest the shallow portion of the Chihshang fault was less locked and had little strain accumulation. On the other hand, the deeper portion of the Chihshan fault was locked and capable for generating earthquakes. The postseismic deformation is mainly concentrated at the shallow depth (Cheng et al. 2009;Hsu et al. 2009a) and dominated by afterslip on the Chihshang fault (Perfettini and Avouac 2004;Savage and Yu 2007;Hsu et al. 2009a).

Postseismic deformation model of the semi-kinematic reference frame
To better characterize postseismic deformation for the maintenance of the Taiwan semi-kinematic reference frame, we examine the implementation of postseismic deformation models based on the concept of earthquakecycle deformation. Particularly, we propose postseismic deformation models in the following three aspects: (1) when the GPS time series are sufficient to differentiate earthquake-cycle deformation (i.e., having interseismic, coseismic and postseismic signals), we examine the efficiency to generate a deformation model with a nonlinear fitting function; (2) when GPS time series lacks pre-earthquake observations, the interseismic velocities are interpolated from surrounding stations to provide a basis for fitting non-linear postseismic decay; and (3) after the pattern of the postseismic decay at each station is fitted, we use a grid model to regionally estimate the amount of postseismic correction for the semi-kinematic reference frame. In the following sections, we presented the methods and results of these three aspects.

GPS time series data used in this study
We used GPS time-series data from the GPS Lab of the Taiwan Earthquake Center (http://gps.earth .sinic a.edu.tw), which provides GPS daily solutions under the ITRF2008 (Altamimi et al. 2011). We selected 14 stations in the source are of the 2003 Chengkung earthquake for this study (Fig. 1). Eleven of 14 stations have observations of at least 2 years before and after the Chengkung earthquake. Six of 11 stations (KNKO, JPIN, CHEN, TUNH, TAPO and S104) are located in the hanging wall area (east of the Chihshang fault) and five stations (JPEI, YULI, TAPE, SHAN, S105) are in the footwall area (west of the Chihshang fault) (Fig. 1). The rest three stations located in the hanging wall (T102 and T103) and footwall area (T105) only have observations of less than 2 years with no preseismic data ( Fig. 1). Because the M w 6.1 Taitung earthquake occurred on April 1, 2006 (e.g., Wu et al. 2006b) around this region, which produced coseismic and postseismic signals at surrounding GPS stations, in this study we examined the postseismic deformation in the period between the 2003 Chengkung and 2006 Taitung earthquakes. Table 1 shows the information and observation periods of all the stations.

Non-linear postseismic fitting
Since we aimed to characterize the non-linear postseismic pattern, one fundamental basis is to isolate the trend of the different earthquake-cycle period. The most common approach is a mathematical fitting of GPS time series, also called trajectory models, with linear, step and non-linear functions for the interseismic, coseismic and postseismic periods, respectively (e.g., Nikolaidis 2002;Bevis and Brown 2014). Because afterslip is the major mechanism of the postseismic deformation of the Chengkung earthquake (Perfettini and Avouac 2004;Savage and Yu 2007;Hsu et al. 2009a), we chose a logarithmic function to describe the non-linear postseismic pattern  (Marone et al. 1991;Scholz 1998). The fitting equation of the GPS time series is as below.
where a is the initial coordinate and b is the steady interseismic velocity. c and d are the amplitude of annual variations, and e and f are the amplitude of semi-annual variations. g is the coefficient of step function H for permanent displacement, such as the antenna changes or coseismic displacements. h is the coefficient of the step function H for velocity change after the earthquake. k is the coefficient of the logarithmic decay and p is the characteristic time. ε is the residual between the observations y and fitting results. Because the characteristic time p could be different for each time series, we used a grid search method to select the optimal value of p for every time series. All the fitting time-series are available in Additional file 1.
To examine the efficiency and stability of the non-linear postseismic fitting with varied data lengths, we tested the fitting method with different postseismic data periods for every 30 days. For most crustal earthquakes without complex postseismic pattern, postseismic deformation (1) can be characterized and fitted with a non-linear function. Conversely, the fitted non-linear function should be able to describe the non-linear postseismic pattern. Therefore, when the postseismic signals are well fitted, the fitted non-linear curve can predict further postseismic pattern if no significant coseismic jump occurs due to a large aftershock. Following this idea, we fitted GPS time series with postseismic observations every 30 days to examine monthly fitting results, and we then computed the root mean square (RMS) values between all fitting results and GPS data. Figure 2 shows two examples of 8 time-series fitting (30, 60, 90, 180, 360, 540, 720, and 840 days) (Fig. 3). This is probably because the fitting functions of semi-annual and annual variations have a strong influence in such a short time series.
In the current realization of the semi-kinematic reference frame in Taiwan, Li et al. (2019) fitted GPS time series with different linear velocities before and after an earthquake following the philosophy of piecewise linear fitting. Therefore, Li et al. (2019) practically incorporated and our non-linear fitting, we used the first 60-day GPS results after the 2003 Chengkung earthquake, which is the time when our non-linear fitting starts to be effective as described above. We selected 5 GPS stations (CHEN, TUNH, TAPO, TAPE, and SHAN) closest to the epicenter, which have larger postseismic deformation.
We computed the misfit for our non-linear fitting result and stepwise linear fitting result following the method of Li et al. (2019). The accuracy is 0.66 mm and 1.73 mm in the east and north components in our method, respectively. However, the accuracy is 4.12 mm and 5.58 mm in the east and north components using the piecewise linear fitting method, respectively. As the results, it shows that the improvement by our method is 70-80% than the method of Li et al. (2019).

Postseismic fitting for incomplete time series
The GPS installation may be conducted in response to postseismic observations after major earthquakes (e.g., Shen et al. 1994;Freed et al. 2006). However, it is more difficult to analyze the postseismic deformation from time series without preseismic observations. Although these stations do not have data before the earthquakes, they provide valuable observations of postseismic deformation that should be used to constrained reference frames. Nevertheless, direct fitting of only postseismic observations with the logarithmic and linear functions simultaneously might result in unrealistic results since the steady interseismic velocity cannot be well constrained by postseismic observations. To better fit the logarithmic function for the time series without preseismic observations, we had to subtract interseismic velocities from the postseismic time series. Because the interseismic velocities result from steady fault creeping at depth, the pattern should be spatially correlated. For example, Yu and Kuo (2001) estimated interseismic velocities in eastern Taiwan, which is the area of this study, and showed that neighboring stations generally have velocities with similar movement directions if they are located at the same side of the fault (i.e., stations on the hanging wall are moving toward northwest while stations on the hanging wall have little motion with respect to the Central Range). Therefore, we interpolated the interseismic velocities for the stations installed after the earthquake from the velocities at the neighboring stations.
We first tested this method for the station with a complete time series to validate the method. We chose the station SHAN, which is located at the center of the study area in the footwall area, and estimated the interseismic velocity according to the surrounding stations also in the footwall area (JPEI, S105, TAPE, and YULI) using the linear interpolation. The interpolated velocities of 8.51 and − 6.34 mm/year are consistent with the velocities of 9.45 ± 0.34 and − 6.09 ± 0.32 mm/ year derived from the time series fitting for the north and east components, respectively (Table 2). Thus, we were able to apply this method to the stations with incomplete time series.
We estimated interpolated velocities for the stations installed after the Chengkung earthquake and fitted the time series. We selected three stations (T102 and T103 on the hanging wall and T105 on the footwall) with apparent non-linear postseismic signals to examine the method for fitting incomplete time series. Since these three stations are in different sides of the fault, we interpolated velocities for the hanging wall and footwall, separately. We used the velocities derived from the fitting of the complete time series (Table 2) for the spatial interpolation. The interpolated velocities are 25.11 and − 18.30 mm/year for T102, 28.44 and − 19.77 mm/year for T103, and 9.79 and − 7.51 mm/year for T105 in the north and east directions, respectively (Table 2). Once the interseismic velocities were generated, we subtracted the velocities from the postseismic observations. As a result, the observations should only have the non-linear postseismic component. We then fitted the time series by the following equation.
Here, we only fitted the time series with the initial coordinate a and the logarithmic pattern. We ignored the seasonal variation because we cannot correctly estimate the (2) y(t i ) = a + n k 1 k j · log 1 + t i − t eq /p + ε annual and semi-annual variations by spatial interpolation. Figure 4 shows the results of the logarithmic fitting for the three stations and the postseismic decay patterns can be characterized. All the results of time-series fitting for the stations with incomplete data are available in Additional file 2.

Postseismic grid model
As we estimated the logarithmic decay functions and interseismic velocities at all stations, we generated a postseismic deformation grid model for the maintenance of semi-kinematic reference frames. According to the GPS time series fitting, the surface deformation after the earthquake should contain at least interseismic velocities and postseismic decays. After properly estimating the linear interseismic velocity and non-linear postseismic decay at each GPS station, we are then able to predict the amount of surface displacements at a given time duration after the earthquake. Since we have the linear velocities derived from either time series fitting or spatial interpolation for interseismic deformation and logarithmic fitting curves for postseismic deformation, we can estimate the accumulative displacements for postseismic deformation and post-earthquake (postseismic + interseismic) deformation (data available in Additional file 3). Moreover, following the method of Ching and Chen (2015) and Li et al. (2019), we used the kriging interpolation to generate a grid model of postseismic deformation and postearthquake deformation within a certain time period using the estimated interseismic and postseismic displacements at every continuously recorded GPS stations. Figure 5 shows the grid models of postseismic deformation and post-earthquake deformation with accumulative displacements in the period between the Chengkung and Taitung earthquakes. We conducted the grid models for the hanging wall and footwall areas separately and integrated the results together to generate regional grid models as shown in Fig. 5.
To validate the grid model, we used independent GPS observations within the study area. The grid model can be validated using other geodetic measurements or surface displacements derived from imaging geodesy. After searching for geodetic data, we found that two datasets of continuous and campaign GPS could be used for the validation. There are three continuous GPS stations in the area, which were not used to construct the grid model, including ERPN, LONT, and PING. These three stations were installed only a year before the Chengkung earthquake so we did not use their time series for our tests. Because the first two stations are located in the southernmost Longitudinal Valley within the fault zone of the southernmost Longitudinal Valley fault (Chuang et al. 2014), their deformation behaviors are complicated. Therefore, we chose to use the observations of PING, which is at least 10 km away from surrounding GPS stations, to validate our grid model (Fig. 1). We computed the total displacement at the station based on our grid model (Fig. 5b) and compared the interpolated displacement to the actual displacement derived from the time series. The result shows that the model displacements are − 20.64 mm and 97.48 mm while the actual displacements are − 17.89 mm and 77.20 mm in the east and north directions, respectively, for the time period of ~ 2 years after the Chengkung earthquake. The difference is well below the accuracy requirement of cadastral surveys of 60 mm. In addition, we used campaign GPS data of postseismic deformation , which was surveyed ~ 4 months after the earthquake, to validate our grid model. We computed a grid model of total displacements for 4 months and compared its result to the campaign data. The accuracy is 11.3 mm and 16.1 mm in the east and north directions, respectively. Through these two accuracy assessments, the accuracy of our grid models is generally below 20 mm, which fulfills the need of cadastral surveys. Further validation may be required using other geodetic observations or results derived from synthetic aperture radar (SAR) images (e.g., Champenois et al. 2012).  Fig. 1. Figure 5a provides amount of postseismic displacements, which can be used to calculate how much postseismic deformation occurred. Figure 5b shows amount of total displacements after the earthquake, which can be used to justify if the total displacements exceed the tolerance of the reference frame

Usage of postseismic grid model to the semi-kinematic reference frame
We proposed the grid models for the postseismic deformation to correct the semi-kinematic reference frame in Taiwan. We created the grid model with the latitude and longitude spacing of 0.01°, which is similar to the baseline of fourth-order control points (original triangulation benchmarks) and cadastral benchmarks in Taiwan. Once the logarithmic patterns are estimated, one can easily estimate accumulative displacements at a specific time and generate the grid model by kriging interpolation. For practical usage, we demonstrated two types of grid models. One has postseismic displacements only (Fig. 5a), which provides the estimates of the excess postseismic displacements in addition to steady interseismic deformation. The other model has the accumulative displacements including interseismic and postseismic movements, which provide an immediate view to evaluate the total displacements after an earthquake.

Discussion
In this study, we explored approaches to accommodate the non-linear pattern of postseismic time series and to generate a deformation model for the maintenance of the semi-kinematic reference frame in Taiwan. The transient characteristics of postseismic deformation not only reflect the dynamic properties of seismogenic faults but also complicate the process to approximate GPS time series.
Following the idea of earthquake-cycle deformation and the convention of GPS time series fitting, we implemented the logarithmic function to fit the postseismic time series. Compared with the method of piecewise linear fitting, our method can fit the non-linear pattern with one single function. Since the 2003 Chengkung earthquake is dominated by afterslip, which is mostly the prominent mechanism for many seismic events especially for crustal faults, the logarithmic function can fit and predict the non-linear pattern well with early observations after the earthquake (Fig. 2). In addition, as no obvious poroelastic deformation of the earthquake, which is another postseismic deformation mechanism generally occurred shortly and near the fault (e.g., Peltzer et al. 1998), we do not need to exclude early observations (Bevis and Brown 2014) and the logarithmic function can work properly. Consequently, the logarithmic function can serve as a useful approximation for the postseismic deformation of the 2003 Chengkung earthquake. Considering the tectonic environment of a fault-and-thrust belt and active orogen, the island of Taiwan is subject to earthquake-related deformation on crustal faults (Shyu et al. 2005(Shyu et al. , 2016 and presumably it is reasonable to simply use the logarithmic function to accommodate non-linear postseismic deformation for earthquakes in the Taiwan area. In terms of the realization of semikinematic reference frames, such a logarithmic approximation can well predict the non-linear trend, providing useful information about the amount of postseismic displacements at a certain time period and helping to estimate the timing to correct the reference frame.
The method of GPS time series fitting with the functions analogous to earthquake-cycle deformation can provide a stable approximation of crustal movement even with incomplete time series. The idea of the time series fitting simply decomposes the observations into interseismic, coseismic and postseismic deformation approximated by linear, step and logarithmic functions, respectively. Although this method is technically a regression for time series data, it follows explainable geophysical mechanisms of fault behaviors. Together with the logarithmic function, the linear function for interseismic deformation can help stably fit postseismic deformation and provide simple velocity reference to determine whether the postseismic deformation is ceased. On the other hand, however, if GPS time series is incomplete, more specifically without observations before the earthquake, the interseismic linear function is not constrained, making the logarithmic fitting unstable. As the installation of continuously recorded GPS stations is conducted by different agencies resulting in nonuniform installation time, the observation periods of the GPS time series in Taiwan are commonly varied and often generate incomplete time series not to capture enough observations associated with interseismic, coseismic and postseismic deformation. Therefore, to better use data at continuous GPS stations with incomplete time series, we estimated interseismic velocities at the locations of those stations by spatial interpolation from surrounding stations with preseismic observations. With the constrained interseismic velocities, postseismic movements at the sites without preseismic observations can easily be approximated.
After the postseismic decays were properly estimated at all GPS sites, we proposed the postseismic deformation grid models to correct the reference frame. Following the ongoing architecture of the semi-kinematic reference frame of Taiwan, which proposed deformation models to correct interseismic deformation (Ching and Chen 2015) and coseismic displacements (Li et al. 2019), we constructed two grid models of accumulative postseismic displacements with given postseismic time period (Fig. 5) using the estimated logarithmic curves at the GPS sites. The two grid models estimate only postseismic deformation and total deformation after the earthquake. The first model with only postseismic deformation can provide clear sense of the amount of accumulative postseismic displacements, which is easier to locate the area of postseismic deformation and track the postseismic decay. On the other hand, because afterslip might generate a similar displacement pattern with coseismic deformation, the total displacements may be canceled if the postseismic displacement direction is opposite to the interseismic displacements especially in the footwall area in this study (Fig. 5b), which consequently will extend the timing to correct the reference frame in such area. Therefore, the latter model with total post-earthquake (postseismic + interseismic) deformation is straightforward for users to determine whether they have to make the correction after an earthquake. Thus, the two proposed grid models can provide useful information to track postseismic deformation and to determine the timing to correct the reference frame in different ways.
In terms of maintaining the reference frame accuracy, it is crucial to characterize time-dependent information of surface deformation. If the information is well constrained, the time-dependent displacements can be used to update official coordinates and to access future deformation patterns. In previous work for examining the realization of the semi-kinematic reference frame in Taiwan, Ching and Chen (2015) explored the usage of interpolation and block models to characterize interseismic deformation, and Li et al. (2019) established a framework to integrate interpolation models for interseismic and coseismic deformation. In this study, we examined the approaches to use complete and incomplete time series to characterize non-linear postseismic deformation for the maintenance of the semi-kinematic reference frame. Once we estimate the amount of postseismic displacements as described above, we can plan and control the timing for updating the reference frame. By incorporating the fitting of non-linear postseismic deformation, we considered all the components of a full earthquake cycle (interseismic, coseismic and postseismic deformation), which could appropriately describe time-dependent displacements. The current examination of the Taiwan semi-kinematic reference frame including this and previous studies (Ching and Chen 2015;Li et al. 2019) focuses on the horizontal component. The realization of the full 3D reference frame will be one future task, and the validation and comparison of the reference frame using other geodetic data are also important for future work.

Conclusion
We explored the methods to implement a non-linear postseismic deformation model for the Taiwan semikinematic reference frame using the 2003 Chengkung earthquake as an example. We adopted a logarithmic function in response to the afterslip-dominant postseismic deformation of the Chengkung earthquake to approximate the non-linear postseismic decays at continuous GPS sites. Following the concept of earthquake-cycle deformation, we fitted the GPS observations with linear, step, and logarithmic functions for the complete time series. For the time series without preseismic observations, we subtracted the interseismic velocities by spatial interpolation from the incomplete time series and fitted the resultant time series with the logarithmic function. After properly estimating postseismic decays at all GPS sites, we conducted two grid models for accumulative postseismic displacements. The first grid model can clearly track postseismic movements, and the latter model can access the total amount of surface displacements. The implementation of the grid models can provide useful prediction for tracking surface movements to access the timing and amount of displacements to correct the semi-kinematic reference frame (Additional files 1, 2, 3).