Fault slip rate estimates along the Altyn Tagh and Kunlun strike-slip faults in northern Tibet vary considerably between short-term geodetic and long-term geologic studies. Here we reanalyze and model all global positioning system (GPS) data from northern Tibet to determine if these differences might be explained by previously unmodeled transient processes associated with the earthquake cycle, which can bias slip-rate estimates from geodetic data. We find that these effects cannot reconcile the geodetic data with the lowest bounds on the geologic slip rates along these faults, even in the presence of low (<1018 Pa s) viscosities within the mid-crust or crust and mantle lithosphere. Surface velocities derived from GPS measurements are best reproduced with models with a high-viscosity (≥1018 Pa s) middle to lower crust and mantle lithosphere.
Fault slip rates inferred from space geodesy sometimes disagree with those determined from geologic studies, which has led some workers to question whether fault slip rates inferred from geodetic surface velocities are representative of those averaged over longer time scales (Mériaux et al., 2004). Nowhere has this discussion been more contentious than in northern Tibet, along the left-lateral Altyn Tagh and Kunlun faults, where slip rate estimates have been used to argue for tectonic models as diverse as eastward rigid-block translation (e.g., Molnar and Tapponnier, 1975) and distributed deformation throughout the lithosphere (Houseman and England, 1986; Royden et al., 1997). Each of these models has support from different data sets: low slip rate estimates of 4–10 mm/yr from geodetic data collected across an ~250 km swath perpendicular to the Altyn Tagh fault (Bendick et al., 2000; Wallace et al., 2004) are cited by some as evidence that deformation is distributed among many faults, while inferences of slip rates of 20–34 mm/yr from geologic data are cited as evidence that since at least 6 ka, eastward translation of rigid Tibetan lithospheric blocks may instead accommodate a large portion of the Indo-Eurasian convergence (Mériaux et al., 2004; Peltzer et al., 1989).
Some workers have sought to attribute the discrepancy in slip rate estimates to systematic errors in geologically determined fault slip rates (e.g., Cowgill, 2007). Alternatively, estimates of fault slip rate from geodetic data may be biased by models that neglect the episodic loading and viscous relaxation of the middle to lower crust and mantle lithosphere that would cause surface deformation rates to vary throughout the earthquake cycle (Hetland and Hager, 2006; Hilley et al., 2005; Johnson and Segall, 2004; Pollitz, 2001; Savage and Prescott, 1978). If transient deformation is neglected in models of geodetic data, inversions could systematically underpredict or overpredict slip rates during different times in the earthquake cycle (e.g., Savage and Prescott, 1978). Previous studies neglected transient deformation processes in models of geodetic data and estimated slip rates by assuming that slip below a locked zone accrues steadily in an elastic half-space over the course of the earthquake cycle. We pose the question, how are surface global positioning system (GPS) velocities, crust and mantle rheology, and long-term fault slip rates within Tibet related through earthquake-cycle effects, and to what extent are previous geodetic estimates of fault slip rates biased by oversimplified models of interseismic deformation? Herein we provide an answer to this question by reanalyzing GPS data from central and northern Tibet and conducting a model exploration of the surface velocity data in the context of time-dependent earthquake-cycle deformation along the Altyn Tagh and Kunlun faults.
DATA ANALYSIS AND MODELING METHODS
Regional (Shen et al., 2000; Wang et al., 2001; Zhang et al., 2004) and local (Bendick et al., 2000; Wallace et al., 2004) GPS networks constrain the far-field and fault-related motions of northern Tibet, but these networks have never been integrated into a consistent reference frame. We have reprocessed data from regional and near-fault central Asian GPS stations into a common reference frame using the GAMIT/GLOBK software (Herring, 2005; King and Bock, 2005) and applied appropriate corrections to the measurements collected since the 2001 Kokoxili earthquake (see the GSA Data Repository1). We find that although station velocities that were formerly only referenced to a local network (Bendick et al., 2000; Wallace et al., 2004) have large uncertainties, they are nonetheless consistent with the regional station velocities nearby (Fig. 1). In our model analysis, we describe all station velocities with respect to the stable interior of the Tarim block. We define the rigid-block rotation of Tarim using GPS stations located within its interior, and subtract this motion from the remainder of the geodetic surface velocity field.
These geodetically measured surface velocities (Fig. 1) record the interseismic accrual of strain, while fault slip reflects the release of this strain during frictional failure. Thus, one must model long-term fault slip rates based on limited snapshots of strain accrual rates. Three general classes of models have thus far been employed for this purpose. (1) The interiors of tectonic blocks far from fault boundaries undergo rigid-body rotation (Thatcher, 2007). (2) Rigid-body block rotations are assumed far from fault boundaries, and steady interseismic elastic strain accumulation near faults is modeled with dislocations in an elastic half-space (Savage and Burford, 1973), with the assumption that the elastic strain accumulation is completely recovered during earthquakes (Meade, 2007). (3) The current deformation field is assumed to result from distributed viscous flow within the lithosphere (Bendick and Flesch, 2007; England and Molnar, 2005). These models represent various approximations of processes that are likely occurring at different depths in the lithosphere: the upper crust deforms primarily by brittle faulting in which elastic strain accumulation along block boundaries is released during earthquakes, while the middle to lower crust and mantle lithosphere probably flow viscoelastically (e.g., Nur and Mavko, 1974; Savage and Prescott, 1978). In this system, the coupled, time-variable deformation of the elastic upper crust and the flowing lithosphere leaves a potentially large time-dependent imprint on the interseismic surface velocities. When the middle to lower crust and upper mantle readily relax under earthquake-generated stresses, far-field surface velocities may be reduced during much of the earthquake cycle: this may cause the appearance of slow strike-slip deformation rates over broad areas when, in reality, long-term slip rates may be rapid. Deep afterslip or localized shear within the middle to lower crust may also contribute to accelerated surface deformation rates early in the earthquake cycle (e.g., Johnson and Segall, 2004). Time-dependent viscous flow at depth has the potential to reduce velocities at large distances from the zone of strike-slip deformation, which is the requisite condition to accommodate both rapid, long-term strike-slip rates and low interseismic surface velocities.
We assess the impact of earthquake-cycle deformation on surface velocities by modeling elastic strain accrual and release in the upper crust coupled to flow in the underlying Maxwell viscoelastic lithosphere. In the model we treat the upper crust as a series of deforming and rotating blocks in an elastic plate. These blocks are spherical caps that rotate about rotation poles. Small components of convergence or divergence along the block boundaries are accommodated using elastic opening-mode dislocations embedded within the upper elastic plate, such that long-term deformation accrues in the crust surrounding faults as convergence between blocks proceeds. Strike-slip motion between blocks is treated by prescribing a set of vertical faults whose long-term strike-slip rates are determined by the relative motion between rotating blocks. Transient deformation associated with steady interseismic locking and sudden coseismic slip events along the bounding faults is incorporated in the model. For simplicity, we assume that earthquakes occur periodically. At the beginning of each earthquake cycle, an earthquake occurs on a fault in the elastic plate of thickness H, with slip magnitude that is equal to the long-term slip rate on the fault multiplied by the earthquake recurrence time (T). These sudden slip events load the visco-elastic lithosphere below, which relaxes stress over a characteristic Maxwell relaxation time (τ) equal to twice the viscosity divided by the shear modulus of the flowing lithosphere (Nur and Mavko, 1974). The flow in the lithosphere coupled to the elastic crust creates a transient postseismic and interseismic phase in surface velocities with duration and magnitude that depends on the ratio of τ/T and the time since the most recent earthquake relative to the recurrence interval (t/T).
We consider three deforming blocks: the Tarim block, bounded on the south by the Altyn Tagh fault, the Qaidam block, bounded on the north and south by the Altyn Tagh fault and Kunlun fault, respectively, and the Songpan block, located south of the Kunlun fault (Fig. 1). Using the measured surface velocities, we determine the best-fitting rotation poles that define the slip rate on bounding faults using a weighted least squares inversion. The reduced chi-square statistic (χr2), defined as the ratio of the weighted residual sum of squares to the number of degrees of freedom in the model (number of each component of observed velocities minus the number of free model parameters), gauges the misfit between modeled and observed velocities in our analysis. Because the time-dependent component of deformation varies nonlinearly with t/T and τ/T, we compute Greens functions for each component of the rotation poles to perform our least squares inversion for each set of t/T and τ/T considered. By systematically varying t/T and τ/T while inverting for the best-fit rotation pole and computing the χr2, we are able to determine how the rheology of the crust and mantle affect the rotation poles, and hence the estimated long-term slip rates along the Altyn Tagh fault and Kunlun fault. Our simple two-layer model additionally assumes that earthquakes along all faults are synchronous and load an infinitely deep viscoelastic material. This provides an end-member scenario in which reduction of far-field velocities due to earthquake-cycle effects are most pronounced.
We varied H, τ/T, and t/T, and inverted for the rotation poles to determine how the layering and rheology of Tibet may affect surface velocities (Fig. 2). Low values of τ/T (corresponding to viscosities <~1017 Pa s for recurrence times of between 300 and 1000 yr; e.g., Washburn et al., 2001) cause viscous relaxation of the middle to lower crust and mantle lithosphere (referred to hereafter as plastosphere, after Scholz, 1988) to broadly distribute and reduce surface velocities relative to their long-term average throughout much of the earthquake cycle (t/T) in areas south of the Kunlun fault (Fig. 2). Depending on the time in the earthquake cycle and the viscosity of the flowing layer (encapsulated by τ/T), flow of the plastosphere may allow model slip rates along the Altyn Tagh fault to reach 12 mm/yr while still matching observed surface velocities. This behavior is generally enhanced by a thicker, elastic crust (larger H) that ruptures during earthquakes (Fig. 2), which loads deeper portions of the flowing plastosphere and forces them to flow more rapidly during the earthquake cycle. As a result, if viscosities of the plastosphere are low, depending on the time in the earthquake cycle, the long-term slip rate along the central Altyn Tagh fault would be systematically underestimated by as much 5 mm/yr in standard slip-rate analyses of the geodetic data (Fig. 2) (e.g., Bendick et al., 2000; Meade, 2007; Savage and Burford, 1973; Wallace et al., 2004).
Should viscosities underneath Tibet be low, our model indicates that slip rates along the Kunlun fault could be as high as ~30 mm/yr, if surface velocities are measured late in the earthquake cycle. Along the Altyn Tagh fault, inferred slip rates approach the lower end of estimates derived from block models as the lower layer's viscosity increases (Fig. 2A). Summed slip rates across the Altyn Tagh fault and Kunlun fault generally increase with decreasing τ/T; however, this increase is accommodated primarily along the Kunlun fault. Nonetheless, in all of the scenarios considered, slip rates along the Altyn Tagh fault are within the broad bounds of 2–12 mm/yr based on our analysis. These geodetically derived slip rates are inconsistent with geologic estimates of 20–34 mm/yr. Geologic estimates of slip rates along the Kunlun fault (9–16 mm/yr; Kidd and Molnar, 1988; Van der Woerd et al., 2002) are much lower than the maximum values permitted by our model. Thus, earthquake-cycle effects cannot reconcile the discrepancy reported between estimates of fault slip rate along these two faults.
The surface velocity field produced by selected models (denoted in Fig. 2) shows how predicted surface velocities compare to those measured for different t/T, τ/T, and H (Fig. DR1). Best-fitting models are those with high viscosities (Fig. 2C), as lower-viscosity scenarios systematically underpredict surface velocities within the Songpan block (Fig. DR1). Also, geologically estimated fault slip rates along the Kunlun fault (9–16 mm/yr) are more consistent with those predicted by high-viscosity scenarios (Fig. 2B). As the thickness of the elastic upper crust increases, observed surface velocities are matched best by high τ/T scenarios (Fig. 2). When all other factors are held constant, increasing t/T requires larger τ/T in the middle to lower crust and upper mantle, while maintaining similar χr2 values (Fig. DR1). Thus, the observed surface velocities generally tend to favor higher values of τ/T and t/T.
These models provide first-order insight into how time-dependent earthquake-cycle deformation involving viscoelastic flow in the interior of the Tibetan Plateau affects the surface deformation and geodetically determined fault slip rates in northern Tibet. We acknowledge that we have made a number of assumptions in our analysis that require further explanation. First, our viscoelastic earthquake-cycle model idealizes the flowing portion of the lithosphere as a Maxwell linear elastic material that is homogeneous, isotropic, and extends to infinite depth. More complicated rheologies have been treated in the context of two-dimensional models (e.g., Hetland and Hager, 2005), and heterogeneities in viscosity within and across lateral, lithospheric-scale boundaries may further affect surface velocities. More complicated models that consider a layered rheology predict that when upper mantle viscosities are low relative to lower crustal viscosities, interseismic velocity gradients near the faults may be pronounced, whereas far-field velocities may be only a fraction of the total transform motion across the boundary late in the earthquake cycle (Johnson et al., 2007, their Fig. 5). This layering could conceivably reconcile high strain rates near the fault with lower far-field velocities that result from the relaxation of the upper mantle. However, layered rheological models of Tibet generally invoke low viscosities in the middle to lower crust relative to the upper mantle (e.g., Klemperer, 2006; Royden et al., 1997). Such a layered viscosity structure produces far-field velocities that are less than would be expected when upper mantle viscosity is used with our simple two-layer model (Johnson et al., 2007, their Fig. 5). A similar reduction in far-field velocities would be expected if a vertical low-viscosity zone extended to great depth beneath the faults; however, these low viscosities would need to persist throughout the middle to lower crust and upper mantle to produce similar effects to a laterally extensive low-viscosity zone.
Our choice of a simple rheology was for expediency and consistency with other earthquake-cycle models (e.g., Pollitz, 2001) and geodynamic models of Tibet. Generally, such simple models are also found to capture much of the observed time-dependent deformation following large recent earthquakes (Bürgmann and Dresen, 2008). With this simple rheology, we are able to treat three-dimensional effects of the northern Tibetan faults in a computationally efficient manner. Our three-dimensional approach captures the first-order relationships between the observed surface velocity field, the fault slip rates that separate rotating blocks, and the flow of the lower portion of the lithosphere during the earthquake cycle.
Our analysis of the GPS data in northern Tibet indicates that earthquake-cycle effects cannot reconcile the reported lower bound on long-term geologic slip rates along the central Altyn Tagh fault with geodetically observed surface velocities. Our models cannot reconcile these data even if plastosphere viscosities are low (<1018 Pa s) when earthquake-cycle effects should be most pronounced. Effective viscosities of 4–10 × 1018 Pa s (rheological scenario shown in Fig. 1 resulting in 5.6 mm/yr and 25.6 mm/yr of slip along the Altyn Tagh fault and Kunlun fault, respectively) were estimated from the observed postseismic deformation transients of the M 7.6 1997 Manyi earthquake near the western end of the Kunlun fault (Ryder et al., 2007), and suggest that the existence of a thin, low-viscosity layer is not required by the geodetic data. In general, misfits between model predictions and observed surface velocities decrease as the viscosity of the flowing lithosphere increases. In addition, cases in which viscosities are low produce slip rates along the Kunlun fault well in excess of those inferred based on surface offsets along its central portion (9–16 mm/yr; Kidd and Molnar, 1988; Van der Woerd et al., 2002). Thus, our models indicate that the observed surface velocity field is most consistent with a scenario in which relatively slow slip rates characterize the boundaries of elastic, rotating blocks, under which a relatively high-viscosity middle to lower crust and mantle are present (Bendick et al., 2000; Meade, 2007; Wallace et al., 2004).