In this paper we present a global model (GSRM-1) of both horizontal velocities on the Earth's surface and horizontal strain rates for almost all deforming plate boundary zones. A model strain rate field is obtained jointly with a global velocity field in the process of solving for a global velocity gradient tensor field. In our model we perform a least-squares fit between model velocities and observed geodetic velocities, as well as between model strain rates and observed geological strain rates. Model velocities and strain rates are interpolated over a spherical Earth using bi-cubic Bessel splines. We include 3000 geodetic velocities from 50 different, mostly published, studies. Geological strain rates are obtained for central Asia only and they are inferred from Quaternary fault slip rates. For all areas where no geological information is included a priori constraints are placed on the style and direction (but not magnitude) of the model strain rate field. These constraints are taken from a seismic strain rate field inferred from (normalized) focal mechanisms of shallow earthquakes. We present a global solution of the second invariant of the model strain rate field as well as strain rate solutions for a few selected plate boundary zones. Generally, the strain rate tensor field is consistent with geological and seismological data. With the assumption of plate rigidity for all areas other than the plate boundary zones we also present relative angular velocities for most plate pairs. We find that in general there is a good agreement between the present-day plate motions we obtain and long-term plate motions, but a few significant differences exist. The rotation rates for the Indian, Arabian and Nubian plates relative to Eurasia are 30, 13 and 50 per cent slower than the NUVEL-1A estimate, respectively, and the rotation rate for the Nazca Plate relative to South America is 17 per cent slower. On the other hand, Caribbean–North America motion is 76 per cent faster than the NUVEL-1A estimate. While crustal blocks in the India–Eurasia collision zone move significantly and self-consistently with respect to bounding plates, only a very small motion is predicted between the Nubian and Somalian plates. By integrating plate boundary zone deformation with the traditional modelling of angular velocities of rigid plates we have obtained a model that has already been proven valuable in, for instance, redefining a no-net-rotation model of surface motions and by confirming a global correlation between seismicity rates and tectonic moment rates along subduction zones and within zones of continental deformation.
Plate tectonics has established itself as the most fundamental conceptual model in modern-day geology and geophysics by explaining most first-order geophysical observations on Earth in a comprehensive and self-consistent manner (e.g. Vine & Matthews 1963; Isacks et. al. 1968; Dewey & Bird 1970; Sclater & Franchateau 1970). Although successful as a theory, the concept of plate tectonics is built on the approximate assumption that plates are rigid and plate boundaries are narrow. Both assumptions are only partially correct (e.g. Gordon 1998). Plate boundaries are narrow (∼1–60 km) along the majority of oceanic ridges and transforms, but they appear much wider within most continental boundary zones and within some oceanic areas. The combined oceanic and continental diffuse plate boundary zones cover ∼15 per cent of the Earth's surface (Gordon & Stein 1992). Only recently, with the vast increase in geodetic observations of crustal motions, has it become possible to measure relative (horizontal) velocities and quantify crustal deformation within (diffuse) plate boundary zones with high accuracy (e.g. Stein 1993). As a result, our understanding of the behaviour of plate boundary zones has increased significantly and enables us, for example, to begin to distinguish between ‘block-like’ (e.g. Tapponnier et. al. 1982; Avouac & Tapponnier 1993; Peltzer & Saucier 1996) versus continuous deformation (e.g. England & McKenzie 1982; England & Houseman 1986; Molnar 1988). Ultimately, inferences on the fundamental distribution of plate boundary deformation as well as the quantification of strain rates and rotation rates (around a vertical axis) within the crust and lithosphere will allow us to improve upon our understanding of the dynamics that govern plate boundary deformation (e.g. Lamb 1994; Thatcher 1995; Jackson 2002). Furthermore, a complete quantification of the velocity gradient tensor field within plate boundary zones will be of great interest in seismic hazard studies (e.g. Ward 1994; Giardini 1999) and will, most likely, become an integrated part in future plate tectonic and plate motion models.
In this paper we present results of a global velocity gradient tensor field model associated with the accommodation of present-day crustal motions. We refer to this model as the global strain rate model, or GSRM-1. Because we solve for the velocity gradient field over almost the entire surface of the Earth, this is the first study to include plate boundary zones into one global model quantifying the complete present-day surface kinematics. Other global kinematic models that have recognized the importance of estimating the kinematics of plates and plate boundaries simultaneously (Drewes & Angermann 2001) are still limited by a relatively small number of geodetic velocities. Our model combines space geodetic data (i.e. global positioning system (GPS), very long baseline interferometry (VLBI) and Doppler orbitography and radiopositioning integrated by satellite (DORIS)) measured on rigid plates as well as within plate boundary zones. Some regional Quaternary fault slip rate data and seismic moment tensor information from shallow earthquakes are included as well. We show strain rate field solutions for several plate boundary zones and, using the assumption of plate rigidity for plates and blocks far away from known deformation zones, give present-day angular velocities for most plate pairs. We also present a new no-net-rotation model, which is a slight revision of the model recently presented by Kreemer & Holt (2001).
The methodology adopted to estimate a global strain rate and velocity model has been presented in Kreemer et. al. (2000a), but a brief overview is given here. We make the assumption that the lithosphere in plate boundary zones behaves as a continuum. This is a reasonable approximation when considering large-scale deformation for areas that have horizontal dimensions several times the thickness of the brittle elastic layer (e.g. England & McKenzie 1982). We adopt the methodology of Haines & Holt (1993) to estimate the horizontal velocity gradient tensor field on a sphere. A bi-cubic Bessel interpolation is used instead of polynomials, however, to expand a model rotation vector function that is obtained by a least-squares minimization between model and geodetic velocities, and model and geological strain rates (Holt & Haines 1995). Geological strain rates are obtained through a summation of Quaternary fault slip rates using a variant of Kostrov's (1974) formula (Shen-Tu et. al. 1999; Holt et. al. 2000a). This type of modelling has been presented in numerous papers concerning regional tectonics in zones of diffuse deformation; e.g. western United States (Shen-Tu et. al. 1999) and central Asia (Holt et. al. 2000a). A comprehensive overview of the methodology can be found in Haines et. al. (1998), Holt et. al. (2000b) and Beavan & Haines (2001).
Our model grid is continuous in longitudinal direction and covers the globe between 80°N and 80°S. Each grid area is 0.6°× 0.5° in dimension. Whether an area is considered to be deforming or not is based primarily on seismicity occurrence (Engdahl et. al. 1998) and the Harvard Centroid Moment Tensor (CMT) catalogue (Dziewonski et. al. 1981; Dziewonski et. al. 2000). All oceanic ridge and transform zones are allowed to deform. Within oceanic and continental regions of diffuse deformation, where seismicity rates are often low, the designation of boundaries between deforming and plate-like regions was often subjective. Therefore, the geometry of deforming regions in this model should be viewed as being approximate. About 24500 grid areas cover the Earth's deforming regions; all other areas are considered to be rigid. The rigid areas mimic 25 independent plates and blocks, including a number of relatively small entities such as the Rivera, Anatolian, Okhotsk, Caroline, Scotia, Sunda, Tarim, Amurian and South China plates and blocks, among others.
To accommodate the fact that not all plate boundary zones or areas within one single plate boundary zone strain with the same magnitude we constrain the magnitude of the a priori strain rate variance to vary globally. In order to do this all plate boundary zones are divided into 218 smaller areas. For each of these areas a value is assigned, depending on the expected magnitude of the strain rate for the area. For conformity the strain rate variance is chosen from a range of four values, depending on the expected strain rate (mainly based on seismicity occurrence); the highest a priori variance value is reserved for the zones that are expected to deform with the highest strain rate. Following this procedure we maintain some control on the distribution of expected strain rates; i.e. by absence of geodetic velocities and Quaternary fault slip rates that could provide constraints on the strain rate distribution, we can constrain some areas to deform at a higher rate than others. The a priori assignment of a heterogeneous strain rate variance distribution is therefore particularly important in obtaining strain rate distributions for areas such as Iran, East African Rift Valley, western Mediterranean, and central America that are consistent with the regional seismicity distribution.
For regions that are not densely sampled with geodetic observations, the interpolation of geodetic velocities can be highly non-unique in describing the regional strain rate field (Kreemer et. al. 2000b; Beavan & Haines 2001). However, the design of the strain rate covariance matrix can place a priori constraints on the style and direction of the model strain rate field. The constraint on the direction of the principal axes of the model strain rate field involves an uncertainty of ±10°. Information concerning the style and direction of expected strain rate is inferred from the principal axes of the seismic strain rate field, which is obtained through a Kostrov (1974) summation of seismic moment tensors in each grid area. Seismic moment tensors are used only from shallow events (≤40 km) in the Harvard CMT catalogue (1977 January–2000 December). For this purpose all events are weighted equally in the Kostrov summation. This approach is adopted here, except for areas where Quaternary fault slip rate data are used to infer geological strain rate estimates (currently only for central Asia). It should be noted that only the style and direction of the model strain rate field is constrained using the direction and style of the seismic strain rate field, not the magnitude. Also, constraints on the style and direction of the model strain rate field do not significantly affect the fit of the model velocities to the observed velocities. Moreover, including constraints on the style and direction of the model strain rate field results in a solution that is as consistent as possible with regional seismotectonics, while providing more stability in the model strain rate field from one grid area to another (Kreemer et. al. 2000b). More specifics on this methodology can be found in Kreemer et. al. (2000a) and Haines et. al. (1998).
One of the main advantages of the methodology used for this model is that an unlimited number of geodetic studies may be combined. The original reference frame of each individual study does not need to be adopted and, in fact, can be left undefined, a priori, in the inversion; the reference frame is solved upon fitting geodetic velocities to one self-consistent velocity gradient tensor field. That is, implicit in our procedure is the assumption that a single rigid-body rotation can be applied to velocity vectors from each individual study (one rotation vector for each study) such that the model velocity field provides a ‘best fit’ to the observed vectors that have been rotated into a single model frame of reference. We acknowledge that the assumption that the reference frames of various studies differ only by a single rigid-body rotation is an oversimplification. In addition to the rigid-body rotation there could be a corresponding translation to describe the reference frame difference (formally there is an additional seventh parameter to describe the reference frame difference, but this scaling parameter can be ignored here because we only use horizontal velocities). For small-scale networks this translation vector will be analogous in effect to a rigid-body rotation, and can be ignored. However, our approach to ignore the translational correction may lead to the possibility that we cannot satisfactory resolve the reference frame differences between large-scale networks. Except for one case (described in the next session), we do not observe significant residuals between velocities at colocated sites for larger networks after we have applied the rigid-body rotation. However, we do acknowledge that future modelling may need to incorporate the translational difference as well.
In this study we perform a least-squares fit to 3000 geodetic velocities from 50 different studies. Site locations are shown in Fig. 1 and references are given in Table 1. The majority of the geodetic data consist of GPS velocities. As part of the GPS velocities we include the GPSVEL (version 0.2) model (Lavallée et. al. 2001). GPSVEL is a global velocity solution determined for 175 global sites that are coordinated by the International GPS Service (IGS) as part of the IGS Densification Project (Zumberge & Liu 1995). The GPSVEL solution is obtained by a rigorous combination of weekly station coordinate estimates from seven global analysis centres and three regional associate analysis centres (Europe, Australia and South America) using a free-network approach (Davies & Blewitt 2000). Only stations with a time-series of at least 4.5 yr are analysed and the methodology of Zhang et. al. (1997) is used for error analysis. Geodetic velocities that are part of this solution will be referred to as GPSVEL velocities in the remainder of the paper. We also include the VLBI global velocity data set (Ma & Ryan 1998) as well as a subset of the global velocity solution obtained from the DORIS system (Crétaux et. al. 1998). We found a consistent global misfit between the DORIS velocities and the velocities obtained by VLBI and GPS. This discrepancy could be caused by the fact that we only applied a rigid-body rotation and not a translation to transform the DORIS velocities from their original reference frame into the model reference frame. We tested this possibility by inverting for both the rotation and translation vector using nine colocated sites of the DORIS and GPSVEL networks. We found indeed a significant translation vector. However, although the misfit between DORIS and GPSVEL became much smaller after we applied the obtained rotation and translation vectors to the DORIS reference frame, a significant misfit between DORIS and GPSVEL velocities was still present, which hints at an additional source behind the reference frame discrepancy. Therefore, for now, we have only included DORIS velocities measured at the Nubian (four sites) and Somalian Plate (one site) and at the Djibouti site location. We observe a good consistency between the DORIS velocities and other geodetic velocities when only this subset of the DORIS global data set is used. Most of the studies used in the model are published, but the most current data from some sources are only available on-line (United States Geological Survey (USGS); , and the Southern California Earthquake Center (SCEC); ). Large geodetic data sets in Asia and for the entire globe have been published very recently by Wang et. al. (2001) and Sella et. al. (2002), respectively, but these studies have not (yet) been included in the results presented here (with the exception of a N–S transect in Tibet that is part of the Wang et. al. (2001) data set and that was given to us by J. Freymueller before publication of that paper).
Formal errors in geodetic velocities often underestimate the ‘true’ uncertainty (Johnson & Agnew 1995). This underestimation is generally owing to the fact that only ‘white noise’ is considered in the data analysis while effects of time-correlated noise are often not included (Mao et. al. 1999). To acknowledge this fact and to obtain a fairly even distribution of velocity uncertainties between the different data sets we have multiplied standard errors in velocity for some of the studies by a factor of 2 or 3 (see Table 1). Alternatively, following Argus & Gordon (1996) a time-dependent velocity variance σnew is defined, where and Δt is in years, such that for sites on the same plate the requirement of plate rigidity is met. We apply this approach to the VLBI velocities (Ma & Ryan 1998) for which we had a good knowledge of the length of the time-series for each site. We use a value of C= 5.5 mm yr−1 as suggested and used by Argus & Gordon (1996) and Larson et. al. (1997).
To date the only Quaternary fault slip rate data used in the global model are from central Asia (England & Molnar 1997; Holt et al. 2000a and references therein). Although an extensive set of Quaternary fault slip rate data is available for the western United States (e.g. Peterson & Wesnousky 1994; Peterson et. al. 1996; Shen-Tu et. al. 1999) these data are not used in this model because of the large amount of available geodetic velocities in this region. That is, in the western United States the geodetic data alone will suffice to accurately describe ongoing present-day deformation at a scale of resolution we are seeking within a global framework.
Owing to the absence of a global data set of slip rate estimates, it is currently not possible to obtain a global strain rate field that reflects the long-term deformation pattern in diffuse plate boundary zones. On the other hand, for some plate boundary zones (e.g. western United States and the eastern Mediterranean) there are sufficient geodetic velocity measurements to constrain the strain rate field associated with ongoing motions reasonably accurately. In most cases, however, strain rate estimates from geodetic velocities will be more laterally distributed than long-term geological strain rates (e.g. Shen-Tu et. al. 1999), due to effects of elastic loading of the lithosphere. Nevertheless, over sufficient length-scales, of the order of the width of plate boundary zones, geodetic velocities are often indistinguishable from relative plate velocities (e.g. Stein 1993; Shen-Tu et. al. 1998; Kreemer et. al. 2000b). We perform a smoothing of the geological strain rates over an appropriate length-scale ∼50–100 km) such that geodetic and geological strain rates compliment each other in reflecting the ongoing horizontal deformation in areas where the lateral dimension is several times larger than the elastic thickness. Geological strain rates add especially useful kinematic constraints in areas where geodetic data are limited or absent. However, as a consequence of combining these different data sets, the resolution for the global strain rate model has a lower limit of ∼50–100 km. Moreover, we do not attempt to solve for slip on individual faults. To avoid modelling of temporal changes in strain rate due to co-seismic and post-seismic relaxation of the crust, we have not included any geodetic velocity that reflects a significant component of co- and post-seismic deformation.
Global Strain Rate Model
In determining a global strain rate model we seek a best fit between model velocities and 3000 observed velocities (the reduced chi-squared is 1.06) and model and geological strain rates inferred from Quaternary fault slip rates. Table 1 lists for each study the obtained angular velocity with which the original velocity vectors are rotated into the model reference frame (Pacific).
In Fig. 2 we present the second invariant of the global model strain rate field, where the second invariant is defined as , where and > are the horizontal components of the strain rate tensor. Fig. 2 also indicates the assigned geometries of the plate boundary zones. It is evident that the highest strain rates reflect spreading along mid-oceanic ridges. Other high strain rate areas are the subduction zones along the Pacific Rim. Significantly high strain rates are apparent in diffuse boundary zones, such as central Asia and the western United States, and deformation in these zones is often heterogeneously distributed. We discuss the strain rate fields in the broad plate boundary zones of central Asia, Indian Ocean and the Aegean in more detail in the sections below describing the appropriate relative plate motions.
Plate Motions and Plate Boundary Deformation
The motions of most major plates are well constrained by a model velocity field obtained in a least-squares fit to the observed geodetic velocities for those plates. Fig. 3 shows the locations of the sites used in this study that are assumed to be located on a stable plate/block. Table 2 shows the reduced chi-squared for each plate for which the motion has been measured at a minimum of two positions. The plates that are defined in the model, but for which motions are not directly measured using geodetic velocities, are the Caroline, Capricorn, Cocos, Juan de Fuca, Rivera and Scotia plates. Consequently, plate motion results for these plates are currently uncertain and not presented. Table 3 displays the angular velocities of all (constrained) plates with respect to the Pacific Plate, whereas angular velocity vectors describing relative motion between selected plate pairs are shown in Table 4. For comparison, Table 4 also contains relative rotation vectors from the NUVEL-1A plate motion model (DeMets et. al. 1994a) and from Sella et. al. (2002) who have presented the most extensive set of geodetically derived plate motions to date.
The motion between the Pacific (PA) and North American (NA) plates has been widely studied (e.g. DeMets et. al. 1987; Argus & Gordon 1990; Ward 1990; Antonelis et. al. 1999; DeMets & Dixon 1999), because of the implications of this motion for plate boundary deformation and seismic hazard mitigation in California and Alaska (e.g. Minster & Jordan 1987; Humphreys & Weldon 1994; Shen-Tu et. al. 1998). The angular velocity found here describing PA–NA motion (50.8°S, 102.2°E, 0.77°Myr−1) is close to, although slightly faster than, the NUVEL-1A result. However, our Euler vector is significantly more western and at a lower rate than the geodetic estimates obtained by Larson et. al. (1997) and Crétaux et. al. (1998). Kreemer et. al. (2000a) speculated that the different pole location found by others could be (partly) explained by the inclusion of the velocity at Fairbanks, Alaska, in the estimation of North America motion, which results in a PA–NA Euler pole more easterly and with a higher rotation rate. Kreemer et. al. (2000a) found that a significantly better fit to observed North American velocities is obtained when Fairbanks is considered to lie within the PA–NA plate boundary zone, which is what we have assumed here. Nome in western Alaska (Ma & Ryan 1998) does not appear to be part of the North American Plate either. A similar conclusion was found by Argus & Gordon (1996). The residual velocity vector we find at Nome with respect to North America is 4.4 mm yr−1 towards S43°E. This difference is significant at the 95 per cent confidence level. When the observed velocity at Nome is taken out of the model, a better fit to the observed velocities in the eastern part of the North American Plate is obtained, but the PA–NA rotation vector remains unaffected. (In fact, the PA–NA angular velocity changes insignificantly when all North American stations west of 104°W are excluded.) Nome is located near the Bering Strait where seismicity levels are relatively high (e.g. Biswas et. al. 1986). Seismicity here possibly delineates the northeastern boundary of the proposed independent ‘Bering block’ (Mackey et. al. 1997). Initial results using the velocity obtained by Kogan et. al. (2000) at Bilibino in eastern Siberia (750–1000 km east of the (diffuse) NA–Eurasia (EU) plate boundary zone Chapman & Solomon 1976) led to a significant residual velocity for this site as well. However, given a new revised 5 yr velocity estimate at Bilibino (M. Kogan, pers. comm. 2002), this station now appears to be part of rigid NA.
PA–NA motion within our model velocity field (Fig. 4) has a magnitude of 53.3 ± 1.0 mm yr−1 towards N56°W ±1° at a point on the Pacific Plate just off-shore of Cabo San Lucas at the southern tip of Baja California (22°N, 111°W) (Baja California itself is assumed to be part of the PA–NA plate boundary), 48.8 ± 1.0 mm yr−1 in a direction of N44°W ± 1° off-shore of southern California (32°N, 118°W), 48.2 ± 1.0 mm yr−1 towards N37°W ±1° near San Francisco (37.5°N, 123°W) and 59.0 ± 1.0 mm yr−1 towards N28°W ± 1° south of Kodiak Island (54°N, 153°W) (here and hereafter, uncertainties represent 95 per cent confidence limits). Our results predict a PA–NA motion along southern California and Baja California that is about 2–3 mm yr−1 faster than NUVEL-1A. This result is consistent with the GPS results found by others (DeMets & Dixon 1999; Beavan et. al. 2002; Sella et. al. 2002).
Along the PA–NA plate boundary we find that model velocities are directed anticlockwise from the NUVEL-1A PA–NA motion direction, with a maximum discrepancy of almost 3° in central and southern California (Fig. 4). This result is consistent with the 2°–3° anticlockwise discrepancy in direction found by Shen-Tu et. al. (1999) and Larson et. al. (1997) based on GPS measurements. Based on the summation of Quaternary fault slip rates over the entire PA–NA plate boundary zone Humphreys & Weldon (1994) and Shen-Tu et. al. (1999) found that the Pacific Plate moves up to 6° anticlockwise from the NUVEL-1A estimate.
Caribbean–North America and Caribbean–South America
We have assumed that San Andres Island, Aves Island and St Croix (US Virgin Islands) are located on the stable Caribbean (CA) Plate (Fig. 5). Site velocities at all three locations are taken from Weber et. al. (2001). The site at St Croix is also measured as part of the VLBI global network (Ma & Ryan 1998), as part of the GPSVEL solution, and as part of a regional GPS network in Puerto Rico (Jansma et. al. 2000). There is an additional measured velocity at Aves Island by Pérez et. al. (2001). Using the least-squares fitting procedure we find a reduced chi-squared of 0.50 (Table 2) for the fit between the seven geodetic velocities and the predicted model velocities (Fig. 5). No apparent evidence is found for possible intraplate deformation as was inferred by Leroy & Mauffret (1996). From our inversion we find that CA–NA motion is best described by an Euler pole at 75.3°N, 136.5°W at a rate of 0.18°Myr−1, close to the result found by DeMets et. al. (2000). With respect to the South-American (SA) Plate the Caribbean Plate moves about a pole at 53.8°N, 71.0°W at a rate of 0.26°Myr−1, consistent with the result by Weber et. al. (2001). The relatively high inferred velocity of the Caribbean Plate compared with the NUVEL-1A estimate (Fig. 5) is consistent with earlier geodetic studies (Dixon et. al. 1998; MacMillan & Ma 1999; DeMets et. al. 2000; Weber et. al. 2001; Sella et. al. 2002). Following Dixon et. al. (1998) we believe that the discrepancy in rate can be explained by a systematic error in the NUVEL-1A estimate of Caribbean Plate motion. The cause of the error is probably the inclusion of some earthquake slip-vectors (especially along the Middle America trench) that may not reflect actual relative plate motion direction (Deng & Sykes 1995). This idea is supported in a study by DeMets (1993) in which he found a significant increase in the rate of CA–NA motion when earthquake data are omitted from the NUVEL-1 model (DeMets et. al. 1990).
South America–North America
Although the separateness of the North and South American plates has been recognized for a long time (e.g. Ball & Harrison 1970), estimates for the relative motion of this plate pair were for a long time solely based on plate circuit closures (e.g. Minster & Jordan 1978; DeMets et. al. 1990). Direct observations using GPS have confirmed that North and South America are distinct plates (Argus & Heflin 1995; Larson et. al. 1997; Dixon & Mao 1997; Crétaux et. al. 1998; Sella et. al. 2002). Geodetic velocities in our model are fitted well on both the North and South American plates (reduced chi-squared is 1.03 and 0.73, respectively, Table 2), including at stations on the South American Plate as far east as Ascension Island (GPSVEL) and as far south as Punta Arenas (Bevis et. al. 1999) and Rio Grande, Argentina (GPSVEL). We obtain an SA–NA angular velocity vector (13.9°S, 126.2°E, 0.15°Myr−1) close to the estimate by Dixon & Mao (1997) and Sella et. al. (2002), which are the only two studies that include a fair number of stations on both plates.
In our model we assume the diffuse SA–NA plate boundary zone between the Lesser Antilles and the Mid-Atlantic Ridge to be ∼325–500 km wide and to include a zone of anomalous ridge and trough features that may be related to SA–NA boundary deformation (Müller & Smith 1993). The boundary zone also includes the 50–20 fracture zone, which has been suggested to be the SA–NA plate boundary (Roest & Collette 1986). The observed seismicity in the boundary zone is low, but significantly higher in comparison with other areas in the Atlantic Ocean (e.g. Bergman & Solomon 1980; Bergman 1986; Wysession et. al. 1995), and for three events a moment tensor solution is available (e.g. Stein et. al. 1982). These three events are located in the western part of the boundary zone and the moment tensors are consistent in style; all three events show roughly N–S-trending P axes, and events are characterized by a combination of thrusting and right-lateral strike-slip faulting along WNW-trending fault planes, subparallel to regional fracture zones (Stein et. al. 1982; Bergman 1986). Our model results places the NA–SA rotation pole within the assumed NA–SA plate boundary zone. Although predicted relative motions in the plate boundary zone have a relatively high uncertainty due to the proximity of the rotation pole, predicted NA–SA relative motions are consistent in direction with the moment tensor solutions, albeit the total relative motion does not exceed 2 mm yr−1. However, due to the lack of additional moment tensor solutions for the (very few) events in the boundary zone and, again, because the Euler pole is very close to the boundary zone, the exact nature of the seismotectonics in this area remains uncertain.
Nazca-South America and Nazca-Pacific
From plate reconstructions (Pardo-Casas & Molnar 1987) and new seafloor data (Somoza 1998), it has been suggested that Nazca (NZ) Plate motion has slowed down gradually since 20–25 Ma. Recently, the decelerating of the Nazca Plate has been confirmed by geodetic measurements (Larson et. al. 1997; Angermann et. al. 1999; Norabuena et. al. 1999). Our angular velocity vectors describing NZ–SA (53.1°N, 92.1°W, 0.61°Myr−1) and NZ–PA (55.2°N, 89.6°W, 1.27°Myr−1) relative motion (Table 4), confirm slower relative motions compared with the NUVEL-1A model as was noted before by others (Angermann et al. 1999; Norabuena et al. 1999; Kreemer et al. 2000; Sella et al. 2002). The NZ–SA Euler pole yields convergence directions along the South American margin that are remarkably indistinguishable from the predicted NUVEL-1A motion directions. Convergence rates, on the other hand, range from 58 ± 2 mm yr−1 in Ecuador to 67 ± 2 mm yr−1 in Central Chile, which is, respectively, 18 and 16 per cent slower than the NUVEL-1A predicted convergence rates. Similarly, we find that predicted spreading rates along the East Pacific Rise between the Nazca and Pacific plates are indistinguishable in direction from the NUVEL-1A model, but ∼9 mm yr−1, or ∼7 per cent, slower.
Because the NZ–PA relative motion is very well constrained in the NUVEL-1A model, it is unlikely that the discrepancy between the geological and geodetic speed of the Nazca Plate is the result of a systematic error in the NUVEL-1A model. Instead, the relatively low present-day speed of the Nazca Plate compared with its speed over the last 3 Myr hints at a true deceleration of the plate (Norabuena et. al. 1999). Norabuena et. al. (1999) suggested that slowing down of the Nazca Plate may be the result of increased drag between the subducting Nazca and overriding South American Plate caused by thickening of South America's leading edge owing to growth of the Andes, resulting in ‘flat-slab’ subduction. Because ‘flat-slab’ subduction is not observed uniformly along the NZ–SA boundary (e.g. Stauder 1975; Barazangi & Isacks 1976) a more likely mechanism for ‘flat-slab’ subduction is the subduction of buoyant ‘aseismic’ ridges (Vogt 1973; Vogt et. al. 1976; Kelleher & McCann 1976). The Nazca Plate seems somewhat outstanding in this respect since it is floored by several ‘aseismic’ ridges that are all being subducted beneath South America. Pilger (1981) pointed out the correlation between zones of low-angle subduction and places where the Nazca and Juan Fernández ridges subduct underneath southern Peru and central Chile, respectively. Gutscher et. al. (1999) determined a similar correlation for Ecuador where the ‘aseismic’ Carnegie Ridge subducts. It has been inferred that subduction of the Carnegie, Nazca and Juan Fernández ridges all initiated 8–9 Myr (Cande 1985; Von Huene et. al. 1997; Gutscher et. al. 1999), which correlates well with the estimated 31 mm yr−1 decrease in NZ–SA convergence rate between 10.8 and 4.9 Ma (Somoza 1998). Similarly, the ∼14 mm yr−1 decrease in convergence rate since 3 Ma, based on a comparison of the NUVEL-1A estimated rate with the GPS measurements (Norabuena et. al. 1999; Angermann et. al. 1999) and our result, could be a reflection of the continuation of the subduction process of these aseismic ridges. In addition, the repeated subduction of ridge crest segments of the active Chile ridge underneath southern Chile since 14 Ma (Cande & Leslie 1986; Ramos & Kay 1992) could have also contributed to the deceleration of the NZ–SA plate motion.
In the NUVEL-1A model the Australia (AU)–PA plate pair is another example for which the relative motion is inferred indirectly through plate circuit closure due to lack of useful geological data along its boundary. A good knowledge of relative motion along the AU–PA boundary is, however, of great interest in studying some of the segments along this boundary; the Tonga–Kermadec subduction zone that is the source of the majority of (deep) global seismicity, the Alpine fault in New Zealand, and the Macquarie ridge along which one of the largest earthquakes (MW= 8.2) in the last couple of decades occurred. Until recently, studies of the PA–AU relative plate motion were hindered by a limited distribution of sites on the Australian and Pacific plates (e.g. Smith et. al. 1994; Argus & Heflin 1995; Larson et. al. 1997; Crétaux et. al. 1998). Fortunately, recent studies by Sella et. al. (2002) and Beavan et. al. (2002), as well as this study, have had access to a much wider distribution of sites on these plates (Fig. 3). Beavan et. al. (2002) found that velocities of stations 5508 and OUSD in eastern South Island, New Zealand, do not represent Pacific motion. We find that the velocities of these two stations, that we have assumed to be located on the Pacific Plate but that we have taken from the Beavan & Haines (2001) study instead of Beavan et. al. 's (2002) paper, are consistent with rigid plate motion. We obtain an angular velocity vector for PA–AU motion of 62.2°S, 175.6°W, 1.06°Myr−1 (Table 4). Our result is within the 95 per cent confidence limit of the NUVEL-1A estimate. This finding is somewhat surprising, because there have been indications that the present-day PA–AU Euler pole is located about 5° south of the NUVEL-1A pole (Larson et. al. 1997; Crétaux et. al. 1998), consistent with the result from plate reconstructions that the pole location has moved southward since at least 49 Ma (Molnar et. al. 1975; Weissel et. al. 1977; Stock & Molnar 1982; Sutherland 1995). However, uncertainties in stage poles have been large and generally overlap for different stages. Furthermore, the GPS-derived Euler pole by Larson et. al. (1997) has been suggested to be incorrect (Beavan et. al. 1999) based on discrepancies between predicted motions of New Zealand's South Island and locally observed GPS velocities. Our result suggests that the PA–AU angular velocity has changed insignificantly since at least 3 Ma. A similar conclusion was reached by Beavan et. al. (2002). Nevertheless, because of the proximity of the PA–AU pole to the PA–AU plate boundary zone, very small variations in the pole location can result in significant variations in predicted velocities along the boundary; we find that our model predicts velocities from Tonga to Macquarie that are 2–5 mm yr−1 faster than the NUVEL-1A prediction. The azimuths of the predicted slip vectors along the PA–AU plate boundary zone are within the 95 per cent confidence level of slip-vectors predicted by the NUVEL-1A model, with the only exception being for the most southern part of the Macquarie Ridge Complex where our model predicts a PA–AU velocity of up to 9° more westward than the NUVEL-1A direction.
Estimating the present-day motion of the Indian (IN) Plate is of great importance for continental deformation studies in central and southeast Asia (e.g. England & Houseman 1986; Holt et. al. 1995; England & Molnar 1997). Because IN–EU relative plate motion cannot be inferred directly from conventional plate motion data, space geodetic velocities are of special interest here. The southernmost part of the Indian subcontinent is particularly well covered by stations (Paul et. al. 2001). Besides a number of velocity estimates at Bangalore from a series of studies, our model also includes velocity estimates at Bhopal and Delhi (JNUC) (Paul et. al. 2001) and at Biratnagar (Larson et. al. 1999) (Fig. 6a). Biratnagar is located ∼100 km south of the active Himalayan thrust belt, but we found no evidence that the velocity at this station contains any significant component of elastic deformation.
We find that IN–EU model velocity is best described by an Euler pole at 25.9°N and 15.0°E with a rotation rate of 0.36°Myr−1. The pole location is within 95 per cent confidence limits of both the NUVEL-1A result (Table 4) and the pole location by Paul et. al. (2001) (25.6°N, 11.1°E, 0.44°Myr−1), but the expected rotation rate is, respectively, 0.15 and 0.08°Myr−1 slower. Our result is consistent with the findings of Sella et. al. (2002) and with previous modelling results using GPS and fault slip rates in central and southeast Asia (Holt et. al. 2000a; Holt et. al. 2000b). Gordon et. al. (1999) investigated a possible systematic error in NUVEL-1A India motion; they re-analysed spreading rate and transform fault data while redefining plate closure circuits and embedding a divided African Plate in their modelling. Gordon et al. (1999) found an angular velocity vector that yields a velocity at Bangalore with respect to Eurasia of 37 mm yr−1 with an azimuth of N31°E, which suggests that India's motion has not slowed down since 3 Ma, but merely that the NUVEL-1A IN–EU estimate is erroneous. Our model velocity at Bangalore with respect to Eurasia is 34.1 ± 1.2 mm yr−1 towards N23°E ± 2°. This model velocity is consistent with results by Holt et al. (2000a), Wang et. al. (2001) and most of the six used GPS vectors at this site (with the exception of the estimate by Shen et. al. 2000). However, in several cases the GPS velocities at Bangalore are only consistent with the model velocity after the GPS velocities have been rotated from their original Eurasian reference frame into the best-fitting model Eurasian reference frame presented here. We find that the discrepancy between our model Eurasian reference frame and the Eurasian reference frame defined by others exists mainly for studies that used the NUVEL-1A no-net-rotation (NNR) model (Argus & Gordon 1991; DeMets et. al. 1994a) to convert from an ITRF to a tectonic reference frame. The current NNR model, however, does not adequately describe the motion of most plates in a present-day NNR reference frame (Kreemer & Holt 2001). We believe that this is also the reason why the IN–EU rotation rate estimate by Paul et. al. (2001) is relatively high, yielding a velocity of 44 mm yr−1 at Bangalore. At Biratnagar, in northern India (26.48°N, 87.26°E), we infer a convergence rate of 35.5 ± 1.3 mm yr−1, and a convergence direction of N17°E ± 2°. Our result implies an IN–EU convergence rate that is 14–15 mm yr−1 lower than predicted from the NUVEL-1A model.
Strain rate and velocity field in central Asia
Central Asia is the only region in the current global model in which we have attempted to not only fit geodetic velocities but also geological strain rates inferred from Quaternary fault slip rates. Fig. 6(a) shows model and observed velocities with respect to model Eurasia. Recently, Wang et. al. (2001) combined many of the same regional studies in China that we have incorporated in this study (their study, however, combined raw data, whereas we combine velocity vectors) and their solution is very similar to the magnitude and directions of the geodetic vectors shown in Fig. 6(a). Fig. 6(b) contains the principal axes of the model strain rate field (black) and the geological strain rate field (white). Although we obtain a solution with a resolution of 0.5 × 0.62, for illustration purposes we show the strain rate field only as averages for larger grid areas (Fig. 6b). These larger grid areas are equal to the areas used by Holt et. al. (2000a). Our strain rate tensor field within central Asia is, in general, similar to the result from the regional study by Holt et. al. (2000a), but differs at two locations as a result of the use of new geodetic data (Chen et. al. 2000; Shen et. al. 2000, 2001; Bendick et. al. 2000; Wang et. al. 2001). Model shear rates along the Sagaing fault are higher in our result compared with that of Holt et. al. (2000a). These high shear strain rates are a result of the SSE motion of southwestern China with respect to Siberia (Chen et. al. 2000) (Fig. 6a). The other new result inferred from our model is the indication that a NNE-trending shear zone exists in north central Tibet (between 85°E and 90°E), subparallel to and south of the Altyn Tagh fault. Recent GPS results (Bendick et. al. 2000; Shen et. al. 2001) have suggested that the slip rate along the Altyn Tagh fault is roughly 9 mm yr−1, significantly lower than geological estimates of 20–30 mm yr−1 (e.g. Peltzer et. al. 1989). The total shear strain rates required to accommodate the 25 mm yr−1 of eastward motion of eastern Tibet with respect to the Tarim Basin (Fig. 6a), as indicated by GPS estimates (Heki 1996; Chen et. al. 2000; Shen et. al. 2000, 2001), significantly exceeds the geodetically observed slip rate, and hence shear strain rates, on the Altyn Tagh fault. The amount of additional shear strain required to accommodate the eastward motion of eastern Tibet with respect to the Tarim Basin seems to localize south of the Altyn Tagh fault on the Kun Lun fault and parallel to the westward extension of the Kun Lun fault. It is interesting to note that the epicentres of the 1997 MW= 7.6 Manyi and 2001 MW= 7.8 Kokoxilli earthquakes are roughly located in this zone of relatively large model shear strain rates. The moment tensor of these events are consistent with the left-lateral predicted shear strain rate direction (the Manyi event ruptured a 170 km long fault with a strike of N76°E Peltzer et. al. 1999).
The separateness of the Australian and Indian plates and the diffuse nature of their mutual boundary have been established for some time (e.g. Stein & Gordon 1984; Wiens et. al. 1985; Gordon et. al. 1990; Royer & Chang 1991). An estimate of the relative motion between these two plates is of significant importance in studies quantifying the deformation within the Indian Ocean and understanding the relative high level of seismicity that occurs there (e.g. Gordon et. al. 1990; Tinnon et. al. 1995). We obtain an AU–IN Euler pole at 5.2°S and 70.5°E with a rotation rate of 0.40°Myr−1 (Table 4). The pole location is within 95 per cent confidence limits of the NUVEL-1A, Gordon et. al. (1990) and DeMets et. al. (1994b) Euler pole estimates, all of which are based on conventional plate motion data and plate circuit closure requirements. However, the rotation rate found here, as well as by Sella et. al. (2002), is 0.09–0.11°Myr−1 faster than these prior estimates.
Strain rate field in the Indian Ocean and the Andaman Sea
For this study we assume a geometry of the IN–AU–Capricorn (CP) plate boundary zone similar to Royer & Gordon (1997), comprising the Chagos Ridge, Central Indian Basin, Ninetyeast Ridge, southern Bay of Bengal and Wharton Basin (Fig. 7a). The model strain rate field (Fig. 7a) indicates ∼3 mm yr−1 of roughly N–S extension near the Chagos Ridge between the Capricorn and Indian plates. In the Central Indian Basin we see pure compression associated with ∼5 mm yr−1 of shortening and the compression direction rotates from NNE–SSW in the western part of the basin to more NNW–SSE in the eastern basin. In the northeastern Indian Ocean, comprising the southern Bay of Bengal, northern Ninetyeast Ridge, and Wharton Basin, the compression direction is NW–SE to NNW–SSE and the estimated relative motion between the Indian and Australian plates is ∼11 mm yr−1. Our inferred AU–IN convergence rate of ∼11 mm yr−1 for the northeastern Indian Ocean is remarkably consistent with the total convergence of 125 ± 28 km since 11 Ma as was inferred by Royer & Gordon (1997) based on plate reconstructions in the Indian Ocean.
Our strain rate results are in general quite similar to expected extension and shortening directions discussed by Royer & Gordon (1997), and consistent with shortening directions inferred from seismic profiles (e.g. Bull & Scrutton 1992; Chamot-Rooke et. al. 1993; Van Orman et. al. 1995), and earthquake P–T axes (Fig. 7a). The consistency between the principal strain rate axes and earthquake focal mechanisms is partly imposed, because of the constraints from seismic moment tensors that are applied, a priori, to the style and direction of the model strain rate field. Tinnon et. al. (1995) derived a seismic strain rate field from an inversion of earthquake focal mechanisms. Although Tinnon et. al. (1995) used a different plate boundary geometry and did not define a Capricorn Plate, they obtained a solution quite similar in style and direction to our solution, with the exception of the Wharton Basin where Tinnon et. al. (1995) predicted more N–S, instead of more NW–SE to NNW–SSE, shortening. The result presented here for the Wharton Basin is more consistent with observed NE–SW-trending undulations in the gravity field, which are thought to have resulted from NW–SE-directed compression (Weissel et. al. 1980; McAdoo & Sandwell 1985; Haxby 1987). Also, the pure strike-slip earthquakes in the Wharton Basin and southern Ninetyeast Ridge region are consistent with the expected NW–SE compressional strain rates, indicating reactivation of N–S-trending fracture zones (Bull & Scrutton 1992; Deplus et. al. 1998). The occurrence of the 2000 June 18, MW= 7.8 strike-slip event in the Wharton Basin near Cocos Island (Robinson et. al. 2001; Abercrombie et. al. 2003) suggests that reactivation of old N–S fracture zones may occur as far east as 97°E in order to accommodate the NW–SE shortening between India and Australia.
We also present the model strain rate field for the Andaman Sea region (Fig. 7b). Our result resolves the partitioning of pure east–west compression along the Burma Trench and extension-dominated strain rates in the Andaman backarc. The predicted direction of extension, as it occurs along the spreading segments, is consistent with earthquake slip-vectors ∼N30°W). However, a significant component of right-lateral shear is also present and is in agreement with seismotectonics indicated by the focal mechanisms (Fig. 7b).
Arabia–Eurasia and Anatolia–Eurasia
The Arabian (AR) Plate motion is defined by a fit to the observed velocities at the GPSVEL station in Bahrain and station KIZI in southern Turkey, located about 100 km south of the Bitlis suture zone (McClusky et. al. 2000). The model angular velocity describing AR–EU motion predicts 0.44°Myr−1 about a pole at 26.2°N and 20.4°E (Table 4), which is significantly different from the NUVEL-1A estimate. It is evident that the plate moves at a lower speed with respect to Eurasia than the NUVEL-1A predicted rate; we find an AR–EU velocity of 22.1 mm yr−1 at Bahrein, whereas NUVEL-1A predicts 30.5 mm yr−1. This implies that AR–EU convergence at Bahrein is 27 per cent slower than the NUVEL-1A prediction.
We find an angular velocity vector describing Anatolia (AT)–EU relative motion (32.0°N, 33.4°E, 1.35°Myr−1) that is, although at a somewhat higher rotation rate, insignificantly different from McClusky et. al. 's (2000) estimate based on a similar set of GPS velocities. If all AT–EU relative motion is accommodated along the North Anatolian fault, slip rates can be inferred from the obtained AT–EU angular velocity. We find a slip rate for this fault of about 22–23 mm yr−1, equal to the slip rate inferred by McClusky et. al. (2000).
Strain rate field in the Aegean
The greater Aegean area is kinematically complex and seismically very active. In order to better understand the driving forces that are at the source of the tectonic complexity and to mitigate seismic hazard in the many populous areas in the region, a well-defined strain rate field is essential. Fig. 8(a) shows the principal axes of the model strain rate field for the Aegean area. The strain rate tensors are plotted as averages for each model grid area, indicating the true resolution of the global model. The strain rate field is similar to previous strain rate tensor calculations using similar sets of GPS velocities (Kahle et. al. 1999, 2000), which are fitted very well by the model velocities (Fig. 8b). As discussed in detail by Kahle et. al. (1999, 2000) many of the characteristics in the strain rate field solution are in good agreement with regional geology (e.g. grabens, troughs, and major strike-slip fault zones) and seismicity distribution (Fig. 8a).
Block motions in central and southeast Asia relative to Eurasia
It has been hypothesized that the northward motion of India with respect to Eurasia has caused the eastward extrusion of crustal material in east and southeast Asia (Molnar & Tapponnier 1975). Geological observations (e.g. Armijo et. al. 1989; Peltzer et. al. 1989; England & Molnar 1997) as well as modelling of earthquake moment tensors in central and east Asia (Holt et. al. 1995) have confirmed the eastward motion of eastern Asia with respect to Eurasia. The model presented here, incorporating both Quaternary fault slip rates and geodetic velocities in the region, builds on the models presented by Holt et. al. (2000a,b). Here, we present the motions of the Sunda, South China and Amurian blocks relative to Eurasia (Table 4, Fig. 9). The results for these blocks are a likely improvement over the results by Holt et. al. (2000a), due to an improved model Eurasia reference frame and the inclusion of more GPS vectors.
The estimation of the Amurian (AM) block angular velocity is particularly difficult because of the expected low rate and because of the poor coverage of GPS sites on this block. For most GPS sites assumed to be on the Amurian block it is in fact arguable whether or not they are located on the stable portion of the block. In the current plate geometry seven GPS velocities are assumed to have been measured on the Amurian block, but they are taken from six different studies (Table 2). We find that the angular velocity vector that best describes AM–EU motion has a rotation pole at 58.8°N, 157.5°E at a rate of 0.03°Myr−1 (Table 4, Fig. 9). Although our result is subject to a large uncertainty (Fig. 9), the obtained angular velocity is close in pole location but 0.04°Myr−1 slower than the result by Holt et. al. (2000a). Furthermore, our result is quite distinct from AM-angular velocities inferred from GPS vectors by Heki et. al. (1999). The expected extension directions across the Baikal Rift are subparallel to the T-axes of normal fault earthquakes. Expected extension rates across the rift, on the other hand, do not exceed ∼2 mm yr−1, which is lower than the ∼4.5 mm yr−1 from a regional GPS study (Calais et. al. 1998) and also lower than the recently inferred 3.2 mm yr−1 Quaternary slip rate for the northern Baikal Rift (San'kov et. al. 2000). We believe that a better constrained estimate of the Amurian Plate motion has to wait until more GPS measurements at more locations and over longer time-spans are available.
The angular velocity vector for South China (SC) with respect to Eurasia is subject to the assumed size and geometry of the South China Block. Owing to the presence of right-lateral faulting along the coast of South China and the occurrence of at least three earthquakes (MW= 5) near Shanghai (Harvard CMT catalogue), we constrain the South China Block to be relatively small and bounded by the Longmen Shan in the west and the Qinling fault in the north. We assume that Wuhan and three other sites, measured by Chen et. al. (2000), are located on this block. We find an SC–EU rotation rate of 0.22°Myr−1 about a pole at 47.3°N and 126.8°E, which yields an anticlockwise SC–EU motion (Fig. 9), consistent with results by Heki et. al. (1999) and Holt et. al. (2000a). On the southern part of the South China Block the predicted SC–EU motion amounts to 12 mm yr−1 in a direction roughly S66°E. We do not include Shanghai on the South China Block; whether or not Shanghai is part of the South China Block is reflected in a less than 1 mm yr−1 difference in model velocity in Shanghai. We estimate a model velocity of 6.5 ± 1.8 mm yr−1 at Shanghai in a direction of S59°E. Our obtained rate at Shanghai is consistent with the VLBI estimate by Molnar & Gipson (1996) and the GPS measurement by Chen et. al. (2000), but 3.5–6 mm yr−1 lower than the VLBI estimate of Heki (1996) and GPS estimates of Kato et. al. (1998), Heki et. al. (1999) and Larson et. al. (1999).
The notion that southeast Asia is moving with respect to Eurasia (Cardwell & Isacks 1978) has recently been confirmed by GPS measurements (Chamot-Rooke & Le Pichon 1999; Rangin et. al. 1999; Simons et. al. 1999; Michel et. al. 2001). The independently moving entity has been named the Sundaland or Sunda (SU) block. Here we will use the latter name. In our model the motion of the Sunda block is inferred from eight GPS velocities from the GEODYSSEA project (Michel et. al. 2001) and an additional velocity from the GPSVEL model in Singapore (NTUS). From the fit to the GPS velocities (Fig. 10) we find an SU–EU angular velocity with a rate of 0.13°Myr−1 about a pole at 26.0°N and 80.4°W (Table 4), which implies clockwise SU–EU motion (Fig. 9). This result is significantly different in location and rate from previous estimates based on regional modelling of earlier GEODYSSEA data (Chamot-Rooke & Le Pichon 1999; Kreemer et. al. 2000b). From our model velocity field we find roughly 10 ± 1 mm yr−1 towards S78°E ± 9° along the northern margin of the Sunda block and about 6 ± 2 mm yr−1 towards S61°E ± 12° along the southern edge of the Sunda block (Fig. 10). The obtained ESE direction of the model velocities is consistent with the observed velocities after they have been rotated into a model Eurasia reference frame; at Tanjung on the southern part of the block the observed velocity is 7.0 mm yr−1 towards S78°E, and at Non Nuoc in the northern part of the block the rotated observed velocity is 11 mm yr−1 towards S72°E. Previous models of Sunda block motion predicted a motion in a roughly NE–E direction (Chamot-Rooke & Le Pichon 1999; Kreemer et. al. 2000b). However, critical in the analysis of Sunda's motion with respect to Eurasia are the constraints placed on the motion of Sunda relative to South China for which the motion with respect to Eurasia is relatively well constrained. Any inferred Sunda–South China motion needs to be consistent with the lack of any significant deformation indicators between the two blocks (Fig. 10). Michel et. al. (2001) found that velocities in northern and southern Sunda with respect to Eurasia equal 14 and 10 mm yr−1, respectively. We find that an expected northern Sunda Block motion that exceeds ∼12 mm yr−1 (the inferred speed of the most southern edge of the South China Block with respect to Eurasia) would require left-lateral motion between the South China and Sunda blocks. Such left-lateral motion would be opposite in sense from the slip inferred from the few known earthquake slip vectors for events near the South China coast (Fig. 10) and is inconsistent with right-lateral motion of the Red River Fault (e.g. LeLoup et. al. 1995). In the model presented here right-lateral motion is predicted on the southern end of the Red River fault at a rate of about 1–2 mm yr−1, which we find to be insignificant. Unfortunately, the actual present-day slip rate on the Red River fault is still under debate (e.g. Zhao 1995; Cong & Feigl 1999).
Nubia—Somalia and Nubia—Eurasia
Whether or not the African Plate is comprised of two distinct entities (Nubia (NU) on the west and Somalia (SO) on the east) has been the source of much debate since the earliest plate model proposals (e.g. McKenzie et al. 1970; Chase 1972; Stein & Gordon 1984; Jestin et al. 1994; Chu & Gordon 1999). This question is of importance for the understanding of the relatively high seismicity occurrence in East Africa. In the case when Nubia and Somalia move separately a correlation is expected between relative plate motion and the seismotectonics in Eastern Africa. On the other hand, when Nubia and Somalia are kinematically one single entity then the East African Rift should be interpreted as a zone of intraplate deformation. The most recent studies that analysed plate motion data favour an independent Nubian and Somalian plate (Jestin et. al. 1994; Chu & Gordon 1999). We know of only of two studies that directly measured the relative motion between Nubia and Somalia geodetically (Crétaux et. al. 1998; Sella et. al. 2002). Crétaux et. al. (1998) inferred that the motions of the plates did not fit a single plate model, but they only used one station at the Somalian Plate (La Reunion). Sella et. al. (2002) used the IGS stations of the Seychelles and near Malindi, Kenya, on the Somalian Plate and derived a pole for NU–SO relative motion near the southern extent of their mutual plate boundary zone. The extent of the plate boundary zone can be seen in Fig. 2 and is mainly determined from distributed seismicity occurrence using the Harvard CMT catalogue. In the study presented here we include two velocities at the Somalian Plate (the DORIS velocity at La Reunion and the GPSVEL estimate at Malindi) and eight stations on the Nubian Plate (Libreville, Arlit, Sainte-Hélène, Dakar (all from Crétaux et. al. 1998), MATR (Egypt) (McClusky et. al. 2000), and Mas Palomas, Gough Island and Sutherland (GPSVEL)). Stations at Djibouti and Hartebeestoek are considered to be situated in the NU–SO plate boundary zone. From this study we find that the motion between the Nubian and Somalian plates is not significant. We find an NU–SO model angular velocity of 0.03°Myr−1 of about 23.3°S and 21.9°E (Table 4). For most of the plate boundary the model predicts extension directions that are consistent with the T-axes of the majority of the focal mechanisms from the Harvard CMT catalogue. However, model velocities across the East African Rift do not exceed 1.6 mm yr−1 of extension. The anomalously low NU–SO relative motion is a surprising result, since it is much lower than any previous reported estimate of extension rate across the East African Rift (about 6 mm yr−1 is reported by Chu & Gordon 1999 and 5–7 mm yr−1 by Sella et al. 2002). It needs to be emphasized though that the uncertainty in our estimated pole location is very large, and we believe that the inclusion of more geodetic velocities in the future on both the Nubian and Somalian plates will place better constraints on this problem.
Of particular interest to the kinematics of the two ‘African’ plates is the observation that both plates move relatively slowly with respect to Eurasia, when compared with previous estimates (DeMets et. al. 1994a; Larson et. al. 1997). Our model predicts that NU–EU convergence amounts to about 4 ± 1 mm yr−1 near the Azores and 5–6 mm yr−1 in the eastern Mediterranean, whereas NUVEL-1A predicts convergence rates ranging from ∼6 to 10 mm yr−1 going from west to east along the Mediterranean Plate boundary. If NU–EU convergence is indeed as slow as suggested by the GPS measurements, it will have implications for both kinematic and dynamic models of the Mediterranean Plate boundary zone and, ultimately, for seismic hazard assessment for this region.
A New No-Net-Rotation Model
We have presented a global kinematic model for almost the entire Earth's surface. From an earlier version of this model Kreemer & Holt (2001) had determined a no-net-rotation reference frame and had calculated the angular velocities of most plates with respect to the NNR reference frame. Here we present a new NNR model that is a small revision from the model by Kreemer & Holt (2001). The small difference between this model and the recent model is twofold. First, the estimates of some relative angular velocities presented in this paper are small revisions of the model results on which Kreemer & Holt (2001) based their calculation. Secondly, in the model on which Kreemer & Holt (2001) based their calculation they did not correct for the motions of the plates that are not constrained by geodetic velocities and which are otherwise likely to be erroneous (particularly in rate). Because the NNR model is based on the velocity field of the entire Earth's surface, it is important to add proper motions of the unconstrained plates to the model for the purpose of determining the NNR model. We have done this by adding the angular velocities of these plates in the calculation of the NNR reference frame. Unfortunately this means that we are adding angular velocities obtained from geological, seismic and seafloor data to the angular velocities derived geodetically in our model. However, this approach would lead to a closer estimate of the Earth's ‘true’ total velocity field in comparison with a case when we do not constrain the motions of the unconstrained plates at all. The angular velocities that we have applied are shown in Table 5. Table 6 shows the NNR angular velocities for all constrained plates (the reader is referred to Kreemer & Holt 2001, for a detailed explanation of the determination of the NNR model). Our new result indicates slightly faster plate motions than the model presented by Kreemer & Holt (2001); we find a total global root-mean-square velocity (Vrms) (i.e. the minimum average velocity of the total Earth's surface) to be 38 mm yr−1, whereas they found a Vrms of 37 mm yr−1. However, for most plates the new angular rotations are not significantly different from the previous result. We find that the effect of including geologically derived angular velocities (in the NNR frame calculation) for the plates that are geodetically unconstrained increases the Vrms for most plates with a few tenths of one mm yr−1. Similar to what Kreemer & Holt (2001) found, the NNR angular velocities presented here are at many places significantly different from the NNR model of Argus & Gordon (1991) (NNR-NUVEL-1A). In NNR-NUVEL-1A the surface of the Earth was assumed to be made up entirely of rigid plates. Regions of plate boundary deformation were defined to be part of rigid plates and the angular velocities were taken from the NUVEL1(A) model (DeMets et. al. 1994a). We now know that the NUVEL-1A motion for several plates (e.g. Nazca, India, Arabia, Nubia) differs significantly from the geodetically determined present-day motion and that plate boundary zones often move at a significant speed compared with their adjacent plates. Kreemer & Holt (2001) found that the differences between present-day motions and NUVEL-1A motions is found to be the largest contributor to the different results between the new NNR model and NNR-NUVEL-1A, but the inclusion of plate boundary zones has in some cases, e.g. Eurasia, a significant effect on the angular velocity as well.
Two of the key assumptions in the concept of plate tectonics are that plates are rigid and boundaries are narrow. While the model presented here incorporates diffuse deformation zones of considerable width, which is an improvement on previous plate motion models, the remainder of the Earth's surface is still assumed to behave rigidly. We obtain a good fit between the model and observed velocities on stable plates (Table 2) suggesting that, to first order, the assumption of plate rigidity is valid. Nevertheless, seismicity can be significant within the interiors of some plates (e.g. the Indian, Australian and North American plates) and, to the extent that intraplate seismicity is an expression of the same tectonic forces that govern seismicity within plate boundary zones, future models should ideally allow for intraplate strain rates to be accommodated as well. The current limitations to include intraplate strain rate calculations in the model presented here lie mainly within the relatively large uncertainties of geodetically measured velocities. Alternative estimates of intraplate strain rates could be obtained by dynamic models using force balance requirements (e.g. Wdowinski 1998; Bird & Liu 1999; Porth 2000). Also, Kreemer et. al. (2002) speculated that perhaps the number of shallow earthquakes above a cut-off magnitude can be used to set bounds on intraplate strain rates. After initial findings by Kagan (1999, 2002) that there appears to be a strong global correlation between tectonic moment rate and the number of events above a cut-off magnitude along subduction zones, Kreemer et. al. (2002) performed a more rigorous analysis using the strain rate model presented here. Kreemer et. al. (2002) found that the global correlation found by Kagan (1999, 2002) appears to be applicable to both subduction zones and most regions of continental deformation. Based on their findings, Kreemer et. al. (2002) argued that a correlation between seismicity rate and tectonic rate may also exist for intraplate settings. Whether this will prove to be correct or not, a combined knowledge of the regional strain rate and the seismicity rate could be very useful for all plate boundary zones to infer regional variations in seismic coupling, maximum moment and seismic hazard in general.
In the model presented here we have assumed that several smaller rigid blocks are present in the central and eastern Asia deformation zone; e.g. Amurian, Tarim, South China, Okhotsk and Sunda. The answer to the question of whether these blocks do in fact behave as independent rigid entities within a wide zone of distributed deformation has started to emerge more urgently with the use of GPS measurements. Although the relative angular velocities for these smaller blocks suffer in general from large uncertainties, the relatively good agreement between model and GPS velocities on these blocks (Table 2) suggests that, within the current uncertainties in GPS velocities, rigidity of these blocks can be assumed. However, we would like to point out that due to both the small number of velocity measurements on some of these blocks and the occasionally low seismicity rates along their boundaries, the presented geometries of these blocks are subject to change. Therefore, to test the robustness of our results, we estimate the angular velocities for the Tarim, Amurian and South China blocks, with the assumption that these regions are low strain rate areas, instead of a priori constraining one single angular velocity to the areas comprising each block. For this test case we average the rotation vector function obtained for the deforming grid areas comprising each block. The angular velocities obtained in this way fall well within the uncertainties of the angular velocities for these blocks presented in Table 4, indicating that, given the uncertainty in our result, a priori rigidity can be assumed for these blocks. In another test we evaluate different models that assumed a priori rigidity the Tarim, Amurian and South China blocks but did not all have geodetic velocities on these blocks included. In this way we could (indirectly) investigate the effect of assuming different block geometries. We find, for instance, that the exclusion of the two sites in the Sea of Japan (believed to be on the rigid Amurian block) changed the AM–EU angular velocity only marginally. Similarly, we find that the exclusion of the velocities at Wuhan (located in the northeastern corner of the defined rigid South China Block) and at two sites located at the northern margin of the Tarim Basin (part of the assumed rigid Tarim block (TA)) does only slightly change the SC–EU and TA–EU angular velocities, respectively, with the result being within the 1σ uncertainty of the originally obtained angular velocities. These results indicate that the obtained angular velocities for the Amurian, South China and Tarim blocks are relatively robust (i.e. insignificantly different) with respect to chosen block geometries, and that an a priori assumption of rigidity of these blocks is appropriate.
For some of the plates and blocks the geometries are subject to uncertainty, and, because they do not contain geodetic measurements, their motions are ill-determined as well. Consequently, the validity of including proposed independent rigid blocks, such as the Caroline Plate (Weissel & Anderson 1978), cannot yet be tested. Similarly, it is also difficult to assess whether proposed entities such as the Gorda Block (Silver 1971) and Sinai subplate (McKenzie et. al. 1970), which have not been included in the model, should be included as independent blocks.
To improve upon the kinematic model estimates for continental regions (in particular, for areas with low strain rates and few to no geodetic observations) additional geological information is welcomed. Currently we have only included Quaternary fault slip rates in central Asia. Concurrently, a large global catalogue of active, or potentially active faults is in the process of being compiled (Machette 2000). Such a data set may provide a tremendous amount of additional kinematic information that could be included in future modelling. Furthermore, the inclusion of regional centroid moment tensor data sets containing focal mechanisms with 5.5 ≥M≥ 4.5 (e.g. Pondrelli et al. 2002; Braunmiller et al. 2002) will provide additional new constraints, particularly for low strain rate areas where moderate- and large-sized events are rare.
We have obtained a global strain rate field and velocity field in a joint fit to 3000 geodetic velocities and geological strain rates, inferred from Quaternary fault slip rates in Asia. Since fault slip rates are only included for central Asia, seismic moment tensors are used to a priori constrain the style and direction (but not magnitude) of the model strain rate field for all regions for which no fault slip rate estimates are included. By solving for the velocity gradient tensor field on almost the entire Earth's surface we have attempted to integrate the principle of rigid-body rotations of plates with the kinematics of plate boundary zones. We have shown strain rate field examples for central Asia, the Indian Ocean and the Aegean, all of which are consistent in style and distribution with existing regional information from geology and seismology (as well as gravity undulations in the Indian Ocean). We have also discussed inferred angular velocities for several plate pairs. We obtain significantly lower rotation rates for the motions of Nazca, Indian, Arabian and Nubian plates relative to Eurasia when compared with the NUVEL-1A model estimates, while a significantly higher rotation rate is found for the Caribbean Plate with respect to its adjacent plates. We find significant motions for the South China, Tarim and Sunda blocks with respect to Eurasia within a global kinematic model. On the other hand, no clearly distinct motion was observed between Nubia and Somalia. The model presented here has already made valuable contributions to the determination of a new present-day no-net-rotation model and has also been used to positively confirm the correlation between tectonic moment rate and the number of events above a cut-off magnitude along subduction zones and within most areas of continental deformation (Kreemer et. al. 2002).
We thank J. Beavan, R. Bennett, E. Calais, J.-F. Crétaux, J. Freymueller, B. Hager, T. Kato, J. Kellogg and T. Sagiya for providing GPS vectors. We are grateful to W. Simons and the rest of the GEODYSSEA team for providing GPS vectors in southeast Asia. We wish to thank all members of the IGS community, and we are particularly grateful to G. Blewitt and D. Lavallée for making the GPSVEL version 0.2 model available. We wish to thank the USGS and SCEC communities for making their GPS solutions available. We thank J. Freymueller for a constructive review. We are grateful to L. Estey and C. Meertens at UNAVCO for setting up and maintaining an interactive website that displays the results of GSRM-1: . WEH is funded by NSF grants EAR-9909621 and EAR-0001160, and NASA grant SENH99-0325-0015. WEH also thanks the International Lithosphere Programme for support of the Global Strain Rate Map Project, ILP II-8. Figures were prepared using GMT 3.0 (Wessel & Smith 1991).