## Abstract

The strength of the lithosphere controls tectonic evolution and seismic cycles, but how rocks deform under stress in their natural settings is usually unclear. We constrain the rheological properties beneath the Taiwan orogenic belt using the stress perturbation following the 1999 Chi-Chi earthquake and fourteen-year postseismic geodetic observations. The evolution of stress and strain rate in the lower crust is best explained by a power-law Burgers rheology with rapid increases in effective viscosities from ~10^{17} to ~10^{19} Pa s within a year. The short-term modulation of the lower-crustal strength during the seismic cycle may alter the energy budget of mountain building. Incorporating the laboratory data and associated uncertainties, inferred thermal gradients suggest an eastward increase from 19.5±2.5°C/km in the Coastal Plain to 32±3°C/km in the Central Range. Geodetic observations may bridge the gap between laboratory and lithospheric scales to investigate crustal rheology and tectonic evolution.

## INTRODUCTION

The rheology of the crust and upper mantle is a strong control over the size and recurrence interval of earthquakes (*1*) and the evolution of plate tectonics processes over geological time scales (*2*, *3*). The strength of the lithosphere is usually inferred from laboratory rock experiments, but the scaling from laboratory samples to plate boundaries (*4*), the uncertainties regarding in situ conditions (such as grain size, water content, and temperature) (*5*), and the limited laboratory constraints on transient creep (*6*–*8*) make it difficult to predict the local rheological behavior. These shortcomings highlight the need for a more direct method to explore rock rheology in natural settings across different spatial and temporal scales.

Postseismic deformation excited by the stress perturbation of large earthquakes provides an opportunity to probe into the rheology of rocks in the deep crust. The mechanisms commonly involved in postseismic transients include viscoelastic flow and continuous, mostly aseismic slip (afterslip) on faults surrounding the epicenter (*9*, *10*). Although afterslip localizes to the fault zone involved in the mainshock and decays rapidly, viscoelastic relaxation usually involves long-lasting, broadly distributed deformation at rates many times higher than of tectonic loading. Therefore, geodetic measurements of postseismic deformation allow us, at least in principle, to investigate the rheology of the deep crust.

Taiwan is a young and active orogenic belt, created by the convergence between the Eurasian and Philippine Sea plates with rates of up to 82 mm/year (*11*). On 21 September 1999, the moment magnitude (*M*_{w}) 7.6 Chi-Chi earthquake ruptured the Chelungpu fault (CLPF) in central Taiwan, generating notable postseismic deformation for more than a decade that provides a chance to study the rheology of the lower crust. Rapid afterslip on the southern segment of the CLPF in regions that exclude the large coseismic rupture area was found to dominate the early postseismic deformation (*12*), whereas the contribution of accelerated viscoelastic flow in the lower crust became gradually more substantial (*13*–*15*). However, the rheological behaviors in the deep crust remain unclear, in part because of the difficulty of probing the rheological parameters with forward modeling methods. To overcome this limitation, we use the recently developed methodology that extends geodetic imaging to include off-fault deformation (*16*, *17*) and constrain the rock rheology in situ.

We use 14 years of Global Positioning System (GPS) measurements following the Chi-Chi earthquake (Fig. 1) and formulate a joint inversion for slip on triangular fault patches (*18*) and strain in volume elements (*19*) to directly invert for the temporal evolution and spatial distribution of strain and strain rate in the lower crust. The method reveals the rheological behavior of rocks in their natural settings. Comparing the temporal evolution of stress and strain rate in the lower crust to various plausible rheologies approximated by spring-dashpot assemblies, we show that accelerated viscoelastic flow can be explained by a power-law Burgers rheology. This suggests that the rheological behavior beneath the Taiwan orogenic belt includes a nonlinear constitutive behavior at the steady state (*20*) and the transient creep effect (*21*, *22*). The shortening that occurs at low strength during the postseismic period may reduce the energy budget of mountain building in the Taiwan orogenic belt.

## RESULTS

### GPS data of postseismic deformation

We investigate the postseismic displacements following the Chi-Chi earthquake by continuous-recording and survey-mode GPS in Taiwan from 1999 to 2013 (see details in Materials and Methods). The dense station coverage in the rupture area provides strong spatial constraints on the near-field postseismic deformation. All GPS position time series are with respect to the site S01R (Fig. 1A). The secular velocity is removed by using GPS measurements from 1993 to 1999 before the Chi-Chi earthquake (fig. S1; see Materials and Methods). We decompose GPS time series at each site into two specified temporal basis functions that describe the short- and long-term postseismic displacements, respectively (Fig. 1B and figs. S2 and S3). Both early and late postseismic transients are well captured by GPS data available before 2001. Moreover, we incorporate continuous-recording GPS sites installed after 2001 in central Taiwan to provide crustal postseismic deformation in great detail (see Materials and Methods). We find that even a decade after the mainshock, velocities at a large number of stations are still 3 to 8 mm/year faster than those before the mainshock (Fig. 1B). The localized afterslip is unlikely to strongly influence the far-field sites in eastern Taiwan, suggesting that activated viscoelastic flow beneath the Taiwan orogenic belt after the Chi-Chi earthquake needs to be considered.

### Geodetic imaging of afterslip and viscoelastic flow

Imaging the postseismic deformation following the Chi-Chi earthquake can be formulated as a linear inverse problem with the localized afterslip on the CLPF and ductile deformation in the viscoelastic lower crust (see details in Materials and Methods). We mesh the CLPF with triangular patches to avoid overlaps and gaps between adjacent patches and divide the lower crust at depths from 20 to 40 km into 30 km × 30 km × 10 km cuboids above the elastic uppermost mantle (Fig. 2, A and B). We construct the deformable volume elements with their horizontal axes aligned approximately parallel and perpendicular to the strike of Taiwan orogenic belt (~N20°) and the direction of plate convergence (Fig. 1A). Additionally, we extend the brittle-ductile transition up to 10-km depth beneath the Central Range (CR) because of the low seismicity (*23*) and high heat flow (*24*, *25*) observed in this region (Fig. 3A), where a lateral contrast of viscosity is often inferred (*14*, *15*).

The coefficients of two basis functions (Fig. 1B and fig. S3) at each GPS site are used to jointly invert for time-dependent afterslip and lower-crustal strain. Afterslip on the fault is the sum of strike- and dip-slip components, and the strain tensor in each deformable volume element includes six components (three axial and three shear strains). A large amount of parameters leads to an underdetermined inverse problem. To avoid overfitting short-wavelength signals in the data, we regularize the fault slip and strain distribution by distance-dependent Laplacian operators. We also penalize the volumetric strain and the strain directions that are perpendicular to the coseismic deviatoric stress in each volume element (see Materials and Methods).

We evaluate our model recovery with a series of checkerboard tests (figs. S4 to S7; see Materials and Methods). We produce synthetic surface displacements from an input model with a unit strain tensor aligned with the coseismic stress tensor at each volume element and a unit thrust slip on the southern portion of the CLPF fault. Uncertainties in synthetic displacements are assigned according to real GPS observational errors. The recovery of models shows good reliability in the rupture area covered by the dense GPS network, whereas the regions with insufficient spatial coverage of GPS sites are unrecoverable. Inevitably, we find smearing between fault slip and viscoelastic flow due to smoothing constraints, yet inverted strain and slip patterns are recovered satisfactorily. Overall, the synthetic tests demonstrate that the inversion approach can image viscoelastic flow in the lower crust reasonably well, with a horizontal resolution of 100 km or finer.

Our preferred model shows a notable cumulative afterslip of ~1 m on the southern segment of the CLPF away from the large coseismic slip area in the north (Fig. 2A), in accordance with results from previous studies (*12*, *14*). The laterally distributed strain in the lower crust generally follows the contours of coseismic stress change (Fig. 2A), with the highest value of ~20 microstrains beneath the middle section of the CR (Fig. 2B). The inferred principal compressional and extensional strain axes beneath the CR gradually rotate counterclockwise from deep to shallow depths (Fig. 2B). This feature is well aligned with a near vertical compressional strain from the surface to 10-km depth indicated by pervasive normal-faulting aftershocks following the Chi-Chi earthquake (Fig. 3A) (*26*). The model produces a satisfactory fit to both horizontal and vertical postseismic displacements (Fig. 2, C and D) with 92% variance reduction for the entire dataset. Our model can explain the 14-year GPS observations (Fig. 2D), although part of the misfit may be related to other secondary features not taken into account in the model, such as poroelastic effects.

### Lower-crustal rheology of the Taiwan orogenic belt

We use rheological spring-dashpot assemblies, which are widely used to approximate rheological behaviors in the crust after the earthquake (*5*, *7*), to fit the inverted time-dependent strain rate in the CR and Coastal Plain (CP) (volumes 1 and 2 in Fig. 3A), where a dense coverage of GPS stations provides satisfactory estimates. The Maxwell rheology and the power-law rheology are analogs to steady-state diffusion creep (stress exponent *n* = 1) and dislocation creep (stress exponent *n* > 1), respectively. When an initial transient creep is involved, the Burgers rheology (*6*) and the power-law Burgers rheology (*8*), respectively, are applicable to model diffusion creep and dislocation creep. These rheologies simulate the temporal evolution of stress and strain rate in the lower crust after a stress perturbation, which can be compared with the results inferred from our geodetic inversion. The rheological properties governing steady-state creep include the shear modulus (*G*_{M}) and a steady-state material constant (*A*_{M}), and the others control transient creep including the work-hardening coefficient (*G*_{K}) and a transient material constant (*A*_{K}) (Fig. 3C; see details in Materials and Methods). We fit the inverted strain rate evolution at volumes 1 and 2 with the theoretical strain rate curve for different rheological models through grid searching (see details in Materials and Methods). The approach enables us to examine the stress and strain rate conditions as a function of time, directly revealing rheological behavior in the lower crust. In contrast, other numerical approaches require prescribing a rheology a priori for the forward model and select rheological properties by trial and error (*13*, *15*). Our approach is similar to the analysis of data obtained by laboratory rock creep experiments, except for the length scale relevant to the tectonics and the time scale of years to decades.

To first order, we find that neither a Maxwell (fig. S8A) nor a power-law rheology (fig. S8B) can explain the stress/strain rate paths. Steady-state rheologies underestimate the strain rate during the early stage of the viscoelastic transient and overestimate it at the late stage. This strongly suggests that steady-state creep alone, whether accommodated by diffusion creep or dislocation creep, is inadequate to explain the 14-year lower-crustal viscoelastic behavior, indicating that an initial transient creep is involved with the viscoelastic relaxation. Accordingly, the Maxwell rheology adopted by previous studies of the Chi-Chi postseismic deformation may only be applicable for a specific time period and is unable to explain observations across the entire postseismic period (*13*–*15*). This is in accordance with similar findings in southern California (*27*) and Turkey (*28*).

Taking transient creep into account, either with a linear Burgers (fig. S8C) or a power-law Burgers formulation (Fig. 3C and fig. S8D), largely improves the fit to the early stage of the viscoelastic behavior. However, the linear Burgers rheology fails to capture the strain rate evolution near the transition from transient to steady-state creep (fig. S8C), implying that the role of linear rheology in the lower crust may be secondary. A power-law Burgers rheology with a stress exponent of *n* = 3 (Fig. 3C) or greater (fig. S8D) can simultaneously approximate the early and late stages of the viscoelastic behavior in a self-consistent manner and best explain our data. Although we could not quantify the stress exponent precisely, our results are fundamentally incompatible with a diffusion creep–dominated rheology in the lower crust. Ambient noise tomography in Taiwan reveals the presence of strong anisotropy in the lower crust (*29*), supporting the implication that rock behavior is likely to be controlled by dislocation creep over the long term. Therefore, we suggest that transient creep followed by a steady-state dislocation creep dominates the behavior of accelerated viscoelastic flow in the lower crust after the Chi-Chi earthquake.

### Implications for thermal gradient

The rheological properties of rocks result from the interactions between stress, mineral content, grain size, fluid, pressure, and temperature (*5*, *8*). In spite of the disparate strain rates applied in laboratory experiments and found in a natural setting, the broad agreement of crustal rheology between inversion and laboratory results encourages us to investigate the plausible range of thermal gradients by fixing the laboratory parameters. Laboratory data may help us to further narrow down the range of these rheological parameters under sensible assumptions. The relationship between our optimal rheological model and flow law of rocks in a steady-state dislocation regime is given by(1)where *A*_{M} is a steady-state material constant in the Maxwell element of the power-law Burgers rheology model, *A* is a pre-exponential factor, *f*_{H2O} is the water fugacity related to confining pressure and temperature with its exponent *r*, *Q* is the activation energy, *p* is the confining pressure, *V* is the activation volume, *R* is the gas constant, and *T* is the absolute temperature.

We adopt the equation of state for pure water (*30*) to isolate the effect of water fugacity and assume confining pressures of 450 and 750 MPa in the middle CR and CP, respectively. We use the flow laws of quartz (*31*–*35*) and feldspar (*36*, *37*) (table S1) from laboratory data to constrain the rheological parameters in the lower crust of Taiwan. We take steady-state material constants (*A*_{M}) obtained by models with stress exponents of either 3 or 4 to reconcile with the laboratory results. We then calculate the temperature in volumes 1 and 2 as the remaining unknown in Eq. 1. Taking into account the uncertainties in the activation energy for each flow law, our estimates give a wide range of temperature (fig. S9). Two of the flow laws for quartz proposed a large SD for the activation energy resulting in large uncertainties in temperature estimates (*31*, *32*), indicating that the rheological parameters vary with different conditions in laboratory experiments. Because of a large distribution of temperature range given by various felsic rocks, we then suggest the plausible temperature range that meets at least four of the flow laws adopted in this study. As a result, we suggest the temperature beneath the middle CR and CP ranging from 437° to 530°C and 423° to 557°C, respectively. These results generally agree with the brittle-ductile transition of 300° to 450°C for quartz and feldspar (*38*) and the distribution of seismicity in Taiwan (*39*), implying the present-day thermal gradients of 29° to 35°C/km above the depth of 15 km in the middle CR and of 17° to 22°C/km above the depth of 25 km in the CP (Fig. 3B).

Overall, the geothermal gradient increases from the foreland to the hinterland of the Taiwan orogenic belt. Although inferred thermal gradients are subject to uncertainties, they are consistent with a number of independent studies. The inferred high thermal gradient of the CR shallower than 15 km (Fig. 3B) is generally compatible with estimates of about 17° to 40°C/km from thermal modeling approaches using either surface heat flow measurements (*24*, *40*, *41*) or Raman spectroscopy of carbonaceous materials (*25*, *42*), 34° to 38°C/km derived from high seismic attenuation (*43*), and about 30°C/km inferred from the *Vp*/*Vs* ratio (*44*). Another estimate of more than 50°C/km based on the Curie point depth from magnetic data (*45*) deviates from other findings largely. We note that our model may predict a higher temperature when the lower crust contains mafic material. Resolving the difference in temperatures inferred from various components of rocks may be informative, but it is beyond the scope of this study.

### Role of transient creep

To better comprehend the strength evolution of lower-crustal rocks, we examine the spatial-temporal pattern of effective viscosity (the ratio of stress and strain rate) during the postseismic period (Fig. 4). To do so, we assume laterally homogeneous rheological properties beneath the brittle crust and extrapolate the best-fit rheological properties in volumes 1 and 2 (Fig. 3C) to the surrounding blocks at the same depth. We then assume a stress exponent of *n* = 3 and calculate the evolution of effective viscosity subject to the coseismic stress change of each volume element by a power-law Burgers rheology. The results show that the initial effective viscosity surrounding the ruptured area is on the order of ~10^{17} Pa s (Fig. 4A), rapidly increases to ~10^{19} Pa s within a year (Fig. 4B), and then reaches a steady viscosity of about 2 × 10^{19} Pa s after a decade (Fig. 4C). The effective viscosities as a function of time estimated from 14-year geodetic observations (Fig. 4, D and E) show that background viscosities in the lower crust are approximately 10^{19} Pa s, about two orders larger than that observed in the first year following the Chi-Chi earthquake. This indicates a strong variability of viscosities over different time scales in the crust. Inferred effective viscosities only vary within an order of magnitude if we use the best-fit rheological properties of each volume element instead of using uniform properties (fig. S10).

A fraction of the postseismic viscoelastic strain is facilitated by transient creep (Fig. 4, F and G). Because of a combination of the nonlinear constitutive behavior and transient creep, about 20% of the shortening that accrues during the 14 years of observation takes place with effective viscosities at least one order of magnitude lower than the long-term strength. We note that strain rate in the Kelvin element shows negative excursions, raising the effective viscosity slightly in the process as a result of dynamic equilibrium. The mechanical coupling between brittle and ductile deformation, most manifest during the postseismic period, implies that the modulation of the lower-crustal strength during the seismic cycle may reduce the mechanical energy required to shorten the lower crust during the Taiwan orogeny.

## DISCUSSION

Postseismic geodetic observations spanning different time periods may be explained by a variety of models with different rheological properties (*17*, *27*, *28*). However, an ideal rheological model should be able to satisfy time-dependent deformation throughout the entire seismic cycle or at least explain data spanning a decade when the strain rate is relatively steady (*28*). Such a model is also expected to match independent geophysical and geological observations such as seismicity, seismic anisotropy, and thermal modeling.

The rheological behavior of the lower crust following the 1999 Chi-Chi earthquake can be explored through a theoretical spring-dashpot model involving transient and steady-state creeps within a dislocation creep regime with a shear strain rate ranging from 10^{−12} to 10^{−14} s^{−1}. (Fig. 3C). The mechanism of accelerated viscoelastic flow seems to switch from a transient creep to steady-state creep approximately after a year beneath the CR and CP regions. An inferred steady-state viscosity of ~10^{19} Pa s and viscosity structure can provide valuable constraints on the geodynamics of the Taiwan orogeny. The high heat flow (*40*), thicker crust (*46*), and the absence of earthquakes (*23*) may be a consequence of a weak lower crust in the CR. Previous models for Taiwan mostly considered the influence of slab geometries, boundary conditions, and layered viscosity structures on strain orientations, surface uplift, and lithosphere dynamics (*47*–*49*). Studies on geodynamic models of the Himalayan orogen have shown promise for the contribution of the crustal strength heterogeneity to the observed topography, deformation pattern, and tectonic evolutions (*3*) and would be interesting avenues for future investigation in Taiwan.

The 14-year-long postseismic deformation following the 1999 Chi-Chi earthquake illuminates the rock rheology in the lower crust beneath the Taiwan orogenic belt. The study demonstrates a potential for geodetic observations to access rock rheological properties under their natural tectonic settings and is applicable to any region that experienced a large earthquake with sufficient geodetic coverage. Our findings highlight that dislocation creep assisted by a transient creep phase accommodated the viscoelastic deformation in the lower crust after the mainshock. While our constraints on crustal rock rheology are compatible with laboratory experiments at the steady state, our results stress the important role of transient creep during the postseismic deformation. Despite the large disparity of spatial and temporal scales between laboratory and nature, we find inverted rheological properties that are in broad agreement with laboratory rock experiments, predicting substantial variations of the thermal gradient between the CP and CR. Our rheological model reconciles geodetic observations, seismic anisotropy, and heterogeneous thermal gradients in the Taiwan orogenic belt. These parameters may be integrated in geodynamic models of the Taiwan orogeny to incorporate the short-term effect of the seismic cycle.

## MATERIALS AND METHODS

### GPS data processing

We collected GPS data from 1999 to 2013, including 89 continuous and 50 campaign-mode sites in Taiwan. GPS data used in this paper are processed with GAMIT 10.42/GLOBK 5.16 software (*50*) developed by the Massachusetts Institute of Technology and the Scripps Institution of Oceanography. The positioning method is based on the double-differenced carrier phase observations with the ionosphere-free linear combination of L1 and L2. To constrain the regional deformation pattern and improve the accuracy of positioning time series in Taiwan, data from International Global Navigation Satellite Systems Service, including IGB08 core stations (*51*), in the Asia-Pacific area were incorporated in the data processing. Daily solutions were aligned to the ITRF2008 reference frame with SDs of 8 and 14 mm in average for horizontal and vertical components, respectively. We eliminated bad daily solutions with SDs larger than 25 and 50 mm in horizontal and vertical components, respectively, which is equivalent to ruling out 1% of the time series data.

### Removing GPS secular velocity

To estimate postseismic displacements, we interpolated interseismic velocities at sites installed after the Chi-Chi earthquake using available GPS velocity data from 1993 to 1999 before the Chi-Chi mainshock (fig. S1) (*52*, *53*). We performed a triangular interpolation by constructing a Delaunay triangular mesh using GPS stations. Given known interseismic velocities *V*(**r**_{1}), *V*(**r**_{2}), and *V*(**r**_{3}) at three GPS stations with Cartesian coordinates **r**_{1}, **r**_{2}, and **r**_{3}, respectively, the velocity *V* of the position **r** inside the triangle is estimated by(2)where λ_{1}, λ_{2}, and λ_{3} are the barycentric coordinates (λ_{1} + λ_{2} + λ_{3} = 1) at the position of the interpolated point. In addition, the SD of the interpolated velocity is also calculated by the error propagation. We removed linear interseismic horizontal motions at each GPS site before investigating postseismic deformation. GPS interseismic vertical velocities before the mainshock were mostly <10 mm/year, except for some sites affected by the groundwater withdraw. Therefore, we removed a few GPS stations affected by severe subsidence associated with groundwater pumping on the alluvial fan (*54*).

### GPS time series analysis

We analyzed GPS position time series from continuous and campaign-mode stations by weighted least squares regression. The equation is given by(3)

Here, *t* is time and *X*_{0} is intercept; *S*_{n} and *C*_{n} are amplitudes for annual and semiannual periodic motions with ω = 2π/*T* (*T* = 1 year). *N*_{off} is the number of offsets due to earthquakes or instrumental changes at time *t*_{i}^{off}; *O*_{i} is the amplitude of the Heaviside function *H*(*t*). *A*_{1} and *A*_{2} are the amplitudes of the Chi-Chi postseismic transient, and *E*_{j} is the postseismic deformation associated to other events; *t*_{eq} = 1999.7208 (in decimal); *b*_{1}(*t*) and *b*_{2}(*t*) are two basis functions for the postseismic transient(4)(5)where *t*_{c} is the characteristic time controlling the time scale of the postseismic transient; *k* is a constant associated with the rate-strengthening frictional behavior, affecting the rate of postseismic transient (*55*). We decided the parameters of two functions (fig. S2) by minimizing the difference between observed and modeled GPS time series at selected GPS sites (i.e., 5936, CHEN, FLNM, I007, S167, SUN1, and YUSN, marked in fig. S3A). The selected GPS sites are operated before 2001 and located at central and eastern Taiwan. The first basis function, *b*_{1}(*t*), with *t*_{c} = 2 years and *k* = 5.5 captures the early transient associated with afterslip or transient creep. The second basis function, *b*_{2}(*t*), with *t*_{c} = 15 years fits the long-term transient deformation, which is mostly associated with the steady-state creep. We used *b*_{2}(*t*) with *t*_{c} = 1 year at time *t*_{j}^{post} to remove postseismic transients following other *M*_{w} > 6 events. Data with residuals exceed three times the root mean square after the regression is eliminated in each GPS time series.

For GPS data available before 2001, we were able to produce a satisfactory fit to the postseismic GPS time series of the Chi-Chi earthquake using two basis functions. For GPS stations installed after 2001, the amplitude of the first basis function is not able to be accurately estimated because of data gaps in the early postseismic period. We estimated the amplitude of the first basis function for stations unavailable before 2001 using sites installed before 2001 and the same interpolation method shown in Eq. 2. Then, we reformulated Eq. 3 into a matrix form(6)where *X* is the position time series of any GPS site and κ is a parameter controlling the weight of the priori amplitude for the first basis function . Our results are shown in fig. S3.

### Joint kinematic inversion for afterslip and viscoelastic flow

We inverted for afterslip on the fault patches and viscoelastic flow in deformable volume elements using the amplitudes of basis functions as a data vector, **d**, and a model vector, **m**, is defined as(7)where **m**_{ss} and **m**_{ds} are the strike- and dip-slip components of fault slip, and **m**_{ij}, *ij* ∈ {11, 12, 13, 22, 23, 33}, are six components of the strain tensor in each volume element. The parameters containing fault slip and strain components are simplified as vectors **m**_{slip} and **m**_{strain}, respectively. Surface displacements contributing from fault slip and strains in volume elements can be written as(8)where ɛ denotes residuals in the data vector **d**, and Green’s matrix *G* is given by(9)

*G*_{ss} and *G*_{ds} are surface displacement results from unit strike- and dip-slip on fault patches, and *G*_{ij}, *ij* ∈ {11, 12, 13, 22, 23, 33}, represent surface displacements caused by a unit strain of each component in the strain tensor. We then simplified Green’s functions related to fault slip and strain into two matrices, *G*_{slip} and *G*_{strain}, respectively. *G*_{slip} is calculated using a triangular dislocation proposed by Nikkhoo and Walter (*18*), and *G*_{strain} is obtained by the analytical approach proposed by Barbot *et al*. (*19*).

We imposed constraints on fault slip and strain in volume elements to avoid overfitting data. For slip on the fault, we expected a gradual spatial change of slip amplitude and azimuth. To do so, we used a distance-dependent discrete Laplacian operator for triangular elements modified from Maerten *et al*. (*56*)(10)where *i*, *j*, *h*_{ij}, and *L*_{i} denote the corresponding triangular element, neighbor elements on three edges, the centroid distance between elements *i* and *j*, and the sum of the centroid distances, respectively. *s*_{X} and *d*_{X}, *X* ∈ {*E*, *N*}, are east and north components of the unit vectors for strike- and dip-slip directions. The modified operator provides a smooth spatial transition for both slip amplitude and azimuth in adjacent elements and is appropriate for a complex fault mesh with different strikes and dips on the triangular element. Using Eq. 10, we regularized the fault slip as(11)where , , , and , respectively, are the corresponding matrices of east and north components for strike- and dip-slip as shown in Eq. 10, and the parameter α_{1} controls the weight of this constraint.

For the strain in each volume element, we regularized each component of strain independently with distance-dependent terms as(12)where *i*, *j*, and *h*_{ij} denote the corresponding volume element, neighbor blocks, and the centroid distance between elements *i* and *j*, individually. We did not include azimuthal constraint because each volume element was orientated in the same direction. The matrix form is given by(13)where *L*_{p} is the regularization matrix shown in Eq. 12 and the parameter α_{2} represents the weight to put on the smoothness of strain distributions. Additionally, we considered viscoelastic deformation as an incompressible flow and penalized isotropic strain by using(14)where ** I** denotes the identity matrix with its size equivalent to the number of volume elements and β is the parameter controlling the weight. Finally, we penalized the strain directions perpendicular to the directions of coseismic deviatoric stress change(15)where

*Pn ij*,

*n*∈ {1, 2, 3, 4, 5}, are the strain components

*ij*of the five strain directions that are orthogonal to a six-dimensional vector formed by the coseismic deviatoric stress in corresponding volume elements. The parameter γ is the overall weight of this penalization. The weight of each volume element is also scaled by the norm of the coseismic deviatoric stress. Five orthogonal vectors span the null space of the preferred proportion of each strain component aligned with coseismic deviatoric stress change and are obtained by the singular value decomposition.

We then formulated an equation with all constraints as(16)where *C*_{d} is a covariance matrix for each set of data vector **d**. We used the amplitudes of two time basis functions as the data vectors in our inversion and inverted model parameters independently. The slip and strain evolutions are the sum of results in two inversions. The fit to the data is quantified from the variance reduction (VR), defined as(17)where *d* is the data with SD σ and *r* is the residual.

### Checkerboard tests

We adopted an input model for the checkerboard tests by taking a normalized strain tensor aligned with the coseismic stress change at each deformable volume element in the lower crust. Because afterslip has a notable contribution to the Chi-Chi postseismic motion, we considered two scenarios for our input models. The first model assumes zero slip, and the second assumes a unit thrust slip on the southern portion of the fault. We first produced a forward model by using a checkerboard sampling of 4-by-4 for the lower-crustal strain without afterslip (fig. S4) and with afterslip on the fault (fig. S5). The variance of synthetic surface displacements was assigned to be identical to the uncertainties in observed GPS data. We used the same regularizations and constraints applied to the model described in the main text and jointly inverted for strain in the lower crust and slip on the fault. Although the model recovery is satisfactory, we found a little smearing between afterslip and viscoelastic flow localized in the ruptured area (figs. S4 and S5). This feature is often inevitable because of the regularizations applied in our model. We also tested another case with a checkerboard sampling of 3-by-3 for the strain in the lower crust (figs. S6 and S7), and the result was comparable to the case with a sampling of 4-by-4. This indicates that our model is able to distinguish viscoelastic flow from afterslip even if a little smearing effect is unavoidable.

### Spring-dashpot models for rheology

We reconstructed the strain rate evolution following the Chi-Chi earthquake in each volume element by(18)where *m*_{1} and *m*_{2} are model parameters for the first and second basis functions, *b*_{1}(*t*) and *b*_{2}(*t*), respectively. We replaced the basis function with their time derivative function(19)where(20)and(21)to obtain the temporal evolution of the strain rate.

We examined four types of spring-dashpot assemblies, including Maxwell, power-law, Burgers, and power-law Burgers rheology, to investigate the rheological behavior of the lower crust after a large coseismic stress perturbation. Maxwell and power-law rheology represent the steady-state creep behavior and can be simulated by an elastic spring in series with a viscous dashpot. For a Maxwell rheology, the relationship between stress and strain rate is linear (constant viscosity), corresponding to a steady-state diffusion creep, and the governing equation is given by(22)where τ is stress and is a background strain rate. *G*_{M} and η_{M} are shear modulus and viscosity in the Maxwell system, respectively. For a power-law rheology, the relationship between stress and strain rate becomes nonlinear (variable viscosity), corresponding to a steady-state dislocation creep, and the governing equation is(23)where *A*_{M} is a coefficient for a nonlinear dashpot with a stress exponent *n* (*n* > 1).

On the other hand, the effect of the transient creep of rocks is well known in the laboratory creep experiments as indicated by a rapid strain change during the early postseismic deformation (*6*–*8*, *27*). The Burgers system, where a Kelvin element and a Maxwell element are in series, is appropriate to incorporate transient creep in rheological modeling. A Burgers rheology represents a combination of steady-state and transient creep in the diffusion creep regime. The set of governing equations is given by(24)where and are the viscous strain rate in the Maxwell element and Kelvin element, respectively. *G*_{K} and η_{K} are the shear modulus and viscosity in the Kelvin system, respectively. Power-law Burgers rheology is similar to the Burgers rheology, but the relationship between stress and strain rate is nonlinear. Power-law Burgers rheology can model a combination of steady-state and transient creep in the dislocation creep regime. The set of governing equations is(25)where *A*_{K} is a coefficient for a nonlinear dashpot with an exponent *n* (*n* > 1) in the Kelvin element.

We search for the optimal parameters of each rheological model by fitting the inverted strain rate with theoretical estimates at our interest areas (CR and CP for volumes 1 and 2, respectively, in Fig. 3A). Pre–strain rate is assumed to be 1.3 × 10^{−14} s^{−1}, inferred from the convergence rate of 82 mm/year over a 200-km-long distance across the Taiwan orogenic belt (*11*). Coseismic stress change of 0.6 and 0.3 MPa at CR and CP, respectively, is calculated from the coseismic slip model (*12*). Shear modulus (*G*_{M}) beneath the CR and CP are 30 and 40 GPa, respectively, according to the estimates from *Vp* and *Vs* tomography (*46*). On the basis of these assumptions, we searched for optimal rheological properties by minimizing the misfit between the inverted and experimental strain rates at volumes 1 and 2 for each candidate rheology. We found that a steady-state creep alone cannot explain the 14-year postseismic strain rate evolution in the lower crust (fig. S8, A and B). A Burgers rheology captures a rapid strain change during the early postseismic deformation but fails to produce the strain rate in the late stage (fig. S8C). A power-law Burgers rheology reconciles the early and late stages of viscoelastic relaxation (Fig. 3C and fig. S8D), indicating that the ductile deformation in the lower crust is probably dominated by a dislocation creep.

The rheological properties inferred from a power-law Burgers rheology include work-hardening coefficient *G*_{K} and material constants for steady-state creep *A*_{M} and transient creep *A*_{K}. The best-fit model with *n* = 3 gives *G*_{K} = 9 × 10^{10} Pa, *A*_{M} = 1.4 × 10^{−31} Pa^{−n} s^{−1}, and *A*_{K} = 2.5 × 10^{−30} Pa^{−n} s^{−1} in the middle part of the CR and *G*_{K} = 1.4 × 10^{11} Pa, *A*_{M} = 1.7 × 10^{−31} Pa^{−n} s^{−1}, and *A*_{K} = 3.3 × 10^{−30} Pa^{−n} s^{−1} in the CP (Fig. 3C). For *n* = 4, our model gives *G*_{K} = 9.0 × 10^{10} Pa, *A*_{M} = 5.0 × 10^{−38} Pa^{−n} s^{−1}, and *A*_{K} = 1.0 × 10^{−36} Pa^{−n} s^{−1} for the middle CR and *G*_{K} = 1.5 × 10^{11} Pa, *A*_{M} = 1.0 × 10^{−37} Pa^{−n} s^{−1}, and *A*_{K} = 2.0 × 10^{−36} Pa^{−n} s^{−1} for the CP (fig. S8D). The best-fit parameter for the Maxwell element, *A*_{M}, can then be used to estimate the temperature in the lower crust by assuming the felsic lower crust in the Taiwan orogenic belt (table S1 and fig. S9).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/2/eaav3287/DC1

Fig. S1. GPS horizontal velocity field from 1993 to 1999 before the Chi-Chi earthquake.

Fig. S2. Searching for optimal temporal basis functions.

Fig. S3. The entire GPS dataset including continuous-recording and survey-mode sites.

Fig. S4. Checkerboard test for 4-by-4 grids without afterslip.

Fig. S5. Checkerboard test for 4-by-4 grids with afterslip.

Fig. S6. Checkerboard test for 3-by-3 grids without afterslip.

Fig. S7. Checkerboard test for 3-by-3 grids with afterslip.

Fig. S8. Rheology and stress versus strain rate.

Fig. S9. Temperature-water fugacity relations and flow lows.

Fig. S10. Effective viscosities derived from a power-law Burgers rheology.

Table S1. Flow laws for quartz (qtz) and feldspar (fsp) used in this study.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank the associate editor, K. Furlong, and two reviewers, F. Capitanio and T. Byrne, for their constructive comments. We also thank F. Deschamps, E. Tan, and B. Kuo for helpful discussions. We are grateful to S.-B. Yu, H.-Y. Chen, H.-H. Su, Y.-C. Tsai, and H.-M. Lee, at the Institute of Earth Sciences, Academia Sinica, who have participated in collecting and processing continuous GPS data. The generous provision of continuous GPS data from the Central Weather Bureau, Central Geological Survey, and Ministry of the Interior, Taiwan, and the international GNSS Service community is appreciated. Generic Mapping Tools (GMT) was used to create several figures (

*57*).

**Funding:**This research was supported by the Institute of Earth Sciences, Academia Sinica (CDA-105-M05) and the Ministry of Science and Technology (MOST 104-2628-M-001-008-MY4). S.B. was supported by the National Research Foundation (NRF) of Singapore under the NRF Fellowship scheme (National Research Fellow Awards Number NRF-NRFF2013-04), the Singapore Ministry of Education (AcRF Tier 1 grant RG100/17), and by the Earth Observatory of Singapore, under the Research Centres of Excellence initiative. This work comprises IESAS 2326 and Earth Observatory of Singapore contribution no. 192.

**Author contributions:**C.-H.T., Y.-J.H., and S.B. designed the study. C.-H.T. conducted the study and wrote the manuscript with contributions from all other coauthors. All authors contributed to the discussion and interpretation of the results.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).