Spatial–temporal properties of afterslip associated with the 2015 Mw 8.3 Illapel earthquake, Chile

In order to characterize the spatial–temporal properties of postseismic slip motions associated with the 2015 Illapel earthquake, the daily position time series of 13 GNSS sites situated at the near-field region are utilized. Firstly, a scheme of postseismic signal extraction and modeling is introduced, which can effectively extract the postseismic signal with consideration of background tectonic movement. Based on the extracted postseismic signal, the spatial–temporal distribution of afterslip is inverted under the layered medium model. Compared with coseismic slip distribution, the afterslip is extended to both deep and two sides, and two peak slip patches are formed on the north and south sides. The afterslip is mainly cumulated at the depth of 10–50 km, and the maximum slip reaches 1.46 m, which is situated at latitude of − 30.50°, longitude of − 71.78°, and depth of 18.94 m. Moreover, the postseismic slip during the time period of 0–30 days after this earthquake is the largest, and the maximum of fault slip and corresponding slip rate reaches 0.62 m and 20.6 mm/day. Whereas, the maximum of fault slip rate during the time period of 180–365 days is only around 1 mm/day. The spatial–temporal evolution of postseismic slip motions suggests that large postseismic slip mainly occurs in the early stage after this earthquake, and the fault tend to be stable as time goes on. Meanwhile, the Coulomb stress change demonstrate that the postseismic slip motions after the Illapel earthquake may be triggered by the stress increase in the deep region induced by coseismic rupture.


Introduction
In general, the postseismic deformations with logarithmic or exponential decay can be observed after many large tectonic earthquakes due to the fault instability and coseismic stress transfer. Numerous studies indicate that the main factors of postseismic deformations contain afterslip, viscoelastic relaxation, and poroelastic rebound (Segall and Davis 1997;Hsu et al. 2002Hsu et al. , 2006Politz et al. 2006;Freed 2007;Bruhat et al. 2011;Hoecher et al. 2012). Among them, the effects of afterslip and viscoelastic relaxation are relatively significant in comparison with poroelastic rebound. Meanwhile, the effects of afterslip are cumulated at the near-field region of the source in the early stage after the earthquake occurrence, while the effects of viscoelastic relaxation are cumulated at the middle-and far-field regions in the long-term period. Characterizing the properties of postseismic slip is helpful for understanding the stability of fault and the rheological properties of earth's medium, which is of great significance for the postseismic disaster risk assessment (Kreemer et al. 2006;Diao et al. 2011Diao et al. , 2018Ozawa et al. 2012;Jiang et al. 2017Jiang et al. , 2018. Compared with other space geodesy technology, GNSS technology has high positioning accuracy (the accuracy of daily solution can reach millimeters) and high time resolution (Shen et al. 2009;Guo et al. 2012;Simons et al. 2011;Diao et al. 2012). Therefore, GNSS has incomparable advantages Open Access *Correspondence: yfxiang181@163.com for exploring the postseismic deformations during the period of months, 1 year, or several years. The continuous GNSS network can accurately detect the postseismic signal, which can provide valuable data sets for exploring the spatial-temporal properties of afterslip after some large earthquakes.
According to United States Geological Survey (USGS), a Mw 8.3 earthquake struck the coastal area of central Chile on 16 September 2015. Numerous studies associated with coseismic and postseismic deformations of this earthquake have been carried out. Chen et al. (2016) utilized GPS/GLO-NASS to determine the real-time coseismic displacements of the Mw 8.3 Illapel earthquake, and revealed that GPS/ GLONASS provided more accurate and robust coseismic displacements in comparison with GPS-only. Meanwhile, the coseismic slip distributions of the Illapel earthquake determined by the GPS/GLONASS observations tended to be shallower and larger. Okuwaki et al. (2016) proposed a seismic source of the Illapel earthquake, and characterized the rupture process with a novel kinematic waveform inversion method. The spatiotemporal properties of coseismic slip revealed a complex rupture propagation mode: two up-dip rupture propagation events, and the second rupture event may be triggered by strong high-frequency radiation events. Heidarzadeh et al. (2016) reconstructed the source model of the Illapel earthquake using teleseismic and tsunami data, and tried to explore the relationship between 2010 Maule and 2015 Illapel earthquakes. Williamson et al. (2017) characterized the coseismic slip model of this earthquake by the combined use of GPS, InSAR, and tsunami waveform data. The coseismic slip showed that most coseismic slip was cumulated at a region immediately offshore, and there was almost no obvious slip near the trench. Melgar et al. (2016) adopted geodetic observations to resolve a kinematic slip model, and found the shallow and deep sections of the megathrust were segmented and had fundamentally different behavior. Ruiz et al. (2016) estimated the coseismic rupture process of the Illapel earthquake based on teleseismic P-waves and strong-motion records, and revealed most of aftershocks focused on the deeper section of the plate interface. Barnhart et al. (2016) characterized the coseismic slip and early afterslip of the Illapel earthquake using GPS and InSAR observations. The result indicated that early afterslip after this earthquake appeared in two regions around the coseismic slip with partial overlap. Huang et al. (2017) explored the spatiotemporal properties of early afterslip and aftershocks using the postseismic GPS and InSAR data. The spatial distribution of aftershocks was consistent with that of afterslip inverted by geodetic observations, which seemed to surround the coseismic slip.
The above studies associated with the 2015 Illapel earthquake have revealed the coseismic slip, coseismic rupture process, early afterslip, and aftershocks based on multi-source observations. However, few studies focus on exploring the spatial-temporal properties of afterslip after this earthquake, as well as the Coulomb stress change. Characterizing the spatial-temporal properties of postseismic slip behavior is of great significance to the assessment of the fault stability and postseismic disaster risk. In addition, considering the influence of the Mw 7.3 aftershock on 25 December 2016, the time span of GNSS position time series after this earthquake is selected to be 1 year. Thus, the GNSS observations of 13 sites during the 1 year after this earthquake are utilized to explore the spatial-temporal properties of postseismic fault slip behavior in this paper. Meanwhile, we also discuss the Coulomb stress change after this earthquake, which will be helpful for understanding the triggering mechanism of postseismic slip motions.

Regional tectonic background and GNSS data
The 2015 Illapel earthquake occurs at the region where the Nazca plate collides with the South American plate, which is the seismic active area in Chile (Shrivastava et al. 2016;Barnhart et al. 2016). The epicenter is located at the latitude and longitude of 31.57 o S and 71.67 o W, and the focal depth is 22.4 km. The Nazca plate subducts to the South American plate at a speed of approximately 74 mm/year along the northeast, and collides beneath the Peru-Chile Trench. Meanwhile, the rupture zone of this earthquake is laterally constrained by two prominent structural features: the Challenger fault zone (CFZ) in the north and the Juan Fernandez ridge (JFR) in the south. Historically, many large thrust-type earthquakes occur in the region, including the 1960 M 9.5 earthquake, the 2010 Mw 8.8 earthquake, and several M 7 + earthquakes.
In order to characterize the postseismic deformation mechanism and spatial-temporal properties of afterslip, the daily position time series of 13 GNSS sites situated at the near-field region are utilized (Fig. 1). In general, the postseismic GNSS position time series also contains plate tectonic movement signal. Thus, the GNSS observations of 13 sites with time span of about 2 years before the earthquake are chosen to estimate the plate tectonic signal. The GNSS position time series of 13 sites are provided by Nevada Geodetic Laboratory (NGL), which are available at http://geode sy.unr.edu/. The data processing is conducted by GIPSY-OASIS-II. NGL data processing center utilizes ionospheric independent combination to remove ionospheric effects, Global Mapping Function (GMF) model to correct tropospheric delays, FES 2004 to correct tidal loading effects, and antenna phase center correction parameters refer to the absolute antenna phase calibration file provided by IGS center.

Postseismic signal extraction scheme
In general, the postseismic GNSS position time series mainly contains long-term trend signal caused by plate tectonic movement, seasonal signal induced by surface mass loading (i.e., atmospheric pressure, non-tidal oceanic change, and hydrological loading), postseismic signal, and noise (e.g., white and colored noise) (Zhou et al. 2020). Among them, the postseismic signal can usually be expressed in logarithmic or exponential form. In order to extract the postseismic signal, we should model and remove the effects of other signals. The preseismic GNSS position time series with time span of about 2 years is utilized to estimate the longterm plate movement rate. Based on the above discussion, we design a scheme for extracting and modeling postseismic signal with considering the influence of plate tectonic motions, surface mass loading, and clear noise. The scheme of postseismic signal extraction mainly contains the following steps: (1) Outlier detection and elimination Due to the influence of abnormal observation environment and other factors, there are usually a small number of outliers in the GNSS observations, which greatly deviate from normal observations. For outlier detection and elimination, we first roughly eliminate the GNSS observations with errors > 20, 20, and 40 mm in the north, east, and up components. After that, the outliers in the GNSS position time series are detected and eliminated by the combined use of inter-quartile range (IQR) and 3 α rules.
(2) Velocity interpolation Due to the influence of background tectonic movement, the horizontal components of GNSS position time series generally have significant linear trend term, which can reach several centimeters per year. Therefore, the linear trend term in GNSS position time series caused by background tectonic movement should be considered for the postseismic signal extraction. Whereas, some key sites used to study the postseismic deformations may have missed preseismic GNSS observations. Considering this problem, this paper mainly utilizes the velocity interpolation to estimate the linear trend term caused by the background tectonic movement for these sites, and the spatial interpolation is utilized to estimate the threedimensional crustal movement rate based on the GNSS velocity field around these sites. For the velocity interpolation, the velocity interpolation for stain rate (VISR) program developed by Shen et al. (2015) is adopted.

(3) Determination of relax time constant
Before parameter estimation, the relaxation time constant must be determined and the non-linear function model should be converted into a linear model, and then the parameter estimation can be accomplished. For the determination of relaxation time constant, the Nelder-Mead simplex method is firstly utilized to roughly estimate the relaxation time constant of each site by means of seeking a local optimal solution. According to the roughly estimated relaxation time constant, the approximate range of relaxation time constant is determined. Then, the trial and error method is utilized to search for the optimal relaxation time constant.

(4) Parameter estimation
Considering the influence of background tectonic movement, it is not reasonable to adopt only logarithmic or exponential function to simulate postseismic deformations (Tobita 2016). When the preseismic GNSS observations are added, the entire position time series can be expressed by the following formula: Among them, t is the epoch time in years, y 0 is the initial position of GNSS site, vt is the linear trend term caused by background plate tectonic movement, is the postseismic signal expressed as logarithmic or exponential form, v l is the crustal movement rate altered by earthquake, and ε is the systemic errors. For parameter estimation of GNSS position time series, least square parameter estimation is a general and robust method. However, the noise properties of GNSS position time series may have changed due to the influence of earthquake. Thus, maximum likelihood estimation (MLE) is utilized to simultaneously estimate the parameters and noise of GNSS position time series, and then extract and model the postseismic signal.

Fault slip inversion method
The dislocation theory is the theoretical basis of fault slip inversion, and Okada elastic dislocation model establishes the elastic response relationship between fault dislocation and crustal deformation. The function model can be expressed as: where x represents parameters related to fault (e.g., strike, dip, rake, length, and width), and ε denotes the errors. It can be seen from the above formula that the crustal deformation can be regarded as the function of earth model parameter, fault geometric parameter, and dislocation parameter. After determining the fault geometric parameter, the inversion problem of dislocation model can be transformed into linear inversion, which can be expressed as: where G is the green function, which denotes the elastic response of surface deformation to the unit fault slip, s (1) denotes the fault slip, and y denotes the surface deformation. The calculation of green function involves the selection of earth model, which can be divided into uniform half-infinite space dislocation model, layered half-infinite space dislocation model, spherical uniform dislocation model, and spherical layered dislocation model.
Okada elastic dislocation theory establishes the spatial response of surface deformation to fault dislocation under the uniform semi-infinite space dislocation model. This method has fast calculation speed and high accuracy, and is widely used in seismic research field. However, due to the inhomogeneity of earth medium, the inversion of fault slip under the uniform semi-infinite space dislocation model maybe not reliable. Wang et al. (2013) extended the theory to the layered half-infinite space dislocation model based on Okada elastic dislocation theory, and proposed the steepest descent method (SDM) for the inversion of fault slip under layered medium model. Thus, the inversion of fault slip can be expressed as: where s is the fault slip, y is the surface deformation, M is the Green function, τ (s) is the stress drop, and β is the smoothing factor. y − Ms 2 indicates the data misfit, and β 2 Hs indicates the model roughness. In general, it is suitable to utilize the trade-off value of data misfit and model roughness as the smoothing factor. In order to improve the efficiency and accuracy of fault slip inversion, the inversion process can be strongly constrained by setting relevant parameters: where s i and ϕ i are the slip and rake of sub-fault i . s max , ϕ min , and ϕ max denote the maximum slip, minimum and maximum rake, respectively.

Postseismic signal extraction
After the 2015 Illapel earthquake, the GNSS position time series of 13 sites located near the epicenter occur significant logarithmic signals. In order to characterize the spatial-temporal properties of afterslip, the scheme of signal extraction and modeling described in Sect. 2.2 is utilized to extract the postseismic signal. Taking the site CMBA as an example, this scheme is utilized to extract the postseismic signal. The time span of the GNSS position time series of site CMBA is 2013.0021-2015.7070, and there are two offsets caused by non-seismic factors. The postseismic signal of site CMBA in three components is attenuated in logarithmic form, and some studies have revealed that the short-term near-field postseismic deformation is mainly caused by the afterslip. Thus, the logarithmic function model is adopted to simulate the postseismic deformations. Firstly, the successive iterative processing of IQR and 3σ is performed until all the outliers in the GNSS position time series are removed. Then, the optimal relaxation time constant of this site can be identified as 8 days by the combined use of the Nelder-Mead simplex method and trial and error method. Finally, MLE is utilized to estimate the parameters of GNSS position time series, and then extract and model the postseismic signal. Figure Table 1 presents the coseismic offsets of the 13 GNSS sites and the cumulative deformations during the 1 year after this earthquake. The above results indicate that the postseismic signal extraction and modeling scheme introduced in this paper can effectively extract postseismic signal and estimate parameters of GNSS position time series with consideration of the influence of background tectonic movement. Meanwhile, this scheme performs parameter estimation of GNSS position time series through the overall modeling, which can effectively reduce the error introduced by the segmented processing.

Coseismic slip inversion
In order to characterize the postseismic fault slip, we should firstly determine the distribution of coseismic fault slip, which can be regarded as a reference to understand the spatial-temporal properties of postseismic slip. USGS conducts the preliminary inversion of coseismic slip using the teleseismic broadband waveform, and confirms that the nodal plane with the strike and dip of 6° and 19° fits well with waveform data. According to the relevant information provided by USGS, we construct a fault model with the length and width of 288 and 160 km along the strike and dip. The fault is divided into 10 × 10 km sub-faults, and the number of sub-faults is 464. Based on the theory of elastic dislocation, SDM is utilized to determine the coseismic slip under the layered medium model. The Green function of dislocation is generated by DGRN code developed by Wang et al. (2003),  Hreinsdottir et al. (2006) and Ding et al. (2015) confirmed that it is reasonable to use the difference of the average value of 4 daily coordinates of GNSS site before and after the earthquake to estimate the offset, which can effectively reduce the noise in GNSS observations and the effect of postseismic deformations. Based on this method, the offsets of 13 GNSS sites are calculated, and the coseismic deformations field is depicted as Fig. 7. It can be seen that large coseismic deformations are mainly cumulated at the near-field region of the epicenter, and the coseismic deformations decreases with the epicentral distance increasing. For the horizontal components, the largest coseismic deformation occurs at site EMAT, reaching According to the fault model and coseismic deformations, the coseismic slip can be characterized under the layered medium model. In order to retrieve the smooth and accurate fault slip model, the trade-off value of data misfit and model roughness is regarded as the optimal smoothing factor. Considering the weight of data misfit and model roughness, the optimal smoothing factor of coseismic slip inversion is identified as 0.1. As depicted as Fig. 7, the coseismic slip is mainly cumulated at the depth of 3-50 km, and the average value of coseismic slip is 1.35 m. The maximum coseismic slip is 7.60 m, which is located at latitude of − 31.13°, longitude of − 72.06°, and depth of 14.59 m. The moment released by this earthquake is equivalent to a moment magnitude of ~ 8.2, the maximum stress drop is 4.34 MPa, and the fitting coefficient of the inversion data and the model reaches 0.9998. The inversion results of the coseismic slip distribution of this earthquake are highly consistent with the results of Heidarzadeh et al. (2016), Okuwaki et al. (2016) and Williamson et al. (2017).

Postseismic slip inversion
Numerous studies have indicated that the effects of afterslip are cumulated at the near-field region of the source in the early stage after the earthquake occurrence, while the effects of viscoelastic relaxation are cumulated at the middle-and far-field regions in the long-term period (Perfettini and Avouac 2014;Feigl and Thatcher 2006;Hsu et al. 2006;Ozawa et al. 2011;Bedford et al. 2013). Meanwhile, Shrivastava et al. (2016) and Barnhart et al. (2016) also confirmed that the postseismic deformations of the 2015 Illapel earthquake in the near-field region is induced by afterslip, and the effect of viscoelastic relaxation can be neglected.
Thus, we attribute the postseismic deformations in the near-field region of this earthquake to afterslip, and then ignore the effects of the viscoelastic relaxation on the near-field region. The fault model similar to coseismic slip inversion is utilized to characterize the postseismic slip. Considering that afterslip is usually generated around coseismic slip region, the fault model is slightly extended to both sides and deep, thus constructing the length and width of 324 and 170 km along the strike and dip. Similarly, the fault model is divided into 10 × 10 km sub-faults, and the number of sub-faults is 544.
According to the method described in Sect. 2.2, the postseismic signals at 13 GNSS sites are extracted, and the cumulative deformations of 1 year after this earthquake are determined (see Fig. 8). Unlike the coseismic deformation field, large postseismic deformations occur at site PFRJ in the east component, which reaches − 240.8 mm. Whereas, the nearest site EMAT do not produce the largest postseismic deformation, suggesting that postseismic slip may present different spatial distribution in comparison with coseismic slip. Meanwhile, the postseismic deformations decay slowly with the increase of epicentral distance, which is different from the coseismic deformation field. The postseismic deformations of 13 sites very from − 240.8 to − 25.0 and − 42.8 to 19.7 mm in the east and north component, while it ranges from − 60.9 to 54.5 mm in the vertical component. In the vertical direction, most GNSS sites appear a downward trend, and the largest deformation reaches − 60.9 mm at site CMBA.
Based on the postseismic deformations and fault model, the postseismic slip can be determined under the layered medium model, and Fig. 8 descries the spatial distribution of afterslip. Similarly, the trade-off value of data misfit and model roughness is regarded as the optimal smoothing factor, and the smoothing factor of postseismic slip inversion is 0.08. Compared with coseismic slip, the afterslip is extended to both deep and two sides, and two peak slip patches are formed on the north and south sides. The afterslip is distributed around the region of coseismic slip and partially overlaps. For the two peak slip areas, the larger slip area is situated at the north side, while the smaller slip area is located near the source. The afterslip is mainly cumulated at the depth of 10-50 km, and the maximum afterslip reaches 1.46 m, which is situated at latitude of − 30.50°, longitude of − 71.78°, and depth of 18.94 m. The released moment of afterslip during the one year after this earthquake is equivalent to a moment magnitude of ~ 7.6, the maximum stress drop is 1.43 MPa, and the fitting coefficient of the inversion data and the model reaches 0.9929. Shrivastava et al. (2016) utilized coseismic deformations and the GNSS observations of 43 days after the 2015 Illapel earthquake to explore the coseismic and early afterslip, and confirmed that the afterslip was distributed around the region of coseismic slip, and forming two peak slip areas on both sides. Barnhart et al. (2016) adopted the data sets of InSAR and GNSS to characterize the coseismic and early (i.e., 38 days) postseismic fault slip for the 2015 Illapel earthquake, and demonstrated that two slip areas P1 and P2 are formed after this earthquake. In terms of spatial distribution, the afterslip characterized by GNSS observations of 1 year after this earthquake is consistent with the conclusions of the above studies. Whereas, the region of postseismic slip is slightly enlarged, and the magnitude of fault slip in each patch is increased accordingly.

Spatial-temporal evolution of afterslip
In order to characterize the spatial-temporal evolution of afterslip, the time span of 1 year after this earthquake is divided into 6 time periods (i.e., 0-30, 0-60, 0-120, 0-180, 0-270, and 0-365 days). Meanwhile, the constructed logarithmic or exponential model of postseismic signal is adopted to determine the cumulative postseismic deformation in each time period. According to the postseismic deformation of each time period and the fault model constructed in postseismic slip inversion, the postseismic slip of each time period is inverted under the layered medium model. As depicted in Fig. 9, two slip patches are formed on both sides in the deep, which are the large area on the north side and small area near the source, respectively. After that, the fault slip is continued along the two patches, and the area and magnitude of fault slip is enlarged accordingly. The spatial distribution of fault slip varies significantly in the first 4 time periods, while there is no obvious variation in the last 2 time periods. The result suggests that the fault activity is more obvious during the 180 days after the earthquake, and the fault tends to be stable during the time span of 270-365 days. Moreover, the postseismic slip during the time period of 0-30 is the largest, and the maximum of fault slip and corresponding slip rate reaches 0.62 m and 20.6 mm/day (see Table 2). Compared with the first time period, the rate of fault slip during the time period of 30-60 is significantly decreased. The maximum of fault slip during the second time period is 0.22 m, and the corresponding slip rate is 7.3 mm/day. For the following 2 time periods (i.e., 60-120 and 120-180 days), the maximum of fault slip is 0.23 and 0.14 m, and the rates of fault slip are continuously decreased, which is 3.8 and 2.3 mm/day, respectively. For the time period of 180-270 and 270-365 days, the maximum of fault slip is 0.14 and 0.11 m, and the corresponding rate of fault slip is 1.6 and 0.9 mm/day, which is relatively small in comparison with other time periods. The above results indicate that postseismic slip occurs during the 1 year after the earthquake, and it reaches largest during the time period of 0-30 days. Meanwhile, the maximum rate of fault slip during the last 2 time periods is only around 1 mm/day, suggesting that large postseismic slip mainly occurs in the early stage, and the fault tends to be stable as time goes on.

Coulomb stress change
Generally, coseismic stress tends to transfer after large earthquake, which may relate to seismic activities around the coseismic rupture area (Steacy 2005;Freed 2005). Thus, investigating the Coulomb stress change induced by coseismic slip helps to understand the triggering mechanism of postseismic fault slip motions and assess the postseismic disaster risk. The Coulomb stress change can be calculated based on the Coulomb failure criterion (King et al. 1994;  where CFS denotes the Coulomb failure stress (CFS) change, �τ s and �σ n are the shear and normal stress change on the receiver fault, and µ ′ is the equivalent friction coefficient. In this paper, the PSGRN/PSCMP developed by Wang et al. (2006) is utilized to estimate the static Coulomb stress change, and a moderate effective coefficient (i.e., 0.4) is selected as the equivalent friction coefficient. Figure 10 presents the spatial distribution of coseismic Coulomb stress change induced by the Illapel earthquake at the depth of 10 and 20 km. The coseismic Coulomb stress change demonstrate that the magnitude of stress increase in the main rupture area induced by the Illapel earthquake reaches several MPa, and the Coulomb stress change present complex characteristics and diversity at different depths. For the depth of 10 km, the magnitude of Coulomb stress change ranges from − 4 to 7 MPa, and the stress shadow covers most region, suggesting that the stress in this region is mainly released by coseismic rupture. As for the depth of 20 km, the Coulomb stress change shows different spatial distribution in comparison with the depth of 10 km. The Coulomb stress change appears in the region surrounding the coseismic slip and present a good correlation with the aftershocks and afterslip activities. The Coulomb stress increase also forms two peak areas on both the north and south sides, and the magnitude reaches about 2 MPa. The Coulomb stress change between the two peak areas is negative, and most aftershocks also gather in this area, implying that the coseismic Coulomb stress in this area is uploaded by these aftershocks. Meanwhile, the result of Sect. 3.3 indicates that the postseismic slip is mainly cumulated at the depth of 10-50 km, and the maximum slip is situated at the depth of 18.94 m. The spatial distribution of afterslip suggests that the afterslip is extended to both deep and two sides, and two peak slip patches are formed on the north and south sides. Thus, on basis of the above analysis, we conclude that the postseismic slip after the Illapel earthquake may be triggered by the stress increase in the deep region induced by coseismic phase.

Conclusion
In this study, the daily GNSS position time series of 13 sites with time span of 1 year after the earthquake is utilized to characterize the postseismic deformation mechanism and spatial-temporal properties of postseismic slip. Firstly, we design a scheme to extract the postseismic signal from GNSS position time series, which can extract and model the postseismic signal with consideration of the background tectonic movement and the influence of noise. Then, the coseismic slip is determined based on the coseismic deformations and fault model, which can be regarded as a reference to understand the postseismic deformation effects. At last, according to the postseismic signals and fault model, the spatial-temporal properties of afterslip is explored. The main conclusions are listed as follows: (1) The postseismic signal extraction and modeling method introduced in this paper can effectively extract the postseismic signal and establish simulation function model with consideration of the influence of background tectonic movement. Meanwhile, the parameters of GNSS position time series are estimated by means of holistic modeling, which can effectively reduce errors introduced by segmented processing.
(2) Comparing with coseismic slip, the afterslip is extended to both deep and two sides, and two peak slip patches are formed on the north and south sides. For the two peak slip areas, the larger slip area is situated at the north side, while the smaller slip area is located near the source. The afterslip is mainly cumulated at the depth of 10-50 km, and the maximum of postseismic slip reaches 1.46 m, which is situated at latitude of − 30.50°, longitude of − 71.78°, and depth of 18.94 m.
(3) The postseismic fault slip occurs during the 1 year after the earthquake, and it reaches largest during the time period of 0-30 days. The maximum rate of fault slip during the last 2 time period is only around 1 mm/day, suggesting that large postseismic slip mainly occurs in the early stage, and the fault tend to be stable as time goes on. Meanwhile, the Coulomb stress change denote that the stress increase in the deep region induced by coseismic phase contribute to the occurrence of afterslip associated with this earthquake.