Kinematics and dynamics of the southeastern margin of the Tibetan Plateau

S U M M A R Y On the southeastern margin of the Tibetan Plateau lies a large region which seismicity and GPS data show to be actively deforming. This paper describes the active faulting in the region, and how it relates to the velocity field observed with GPS. In places the velocity field is accommodated by rotations about vertical axes, and most or all of the strain at the surface in the region appears to be released seismically. GPS velocities are then compared to velocities calculated using a model for deformation driven by gravitational driving forces. Using rheologies estimated from experimentally derived mineral flow laws, the model provides velocities that are in good agreement with observed GPS velocities. It is not possible to uniquely determine the rheology or flow velocity at depth, and there are two forms of model solution which match the observed horizontal surface velocities. In one of these, vertical planes deform by pure shear, and in the other vertical gradients of horizontal velocity are present within the crust. Two distinct regions of normal-faulting earthquakes are present in the region, and have mechanisms which are most easily explained by gravity-driven deformation.


I N T RO D U C T I O N
Seismicity and GPS data show that the southeastern margin of the Tibetan Plateau is actively deforming. The geomorphological expression of active faulting is often clearly visible in satellite imagery and the topography, numerous earthquakes have occurred during the instrumental period, and the surface velocity-field is well known (Shen et al. 2005;Socquet et al. 2006). It is therefore possible to compare these data sources to see how the active faulting is accommodating the observed velocity field. The dense measurements of the spatially variable velocity field also present an opportunity to test the hypothesis, dating back over 25 yr (e.g. Bird & Piper 1980;England & M c Kenzie 1982), that the long-wavelength deformation of continental lithosphere can be modelled as that of a continuous fluid deforming in response to gravitational driving forces. Advances in computing power mean that it is now possible to solve the relevant equations in a more complete form than was previously possible. These two topics are the subject of this paper.
In Section 2, the active faulting in the area is described, it is suggested how the faults accommodate the velocity field, and the style and rates of deformation inferred from GPS and seismological data are compared. In Section 3, a numerical model is described which has been used to recreate the major features of the GPS velocity field.

K I N E M AT I C S
The southeastern margin of the Tibetan plateau was the subject of a GPS campaign by Shen et al. (2005), and the motions of the surrounding regions have been studied by Socquet et al. (2006). The velocities obtained from these surveys are shown relative to Eurasia in Fig. 1, along with those of Zhang et al. (2004). Zhang et al. (2004) and Shen et al. (2005) both give velocities for sites on the eastern and southeastern margin of the Tibetan Plateau from the Crustal Motion Observation Network of China. In this paper the velocities of Shen et al. (2005) are used, as they are calculated using data from more occupations than those of Zhang et al. (2004). Socquet et al. (2006) give their velocities in the ITRF2000 reference frame, and these have been converted into a Eurasia-fixed reference frame using the pole of rotation of Altamimi et al. (2002). These converted velocities match the Shen et al. (2005) velocities well where stations appear in both networks [the three stations which appear in both networks have velocity differences between the Shen et al. (2005) and converted Socquet et al. (2006) velocities of 0.2 (WUHN), 0.5 (SHAO) and 1.1 (KUNM) mm yr −1 , and azimuth differences of 10 • (WUHN), 12 • (SHAO) and 8 • (KUNM)]. It therefore seems that velocity differences between stations in the network of Shen et al. (2005) and the surrounding stations of the Socquet et al. (2006) network can be treated as 'real', and interpreted in terms of tectonics, so long as the velocity differences being considered  Shen et al. (2005), those with white bases are from Socquet et al. (2006), and those with blue bases are from Zhang et al. (2004). To aid clarity, the Zhang et al. (2004)  are relatively large. Fig. 2 shows a more detailed view of the GPS velocities in the area considered in this paper. Fig. 3 shows the focal mechanisms of earthquakes which have occurred on the southeastern margin of the Tibetan Plateau. Except for the earthquakes shown in green, which are in the slab subducting eastwards under the central lowlands of Burma (e.g. Ni et al. 1989), all of the earthquakes in the area which have well-determined depths (those with more than 10 depth phase observations and a moment magnitude greater than 5.5 in the catalogue of Engdahl et al. 1998) are shallower than ∼15 km. The majority of the earthquakes which have occurred in the area have strike-slip mechanisms, although some normal faulting earthquakes (shown in red and pink) have also occurred. Upper crustal shortening is virtually absent except on the margins of the Eastern Syntaxis and the Szechwan Basin. Fig. 4 shows the known active faults in the area. The locations of these faults are known from Allen et al. (1984), Wang & Burchfiel (1997), Lacassin et al. (1998), Wang et al. (1998), Wang & Burchfiel (2000), Burchfiel & Wang (2003), and my own observations of topographic data and satellite images. The lowlands on the west of Fig. 4 are the central lowlands of Myanmar (Burma), where a thick pile of sediments overlies crust of unknown nature, although Curray et al. (1978) suggest it is oceanic. This lowland basin is cut by the Sagiang fault, a right-lateral strike-slip fault with a slip rate of ∼18 mm yr −1 (Vigny et al. 2003), which accommodates some of the northwards motion of India relative to southeast Asia. The eastern boundary of the deforming zone is the 'South China Block', between the eastern margin of the Tibetan Plateau and the eastern coast of China (the grey shaded area on Fig. 4), which Shen et al. (2005) show to be undeforming at the 2 mm yr −1 level. The Szechwan basin is part of this undeforming region. Active faults are likely to be present in the deforming area shown on Fig. 4 which have not been recognized (due, for example, to erosion or to a recent onset of faulting). However, the faults shown are likely to give a good indication of the type of faulting occurring in each area, and are sufficient for the purposes of what follows.

The area between the Szechwan Basin and the Eastern Syntaxis
The area between the Szechwan Basin and the Eastern Syntaxis is characterized by high topography dissected by deeply eroding rivers which, in places, obscure the geomorphological signs of active faulting. In the east (near the margin of the Szechwan Basin) the velocity gradient observed in the GPS data is equivalent to leftlateral strike-slip on N-NW striking planes, and is accommodated by N-NNW left-lateral strike-slip on one or more parallel faults, which are continuations from the Xianshuihe fault.
Further west, on the margin of the Burmese lowlands, the deformation is equivalent to right-lateral strike slip on N to NW striking planes. The most striking geomorphological features in this area are two long, straight, river valleys (at 25-28 • N 99 • E, Fig. 3). At least one of these valleys is thought to contain an active right-lateral strike-slip fault (e.g. Wu 1991;Wang & Burchfiel 1997), which is oriented to accommodate the velocity gradient by simple strike-slip motion.
Normal faulting earthquakes with ∼E-W fault planes have occurred at latitudes of ∼29-31 • N (coloured red and pink on Fig. 3). There are a number of relatively short (10-40 km) ∼E-W scarps bounding low relief basins in the area shaded green on Fig. 4, which the seismicity suggests are likely to be normal faults. However, most of the normal faulting earthquakes have occurred in areas where no faults are visible in the geomorphology, suggesting that either erosion is the dominant geomorphological process, or alternatively that the normal faulting is recent or highly distributed. The likely cause of the normal faulting is more easily explained after discussion of the dynamics of the region, so these faults will be considered further in Section 4.5.

The lowlands around the southeast border of the Tibetan Plateau
To the south of the high topography between the Szechwan Basin and the Eastern Syntaxis (discussed above), is a large area of active  Shen et al. (2005), those from white bases are from Socquet et al. (2006). faulting and seismicity (Fig. 3). The faults to the SW of the Red River fault will be described first (see Fig. 4 for a location map), then the area to the NE of the Red River Fault, then the normal faults in the area, and finally the Red River Fault itself.

NE striking left-lateral strike-slip faults (SW of Red River Fault)
To the south and west of the Red River Fault are numerous NE striking left-lateral strike-slip faults (Fig. 4). These faults have clear geomorphological expressions (Fig. 5a) and have produced a number of earthquakes (Fig. 3). The faults have previously been described by Lacassin et al. (1998), who suggested the sense of slip changed from right-lateral to left-lateral during the late Cenozoic. The western boundary of the zone of left-lateral strike-slip faulting is the eastern margin of the Burmese lowlands, which is approximately stationary in a Eurasia-fixed reference frame (Fig. 2). The eastern edge of the zone of left-lateral strike-slip faulting is near the Red River Fault, and is moving south in the Eurasia frame at ∼5-10 mm yr −1 . Therefore, the overall velocity difference across the zone of strike-slip faulting is equivalent to N-S right-lateral strikeslip. The faults, however, are left lateral and strike ∼NE, so must ro-tate clockwise to accommodate the N-S right-lateral shear (Figs 5b and c), in a similar manner to that described in eastern Iran by Walker et al. (2004). The geometry shown in Figs 5(b) and (c) is obviously simplistic, and takes no account of the deformation which must occur at the extremities of the blocks, but illustrates the general mode of deformation in the area. The suggestion that the velocity field is accommodated by rotations about vertical axes is in agreement with the findings of Wang & Burchfiel (1997), who used geological observations to suggest that Cenozoic clockwise rotations have occurred in the region of the strike-slip faults.

N-S striking left-lateral strike-slip faults (NE of Red River Fault)
NE of the Red River Fault, on the margin of the undeforming area of south China, the strain is equivalent to left-lateral strike-slip on north-south planes, and is accommodated by a number of subparallel N-S striking left-lateral strike-slip faults. These faults are the southern continuation of those described above which accommodate N-NW left-lateral shear along the western margin of the Szechwan Basin. Wang et al. (1998) describe ∼NE-SW striking right-lateral strike-slip faults at the southern end of these left-lateral faults  Zhou et al. (1983a,b), Jones et al. (1984), Ekström & England (1989), Molnar & Lyon-Caen (1989) and Holt et al. (1995) are shown in dark red if the rake is within 25 • of pure normal, and otherwise are shown in black. Earthquakes from the CMT catalogue with more than 70 per cent double couple are shown in pink if the rake is within 25 • of pure normal, are shown in green if the depth of the event in the catalogue of Engdahl et al. (1998) is greater than 50 km, and otherwise are shown in grey. The percentage double couple of a focal mechanism is a measure of how closely the mechanism represents the double couple source expected for an earthquake [see Jackson et al. (2002) for further description].
(including the fault which ruptured in the 1970 Tonghai earthquake (Zhou et al. 1983a)). These right-lateral faults are shorter than the left-lateral faults, and do not offset them. Wang et al. (1998) attribute the right-lateral faulting to the accommodation of N-S left-lateral strike-slip by anticlockwise rotations about vertical axes.

Normal-faulting earthquakes with N-S fault planes (Area B on Fig. 4)
Normal faulting earthquakes with N-S fault planes have occurred in the area near 100 • E 26 • N (coloured red and pink on Fig. 3). Faulting with a vertical component is clearly visible in the geomorphology (Figs 6 and 7). These normal faults represent extension approximately perpendicular to the direction of motion relative to Eurasia shown by GPS in the area between the Szechwan Basin and the Eastern Syntaxis. The role of these normal faults is most easily explained with reference to the dynamics of the region, so the possible causes of the normal faulting will be discussed further in Section 4.5. Shen et al. (2005) constructed profiles through GPS measurements perpendicular to the major faults in the area. They suggest a rightlateral slip rate for the Red River Fault of 2 ± 2 mm yr −1 at the northwestern end of the fault, and 1 ± 2 mm yr −1 at the central portion of the fault. Fieldwork by Allen et al. (1984) shows the fault to be active, and Weldon et al. (1994) suggest a slip rate of 1-4 mm yr −1 from a trenching study. However other faults in the area have slip rates in excess of 5 mm yr −1 , and it therefore seems that although the Red River Fault forms a striking geomorphological feature, it does not dominate the present-day kinematics of the area. The GPS measurements of Shen et al. (2005) (Fig. 2) show that around 7-10 mm yr −1 of north-south shortening at the surface must occur between the latitudes of ∼27 • N and the southern margin of their GPS survey at ∼22 • N. The lack of thrusting or reverse faulting earthquakes suggests that this shortening is not accommodated by crustal thickening, so must be accommodated by lateral transport away from the area of shortening. The Red River Fault may help to achieve such deformation by accommodating some of the east/southeastwards motion of South China with respect to the Central Lowlands of Burma.

Comparison of earthquakes and geodetically observed strain
Cubic splines have been used to create a continuous velocity field from the GPS data, which has then been smoothed using a radial gaussian filter with a full width of 400 km, to remove small-scale irregularities related to the errors of individual GPS velocity estimates. The smoothed velocity field can then be used to calculate the horizontal velocity gradient tensor, which can be split into a symmetric part (the strain rate tensor) and an antisymmetric part (the vorticity tensor) (e.g. Malvern 1969). Fig. 8 shows the magnitude of the second invariant of the horizontal strain rate tensor, and also the direction and magnitude of the principal axes. The calculated vorticity is shown in Fig. 9. The style of strain derived from the GPS data is clearly in agreement with the focal mechanisms of earthquakes (Fig. 3). Left-lateral strike-slip on north-south planes is present along the margin of undeforming south China, ∼E-W extension occurs in the vicinity of the N-S striking normal faults (near 100 • E 26 • N), and ∼N-S extension occurs near the area of E-W striking normal faults (at 29-31 • N). Near the western margin of the Szechwan Basin the strain rate field shows ∼E-W compression and ∼N-S extension, which is the result of left-lateral strike-slip on the northwest striking Xianshuihe fault (Fig. 4). The GPS data extends as far south as ∼22 • N, within the northern part of the zone of ∼NE striking left-lateral strike-slip faults which appear to be rotating clockwise about vertical axes (Section 2.2.1). It is interesting to consider what signal strike-slip faulting undergoing vertical axis rotations would be expected produce in GPS data. At one extreme, if the GPS receivers in an area of rotating blocks were all positioned on the same block, the velocity field would record no strain, but would accurately measure the rotation of the individual block. At the other extreme, if the receivers were all placed close to (and on both sides of) one of the block bounding faults, the velocity  Lacassin et al. (1998) suggest that the stretch of river between the points marked C and D represents a recent left-lateral offset, and that the stretch between B and C represents an older right-lateral offset. Panels (b) and (c) show schematically how north-south right-lateral shear can be achieved by anticlockwise rotation of E-W left-lateral strike-slip faults. Panel (b) shows the initial situation, which will evolve through time to the geometry shown in panel (c). field would record the strain due to deformation at the fault, and the calculated vorticity would appear to be in the opposite sense to that produced by the block rotations. In reality the situation will be somewhere between these two extremes. The strain rate tensor will provide information on the deformation at the faults, and the vorticity will show a combination of the effects of the block rotation and the faulting between the blocks. Appendix A describes some calculations which examine what signal a series of rotating faults  Fig. 6. The photograph was taken looking to the west. The surface of the alluvial fan to the right of the stream has been uplifted relative to the rest of the fan along the line marked with white arrows. Photographs courtesy of James Jackson. with similar sizes and rotation rates as the NE striking left-lateral strike-slip faults would be expected to produce in GPS data. These calculations suggest that, with a similar density of GPS receivers to southeastern Tibet, the vorticity field should be relatively smooth, and be of the same sign as, but smaller in magnitude than, the vorticity expected due to the block rotations alone. The GPS vorticity field in the region of the rotating left-lateral strike-slip faults is in agreement with these calculations. The principal axes of the strain rate tensor in area of the strike-slip faults show deformation consistent with left-lateral strike-slip on ∼NE striking planes.
It is possible to compare the GPS-derived strain rate with the seismic strain rate using the method of Kostrov (1974). This comparison is described in detail in Appendix B and, to summarize, it appears that in areas where comparison of the two quantities has been possible, the geodetic strain is accommodated mostly, or entirely, by seismic slip on faults. Therefore, earthquake focal mechanisms are likely to give a good indication of the style of strain occurring in a given area, and any aseismic deformation which occurs is not the dominant component of the strain.

Summary and discussion of the kinematics of the region
To summarize, the strike-slip earthquakes on the southeastern margin of the Tibetan Plateau are associated with faults visible in the geomorphology, and can be simply related to the velocity field observed with GPS. The NE-SW striking left-lateral strike-slip faults to the southwest of the Red River Fault accommodate the velocity field by rotating clockwise about vertical axes. The roles of the two distinct regions of normal faults are less obvious, and will become clearer after discussion of the dynamics of the region (below). It is likely that, at the surface, most of the strain in the areas discussed above is accommodated by seismic slip on faults, which means that earthquake focal mechanisms provide an accurate indication of the style of strain.
M c Kenzie & Jackson (1983) used simple kinematic models to show that in zones of distributed deformation, strike-slip faults can only exist in a stable (i.e. time invariant) arrangement if they are on the margins of the deforming zone, and in the centre of a zone of distributed deformation strike-slip faults will be rotated about vertical axes. The arrangement of faulting described above agrees with the suggestions of M c Kenzie & Jackson (1983); the left-lateral strike-slip faults on the margin of undeforming south China show no evidence of rotation, but the NE-SW left-lateral strike-slip faults to the SW of the Red River Fault, which are entirely within the deforming zone, are rotating.

DY N A M I C S
This section describes a model for the deformation of continental lithosphere in response to the horizontal pressure gradients which arise from topographic contrasts and the velocities on the boundaries of the model domain. The lithosphere is treated as a continuous fluid. The model calculates instantaneous velocities, having input the present-day topography and a rheology, which are then compared to the GPS measurements of Shen et al. (2005). Viscosities are estimated from experimentally derived mineral flow laws, taking into account both the temperature and strain rate dependence of viscosity. The model has a number of poorly known free parameters, for instance the rheology of the crust, which is estimated from a calculated temperature structure and an assumed mineralogy. These parameters have been varied, and it has been found that there are two forms of solution which satisfy the velocity field. First, the numerical model will be described, and then the two types of solution. The solutions presented below should not be taken too literally: there are a large number of rheological combinations which can produce velocity fields which resemble the GPS data, and the purpose of this modelling is not to suggest the values of the free parameters in the model, but is simply to present the different types of solution which are consistent with the constraints imposed on the model by the surface velocities, the Moho depth, and the rheology used for the mantle lithosphere.

Model description
The equations for fluid flow in which inertial terms are negligible are solved using the method described in Appendix C. The equations are solved for a layer with a stress-free upper surface and lithostatic pressure and zero shear stress on the lower surface. The upper boundary represents the Earth's surface, whilst the lower boundary represents the base of the lithosphere. The model surface velocities are relatively insensitive to the boundary condition used on the base of the lithosphere; if the velocity of Eurasia in the 'hotspot reference frame' given by Gripp & Gordon (1990) is imposed, a shear zone forms at the base of the lithospheric mantle but the surface velocities are relatively unaffected. The driving forces for the deformation are the horizontal pressure gradients in the crust resulting from the surface topography, and the velocity boundary conditions (described below). M c Kenzie et al. (2000) showed that the timescale to attain isostatic equilibrium is always much shorter than that related to the removal of crustal thickness contrasts, so it is assumed that surface elevation contrasts are isostatically compensated, which appears to be confirmed by the crustal thickness variations revealed by Figure 8. GPS-derived strain rate and principal axes of the horizontal strain rate tensor on the southeastern margin of the Tibetan Plateau. The colour represents the magnitude of the second invariant of the horizontal strain rate tensor (defined as ˙ i j˙ i j , where˙ i j is the strain rate tensor and summation over repeated subscripts is assumed). The black bars represent the directions and magnitudes of the principal axes of the horizontal strain rate tensor. Thick bars represent contraction, and thin bars represent extension. The white circles show the locations of the GPS sites used to create the continuous velocity field, and so the strain rates should only be considered reliable in the vicinity of these circles. The pink line shows the 1000 m topographic contour.
the receiver function study of Xu et al. (2007). Compensation is assumed to occur at the base of the crust, so the mantle part of the lithosphere deforms only in response to the velocities on the boundaries and the velocities within the lower crust.
The viscosities used in the model are calculated from experimentally derived mineral flow laws. The model uses a rheology for quartz deforming by dislocation creep for the upper crust (Carter & Tsenn 1987), a rheology for anorthite and clinopyroxene aggregates deforming by dislocation creep for the lower crust (Dimanov et al. 2005), and an olivine rheology for the mantle (Hirth & Kohlstedt 2003). Effective viscosities for olivine deforming by both dislocation creep and diffusion creep are calculated, and whichever gives the lowest viscosity at a given point is used. The geochemistry of igneous rocks in the area has led Wang et al. (2001), Guo et al. (2005) and Jiang et al. (2006) to suggest that the mantle lithosphere is enriched with subduction-derived fluids. The 'wet' versions of the olivine and diopside-anorthite aggregate flow laws have therefore been used for most of the area. The Szechwan Basin is thought to be underlain by Precambrian cratonic material (Burchfiel et al. 1995), so the 'dry' versions of the anorthite-diopside aggregate and olivine flow laws have been used in the area of the basin. Curray et al. (1978) suggest that the Burmese lowlands are underlain by oceanic crust, so the wet olivine flow laws have been used for the entire lithosphere in these lowlands. In addition to the use of mineral flow laws, a pressure-dependent failure criterion is used to model the brittle part of the crust, simulating Byerlee's law (e.g. Brace & Kohlstedt 1980), following the approach of Fullsack (1995) and Cattin & Avouac (2000). The maximum deviatoric stress which can be supported is calculated for each point, and the strain rate calculated from a previous iteration (see Appendix C) is then used to calculate an effective viscosity. The deformation can then be treated using a fluid dynamic approach. The effective viscosities calculated using this approximation to Byerlee's law are used wherever they result in a viscosity lower than that calculated from the mineral flow laws. Byerlee's law was chosen because the depth of the boundary between the upper layer (where the viscosity is given by the failure criterion), and the lower parts of the lithosphere (which have viscosities calculated from mineral flow laws), is roughly in agreement with the depth of the deepest earthquakes in the area (10-15 km). Because the upper crust is treated as a continuum, and individual faults are not specified, the model is most applicable in areas containing many spatially distributed faults, which is the case for most of the southeastern margin of the Tibetan Plateau (Fig. 4). . u x refers to eastwards velocity, and u y northwards velocity. Where the perimeter is marked 'T', the velocity on the boundary is tapered between the boundary conditions shown to either side. The cross-hatched area in the southeastern corner of the model domain is given the velocity of undeforming south China, as suggested by the GPS velocities, which are shown in grey. Arrows from red bases are from Shen et al. (2005), those with white bases are from Socquet et al. (2006), and those with pale blue bases are from Zhang et al. (2004).
The boundary conditions on the edges of the model domain are fixed velocities, taken from GPS measurements, and shown in Fig. 10. The eastern boundary is given the velocity of the undeforming area of south China. The western margin is given zero velocity in the lowland area, and an eastwards velocity of 15 mm yr −1 in the area of high topography, with the velocity tapered between these conditions. The velocity of the southern margin is poorly known, so the boundary condition used is a linear velocity profile between the southeastern and southwestern corners, the velocities of which have been set by the conditions on the eastern and western margins. As the GPS velocities vary along the northern margin (within the Tibetan Plateau), for simplicity the velocity boundary condition used is a linear profile between the velocities at the northeastern and northwestern corners. The cross-hatched area in the southeastern corner of the model domain is given the velocity of undeforming south China, as suggested by the GPS velocities (Fig. 10). This undeforming area has been included in the model so that the eastern boundary of the model domain is at a longitude which also includes the Szechwan Basin, so the effects of using alternative rheologies for the basin can be investigated (see below).
The temperature field, which affects the viscosity, is calculated using a Lagrangian scheme (i.e. the reference frame moves with an individual column of rock, so that if there are no vertical gradients in horizontal velocity then horizontal advection can be ignored). Columns of rock on the high part of the Tibetan plateau, in the northern part of the model domain, have a temperature structure calculated following the method of M c Kenzie & Priestley (2007). In their method, a lithospheric column is suddenly thickened, and the temperature then evolves by vertical diffusion and radiogenic heating. The time which has elapsed since thickening (i.e. the age of the plateau) has been varied in different model runs. On the plateau margin, at surface elevations above 2000 m, lithospheric columns are assumed to have been derived from the high plateau, and to have thinned because of the motions to the southeast that are studied here. Low temperature thermochronology (Clark et al. 2005b), and the re-arrangement of major sediment catchment areas (Hall & Morley 2004) suggest that the southeastern margin of the plateau became elevated at around 10-15 Ma, at a time similar to the cessation of motion on left-lateral Ailao Shan shear zone (a major ductile shear zone which parallels the Red River Fault) at ∼17 Ma (Leloup et al. 2001). 10-15 Ma is short enough compared with the timescale of thermal relaxation for the lithosphere that it is possible to approximate the resulting temperature structure by simply vertically compressing a geotherm from the high plateau to fit with the current crustal thickness. Lithospheric columns in the relatively low-lying areas around the plateau margins are assumed not to have been derived from the high plateau, and a steadystate geotherm is calculated for these areas (using the method of M c Kenzie et al. 2005). Using the approach described above neglects the influence of advection into and out of vertical columns that are moving with the fluid, so assumes that the lithosphere deforms by pure shear. This assumption may or may not be accurate (see below), but calculating the temperatures in the manner described provides a starting point from which to investigate the behaviour of the system.
The model used in this paper differs from others used to study deformation in the India-Asia collision zone in a number of respects, as will be described below. In this section 'rigid' is taken to mean that the horizontal velocity on a boundary is zero and 'stress free' means that the shear stress on the boundary is zero, and there is no restriction on the horizontal velocity. 'Undeformable' means that the vertical velocity of a boundary is zero, and 'deformable' means that the normal stress is continuous. England & M c Kenzie (1982) developed the 'thin-viscous-sheet' approach, in which a power-law fluid is assumed to have stress-free upper and lower boundary conditions, and the horizontal velocity does not vary with depth. Copley & McKenzie (2007) made the same assumptions when studying the deformation in a transect between the Szechwan Basin and the Eastern Syntaxis. Ellis (1996) used a model in which the velocity could vary with depth, but the viscosity was constant. The model described here differs from these models by allowing both the physical characteristics of the fluid, and the velocity, to vary with depth. Royden (1996) and Shen et al. (2001) also created models in which the viscosity and velocity can vary with depth. The model presented here differs from theirs by having a viscosity that depends on temperature and strain rate, rather than on depth. In addition, no velocity constraints are imposed on the lithospheric mantle, rather than it imposing a rigid and deformable boundary condition on the base of the crust. Clark & Royden (2000) developed a model in which flow is confined to a fixed-width channel with rigid and undeformable upper and lower boundaries, an approach also used by Clark et al. (2005a), who extended the model to include rigid areas within the channel. The model presented here differs from these models by allowing flow to occur throughout the lithosphere wherever the viscosity is low enough, rather than being confined to a channel, and because it uses deformable upper and lower boundaries.  Shen et al. (2005). Panels (d) and (e) show how the speed and viscosity vary along the line x-y. The contours are labelled with speed in mm yr −1 and log 10 (viscosity in Pa s), respectively. The red line shows the Moho. In model results pictured the temperature in the high plateau was allowed to evolve for 40 Ma after thickening, and a quartz rheology was used to a depth of 28 km.

Model results
There are two forms of solution which can reproduce the surface velocities in southeastern Tibet. Examples of each of these types of solution are described below, and there is a continuum of behaviour between them. It should be noted that the two types of solution are indistinguishable on the basis of horizontal surface velocities, so it is not possible to infer the rheology or flow velocity at depth using only the horizontal surface velocities, and that in both solutions the motions are governed by flow in response to gravitational driving forces. If the gravitational driving forces are neglected, and the deformation occurs only as a response to the velocity boundary conditions, the velocities in the interior of the model domain vary smoothly between those imposed on the boundaries, and the large southwards component of velocity present in the GPS data is absent.
The first type of solution has relatively minor vertical gradients of horizontal velocity (Fig. 11). In the model results shown in the figure the geotherm for the high plateau was constructed by letting the temperature evolve for 40 Ma after thickening, and a quartz rheology was used to a depth of 28 km. The deformation of the lithosphere is close to pure shear, and so is similar to the 'thin-viscous-sheet' model of England & M c Kenzie (1982). An interesting feature is that pure shear deformation can occur even where there is a pronounced viscosity minimum in the mid to lower crust. In the case shown in Fig. 11, the vertical gradients in horizontal velocity are small because the low viscosity layer is too thin to flow as a high-velocity channel, and the lower crust has a high enough viscosity that it flows at a similar rate to the upper crust. It is important to note that low viscosities in the mid to lower crust do not necessarily lead to increased velocities.
The second type of solution has significant vertical gradients of horizontal velocity, and the lower crust and lithospheric mantle flow faster than the upper crust. The only difference between the model shown in Fig. 12 and that described above is that a quartz rheology was used to a depth of 35 km, rather than 28 km, which has the effect of weakening the middle and lower crust, so a horizontal shear zone develops in the middle crust and the lower crust and mantle flow significantly faster than the upper crust. High velocities in the lower lithosphere are present under the most steeply sloping part of the plateau margin (between the Eastern Syntaxis and the Szechwan Basin), but not under the flat surface of the high plateau or the lowlands in the south of the model domain. Situations where the flow at depth is faster than the flow at the surface, and where the surface velocities match the GPS velocities, only occur when a low Figure 12. As Fig. 11, but for a model result when the deformation involves significant vertical gradients of horizontal velocity. In the model results pictured all parameters are the same as the model shown in Fig. 11, except that a quartz rheology was used to a depth of 35 km rather than 28 km. viscosity zone exists in the mid crust, and when the lower crust and upper mantle are significantly less viscous than the upper crust.
The thermal state and mineralogy of the deeper layers of the lithosphere in southeastern Tibet are not known, so it is not possible to suggest which of the types of solution presented above is more likely to exist in the area. However it is clear that, for likely rheologies and geotherms, the horizontal surface velocities on the southeastern margin of the Tibetan plateau resemble those that would be expected for fluid flow in response to gravitational driving forces. These model results demonstrate that where the surface velocities are governed by horizontal gradients of horizontal velocity, the horizontal surface velocities allow us to place constraints on the effective viscosity of the shallow parts of the lithosphere, but give no information on the viscosity or flow velocity at depth.
It should be noted that the vertical velocities in the two types of solution are significantly different. In models where there is rapid flow in the lower crust and upper mantle there is a large flux of mass from beneath the high Tibetan Plateau, which means that relatively rapid crustal thinning occurs. In the model shown in Fig. 12 vertical thinning occurs at rates of up to 5 mm yr −1 . The mass flux at depth results in an increase in crustal thickness in the lowlands beside the plateau, with a maximum rate of thickening similar to the rate of thinning under the high plateau. Where the deformation is by pure shear (Fig. 11) the rate of thinning in the high plateau is simply a result of the divergence of the horizontal velocity, and is significantly smaller (∼1.5 mm yr −1 and under). It is therefore possible, in principle, to distinguish between these two types of behaviour on the basis of vertical velocities, which unfortunately are not known for southeastern Tibet at present. Burchfiel et al. (1995) found little evidence for upper crustal shortening on the southeast margin of the Tibetan Plateau in recent times, so Clark & Royden (2000) suggested that the plateau margin may have become elevated by crustal thickening caused by flow in the lower crust. If these observations are correct, then deformation of the type shown in Fig. 12 is most likely to be occurring.

Comparison of model and actual strain rate and vorticity fields
In this section the model strain rate and vorticity fields will be compared with the geodetic strain rate and vorticity fields, and the information on the style of strain given by the focal mechanisms of earthquakes and the geomorphological expression of active faulting. The calculation of the geodetic strain rate and vorticity fields is described in Section 2.3, and the model equivalents are calculated from the model surface velocities. The strain rate field calculated from the model results in Fig. 12 is shown in Fig. 13, and the vorticity field is shown in Fig. 14  Basin, and strain equivalent to strike-slip deformation (with P-axes directed towards the area between the Szechwan Basin and the Eastern Syntaxis) in the lowlands to the south of the Tibetan Plateau. The model strain rate field agrees well with the style of strain indicated by earthquake focal mechanisms and active faulting (Figs 3  and 4). One difference between the calculated and GPS velocities is that, in the model strain rate field, strike-slip deformation is not concentrated into a narrow zone in the north of the model domain, as it is in the GPS data (at the Xianshuihe Fault). The other main difference is that, in the region of N-S striking normal faults at 100-102 • E 26-28 • N, ∼E-W extension is clearly visible in the geodetic strain rate field but not as visible in the model results. The differ- ence is due to the velocity gradient corresponding to E-W extension being distributed over a greater distance in the model results than in the GPS data, although same overall difference in eastwards velocity is present. A possible reason for the two differences between the model results and the GPS data is that strength heterogeneities or anisotropies, which are not included in the model, may have focused the deformation. Alternatively, M c Kenzie & Jackson (1983) suggested that in zones of distributed deformation the only stable arrangements of faulting are those in which the motion on the faults is pure thrust, normal, or strike-slip. Therefore, another possible reason why the E-W extension is concentrated into the area near 100 • E 26 • N is that the faulting can be pure normal faulting (Fig. 3), rather than a mix of normal and strike-slip motion, as would occur if the extension occurred further west (as is seen in the model velocity field), in the area of N-S right-lateral strike-slip faulting (Fig. 4).
The model and geodetic vorticity fields are in good agreement. Both show positive vorticity (equivalent to anticlockwise rotations) in the east of the area, and negative vorticity in the west. The change of sign occurs at 100-102 • E in both the model and geodetic fields. The magnitudes are higher in the model, probably because the smoothing which was used during the construction of the geodetic strain rate and vorticity fields will smear the signal over a greater area and reduce the amplitude.

D I S C U S S I O N
I have described how active faulting on the southeastern margin of the Tibetan plateau relates to the velocity field observed with GPS. In places the velocity field is accommodated by strike-slip faulting and rotations about vertical axes. I have also shown that a model for gravity-driven flow can explain the main features of the observed velocity field. It is not possible to uniquely determine the rheology of the lithosphere, as a range of possible rheologies can produce velocity fields which match the GPS measurements of horizontal velocity. However, all of the modelled velocity fields which resemble the GPS data of Shen et al. (2005) are governed by flow in response to gravitational driving forces. It would be possible to distinguish between the two types of flow described above (those with and without vertical gradients of horizontal velocity) with knowledge of the vertical surface velocity. Where the lower crust and lithospheric mantle flow rapidly relative to the upper crust, the rate of thinning on the margin of the high plateau is significantly higher than if vertical planes deform by pure shear, and thickening is present in the lowlands bordering the plateau. However, the vertical velocities are not known for the southeastern Tibetan Plateau.
The two types of solution presented above are indistinguishable from each other on the basis of observations of the horizontal surface velocity, and a general result of this work is that in areas where surface motions are governed by horizontal gradients of horizontal velocity, modelling of the horizontal surface velocities will give information regarding the effective viscosity of the upper crust, but not of the deeper parts of lithosphere. In comparison, in areas where deformation occurs by simple shear above a rigid base, such as where the Indian lithosphere underthrusts southern Tibet, the surface velocities will be sensitive to the viscosity at all depths above the rigid base. Therefore, although Copley & McKenzie (2007) suggested a viscosity of 10 22 Pa s for the area between the Eastern Syntaxis and Szechwan Basin, and 10 20 Pa s for southern Tibet, the rheology of the upper part of the lithosphere (at depths above the level of the underthrust Indian lithosphere in southern Tibet) need not be different in the two places. The models presented here suggest that in the area between the Eastern Syntaxis and Szechwan Basin the modelling of Copley & McKenzie (2007) will have been sensitive to the effective viscosity of shallowest part of the lithosphere, whilst in southern Tibet the viscosity they estimated would depend on the rheology of everything above the depth of the underthrust Indian lithosphere, which is at a similar depth to the low viscosity layers in the mid crust in the models presented above.

Finite deformation
The finite deformation resulting from the two types of model solution presented in Section 3 should be similar. Where deformation of vertical planes is by simple shear, and the relief of the flow is large compared with the thickness of the flowing layer, a sharp topographic front is produced (e.g. Huppert 1982;M c Kenzie et al. 2000). Copley & McKenzie (2007) suggested the steep slopes of the Himalayan margin of the Tibetan Plateau were a result of such simple shear deformation above the underthrust Indian lithosphere. In the first model presented above (Fig. 11) vertical planes deform by pure shear. In the second model (Fig. 12) there is a strong zone of vertical simple shear in the mid crust, but the thickness of the flow is large compared with the relief on the upper and lower boundaries. Both of these flows would, therefore, be expected to form a relatively gently sloping margin, as is seen in the present day topography. The differences between the large-scale morphology of the plateau margin in the area of southeastern Tibet described in this paper, and the area further west which is underthrust by Indian lithosphere, may therefore be a result of differences in the nature and depth of the lower boundary condition.
It is possible to integrate rates of crustal thickening and thinning along flowlines in the model velocity fields, and so to calculate the topography which would result if the velocities did not vary through time. The topography calculated in this manner does not resemble the present day topography, which suggests that the velocities on the southeastern margin of the Tibetan Plateau have varied through time.

The Szechwan Basin and Xianshuihe Fault
The relatively rigid behaviour of the Szechwan basin requires it to have a different rheology to the surrounding material; if the same rheology is used as for the remainder of the model domain, large downslope velocities (up to ∼20 mm yr −1 ) are produced on the margins of the basin, and the interior of the basin deforms. The use of anhydrous flow laws for the Precambrian rocks of the basin is a simple method of producing relatively rigid behaviour, despite the large slopes on the margins (a ∼4000 m elevation change over ∼100 km). The ability of anhydrous minerals to withstand the forces generated by the steep slopes on the margins of the basin suggests areas of anhydrous lithosphere are likely to remain tectonically dormant for large periods of geological time, as was suggested by Jackson et al. (2004), who attribute the apparent strength of the Indian lithosphere underthrusting southern Tibet to the presence of anhydrous metastable granulite. The presence of anhydrous rocks may be why the Szechwan Basin has remained undeformed by orogenic events throughout Mesozoic and Cenozoic time.
The Xianshuihe fault extends northwestwards from the western corner of the Szechwan Basin (Fig. 4). If slip on the fault is purely strike-slip then material to the south of the fault will, if the motions are extrapolated forward through time, become part of the southeastwards flow between the Szechwan Basin and the Eastern Syntaxis. Material which is currently north of the fault will move slower and approach the plateau margin on the northwestern margin of the Szechwan Basin. The Xianshuihe fault can therefore be seen as the divide between the higher velocity flow being driven by downhill motions between the Eastern Syntaxis and the Szechwan Basin, and the slower movement of material towards the northernmost part of the Szechwan Basin. It is likely that the velocities to the north of the fault are slower than those to the south because the Szechwan Basin is inhibiting the motions. Given the lithosphere of the basin appears to be rigid enough to resist deformation, material can only move to the southeast by flow over the top of the basin (and thrusting in the brittle upper crust). Such motion will occur slowly because of the small thickness of material overlying a rigid base, which will provide large amounts of resistance to the flow. In comparison, for the flow between the Szechwan Basin and the Eastern Syntaxis the nearest sources of resistance are the relatively distant lowlands on the margins of the flow, so the motions are relatively rapid. The Szechwan basin does not simply move with the flow, as a rigid inclusion, so there must be forces on the eastern boundary of the basin resisting the motion, which may suggest that the lithosphere further east (in the region of the South China Craton) is also strong enough to behave in a rigid manner.

Comparison with other models
Other models of the southeastern margin of the Tibetan plateau (e.g. Clark & Royden 2000) consider flow in a channel confined to the lower crust. The calculations presented here suggest that if the mantle lithosphere is hydrated, as is suggested by the composition of igneous rocks (Wang et al. 2001;Guo et al. 2005;Jiang et al. 2006), it does not have a high enough viscosity to provide the relatively rigid lower boundary required for flow in a channel by simple shear. Rather, if increased flow velocities exist at depth, then the lithospheric mantle also flows rapidly, as illustrated in Fig. 12, and a channel does not form. The numerical experiments described above therefore suggest that in areas where the mantle is hot and hydrated, the lithospheric mantle is too weak to act as a rigid lower boundary to deformation within the crust.

Clockwise rotations about the Eastern Syntaxis
The first order feature of the velocity field is that, in the Eurasian reference frame, the northern portion of the modelled domain is characterized by eastwards velocities but further south the velocities are directed to the south, resulting in part of the clockwise rotation of velocities about the eastern syntaxis. In southern Tibet, in the area underthrust by rigid Indian lithosphere, the surface velocities are, in an India-relative reference frame, radial to the Himalayan arc and ∼15-20 mm yr −1 in magnitude (e.g. Jade et al. 2004). Although the surface of southern Tibet moves radially to the Himalayan arc in an India-relative frame because of gravity-driven flow (e.g. Copley & McKenzie 2007), the motion of India relative to Eurasia is northeastwards at ∼40 mm yr −1 , so the surface of southern Tibet moves north to northeastwards at 20-25 mm yr −1 in the Eurasian reference frame. Further east, in the area discussed here, there is no such northwards push by India so that, in the Eurasian reference frame, material is able to flow southwards, down the topographic slope. Taken together, these effects produce the rotation of velocity vectors about the Eastern Syntaxis. Copley & McKenzie (2007) suggested that the normal faulting earthquakes with east-west fault planes which have occurred in an east-west band at ∼30 • N were caused by extension in the downslope direction caused by increased flow velocities resulting from increased surface gradients. The concentration of extensional strain near the edge of the flat top of the plateau can be clearly seen in the model strain rate field shown in Fig. 13. The normal faulting earthquakes with east-west nodal planes only occur in the area where the model shows north-south extension, and the direction of maximum extension shown by the modelled strain rate and the normal fault focal mechanisms are in agreement. Fig. 15 shows that the locations of the normal faulting earthquakes agree well with the area of maximum extension seen in the model strain rate field. N-S extension shows that the flow of material between the Szechwan Basin and  Fig. 12, that is, the difference between the absolute magnitudes of the principal axes. This quantity is zero for strike-slip deformation, and is high in areas where the deformation involves thickening or thinning. As can be seen, the locations of the normal faulting earthquakes agree well with the area of maximum extension seen in the model strain rate field. the Eastern Syntaxis occurs because of gravitational driving forces arising from the surface slope, rather than because of forcing from the north.

Explanation for the normal faulting
A second population of normal faulting earthquakes have occurred with north-south fault planes at ∼27 • N 100 • E. Because the motions are driven by gravitational driving forces, south of the latitude of the southern margin of the Szechwan Basin the shape of the topography causes the velocities to fan outwards. On the western margin of the deforming zone the surface slopes due south, and on the eastern margin the slope is to the southeast. Because the Figure 16. A cartoon to schematically illustrate some aspects of the deformation discussed in this paper. The diagram on the right shows normal faulting due to increased downslope velocities on the edge of the plateau, and also because of the divergence in flow directions caused by the shape of the topography (see text). The diagrams A and B schematically show the deformation of an initially cuboid-shaped volume in the lower crust. Diagram A shows the deformation if the flow has only minor vertical gradients in horizontal velocity (Fig. 11). Thinning occurs under the high part of the plateau because of increased downslope velocities, and also lower down the plateau margin because the shape of the topography causes the flow to diverge. B shows the deformation if large vertical gradients of horizontal velocity are present, with the lower lithosphere flowing faster than the upper crust (Fig. 12). The thinning under the high plateau is greater than in A, because of the more rapid downslope velocities at depth. In the lowlands, the volume extends perpendicular to the flow direction, because of the same divergence as in diagram A. In addition, the volume also thickens vertically because the velocity at depth decreases as the downslope limit of the zone of high flow velocity is approached, which leads to horizontal compression and crustal thickening. gravitational driving forces act to move material in the direction down the surface slope, the variation in downslope direction leads to extension in the direction roughly parallel to the topographic contours, and is manifested at the surface as the normal faulting earthquakes with N-S nodal planes. Fig. 16 schematically shows the causes of the two regions of normal faulting, and the deformation resulting from the two types of flow discussed in this paper (those with and without vertical gradients in horizontal velocity).

C O N C L U S I O N S
It has been described how active faulting accommodates the velocity field on the southeastern margin of the Tibetan Plateau. As would be expected from simple models of faulting in zones of distributed deformation, the strike-slip faults on the margin of undeforming south China appear to have a stable arrangement, but the NE striking strike-slip faults which are entirely within the deforming zone are rotating about vertical axes. The surface deformation in the region appears to be mostly or entirely accommodated by seismic slip on faults. A numerical model for deformation in response to gravitational driving forces has been constructed, and the velocities from this model resemble those observed with GPS. Using the horizontal surface velocities alone it is not possible to infer the rheology at depth, and model calculations which predict large differences in the flow velocity at depth can match the surface velocities equally well. Knowledge of the vertical velocities would allow the flow velocity in the lower crust to be inferred. Two populations of normal faults on the southeastern margin of the Tibetan Plateau have mechanisms which are most easily explained as deformation resulting from the gravitational driving forces arising from the surface slope.

A C K N O W L E D G M E N T S
The author wishes to thank Dan M c Kenzie and James Jackson for many useful discussions and comments on the manuscript. John Beavan, Susan Ellis and Leigh Royden provided useful reviews. This work was supported by NERC funding to COMET (http:// comet.nerc.ac.uk).

A P P E N D I X A : M E A S U R E M E N T S O F S T R A I N R AT E A N D V O RT I C I T Y I N Z O N E S O F RO TAT I N G C RU S TA L B L O C K S
This appendix describes some calculations which simulate the strain rate and vorticity which would be recorded by a hypothetical GPS campaign conducted in a region of rotating crustal blocks with a similar aspect ratio and rotation rate to the NE striking left-lateral strike-slip faults on the southeastern margin of the Tibetan Plateau. The velocities are calculated in a region composed of three crustal blocks, each rotating clockwise about a point in the centre of the eastern margin of the block. A random number generator is then used to determine the points at which the velocity will be 'measured' by the hypothetical GPS campaign. The strain rate tensor and vorticity tensor are then calculated from these velocities using the method described in Section 2.3, except the continuous velocity field constructed using splines is not smoothed during these calculations. The results of three of these calculations (in which the velocity is sampled at different densities) are shown in Figs A1, A2 and A3. The calculations shown in the Figs are representative of calculations with other distributions of 'GPS sites'. In the calculation where the velocity is only known at 10 points (Fig. A1) there is only a slight variation in the strain rate to show the location of one of the faults, and the vorticity is relatively constant across the zone. Because of the competing influence of the rotation of the blocks and the deformation on the faults, the calculated vorticity has the same sign as, but is smaller in magnitude than would be expected simply because of the block rotations (− 2.5 × 10 −8 ). The calculation in which the velocity is known at 20 points (Fig. A2) shows more structure in both the strain rate and vorticity fields, but the faults are still not clearly visible and the vorticity is still lower than the actual vorticity in the block interiors. The calculation using 100 velocity vectors (Fig. A3) clearly shows the locations of the faults in the strain rate field, and the vorticity field shows the difference between the vorticity related to block rotation and that related to deformation on the faults. These calculations are obviously simplistic, but provide an indication as to how zones of rotating strike-slip faults would be expected to appear in GPS data, which is useful when examining the NE striking left-lateral strike-slip faults on the southeastern margin of the Tibetan Plateau (where the GPS station density is similar to that shown in Figs A1 and A2).

A P P E N D I X B : C O M PA R I S O N O F G E O D E T I C A N D S E I S M I C R AT E S O F S T R A I N
This appendix describes a comparison between the geodetic and seismic rates of strain for the fault systems on the southeastern margin of the Tibetan Plateau. Kostrov (1974) showed that the seismically released rate of strain in a volume (V ) can be estimated by summing the moment tensors of earthquakes in that volume which Figure A2. As Fig. A1, but the velocity has been sampled at 20 locations. occurred over a time T, that is,˙ i j = (1/2μV T ) M i j , where˙ i j is the strain rate tensor, μ is the shear modulus and M i j are the components of the moment tensors of the earthquakes. The following sections will compare the strain rate estimated from the geodetic data and that estimated from earthquakes for all of the areas of faulting which have GPS coverage. The earthquakes used to calculate the rate of moment release are taken from the CMT catalogue and the compilation of Holt et al. (1995) (which in places can be used to extend the record back to the early 1900s). The seismogenic thickness was taken to be 10 km throughout, based on the depths of earthquakes in the catalogue of Holt et al. (1995).

B1 E-W striking normal faults
The E-W striking normal faults at ∼29-31 • N accommodate northsouth extension. The largest component of the horizontal geodetic strain rate tensor can therefore be estimated by looking at northsouth gradients in southwards velocity, and is ∼2.5 × 10 −8 yr −1 . The seismic strain rate can be estimated by summing earthquake moment tensors from the catalogue of Holt et al. (1995) and the CMT catalogue. The oldest earthquake in the region of the E-W normal faults from either catalogue occurred in 1986, but the CMT catalogue should be complete since 1976 (for earthquakes large enough to make a significant contribution to the strain rate), so the time span used in the calculations is 1976 to the present. The calculated seismic strain rate corresponding to N-S extension is ∼6 × 10 −9 yr −1 , which represents ∼25 per cent of the geodetic strain rate. The difference between the seismic and geodetic strain rates is equivalent to an M w 7 earthquake every ∼30 yr or an M w 7.25 earthquake every ∼80 yr. Molnar & Deng (1984) describe an earthquake of M s 7.25 which occurred on 1948 May 25, in the area of the recent normal faulting earthquakes. The historical catalogue of Luo (1980) shows an event of magnitude 7.25 occurred in 1870, also in the area of recent normal faulting. The 1948 and 1870 earthquakes are large enough that, if they occurred on E-W normal faults, they would remove most or all of the difference between the estimates of seismic and geodetic strain rate, assuming that the seismicity recorded since 1976 is representative of the seismicity occurring over similar timeframes before the instrumental period in addition to the rarer large earthquakes. It therefore seems that the surface deformation in the area of E-W normal faulting is likely to be mostly or entirely achieved by seismic slip on faults.

B2 N-S striking normal faults
The main component of the horizontal strain rate tensor in the area of N-S normal faulting is ∼E-W extension, which can be estimated from a transect through the GPS data, and is ∼3 × 10 −8 yr −1 . Summing the moment tensors in the CMT catalogue and the catalogue of Holt et al. (1995) (in which the oldest earthquake in the region of the N-S normal faults occurred in 1966) provides an estimate of the seismic strain rate corresponding to E-W extension of ∼1.6 × 10 −8 yr −1 , which accounts for ∼55 per cent of the geodetic strain rate. The difference between the seismic and geodetic strain rates is equivalent to an M w 7 earthquake every ∼100 yr or an M w 7.5 earthquake every ∼600 yr. Luo (1980) shows that an M s 7 earthquake occurred near the town of Dali (Fig. 6) in 1925, and that in 1515 an event with a magnitude thought to be greater than 7.5 (based on historical records) occurred near the northern end of the north-south striking normal faults. If these two earthquakes both occurred on normal faults they could account for most or all of the difference between the seismic and geodetic strain rates, provided that the recently recorded seismicity is similar to that which occurred before the instrumental period in addition to the historically recorded large earthquakes.

B3 N-S left-lateral strike slip faults
The N-S striking left-lateral strike slip faults along the margin of undeforming south China have produced a number of earthquakes which have been recorded in the historical catalogue of Luo (1980). Shen et al. (2005) estimate the left-lateral slip rate along the margin of undeforming south China to be ∼7 mm yr −1 . The rate of moment release that would be expected along the length of the ∼700 km long left-lateral strike-slip fault system if the deformation was entirely seismic can be calculated asṀ 0 = μldu, where μ is the shear modulus, l is the length of the fault system, d is the seismogenic thickness, andu is the slip rate. If the seismogenic thickness is 10 km, the estimated moment rate is 1.5 × 10 18 N m yr −1 . The earthquakes in the catalogue of Luo (1980) are given surface wave magnitudes, which are estimated from historical records, and the global M s − M 0 relation of Ekström & Dziewonski (1988) is used to estimate the seismic moments of these events. The seismic moments of the 7 historically recorded earthquakes on this fault system can then be summed, and divided by the time since the occurrence of the oldest earthquake (1500 AD) to give an estimate of the moment rate of 5.2 × 10 18 N m yr −1 . There are obvious uncertainties involved with estimating the magnitudes of historical earthquakes which may in part explain why the rate of seismic moment release is greater than the rate estimated geodetically. Another possible reason why the seismic strain rate is higher than the geodetic strain rate, however, is the occurrence of a magnitude 8 earthquake on the fault system in 1833, which accounts for roughly half of the seismic moment release. There is only one magnitude 8 earthquake in the historical catalogue of southeastern margin of the Tibetan Plateau, and earthquakes of such magnitude are likely to have very long repeat times if they occur on faults with relatively low slip rates, so it is likely that the rate of moment release has been overestimated because our period of observation includes such a large, rare, event. Because of the numerous sources of potential error in the above calculations, all that can be concluded is that plentiful large earthquakes have occurred on the system of N-S left-lateral strike-slip faults in historical times, and that these earthquakes could account for the geodetically determined deformation rate, but the errors in the estimation are large.

B4 NE striking left-lateral strike-slip faults (SW of the Red River Fault)
Strike-slip earthquakes which have occurred to the SW of the Red River Fault are present in the catalogue of Holt et al. (1995) dating back to 1923. The summation of earthquake moment tensors from the catalogue of Holt et al. (1995) and the CMT catalogue suggests that the component of the seismic strain rate tensor corresponding to north-south gradients in northwards velocities has a magnitude of −2.4 × 10 −8 yr −1 (the negative sign indicates contraction). The component corresponding to east-west gradients of eastwards velocity has a magnitude of 2.4 × 10 −8 yr −1 , and the off-diagonal component of the horizontal strain rate tensor has a value of 1.5 × 10 −8 yr −1 . The geodetic strain rate in the area of the NE striking strike-slip faults is difficult to estimate, as the faults cross the southern boundary of the GPS network of Shen et al. (2005), so the velocity field is not as well constrained as in the areas discussed above. The ranges of the possible values of the components of the horizontal strain rate tensor are ∼ −1 to −2.5 × 10 −8 yr −1 for the component corresponding to north-south gradients in northwards velocities, ∼1-3 × 10 −8 yr −1 for the component corresponding to east-west gradients of eastwards velocity, and 0-1.5 × 10 −8 yr −1 for the off-diagonal component. Therefore, in the area of NE striking left-lateral strike-slip faults, the seismic strain rate appears to account for most or all of the geodetic strain rate.
The estimates of seismic strain rate presented in this paper should all be treated with caution. The southeastern margin of the Tibetan plateau contains many long faults which may rupture in large earthquakes and have long repeat times. Because the repeat times may be significantly longer than the seismic record, large earthquakes probably occur which are not taken into account in the above calculations, and so the estimates of the rate of moment release may be an underestimate. It is also possible that the instrumental period has seen an abnormally large number of earthquakes and the rate of moment release may have been overestimated. The incomplete and spatially variable nature of the GPS coverage results in uncertainties in the estimates of geodetic strain rates. Uncertainty is also introduced by the choice of seismogenic thickness of 10 km. If the value used for the seismogenic thickness is altered, the estimates of the rate of seismic strain will change in proportion to the size of the alteration. However, given these caveats, it appears that the geodetically measured strain on the southeastern margin of the Tibetan Plateau is likely to be accommodated mostly or entirely seismically in the areas where comparison of the two quantities has been possible. Any aseismic deformation which occurs is, therefore, not the dominant component of the strain, which means that earthquake focal mechanisms give a good indication of the style of strain in a given area.

A P P E N D I X C : D E S C R I P T I O N O F T H E N U M E R I C A L M O D E L
This appendix describes the numerical model used in Section 3. This model calculates the velocities and temperatures in a thermomechanically coupled fluid under the influence of gravity and velocity boundary conditions. The model will be described in two parts. First, the calculation of the velocity field for a given temperature distribution will be explained, and then the calculation of the temperature distribution.

C1 The velocity calculation
At low Reynolds number the equations for Stokes flow are satisfied, that is, where summation over repeated subscripts is assumed and σ ij are the elements of the stress tensor, x 1 , x 2 , and x 3 are the coordinate axes in Cartesian space and f i are the driving forces. The first assumption made is the 'hydrostatic approximation' in the vertical, so that the z-equation reduces to ∂ σ zz /∂ z = ρ g where ρ is density and g is the acceleration due to gravity. The deviatoric stress (τ ) is related to the strain rate (˙ ) as τ i j = 2η˙ i j where η is the viscosity and is given by where n is a constant known as the stress exponent,Ė is the second invariant of the strain rate tensor and the value of B depends, amongst other things, on temperature and grain size. n is equal to 1 if the deformation is by diffusion creep and 3 or more if the deformation is by dislocation creep. Following experimentally and theoretically derived mineral flow laws the temperature and grain size dependence of B for a given mineral deforming by a given creep mechanism is taken to be B = Cd − p exp [− (Q + PV )/RT] where C is a constant, d is the grain size, p is a constant known as the grain size exponent, Q is the activation energy, P is the pressure, V is the activation volume, R is the gas constant and T is the temperature. The value of p is 3 if the deformation is by diffusion creep and 0 if deformation is by dislocation creep. The second assumption which is made to simplify the solution of eq. (C1) is that the horizontal gradients of vertical velocity are neglected. This assumption simplifies the (1, 3) and (2, 3) components of the strain rate tensor (and their symmetrical opposites) and means that the x-component of eq. (C1) can be written where s is the surface elevation of the flow, and u = (u, v, w) is the velocity. To aid the numerical solution of this equation a co-ordinate transformation is undertaken. The model domain is transformed in the vertical so that the top and bottom of the model domain correspond to the top and bottom of the flow to be modelled. Although this transformation makes the co-ordinate system non-Cartesian and the equations themselves more complex, the numerical solution is made considerably easier. If the viscosity is known eq. (C3) and the other horizontal equation, once written in finite difference notation, form a system of linear equations for the horizontal components of velocity, which can be solved using the minimum-residual method of Saad & Malevsky (1995). However, the viscosity is not known.
Although the temperature dependence of the viscosity is known from experimental studies which have calculated the values of the parameters in B in eq. (C2), where deformation is by dislocation creep the viscosity also depends on strain rate. Therefore, an initial viscosity structure is assumed. A velocity field is then calculated using this viscosity structure. The strain rate can be calculated from this velocity field. The viscosity at each point can then be calculated taking strain rate into account. A new velocity field is then calculated using the updated viscosity structure. These new velocities are then compared to those of the previous iteration, and if the difference is above a certain threshold the process is repeated until the threshold is reached. Once the horizontal velocities are known the vertical velocity can be calculated by assuming incompressibility.

C2 The temperature calculation
The temperature structure for the high plateau which is input into the velocity calculation is estimated from a simple model (used by M c Kenzie & Priestley 2007) which takes into account vertical diffusion and radiogenic heating, that is, where T is temperature, t is time, C p is heat capacity, k is the thermal conductivity, and H is the rate of radiogenic heating. This equation is solved using the finite difference method with a Crank-Nicholson (joint implicit and explicit) scheme. Temperature-dependent thermal conductivity is used for the mantle lithosphere, given by the relation of M c Kenzie et al. (2005). The thermal conductivity of the crust is taken to be 2.5 W K −1 m −1 . The upper 75 per cent of the crust is assumed to have a radiogenic heat production of 2 μW m −3 , and the remainder of the crust of 0.4 μW m −3 . The base of the lithosphere has a constant temperature, which is equal to the isentropic temperature at that depth (with a potential temperature of 1315 • C). Because the crust is thick, the temperature increases due to radiogenic heating. In different model runs the length of time over which this heating has occurred has been varied. As described in the text, the lowlands around the plateau margins are assumed to be in thermal steady-state, and the method of M c Kenzie et al. (2005) is used to calculate the geotherms in these areas.