Constraints from GPS measurements on the dynamics of the zone of convergence between Arabia and Eurasia

We investigate the dynamics of deformation in the zone of convergence between the Arabian and Eurasian plates using a physical model that treats the lithosphere as a thin fluid sheet deforming in response to lateral variations in gravitational potential energy (GPE). This model has two free material parameters, the power law exponent, n, of the vertically averaged rheology of the lithosphere and the Argand number, Ar, which expresses the relative importance of GPE and stresses required to deform the lithosphere. With boundary conditions described by three additional free parameters, the model fits the observed deformation, as measured by 367 GPS velocities, with a root‐mean‐square residual of <2.4 mm/yr. We find negligible improvement when variations in material properties are introduced that represent increases or decreases in the lithospheric strength of the Central Iranian and Turkish‐Iranian Plateaux and the Zagros Mountains. Effective viscosity of the lithosphere ranges from 5 × 1022 Pa s at 10 nanostrain/yr to 1022 Pa s at 100 nanostrain/yr. As well as matching the decadal time scale strain rates determined from GPS, the computed distribution of strain rates is also consistent with the distribution and types of earthquake focal mechanism. Significant historical earthquake activity is seen in regions with strain rates lower than 20 nanostrain/yr, implying that it is prudent to base assessments of seismic hazard on regional strain rates.


Introduction
In the longitude range 40-62 ∘ E, oblique convergence between the Arabian and Eurasian plates at 20-30 mm/yr [Reilinger et al., 2006] is accommodated within a region approximately 2000 km long (NW-SE) and between 500 and 1000 km wide which covers Iran, Azerbaijan, and parts of Turkey and Armenia [e.g., Engdahl et al., 2006;Hatzfeld and Molnar, 2010;Jackson and McKenzie, 1988;Tatar et al., 2004;Walker and Jackson, 2004] (Figure 1). The surface manifestation of this deformation is a combination of reverse and strike-slip faulting, seen both in the geomorphic expression of the late Quaternary faulting and in the focal mechanisms of earthquakes [e.g., Baker et al., 1993;Berberian and Yeats, 1999;Copley and Jackson, 2006;Engdahl et al., 2006;Fattahi et al., 2006;Foroutan et al., 2012;Hollingsworth et al., 2006Hollingsworth et al., , 2010Jackson and McKenzie, 1988;Jackson et al., 2002;Javidfakhr et al., 2011;Le Dortz et al., 2009McKenzie, 1972;Meyer and Le Dortz, 2007;Meyer et al., 2006;Nazari et al., 2009;Shabanian et al., 2009Shabanian et al., , 2010Regard et al., 2005;Talebian and Jackson, 2002;Walker et al., 2009;Walker and Jackson, 2004;Walker et al., 2010] (see Figure 1). West of ∼57 ∘ E, this convergence is accommodated within the continental lithosphere of Eurasia by a combination of reverse and right-lateral strike-slip faulting in the Zagros, in the mountain ranges of the Alborz and Caucasus, and by convergence across the Apsheron-Balkhan Sill [e.g., Jackson et al., 2002]. East of 57 ∘ E, approximately three quarters of the convergence is accommodated by the Makran subduction zone [Vernant et al., 2004], with the remaining shortening taken up within the Kopeh Dagh and other mountain ranges in NE Iran. Right-lateral shear between central Iran and Afghanistan is accommodated by N-S trending strike-slip faults in eastern Iran.  [Walker et al., 2009] and seismicity from the Global Centroid Moment Tensor Catalog [Dziewonski et al., 1981;Ekström et al., 2012]. Blue, black, and red focal mechanisms show thrust-, strike-slip-, and normal-faulting events, respectively. Labels show locations of places and features referred to in the text. SCB refers to the South Caspian Block.
[e.g., Jackson and McKenzie, 1988]. While subsets of the GPS velocities in Iran have been used to describe the kinematics by the relative motions of rigid blocks [e.g., Reilinger et al., 2006;Walpersdorf et al., 2014], the velocity measurements are too sparse except in few small areas (Figure 2b), to permit more than ∼5-10 velocity measurements per block, so it is hard to test such models rigorously.
The occurrence of M ≳ 6 earthquakes within central Iran suggests that the plateau may not be rigid [Jackson and McKenzie, 1984], a suggestion that is supported by field evidence of faults that have produced earthquakes in Holocene time [e.g., Foroutan et al., 2012;Le Dortz et al., 2009Meyer and Le Dortz, 2007;Meyer et al., 2006;Nazari et al., 2009] and by recent GPS studies that show measurable strain rates in the region [e.g., Masson et al., 2005Masson et al., , 2007Masson et al., , 2014]. An alternative interpretation of the low rates of strain in central Iran is that they represent a state of low deviatoric stress in which the compressional stresses applied to the region are balanced by its higher gravitational potential energy [e.g., Hatzfeld and Molnar, 2010].
We use the GPS determinations of surface velocities in Iran and neighboring regions to constrain a set of numerical experiments that investigate the extent to which lateral variations in gravitational potential energy and in physical properties of the lithosphere control the distribution of faulting, seismicity, and seismic hazard in Iran.  Kreemer et al., 2014] that is used in this paper (see section 2.4). (b) The box contains the sites within 200 km of the profile X 1 -X 2 shown in Figure 2b, whose velocities are plotted in Figure 4. The polygon shows the boundary of the calculation domain. (b) Orientations and magnitudes of principal axes of strain rate calculated from the velocities shown in Figure 2a using the method of Shen et al. [2015] (Appendix A). Extensional principal axes are shown by red bars, contraction by blue bars (scale at top right). Background colors (scale at bottom) show the distribution of smoothing distances D j using a weighting (see Appendix A) that gives a median smoothing distance of 160 km. Line X 1 -X 2 shows the profile analyzed in section 3 and shown in Figure 4.

Physical Model
The physical model treats the continental lithosphere as a sheet of material whose deformation is determined by the vertical averages of its mechanical properties and of the forces acting upon it-the thin sheet approximation, in which horizontal gradients of the deviatoric stress associated with deformation of the lithosphere are balanced by gradients of the gravitational potential energy (GPE) per unit area of lithosphere [Bird and Piper, 1980;England and McKenzie, 1983;Houseman and England, 1986]. Here we outline the principal features of the model; details and means of solution are to be found in Appendix A1 and in England et al. [2016].

Gravitational Potential Energy
It has often been suggested that deformation of the continental lithosphere represents a balance between lateral variations in gravitational potential energy and the stresses required to deform the lithosphere [e.g., Dalmayrac and Molnar, 1981;England and McKenzie, 1982;England and Molnar, 1997;England and Houseman, 1989;Frank, 1972;Houseman and England, 1986;Jones et al., 1996;McKenzie, 1970McKenzie, , 1972Molnar, 1988;Molnar et al., 1993;Molnar and Lyon-Caen, 1988;Molnar and Tapponnier, 1978]. The gravitational potential energy per unit area (GPE) of a column of lithosphere, Γ, is where h is surface height, z is positive upward, and the level z = 0 is taken to be at a depth L, the nominal lithospheric thickness, below the geoid. If the column is in a state of local isostatic balance, the GPE is proportional to the depth-averaged vertical stress [e.g., Parsons and Richter, 1980]. We calculate the GPE in two ways, in each case assuming that all the columns of lithosphere are in isostatic balance with a reference column of lithosphere, with GPE Γ 0 , crustal thickness S 0 , and surface elevation h 0 (Table 1). One method assumes Airy isostatic balance, using the distribution of surface height shown in Figure 1. We define a reference column of lithosphere, with crustal thickness S 0 , and surface elevation h 0 (Table 1). Surface height is related to crustal thickness, S, under the assumption of Airy isostatic balance

WALTERS ET AL. DYNAMICS OF DEFORMATION IN IRAN
where c and m are the average densities of crust and mantle. For uniform crust overlying uniform mantle, the difference in GPE between the reference column, of GPE Γ 0 , and a column with crust of thickness S is where g is the acceleration due to gravity [Haxby and Turcotte, 1978].
The presence of widespread volcanism across the Turkish-Iranian Plateau suggests, however, that there may be significant lateral variation in the density of the mantle beneath the region, and evidence from seismic receiver functions and surface wave tomography also indicate that the elevation of the plateau is partly supported by relatively low-density lithospheric mantle [Al-Lazki et al., 2003, 2004, 2014Alinaghi et al., 2007;Amini et al., 2012;Innocenti et al., 1976;Kaviani et al., 2007;Maggi and Priestley, 2005;Mutlu and Karabulut, 2011;Sengör et al., 2003;Zor et al., 2003]. We therefore also estimate differences in GPE from differences in geoid height ΔN [Haxby and Turcotte, 1978] (see equation (1)), where G is the gravitational constant and Δ is the difference between the density and density at the same level in the reference column.   Figure 3a shows the GPE estimated from the Shuttle Radar Topography Mission 3 arc sec topography [Farr et al., 2007], smoothed with a Gaussian filter of width 100 km. Figure 3b shows GPE calculated from the EGM08 geoid model [Pavlis et al., 2012]. To remove long-wavelength contributions to the geoid, the model is truncated below degree and order 8, and a tapered filter is applied between degrees and orders 8 and 12; the taper, which was chosen so that the geoid height along the coast of Iran remains approximately constant, is similar to that used by Chase et al. [2002] to calculate the lithospheric contribution to the geoid. The main differences between the two estimates of GPE are in the Turkish-Iranian Plateau and at the southeast end of the Central Iranian Plateau, where the GPE estimated from geoid variations is higher than that estimated from Airy isostasy, and in the Zagros, where the GPE estimated from geoid variations is lower than that from isostasy.

Mechanical Properties
The mechanical properties of the sheet are described by a power law relation between the vertically averaged components of deviatoric stress and strain rate wherėi j is the ijth component of the strain rate tensor, assumed constant with depth through the sheet, ij is the ijth component of the deviatoric stress tensor, with the overbar denoting averages over the thickness, L, of the lithosphere, andĖ is the second invariant of the strain rate tensor, with the convention of summation over repeated indices. The fluid is assumed to be incompressible (̇k k = 0).
The viscosity coefficient, B, and the power law exponent in the rheology, n, define the physical properties of the lithosphere. The quantity 10.1002/2016JB013370 where T is the second invariant of the vertically averaged deviatoric stress tensor, ( kl kl ) 1 2 , is commonly referred to as the effective viscosity, by analogy with the viscosity of a Newtonian fluid in the case that n = 1. Note that for n ≠ 1 the effective viscosity depends on the deviatoric stress, or equivalently upon strain rate, so that lithosphere of spatially invariant viscosity coefficient, B, will in general exhibit lateral variation in its effective viscosity.
The relative importances of deviatoric stresses required to deform the lithosphere and GPE are expressed by the Argand number [England and McKenzie, 1982], where g is the acceleration due to gravity, c and m are the average densities of crust and mantle, respectively, and U 0 is a scale velocity. In what follows, the results of dimensionless calculations are expressed in dimensional terms using the values of parameters given in Table 1. The spatial variations of strain rate and velocity depend, however, only upon the value of Ar, not upon the individual quantities within it (see Appendix A1).

Region of Interest and Boundary Conditions
We divide the border to the calculation domain into seven sections whose lengths and locations are dictated by their geological settings ( Figure 3). The choice of a small number of long boundary segments is justified by a fundamental aspect of the governing equations. When, as is the case in this study, the boundary conditions are of normal velocity or traction, they influence strain rates over a distance from that boundary that is ∼D∕ √ n, where n is the exponent in the power law rheology (equation (7)) and D is the length of the boundary [England et al., 1985]. It follows that any variation in boundary condition over a scale that is small compared with the length scale of the calculation domain has only a local influence and does not significantly affect the distribution of stresses and velocities in the region as a whole.
We use a reference frame fixed to stable Eurasia. Where the calculation domain is directly in contact with a plate or (as in the case of the South Caspian) a smaller block that appears to be rigid, we assign a velocity boundary condition (segments EU, AR, and SC, Figure 3). The northern boundary is placed along the sharp topographic contrasts between Iran and the Caucasus, the Caspian Sea, and Eurasia. We treat this boundary (segment EU in the north of Figure 3) as lying within stable Eurasia and therefore set the velocity on this boundary to zero. We assign zero velocity to the eastern boundary, which follows approximately the border between Iran and Afghanistan (segment EU in the east of Figure 3) because estimates of seismic moment release and GPS measurements show that Afghanistan is moving slowly with respect to Eurasia [Ambraseys and Bilham, 2003;Masson et al., 2005]. Apart from the subduction zone of the Makran, the southern boundary of our calculation domain follows the Iranian coast and the marked step in topography between Iraq and the Turkish-Iranian Plateau (segment AR, Figure 3). These sharp topographic features approximate the boundary between the Arabian plate and deforming Iran; we fix the velocity on this segment using an angular velocity estimated from GPS sites within Arabia that lie close to the boundary. This angular velocity has a pole at 18.1 ∘ E, 26.5 ∘ N and a rotation rate of 0.38 ∘ /Myr, which differs slightly from the Arabia-Eurasia angular velocity of Reilinger et al. [2006] because of internal deformation of Arabia near this boundary.
The southern Caspian Sea moves relative to Eurasia [Jackson and McKenzie, 1984;Jackson et al., 2002;Copley and Jackson, 2006] at approximately 11 mm/yr at an azimuth of 335 ∘ [Copley and Jackson, 2006]. We apply that velocity condition to the boundary segment corresponding to the southeast and south coasts of the Caspian Sea (segment SC in Figure 3); the fits to the observed velocities are almost identical if we use the slightly different definition of the South Caspian motion of Djamour et al. [2010]. The relative velocities across three boundaries (segments M, SWC, and T, Figure 3) are poorly constrained. The section of the southern boundary east of 57.2 ∘ E (segment M, Figure 3) lies to the north of the Makran subduction zone, which accommodates much of the Arabia-Eurasia convergence at this longitude; the choice of a velocity condition for this boundary is, therefore, uncertain. A similar uncertainty arises at the border of the southwest Caspian Sea (segment SWC in Figure 3); the occurrence of reverse-faulting earthquakes on shallowly dipping planes suggests that the South Caspian is being thrust beneath the bordering continental lithosphere in this region [Jackson and McKenzie, 1984;Jackson et al., 2002;Copley and Jackson, 2006]. The third segment (T, Figure 3) represents the transition between predominantly contractional deformation in Iran and neighboring regions and the zone 10.1002/2016JB013370 of horizontal simple shear and extension in central Anatolia and the Aegean, which is the subject of a separate study [England et al., 2016]. At these boundaries to the calculation domain we apply normal forces per unit length, respectively, Γ M , Γ C , and Γ T , and allow these forces to be free parameters in the search described in section 4.

Observations of Crustal Velocity
We constrain our analysis with crustal velocities determined using the Global Positioning System (GPS). Geodetic measurements of the deformation of the lithosphere are taken at the top of the upper crust, which in general deforms elastically. We take the view that the strain of this thin elastic layer follows closely the ductile strain of the layer beneath [e.g., Bird and Piper, 1980;Bourne et al., 1998;Houseman and England, 1986]. Furthermore, except during the postseismic intervals of nearby earthquakes, the ductile strain rate may be considered time invariant on the time scale of geodetic observation. This view contrasts with interpretations of velocity fields that treat the lithosphere as an elastic medium whose surface deformation is controlled by the slip of dislocations buried within it [e.g., Hammond et al., 2011;Meade et al., 2002;Meade and Hager, 2005;Reilinger et al., 2006;Wallace et al., 2007]. Kreemer et al. [2014] combined published measurements of Aktug et al. [2013aAktug et al. [ , 2013b, Djamour et al. [2010Djamour et al. [ , 2011, Hessami et al. [2006], Kadirov et al. [2012], Karakhanyan et al. [2013], Mousavi et al. [2013], Peyret et al. [2009], Reilinger and McClusky [2011], Reilinger et al. [2006], Tatar et al. [2012], Tavakoli et al. [2008], Walpersdorf et al. [2006], and Walpersdorf et al. [2014] and their own processing of continuous GPS sites. They report 630 determinations of velocity from within our area of interest, but 195 of these are duplicated between different studies. We removed all duplicates, retaining only the most recently published velocity for each site, presuming those velocities to be based on the greatest number of reliable occupations for each site and on the most modern processing strategies. We also eliminate a further 68 sites that have a formal error greater than 2 mm/yr in their velocity, leaving 367 measurements of horizontal velocity within the region shown in Figure 2a.
Interpretation of the dynamics of the deformation depends on knowledge of strain rates, which we calculate from the GPS observations using the method of Shen et al. [2015] (Appendix A). The principal horizontal axes of strain rate are displayed in Figure 2b on an equidimensional surface grid with points separated by angular distance of 1 ∘ .

Preliminary Analysis
We first ask whether the observed distribution of velocities shows any influence from lateral variations in GPE, following an analysis suggested by D' Agostino et al. [2014], Dalmayrac and Molnar [1981], and England and Molnar [2015]. In the central part of the orogen, between about 50 ∘ and 54 ∘ E, the velocities and principal horizontal contractional strain rate axes are both aligned approximately north-south ( Figure 2), which allows us to use an approximate analytical solution to the thin sheet equations (equation (A1)). We analyze the velocities within 200 km of profile X -X ′ in Figure 2. Neglecting the component of velocity in the direction perpendicular to the profile and integrating the thin sheet equation in the direction parallel to the profile yields an expression that relates the integral of GPE along the profile to the velocity, u x , in the same direction where W is the length of the profile; U 1 and U 2 are the values of velocity at the south (x = 0) and north (x = W) ends of the profile; eff is the effective viscosity of the lithosphere, assumed constant along the profile; and ΔΓ is the difference between the GPE of the lithosphere at a point on the profile and that at x = 0 (Appendix A2).
Figures 4a and 4b show the variation of topography and of integrated GPE calculated from topography using equation (5), along the profile X 1 -X 2 of Figure 2. Figure 4c shows the variation of u x measured at GPS sites within 200 km of the profile and the linear least squares fit to these velocities, which has an RMS misfit of 1.5 mm/yr. The values of the parameters to this fit are as follows: U 1 , U 2 = 23 and 7 mm/yr and 4 eff L = 2 × 10 27 N m −1 s which, for L ∼ 100 km, is equivalent to an effective viscosity for the lithosphere of ∼5 × 10 21 Pa s. The principal features of the velocity variation are explained by the variation in GPE along the profile. Contractional strain (northward decrease of northerly velocity) occurs at the south of the profile where the topography (and GPE: equation (5)) is low; the contractional strain rate drops close  (11)); dashed line shows the linear part of this fit (first two terms in equation (11)).
to zero over the Central Iranian Plateau (∼300-700 km on the profile) then increases again to the north (∼800-1100 km on the profile), as the surface height decreases. This analysis suggests that the lithosphere of the region deforms in response to lateral variations in its gravitational potential energy. The compressional stresses at the southern border of Iran are transmitted to northern Iran via the higher lithostatic stresses beneath the plateau, in a fashion analogous to that described by Molnar and Tapponnier [1978] for Tibet [see also Hatzfeld and Molnar, 2010].

Estimation of Physical Parameters
Turning now to the two-dimensional problem constrained by the boundary conditions and GPE distribution indicated in Figure 3, we seek values for the five parameters that will allow an adequate fit to the observed GPS velocity field. Three parameters are the forces per unit length Γ M , Γ C , and Γ T applied at the boundary segments corresponding to the Makran, the South Caspian, and eastern Turkey. The other two define the physical properties of the lithosphere: the exponent n in the power law rheology (equation (7)) and the Argand number, Ar (equation (10)). We assume for now that the rheological parameters n and B are spatially invariant.

Parameter Search
The governing equations are solved subject to the boundary conditions described in section 2.3, and with the distribution of GPE derived from both topography and the EGM08 geoid model (section 2.1), using a finite element approximation to the equations for deformation on the surface of a sphere (see Appendix A1). The finite element mesh contains 5400 nodal points, 367 of which are placed at the locations of the GPS sites. The mean spacing between nodal points is approximately 20 km, significantly smaller than the mean spacing between GPS sites (Figure 2a).  (12)) mm yr −1 RMS Misfit (equation (13) a For a Newtonian fluid, the viscosity, , is equal to half the viscosity coefficient, B (equation (9)).
For each combination of parameters, we quantify the misfit between our calculated velocities, u C i , and the GPS velocities, u GPS i , by the weighted root-mean-square misfit, M where V i are the Voronoi areas of the N GPS sites (Appendix A). We use weighted RMS misfits in the parameter search in order to avoid biasing our fit toward regions with the highest density of GPS observations but, in discussing the quality of the fits, we quote the RMS misfit to all the observations weighted equally as follows: We checked the resolution of the calculations using meshes of 10,000 points and found that neither measure of misfit differed by more than ∼1% between the two meshes.
The Makran boundary segment is about 400 km long and lies farther than this distance from all but about 10 GPS observations ( Figure 2). By the argument of section 2.3, the condition applied at this boundary has little influence on the rest of the solution. We therefore estimate the parameter Γ M independently of the others. For n = 1, 3, 5, and 9, and with GPE estimated from the EGM08 geoid (section 2.1 and Figure 3b), we varied Ar and Γ M , finding the value of Γ M that produces the minimum misfits to velocities. These values were −0.625 TN m −1 , (n = 1); −3.0 TN m −1 , (n = 3); −2.75 TN m −1 , (n = 5); and −4.0 TN m −1 , (n = 9). Each of these corresponds to a compressional force per unit length. For subsequent calculations with n = 3, 5, and 9 we fix Γ M to −3.0 TN m −1 , the best-fitting value for n = 3, rather than imposing a different boundary condition for each value of n. For n = 5 and 9 this value fits the observations almost as well (difference of less than 0.01 mm/yr in M) as the bestfitting value. For n = 1 we use Γ M = −0.625 TN m −1 .
We then systematically searched the four-parameter space (n, Ar, Γ T , and Γ C ), finding the values of the parameters that best fit the observations with GPE estimated from both topography and the EGM08 geoid (section 2.1 and Figure 3). For n = 1, 3, 5, and 9, we varied Ar from 0 to 10 in steps of 0.25, and from 10 to 30 in steps of 0.5, and varied Γ T in steps of 0.125 TN m −1 from −6 to 1 TN m −1 and Γ C from −4 to 4 TN m −1 . The best-fitting values are listed, with their associated misfits, in Table 2 and shown in Figure 5. The best-fitting calculations with GPE calculated from topography have RMS misfits (equation (13)) in the range 2.8 to 3.2 mm/yr. The calculations using GPE derived from the geoid fit the data better (Figures 5b-5f and   (12)) obtained for calculations using GPE derived from an Airy isostatic model of topographic support (Figure 3a), varying the boundary forces Γ T and Γ C as described in section 4.1. Curves show, for each value of n, the value of Ar at which the minimum is found; numbers on each curve are at the minima and denote the value of n. (b) As in Figure 5a but for calculations with the GPE calculated from the geoid (Figure 3b). (c-f ) Distribution of M for calculations using GPE derived from the geoid; colors show the lowest value of M for each combination of Γ T and Γ C ; white circles correspond to the minima in Figure 5b. inside the boundary is between 2 and 2.5 TN m −1 (Figure 3b), consistent with the middle of the range of bestfitting Γ T . The best-fitting force per unit length on the boundary segment corresponding to the southwest coast of the Caspian, Γ C , differs from zero, for n > 1, by less than 0.5 TN m −1 .

Lateral Heterogeneity in Viscosity Coefficient
The solutions presented so far have constant material properties. The viscosity coefficient, B, is spatially invariant, and the effective viscosity of the lithosphere varies solely in response to variations in strain rate (equation (9)). It is reasonable, however, to investigate whether there may be large-scale variations in B. For example, it has been suggested that the Central Iranian Plateau is a quasi-rigid microplate or at least a region of intrinsically greater strength than its surroundings [e.g., Jackson and McKenzie, 1988;Sobouti and Arkani-Hamed, 1996]. In contrast, seismic wave speeds in the upper mantle are lower beneath the plateaux of central and northwestern Iran than beneath surrounding regions, suggesting that these regions are underlain by hot, weak uppermost mantle and may therefore be intrinsically weaker than their surroundings. Wave speeds below the Zagros Mountains are higher, suggesting that this region may have stronger upper mantle [e.g., Al-Lazki et al., 2003, 2004, 2014Alinaghi et al., 2007;Amini et al., 2012;Kaviani et al., 2007;Maggi and Priestley, 2005;Mutlu and Karabulut, 2011;Priestley et al., 2012].
We investigate this issue by defining three regions in which we allow the value of B to be a multiple, B ′ , of the value for the rest of the domain of calculation. The first of these regions corresponds to the region of

10.1002/2016JB013370
Holocene volcanic activity in the northwestern end of the Turkish-Iranian Plateau (TIP in Figure 6a), the second corresponds to the Central Iranian Plateau (CIP in Figure 6a) and the third to the Zagros Mountains (ZAG in Figure 6a). We first allow B ′ to vary separately in each region between 0.1 and 10, holding n = 3. Γ T and Γ C were fixed to their best-fitting values for the calculations with spatially invariant material properties ( Table 2). The variations of misfit with B ′ and Ar are shown in Figure 6b (Table 3). The reduction in misfit, in comparison with the best-fitting solutions for the homogeneous sheet (M = 2.7 mm/yr, Table 2), is negligible.
We next allowed the value of B ′ to vary in all three regions simultaneously, finding the best-fitting combination of Ar, B ′ TIP , B ′ CIP , and B ′ ZAG using the downhill simplex minimization of Nelder and Mead [Press et al., 1992]. The best-fitting solution gave B ′ TIP = 1.2, B ′ CIP = 1.0, and B ′ ZAG = 0.83 with M = 2.6 mm/yr and an RMS misfit of 2.2 mm/yr (Table 3). To check whether this procedure had settled on a local minimum, we calculated M for 100,000 random combinations of the parameters within the ranges shown in Figure 7, finding no combination of parameters that produced a lower misfit than that obtained in the inversion. Of these 100,000 solutions, those in which B ′ within the Turkish-Iranian Plateau lies outside the range 0.7 to 1.8 fit the data worse than calculations with a homogeneous lithosphere; the equivalent range of B ′ for the Zagros is 0.7 to 1.0. The lower bound for B ′ in the Central Iranian Plateau is 0.6, but an arbitrarily high B ′ for this region is permitted by the observations (Figure 7). There appears to be negligible trade-off between the parameters, and the ranges of acceptable values are very similar to those obtained when the B ′ s are varied separately ( Figure 6).

Comparison With Observed Velocities and Strain Rates
The RMS magnitude of the velocities measured with GPS in the convergence zone between Arabia and Eurasia ( Figure 2) is 11.6 mm/yr; the best-fitting calculation for the homogeneous sheet with GPE calculated from the geoid fits the observations with an RMS misfit of ≤2.4 mm/yr, a reduction in the original variance of 95% (Table 2). Observed velocities are compared with those from the homogeneous sheet with n = 3 in Figures 8a and 8b; of the 367 GPS velocities, 256 are fit to within 3 mm/yr. When B ′ is allowed to vary simultaneously in all three regions, the best-fitting calculation fits 267 GPS velocities to within 3 mm/yr (Figure 8c). For each calculation a set of observations in the Zagros, marked by black circles in Figures 8b and 8c, are misfit by ∼5 mm/yr whereas nearby sites are fit to better than 3 mm/yr. These sites are all from the study of Hessami et al. [2006], who report that their limited number of receivers did not allow the complete network to be measured in one session; it is possible that these larger misfits result from a problem with reference frames in that study. A similar consideration may apply to parts of the southeast Caspian coast and to the region around 45 ∘ E, 39 ∘ N, where few sites have large misfits, among larger groups with small misfits.
There are clusters of misfitting vectors in the Zagros-Makran transition zone, near the eastern and western coasts of the South Caspian, and at the boundary with Turkey at about 40 ∘ E, 40 ∘ N. The misfits in Zagros-Makran transition zone are confined to a region about 200 km on a side, in which the boundary conditions are probably more complicated than those we use here, and we do not consider these misfits further. The areas of misfit around the South Caspian and at the boundary with Turkey also have length scales of a couple of hundreds of kilometers and lie close to the boundaries of our calculation domain. By the arguments of section 2.3, these misfits could be reduced by adjustments to the boundary conditions in their immediate neighborhood without influencing conclusions about the large-scale deformation. For example, the observed velocities near the southwest Caspian coast (segment SWC, Figure 3) have a larger northerly component than do the calculated velocities; these misfits would be reduced by a corresponding increase in the northerly component of velocity on this section of the boundary. England et al. [2016] demonstrated the viability of reducing residuals by local adjustments to the boundary conditions, and the insensitivity of the overall solution to such adjustments, but we do not repeat that exercise here. Figure 9a shows the horizontal principal axes of strain rate derived from the GPS observations, at the 155 geographical locations of the grid shown in Figure 2b, superimposed on the second invariant of the corresponding strain rate tensor. Figures 9b, 9d, and 9e compare those observations with the results of the bestfitting calculations for the homogeneous sheet (Figures 8a and 9b). Observed strain rates range from 5 to over ∼60 nanostrain per year; at 141 out of the 155 locations the calculated and observed rates agree to within a factor of 2 (Figure 9d). In 123 out of the 155 locations the calculated and observed azimuths agree to within 20 ∘ (Figure 9e). Figures 9c, 9f, and 9g make the same comparisons for the calculations with lateral variations in strength (section 4.2 and Figure 8c). As with the fit to velocities, the influence of permitting lateral variations in B is negligible: the number of locations in which the calculated and observed rates agree to within a factor of 2 increases to 143, while the number of locations in which the calculated and observed azimuths agree to within 20 ∘ decreases to 121.
The RMS misfit of 2.4 mm/yr between calculated and observed velocities suggests that the thin sheet model accounts for the principal features of the deformation of Iran, as measured by the current geodetic data. The remaining misfits between observation and model may be attributed partly to data inconsistencies (see, for example, discussion of misfits at the beginning of this subsection) and partly to local deviations from the continuous behavior expected of a viscous lithosphere, due to shallow tectonic movements, to nontectonic movements, or to monument instability. The agreement between observed and calculated strain rate fields, at the scale of 100 km (Figure 9), allows us to investigate the implications of this model for the rheology of the lithosphere, for active faulting, and for seismic hazard.

Rheology of the Lithosphere
The viscosity coefficient for the lithosphere obtained from the calculations with spatially invariant material properties and n = 3 is B = 2.2 × 10 12 Pa s 1∕3 , which yields effective viscosities in the range 5 × 10 22 − 10 22 Pa s at strain rates of 10 to 100 nanostrain per year (Figure 10), consistent with the preliminary analysis (section 3). The range is similar for n = 5 and 9. This range of effective viscosity is comparable with estimates made for other regions of active continental deformation. Using the approach employed here, England et al. [2016] estimated the viscosity coefficient, B, of the lithosphere in Anatolia to be 6.0×10 11 Pa s 1∕3 for a homogeneous sheet with n = 3, which corresponds to effective viscosity of the lithosphere in the range 10 22 to 3 × 10 21 Pa s at strain rates of 10 to 100 nanostrain per year. Allowing the material properties of the lithosphere to vary spatially, Özeren and Holt [2010] estimated the viscosity to lie in the range ∼10 21 to 10 22 Pa s for most of that region. Estimates of the effective viscosity of the Tibetan lithosphere lie in the range 5 × 10 22 to 5 × 10 21 Pa s at strain rates of 3 to 30 nanostrain per year [England and Molnar, 1997;Flesch et al., 2001]. Estimates of the viscosity of the lithosphere in western North America lie in the range ∼10 21 to 10 22 Pa s [e.g., Flesch et al., 2000;Jones et al., 1996;Townsend and Sonder, 2001]. D' Agostino et al. [2014] obtained an estimate for the Apennines of 1.5 to 3 × 10 21 Pa s at a strain rate of ∼50 nanostrain per year, and England and Molnar [2015] estimated the effective viscosity of the Tien Shan to be 1 to 4 × 10 22 at a strain rate of ∼30 nanostrain per year.
The experiments summarized in Figures 6 and 7 place limits on the variations of the viscosity coefficient B within the Arabia-Eurasia convergence zone. Lateral variations in the effective viscosity by a factor of 5 or more igure 9. Comparison between observed and calculated strain rates. (a) Distribution of the second invariant of the strain rate tensor calculated from GPS velocities (Figure 2b and Appendix A). Superimposed crosses show the orientations and magnitudes of the principal horizontal axes of strain rate. (b) As in Figure 9a but for the best-fitting calculation with laterally invariant material properties (Figures 8a and 8b). (c) As in Figure 9a but for the best-fitting calculation with B ′ allowed to vary (Figure 8c). (d) Comparison between observed and calculated values of the second invariant of the strain rate tensor (equation (8) (9)) for the best-fitting calculation with laterally invariant material properties (Figure 8). arise within lithosphere with constant B (Figure 10), simply from the dependence of effective viscosity on deviatoric stress (equation (7)). These variations comfortably exceed the permitted variation of B ′ (section 4.2), with the possible exception that the upper bound on B ′ for the Central Iranian Plateau is effectively unconstrained.
The small permitted diminution in strength within the Turkish-Iranian Plateau (B ′ > 0.7) may seem surprising, in view of the suggestion that the volcanism of the region represents loss of the lower lithosphere and its replacement by asthenosphere [Ṡengör et al., 2003]. That process is envisaged, however, to involve only the weak, lower, part of the lithosphere [e.g., Houseman et al., 1981;Houseman and Molnar, 1997;Molnar et al., 1998] so, while loss of the lower lithosphere may appreciably increase the GPE of the remainder, it may not greatly alter its integrated strength [e.g., England and Houseman, 1989]. In contrast, the process of delamination is envisaged to remove the strong portion of the lithosphere [e.g., Bird, 1979;Bird and Baumgardner, 1981].
The tight constraint on B ′ for the Zagros (Figures 6 and 7) excludes the possibility that this region is significantly stronger than its surroundings, as might have been inferred from the high upper mantle wave speeds observed beneath the region. Priestley et al. [2012] argue, however, that these high wave speeds represent thickening of the lithosphere as a result of the convergence between Arabia and Eurasia, so they need not imply that the region is strong. Indeed, the 1-D analytical approximation (Figure 4) and the 2-D numerical calculations (Figures 6 to 9) each require that the lithosphere of the Zagros must be weak enough to deform under the observed gradients of GPE.

10.1002/2016JB013370
B ′ for Central Iranian Plateau may be as low as 0.6, consistent with the observations of low wave speeds in the uppermost mantle of the region. An arbitrarily strong plateau is, however, permitted by the observations (Figures 6b and 7) although the misfits for such solutions are negligibly different from those for spatially invariant physical properties. High values of B ′ are permitted because the deviatoric stress in this region is so low that its effective viscosity is already high (Figure 10), and further increasing B has very little influence on the velocity field. We therefore cannot rule out suggestions that the plateau is intrinsically strong (B ′ ≫ 1) [e.g., Jackson and McKenzie, 1988;Sobouti and Arkani-Hamed, 1996] but in view of the fact that it has undergone horizontal shortening during Tertiary time and is likely to be underlain by hot, weak uppermost mantle, it seems more reasonable to suggest that the plateau's present state of low strain rate is simply explained by the fact that its GPE is high enough to resist the compressional stresses applied by its surroundings [Hatzfeld and Molnar, 2010].

Comparison With Focal Mechanisms of Earthquakes
The agreement between the physical model and the strain rates determined from GPS velocities (Figure 9) suggests that the deformation in the convergence zone between Arabia and Eurasia is determined by the balance between gradients of gravitational potential energy and gradients of stress required to deform the lithosphere (equation (A6)). The geodetic strain rates are derived from measurements on the decadal time scale and reflect primarily interseismic deformation of the upper crust; the question arises, therefore, as to whether the same physical model can explain the longer-term deformation of the upper crust, which is primarily accommodated by earthquakes.
In testing the thin sheet model, we have assumed that the style of deformation in the upper crust follows the ductile strain of the lower lithosphere. Under that assumption the expected style of faulting may be calculated from the principal axes of horizontal strain rate in the thin sheet calculations [Gordon and Houseman, 2015;Houseman and England, 1986]. At any location the triaxial strain rate field of the ductile layer is compatible with a pair of double couple mechanisms acting in the brittle layer. The quantity indicates the expected style of surficial faulting, wherė1 anḋ2 are principal horizontal strain rates anḋ The thin sheet calculations show that faulting should lie mainly in the ranges of reverse faulting with subsidiary strike slip (RS), or of strike slip with subsidiary reverse (SR), with a small region of pure reverse faulting (RR) (Figure 11). We use earthquakes from the GCMT catalog [Dziewonski et al., 1981;Ekström et al., 2012] and separate them into reverse-, strike-slip-, or normal-faulting events according to whether the smallest, intermediate, or greatest principal axis of their moment tensor is closest to vertical. The regions corresponding to RR contain seven earthquakes, of which five have reverse-faulting mechanisms and two have strike-slip mechanisms. In the regions where reverse faulting with subsidiary strike slip is expected (RS), there are 141 reverse-faulting earthquakes, 58 strike slip, and 3 with normal-faulting mechanisms. In the regions where strike slip with subsidiary reverse faulting is expected (SR), there are 113 strike-slip earthquakes, 89 reverse faulting, and 9 with normal-faulting mechanisms. The normal-faulting earthquakes in the RS and SR regions account for less than 0.1% and 1% of the moment release in these regions respectively and probably represent stress heterogeneity associated with larger nearby earthquakes. Otherwise, the distribution of earthquake types is consistent with that predicted by the thin sheet calculations, with reverse faulting dominating in RS regions and strike-slip faulting dominating in SR regions.
The agreement between the style of deformation shown by our calculations and the distribution of earthquakes suggests that strains occurring in the ductile lithosphere are compatible with deformation in the seismogenic upper crust. This observation validates a key assumption of the thin viscous sheet model: that vertical gradients of horizontal displacement may be neglected when analyzing the deformation of the lithosphere at the scale of tens to hundreds of kilometers. Many of the earthquakes in our region of study that have magnitudes greater than about 6 are associated with faults that show displacements of hundreds of Figure 11. Predicted distributions of fault types for the best-fitting calculation with laterally invariant material properties and n=3 (Figure 8) compared with the observed distribution of earthquake focal mechanisms. Colors show the predicted distribution of fault types from the model strain rate field (equation (14), section 5.2); in the two-letter designations, N, S, and R, refer to normal, strike slip, and reverse faulting with the first letter referring to the type of faulting that accommodates the greater fraction of strain rate. Focal mechanisms are from the GCMT catalog [Dziewonski et al., 1981;Ekström et al., 2012]; reverse-, strike-slip, and normal-faulting earthquakes are shown with blue, gold, and red compressional quadrants, respectively. meters to kilometers over Quaternary time [e.g., Berberian and Yeats, 1999;Fattahi et al., 2006;Foroutan et al., 2012;Le Dortz et al., 2009;Meyer and Le Dortz, 2007;Meyer et al., 2006;Nazari et al., 2009;Shabanian et al., 2010;Regard et al., 2005;Talebian and Jackson, 2002;Walker and Jackson, 2004;Walker et al., 2010]. We therefore suggest that the Quaternary deformation of this convergent zone is characterized by the same balance between gradients of GPE and deviatoric stress that explains the decadal-scale deformation [see Hatzfeld and Molnar, 2010].

Implications for Seismic Hazard
Although the presence of large fault systems within the convergence zone ( Figure 1) encourages a description of the kinematics in terms of the relative motions of rigid blocks that are separated by those fault systems [e.g., Jackson and McKenzie, 1988;Reilinger et al., 2006;Walpersdorf et al., 2014], the assumption of block-like behavior carries the risk of serious error in the analysis of seismic hazard. If all the fault zones capable of delivering damaging earthquakes are identified and correctly assigned to block boundaries, then block models may provide useful constraints on slip rates and hence on the distribution of hazard. If, however, a block model  [Dziewonski et al., 1981;Ekström et al., 2012]; red dots show epicenters of earthquakes between 1900 and 1976 [Engdahl and Villaseñor, 2002]; yellow dots show the epicenters of earthquakes before 1900 [Ambraseys and Melville, 1982]. (b) Black line shows cumulative moment release for the earthquakes in the GCMT catalog allocated to 100 × 100 km 2 boxes, arranged in order of increasing strain rate of those boxes. Grey vertical lines show the number of earthquakes [Ambraseys and Melville, 1982;Dziewonski et al., 1981;Ekström et al., 2012;Engdahl and Villaseñor, 2002] that occurred within each box.
were to be based upon an incomplete catalog of faulting, then its use in seismic hazard estimation would be seriously misleading: hazard would be overestimated near the boundaries of the blocks and underestimated in their interiors.
The analysis of the previous two subsections indicates that, while some regions deform more slowly than others, there is little evidence for rigid behavior. Though the strain rates may be low, our inference that the brittle upper crust follows the ductile creep of the deforming lithosphere (section 5.3) implies that the brittle upper layer will eventually attain stress differences that can cause earthquakes. A more prudent approach to the seismic hazard of the region may therefore be to seek a relationship between seismicity and strain rate [e.g., Bird et al., 2010Bird et al., , 2015. To investigate this suggestion, we compared the seismic moment release since A.D. 1900 with the distribution of strain rate calculated from GPS measurements. Figure 12 shows the epicenters of earthquakes [Ambraseys and Melville, 1982;Dziewonski et al., 1981;Ekström et al., 2012;Engdahl and Villaseñor, 2002] superimposed on the second invariant of the strain rate calculated from GPS. Although the epicenters are more abundant in regions in which the strain rate is higher than 20 nanostrain/yr, dozens of earthquakes occurred in regions having strain rates lower than 10 nanostrain per year. Only the southwestern Central Iranian Plateau (around 54 ∘ E, 31 ∘ N), the largest region in which strain rate is lower than ∼10 nanostrain/yr, is aseismic on the historical time scale.
We next calculated the second invariant of the strain rate from the GPS observations (Appendix A) in 100 × 100 km boxes, and summed the scalar moments of the earthquakes within each box, using the GCMT catalog [Dziewonski et al., 1981;Ekström et al., 2012] for the interval from 1976 to 2014 and the catalog of Engdahl and Villaseñor [2002] for the interval 1900 to 1975. Figure 12b shows the cumulative seismic moment in these boxes, plotted in ascending order of strain rate in the box. Approximately half of the seismic moment release took place within areas where the calculated strain rate is between 5 and 20 nanostrain/year, with the rest taking place in regions where the strain rates lie between 20 and ∼50 nanostrain/yr (Figure 12b).
Although the Central Iranian Plateau has been relatively aseismic in recent time (Figure 12a), there were significant historical earthquakes in the Plateau [Ambraseys and Melville, 1982]. Furthermore, the region is cut by several faults over 100 km in length which show clear evidence of earthquakes in Holocene time, despite generally having time-averaged slip rates of lower than 1 mm/yr [e.g., Foroutan et al., 2012;Le Dortz et al., 2009Meyer and Le Dortz, 2007;Meyer et al., 2006;Nazari et al., 2009]. Similarly, parts of northeastern Iran (around 57 ∘ E, 37 ∘ E) and northwestern Iran (around 48 ∘ E, 40 ∘ E) that are 200-300 km on a side had few earthquakes since the beginning of the twentieth century and have strain rates lower than about 20 nanostrain per year. The historical catalog [Ambraseys and Melville, 1982] shows, however, abundant earthquakes in those regions ( Figure 12a). These observations suggest that the entire region may be susceptible to seismic hazard on faults that have low slip rates and long recurrence intervals, a situation that is common in regions of continental deformation [e.g., England and Jackson, 2011;Zhang, 2013]. The construction of a denser field of crustal velocities, either from GPS or InSAR, is essential to the better assessment of seismic hazard in Iran and neighboring regions.

Conclusion
A model that treats the lithosphere as a thin sheet with spatially invariant viscosity coefficient fits the observed velocity distribution with a RMS misfit of 2.4 mm/yr, a 95% reduction in variance (section 4.1). Allowing lateral heterogeneity in material properties produces negligible reduction in misfits (section 4.2, Figure 6, and Table 3). While the hypothesis that the Central Iranian Plateau is considerably stronger than its surroundings cannot be ruled out, this condition is not required by the available GPS measurements (section 4.2). It seems more likely that the low strain rate of the plateau is explained by a state of low deviatoric stress as its high GPE balances the compressional stresses applied to the region (see discussion of Figures 4,9,and 12) [Hatzfeld and Molnar, 2010]. A denser and more precise field of velocity is required to resolve whether more complex configurations of rheology are demanded by the observations. We conclude that the balance between gravitational potential energy and the stress required to deform the lithosphere shapes the deformation field of Iran and controls the distribution of active faults, seismicity, and seismic hazard across the region. Figure A1. Comparison between strain rate fields calculated using the method of Shen et al. [2015] with different degrees of smoothing. (a) Orientations and magnitudes of principal axes of strain rate calculated from the velocities shown in Figure 2a. Extensional principal axes are shown by red bars, contraction by blue bars. Background colors show the distribution of smoothing distances D j calculated using a weight of W T = 20, which gives a median smoothing distance of 120 km. (b) As in Figure A1a but with a weight of W T = 47.5, which gives a median smoothing distance of 200 km. (c) Comparison between values of the second invariant of the strain rate tensor (equation (8)) calculated at the location of each of the cross symbols in Figure A1a with W T = 20 and with W T = 35, which is the value adopted in the body of the text. (d) Comparison between the azimuths of the principal horizontal axis of contraction calculated with W T = 20 and 35. Symbols are colored according to the magnitude of the smoothing distance in Figure A1a. (e, f ) As in Figures A1c and A1d but comparing strain rate fields calculated with W T = 47.5 and 35. Symbols are colored according to the magnitude of the smoothing distance in Figure A1b. least squares inversion. The weight on each observation is the product of the area of the Voronoi cell occupied by the site and a Gaussian weight, exp , where R ij is the distance of site i from grid point j. The smoothing distances, D j , at the grid points are chosen so that the total weight of observations (W T in the notation of Shen et al. [2015]) is the same for each grid point. The significance of W T lies in the distribution of D j that it produces for a given data set. In the body of the text we set that weight to 35, which translates into a median smoothing distance of 160 km (Figures 2 and 9-12). Figure A1 shows the influence of changing the choice of weight upon the distribution of smoothing distances. A weight of 20 translates into a median smoothing distance of 120 km and produces strain rates in regions of denser GPS coverage that are up to 50% higher than in the fields with weight equal to 35 that we use in the body of the paper ( Figure A1c). A weight of 47.5 translates into a median smoothing distance of 200 km and produces strain rates in regions 10.1002/2016JB013370 of denser GPS coverage that are up to 25% lower than in the fields with weight equal to 35 ( Figure A1e). The azimuths of the axes of principal horizontal compression generally vary by less than 10 ∘ for the different choices of weight (Figures A1d and A1f.)

A1. Numerical Solutions to Thin Sheet Equation
The thin sheet approximation assumes that vertical tractions on vertical planes may be neglected, so that density variations within the lithosphere are in local isostatic balance. Under this condition, the stress balance equation for creeping flow in spherical coordinates is as follows: where and are, respectively, the colatitude and the longitude (in radians), r is the radial distance from the center of the Earth, L is the thickness of the lithosphere, and Γ is the gravitational potential energy per unit area of lithosphere (GPE) [Bird and Piper, 1980;England and McKenzie, 1983;Gordon and Houseman, 2015;Houseman and England, 1986]. These equations are nondimensionalized using the lithospheric thickness, L, as a scale length and using a scale velocity, U 0 ,̇′ and all variables are identified in Table 1.
This nondimensionalization allows us to express the rheology of the thin sheet by two parameters, n the power law exponent and the Argand number Ar = g c (1 − c ∕ m )L B(U 0 ∕L) 1∕n (A5) England and McKenzie [1982]. In presenting our calculations in dimensional form we use the values of g, c , m , U 0 , and L given in Table 1. The value of U 0 is set here to 30 mm/yr, but the best fit value of Ar varies inversely with U 0 , so that the implied physical value of B is independent of the choice of U 0 .

A2. Analytical Solution in One Dimension
In section 3 we consider a cross section through Iran, along which it is reasonable to treat the velocity as varying only in the along-strike (x) direction. Under this condition, the equation of stress balance reduces to Dalmayrac and Molnar, 1981]. The deviatoric stress is related to horizontal strain rate by We assume that the effective viscosity is constant along the profile; D'Agostino et al. [2014] show that the differences between solutions obtained for a Newtonian thin sheet ( eff constant) and those for a power law thin sheet with n > 1 ( eff varying with strain rate, equation (9)) are small, because the variations in strain rate occur over length scales that are primarily controlled by the variations in gravitational potential energy, Γ. For lithosphere of spatially invariant physical properties, equation (A6) may be integrated twice to yield a relation between the velocity in the direction perpendicular to the chain and the integral, in the same direction, of the GPE [D'Agostino et al., 2014]. We consider a mountain belt of width W and constant effective viscosity, eff , with the velocities at either end of the profiles being U 1 and U 2 . Equation (A6), when combined with equations (A7) and (A8), becomes Integrating twice, we obtain where ∫ x 0 ΔΓ(x ′ )dx ′ is the integral along the profile of the difference between the GPE of the lithosphere and that at x = 0. The boundary conditions U(0) = U 1 ; U(W) = U 2 give [England and Molnar, 2015]. Lynn Evans has made numerous contributions to the development and maintenance of the finite element code basil used to make the thin sheet calculations in this paper. This work was supported by the National Environmental Research Council and the Economic and Social Research Council through grant NE/J02001X/1 and by the Natural Environment Research Council (NERC) through the Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET). All figures were prepared using the GMT package [Wessel and Smith, 2013]. The data used in this analysis can be found in Kreemer et al. [2014] and references therein. We are grateful to Z.K. Shen for providing the code used for calculating strain rates from GPS measurements. We thank two anonymous reviewers for their helpful comments.