From Decades to Epochs : Spanning the Gap Between Geodesy and Structural Geology of Active Mountain Belts

7 Geodetic data from the Global Navigation Satellite System (GNSS), and from sa8 tellite interferometric radar (InSAR) are revolutionizing how we look at instantaneou s 9 tectonic deformation, but the significance for long-term finite strain in orogenic belts is 10 less clear. We review two d ifferent ways of analyzing geodetic data: velocity gradient 11 fields from which one can extract strain, d ilatation, and rotation rate, an d elastic block 12 modeling, which assumes that deformation is not continuous but occurs primarily on 13 networks of interconnected faults separating quasi-rigid blocks. These methods are 14 complementary: velocity gradients are purely kinematic and yield information about 15 regional deformation; the calculation does not take into account either faults or rigid 16 blocks but, where GNSS data are dense enough, active fault zones and stable blocks 17 emerge naturally in the solution. Block modeling integrates known structural geometry 18 with idealized earthquake cycle models to predict slip rates on active faults. Future 19 technological advances should overcome many of today’s uncertainties and provide 20 rich new data to mine by provid ing denser, more uniform, and temporally continu ous 21 observations. 22 Manuscript Click here to view linked References


Introduction
Structural geologists have always wanted to study the growth of structures in real time: how do limbs of folds rotate? How do fault damage zones evolve? How do extensional detachments and thrust belt dé collements work? A myriad of additional questions could be asked about different structural processes. Most structural geology involves the interpretation of process from the finite strain at the end point of the deformation. Clever methods, ranging from growth strata geometries to curved fibrous minerals in pressure shadows and veins, and sophisticated numerical mechanical and kinematic models have been developed to try to extract interpretations of incremental strains and insight into processes.
Today, we can sample in real time the surface effects of on-going structural processes via space-based geodesy. This ability, seemingly, should bring uniformitarianism to structural geology: the processes that we observe in action today should be the same as those that were responsible for the older structures that we see exposed in mountain belts. However, despite two decades of the widespread availability of Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR) geodetic data, such information is still relatively little used in the structural geology community. Two fundamental reasons underlie the gulf between traditional structural geology and geodesy, both related to the fact that much of the geodetic signal we observe is related to the earthquake cycle: first, a significant part of that deformation is elastic. We still need to understand the relationship between this non-permanent, infinitesimal strain and permanent, finite strain with which structural geologists are more familiar. Second, we have just two decades (or less) of space geodesy, but the earthquake cycle on plate boundaries lasts on the order of 100-200 years and in continental interiors can be more than 1000 years. Thus, we are trying to understand a finite processdultimately, mountain buildingdby sampling just 1/10th to 1/100th of the smallest significant unit of that process (i.e., the earthquake cycle).
In this paper, we review what has been learned to date about the relationship between geodetic and structural geology observations. We examine the most common analytical methods currently used (calculation of velocity gradient fields and elastic block modeling) and the issues and artifacts surrounding those methods. For both of these methods, questions about spatial scaling and continuous (or discontinuous) nature of surface deformation are key. Current advances in understanding of earthquake and volcanic processes, as well as regional tectonics, underscore the utility of working with geodetic data from a geological viewpoint. Finally, we examine how new technology will further enhance our ability to use these data.

The nature of space geodesy data
For more than 100 years, repeat surveys of geodetic monuments have revealed the deformation of the Earth's surface from tectonic processes (e.g., Yeats et al., 1997). The advent of space-based geodetic measurements has allowed these surveys to be done more frequently, more densely, more precisely, and over larger areas of the world.

Global navigation satellite system
The most commonly and widely used space-based technique (at least in terms of number of papers published) relies upon the establishment of a geodetic station and antenna that can receive signals simultaneously from several satellites in a global constellation. The constellation is commonly referred to as the Global Positioning System (GPS), but is now more generally termed the Global Navigation Satellite System (GNSS) because it includes the original GPS satellites as well as satellites from the European Galileo and the Russian GLONASS constellations, and may include others in the future (Hofmann-Wellenhof et al., 2008).
The GNSS observations are often divided into two end members: campaign observations that involve occupation of a geodetic benchmark for a few hours or days every year or so, and continuous stations that are not moved around and take a measurement as frequently as 10 times a second (e.g., Larson et al., 2003). The majority of published results to date report campaign data in which errors in vertical velocitydtypically three to five times greater than errors in the horizontaldare so large as to render them unusable (Segall and Davis, 1997). Thus campaign data usually provide only a two-dimensional velocity field and, because of their intermittent nature, cannot capture short-term transient deformation events nor fully describe the potentially significant annual or quasiseasonal component of the deformation field, points to which we return in a later part of the paper. Continuous GNSS observations avoid both of these drawbacks yet still reflect velocities determined at a relatively small number of discrete stations rather than more spatially continuous velocity fields.

Interferometric synthetic aperture radar
Satellite-based Interferometric Synthetic Aperture Radar (InSAR) has greatly expanded the spatial and temporal coverage of ground deformation (e.g., Bü rgmann et al., 2000;Rosen et al., 2000). InSAR is capable of measuring deformation of the Earth's surface with a pixel spacing of order 1-10 m over hundreds of kilometers. Accuracy on the order of 1 mm/yr can be obtained when many overlapping observations are combined in order to create a time series, or when stable ground points can be identified (e.g., persistent or permanent scatters or PSInSAR, Berardino et al., 2002;Ferretti et al., 2001;Hooper et al., 2004). An additional advantage of InSAR is that measurements can be made during every satellite overflight (weeks to months apart) without laborious ground surveys. In contrast to GNSS data, InSAR can provide more spatially complete velocity fields, however the measurement obtained is only of the component of the total velocity vector in the look direction of the satellite. That is, a single interferogram measures ground displacements that occur along the line of sight between the satellite and the grounddmeaning that multiple interferograms with different lines of sight are necessary to reconstruct the three-dimensional ground motion (e.g., Fialko et al., 2005). Furthermore, the frame of reference is local and corresponds only to the region over which the radar signal is coherent.

Effects of incomplete temporal sampling
Space geodetic observations may be affected by a variety of nontectonic processes including changes in the atmosphere the electromagnetic signals travel through, monument instability (due to thermal effects, soil creep, etc.) and real ground movements from hydrological or anthropogenic sources (e.g., Segall and Davis, 1997). Often these later types of deformation vary seasonally, for example, subsidence during a dry season when sub-surface groundwater reservoirs are depleted or recharged during a wet season (Reilinger and Brown, 1981;Schmidt and Bü rgmann, 2003) or downward flexure of the crust under a snow load during the winter (e.g., Heki, 2001;Blewitt et al., 2001). The deformation pattern may vary in amplitude from year to year and does not usually follow a simple seasonal (or sinusoidal) pattern due to other contributions, and so is often called quasi-seasonal. Such quasi-seasonal patterns are routinely removed from continuous GPS observations or when InSAR data or other ancillary information is available (e.g., Bawden et al., 2001;Argus et al., 2005).

Strain at a point from three or more stations
The simplest way to analyze deformation from inherently kinematic geodetic data, and that which is most comfortable to many structural geologists, is to calculate gradients in the velocity field (Allmendinger et al., 2007;Cardozo and Allmendinger, 2009). The pertinent equation is well known to structural geologists, though many are more familiar with the non-time derivative format in which D ij would be the displacement gradient tensor: where u i is the velocity at a station, D ij is the asymmetric velocity gradient tensor, x j is the coordinates of the station, and t i is a constant of integration representing the velocity of a point at the origin of the coordinate system. Because geodetic strain is most assuredly infinitesimal, it does not matter whether the initial or final station coordinates are used; technically the symmetric part of D ij is d ¼ 2 6 6 6 6 6 6 6 6 6 6 4 1 u 1 1 u 2 2 u 1 2 u 2 . . n u 1 n u 1 3 7 7 7 7 7 7 7 7 7 7 5 ¼ 2 6 6 6 6 6 6 6 6 6 6 4 the rate of deformation tensor with respect to the final coordinates and the time derivative of the strain tensor if initial coordinates are used (e.g., Malvern, 1969). In three dimensions, equation (1) represents a system of three equations (i, j ¼ 1-3) with 12 unknowns: the 9 components of D ij and the three components of t i . Because there are three equations for each station (corresponding to the east, north, and elevation components of the velocity vector and station coordinates), the solution for the unknowns is exactly constrained for four non-coplanar stations; for two-dimensional data, more commonly the case with GPS data, 3 non-colinear stations are required to solve for the four components of D ij and two of t i (as shown graphically in Fig. 1). We rewrite the system of equations to isolate the unknowns (in two dimensions):We can solve for the matrix of unknowns, m (the t i plus the components of the velocity gradient tensor), in the over-constrained case (i.e., n > 3) by standard least squares matrix inversion (Menke, 1984): Alternatively, equation (2) can be solved for m with a weighted least squares inversion (Menke, 1984): where W is the diagonalized matrix of weighting values, W. Equation (4) allows us to weight the contribution of each station used according to its distance from the point where the calculation is made (Shen et al., 1996), where a is a smoothing constant that specifies how the effect of a particular station decays with distance from the point of the calculation.
Once the velocity gradient tensor has been calculated, it may be additively decomposed into a symmetric strain rate and antisymmetric rotation rate tensors. For infinitesimal deformation, the first invariant of the strain rate tensor is the dilatation, or volume strain, rate of the region. Where GNSS data are two-dimensional, the full dilatation rate cannot be determined, but the 2D dilatation rate contains considerable information of interest (Allmendinger et al., 2007). Likewise, only a vertical axis rotation rate can be determined in two dimensions.

Spatial variation in strain rate
Of course, calculating a single strain rate from an arbitrary selection of stations assumes that strain is homogenous and is less useful than seeing how strain rate varies spatially across large areas and especially across structures of interest. In general, a rectangular grid or Delaunay triangulation using the actual GPS stations as nodes is constructed across the region of interest. If one starts with velocities and calculates strain rates, there is always a realistic solution, especially where the displacements over the time period of measurement are small compared to the station spacing. However, if one starts with strain rate data and calculates velocities, then St.-Venant's compatibility equations must be satisfied to make sure that the strain is compatible (Malvern, 1969;Haines and Holt, 1993;Holt et al., 2000). The Global Strain Rate Map (Kreemer et al., 2003) is an example of the application of this theory utilizing strain rates inferred from earthquake moment tensor solutions, and Quaternary fault-slip rate data combined with GNSS data.
Once a grid has been defined, there are two simple approaches to calculating strain and rotation rate at each node. In the nearest neighbors method, only the n stations nearest to the grid node are used in the calculation of strain at the node. Because stations are not typically distributed uniformly across the region of interest, this approach has the effect of a spatially varying length scale, as described below. Alternatively, one can use the distance weighted method of Shen et al. (1996) where the strain rate at each node is calculated from all of the GNSS vectors in the data set, but the contribution of each vector is weighted by its distance from the node. The weighting factor, a, in equation (5) provides a means of smoothing the strain rate field at different wavelengths.

Length scale and irregular station spacing/distribution
The pitfalls of inverting GNSS vectors for strain rate should be familiar to structural geologists: they involve the inherent length scale dependence of strain and artifacts that can arise from the fact that our observations are not regularly spaced but are clustered in some areas relative to others. None of this would matter if the strain rate was homogeneous across the region of interest, but it is profoundly heterogeneous at both the finite deformation scale and the geodetically infinitesimal scale (Fig. 1). When calculating the shear strain across a broad shear zone, the result depends on where the measurements are made (Ramsay and Graham, 1970) and two transects across the shear zone will yield quite different results if the measurements on one transect are clustered near the center of the zone and in the other are taken far from the center. In the example shown in Fig. 1, the magnitude of strain and orientation of principal axes would be different if the three stations were closer to the fault (higher strain) or farther away from the fault. The same is true for GNSS data and we illustrate this with an interesting and potentially significant case from the Himalayan front.
The Himalayan front is well known for generating a number of large, damaging earthquakes (Bilham et al., 2001;Bilham, 2006). In northwestern India, the w150 km on either side of the Dehra Dun recess in the main boundary thrust have not slipped since 1400 A.D. (Banerjee and Bü rgmann, 2002;Kumar et al., 2006). Except for a short stretch in Bhutan, this segment is the longest-lived seismic gap along the entire Himalayan front, and thus presents considerable potential seismic hazard. Analysis of the composite GPS data set from Zhang et al. (2004) using a nearest neighbor analysis produces a pronounced shortening rate anomaly in the middle of the seismic gap ( Fig. 2a): the principal horizontal shortening rate (which has a negative sign) there is about twice that on either side of the gap. Because our temporal sampling interval is much less than the earthquake cycle, a common strategy is to use along strike variation as a proxy for time; thus, it is tempting to interpret this anomaly as possible evidence of impending rupture on this segment, the anomaly is more likely due to the dense cluster of GPS stations (Banerjee and Bü rgmann, 2002) in the recess. If we subsample the GPS stations down to a density and distribution more nearly similar to that elsewhere along the front, the anomaly disappears ( Fig. 2b); it also disappears when the distance-weighted algorithm is used (Allmendinger et al., 2007).
The conclusion that the shortening rate anomaly is probably an artifact of station density does not necessarily mean that this seismic gap is quiescent and unlikely to experience a large earthquake in the near future. The rate of strain accumulation on a locked fault segment may be relatively constant throughout the interseismic period and should be assessed from the point of view of dislocation modeling. In the much more densely sampled eastern margin of Tibet, no strain rate anomaly was observed in data collected in the decade prior to the Sichuan earthquake of 2008 ( Fig. 3). However, in that area, the network consists of campaign stations and, continuous observations might be necessary to capture any pre-seismic anomalies, if they existed at all.
We illustrate the effect of length scale by two different calculations of two-dimensional dilatation rate of continuous GNSS data from the Plate Boundary Observatory (PBO) for the western United States (UNAVCO, 2008). If deformation is constant volume, then the 2D dilatation rate should be positive where excess horizontal extension rate occurs, implying a region of crustal thinning, and negative where there is an excess of horizontal shortening rate, suggesting a region of crustal thickening (Allmendinger et al., 2007). These regional patterns are very well developed when the PBO data are analyzed using the distance weighted algorithm with a smoothing distance, a ¼ 100 km, appropriate to the scale of the tectonic province (Fig. 4a). The choice of a depends both on the spacing of the GNSS network and the scale of the problem of interest and thus will vary from province to province and network to network. The Basin and Range is a province of crustal thinning on both finite and infinitesimal scales, whereas Cascadia, dominated by an active subduction zone, is an area of crustal thickening produced by interseismic locking of the plate boundary.
A completely different, but equally valid, picture of the western United States emerges when a nearest neighbor analysis is applied (Fig. 4b). As emphasized above, clusters of stations can produce apparent strain rate anomalies, but they can also highlight real zones of rapid strain rate. When the PBO data are analyzed using twelve nearest neighbors, the areas that emerge are largely ones of active volcanic deformationdMt. St. Helens, the Long Valley Caldera, and Yellowstonedand to a lesser extent the San Andreas fault and Mendocino Triple Junction. For the volcanic centers, the deformation reflects inflation or deflation of the underlying magma chambers so the assumption of constant volume does not hold. Notice, too, that there is a factor of 40 difference in the magnitudes of strain that are calculated from exactly the same data set. The difference, fundamentally, is the length scale over which the calculation is made.

Block modeling
In addition to assessing the deformation recorded by geodetic data using the velocity gradient method, several classes of quantitative models have been used to interpret deformation measured by geodesy, namely continuum, microplate, earthquake cycle, and block models. These methods have been described and compared in previous review papers (e.g., Thatcher, 1995Thatcher, , 2003Thatcher, , 2009Meade and Loveless, in press) but we comment here specifically on the bearing these models have on geodetic and geologic observations. Continuum models assume that crustal deformation can be effectively represented as smoothly varying (e.g., England and Molnar, 2005;Haines and Holt, 1993). The ''thin viscous sheet'' parameterization of bulk lithospheric rheology allows for the linear estimation of crustal strain rates resulting from applied forces (e.g., England and McKenzie, 1982). Microplate models of continental deformation (e.g., Avouac and Tapponnier, 1993;Thatcher, 2007) adopt the kinematic description of motion central to plate tectonic theory, describing geodetically determined velocities as resulting exclusively from the rotation of crustal blocks about Euler poles. By nature, such motion is discontinuous across block boundaries. Continuum and microplate models have been Tibet based on data from Shen et al. (2005). The distance weighting factor a ¼ 50 km and the magnitudes are in nstrain/yr (10 À9 yr À1 ). The location (red star) and focal mechanism for the 2008 Sichuan earthquake, which post-dates the GNSS data, are shown. Note the lack of strain anomaly in the vicinity of the earthquake, despite reasonably dense station spacing. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)  Zhang et al. (2004). More negative shortening rates are larger. The irregular ellipses with four digit numbers show the extent of historic seismic rupture segments of the Himalayan front (Bilham, 2006). (a) Analysis calculated from the 15 nearest neighbors using all available stations. Note the station cluster and large negative shortening rate centered in the region of the 1400 AD rupture. (b) Stations in the vicinity of the cluster have been randomly subsampled to produce an average density of stations more nearly similar to elsewhere along the Himalayan front 1 nstrain is a strain of 1 Â 10 À9 . considered to lie at opposite ends of a spectrum of crustal deformation hypotheses (Thatcher, 1995), with the latter assuming that measured deformation is accommodated on a relatively small number of major crustal structures. The continuum limit approximates a pervasively fractured upper crust with deformation accommodated on an infinite number of faults. Therefore, continuum models do not provide information regarding slip rates on specific geologic structures. Earthquake cycle models (Savage and Burford, 1973;Savage, 1983), which are commonly used to model inter-, co-, and/or post-seismic deformation around individual or small groups of structures, rely on dislocation theory to describe the accumulation and release of elastic strain near faults. Block theory combines microplate and earthquake cycle models, describing interseismic geodetic velocities as the summed effects of crustal block rotation and the earthquake cycle processes that occur on the faults defining block boundaries (Matsu'ura et al., 1986;McCaffrey, 2002;Meade and Hager, 2005;Meade and Loveless, in press). By integrating seismic cycle models with microplate rotations, block models provide information about relative motion across discrete structures as well as the smooth velocity gradients produced by the accumulation of elastic strain on these structures.
In microplate and block models, the geometry of the crustal blocks is generally inferred from maps of plate boundaries and active faults. The entire area of interest is discretized into a network of contiguous faults that intersect to define the boundaries of blocks. Generating a block geometry from an active fault map thus requires connection of structures that have been mapped as, and may in fact be, discontinuous (e.g., en echelon faults, tip line folds, etc.) and assumptions regarding fault intersections. Additionally, it may be impractical to incorporate all active faults into a block model given the distribution of GPS stations, as closely spaced faults may produce deformation signals that are difficult to identify uniquely. These compromises are well understood by practitioners as artifacts of the current state of the modeling procedure. For example, fault terminations within a block result in unbounded stresses at the tip line in most dislocation theory (leading some to explore fractal approaches, King, 1983), even though fault tip lines certainly exist in nature. Dense InSAR observations could be combined with the GPS observations to test the existence of closely spaced faults in a block model, but this has not yet been done (although see Fialko, 2006 for an example of how InSAR contributes to estimates of fault slip rates even in areas with dense GPS).
Quantitatively, the frequency distribution of both geodetically and geologically constrained fault slip rates in southern California suggest that 97% of the deformation between the Pacific and North American plates is accommodated on faults slipping !1 mm/yr (Meade, 2007). Considering only these faults significant in a block model means that the accumulated effects of faults slipping more slowly are aliased or mapped onto the faster-slipping, modeled structures (Meade, 2007). Thus the general question ''How many blocks are necessary?'' can be quantitatively stated as ''How many blocks are necessary to model deformation in a particular region at a given level of precision?'' which can be determined through the testing of multiple fault system geometries. Structural geologists make much the same supposition when they construct a balanced cross section because the implicit assumption is that that the deformation due to faulting, and to a much lesser extent folding, can be largely accounted for by a relatively small number of thrust faults (Marrett and Allmendinger, 1990, 1991, 1992.
Basing a block model geometry on an active fault map facilitates comparison between the geodetically constrained slip estimates from the model and information about the long-term motion on faults generally provided by paleoseismological or geological data. Of course, temporal variations in fault activity can hinder the comparison of geologic and geodetic slip rates, because the latter can be estimated only on those faults that are active, literally, during the deployment of the geodetic network. However, these discrepancies provide an opportunity to understand better temporal variations in fault system behavior. Even if only the sense, but not rate, of long-term slip is known for a particular structure, it provides valuable information against which the block model results can be compared. At the same time, block models are capable of making predictions about the structural characteristics everywhere in a plate boundary zone, including regions where the model may be a poor description of faulting. In particular, localized groups of coherent residual velocity vectors may suggest the presence of a previously unknown structure, or anomalous behavior, such as interseismic creep or partial coupling. Block models thus can be particularly useful in remote regions, such as rugged and/or forested mountains and areas offshore, where the distribution of active faults is likely to be less well known. Additionally, the coupling of fault slip rates and rotational motions in the classic block model formulation implicitly predicts the propagated consequences of changing a single slip rate throughout the entire fault network.

Lessons learned
In the last decade, several significant earthquakes have been captured by GNSS and/or InSAR data (e.g., Fielding et al., 2005;Klotz et al., 1999;Pritchard et al., 2002;Reilinger et al., 2000;Simons et al., 2002). All of these studies demonstrate unequivocally that elastic strain accumulation and release is the dominant signal in geodetic data.

Interseismic deformation and permanent strain
At a regional scale within continents, interseismic deformation is most nearly similar to regional late Cenozoic tectonic deformation (compare the strain in Fig. 1b and c). This is particularly true for orientation but also, somewhat surprisingly, for magnitude. Analyses of strain rates in the Tibetan and Anatolian plateaus (Allmendinger et al., 2007) demonstrated that infinitesimal maximum shear strain planes are parallel to known, active, longlived strike-slip faults (Fig. 5), and that principal axes are orthogonal to active, long-lived thrust and normal faults. Additionally, interseismic vertical axis rotation is consistent with long-term rotation observed in eastern Tibet, which is essentially a rigid body rotation (1.8 AE 0.5 /m.y.) bounded by the Ganze-Xianshuihe-Xiaojiang fault system, and across the Bolivian orocline in the Central Andes, a rotation due to distributed simple shear. In the latter area, the magnitude of the instantaneous vertical axis rotation, integrated over 25 million years is the same as the long-term vertical axis rotation observed by paleomagnetic studies . Likewise, in most but not all areas, 2D interseismic dilatation rate from GNSS data commonly reflects the long-term nature of the tectonic province: for example, the Aegean is an area of crustal thinning, the Altiplano is one of crustal thickening, etc. (Allmendinger et al., 2007).
The Basin and Range of the western United States is a particularly good area to explore some of these questions (Fig. 4). At a regional scale, the areas of positive 2D dilatation rate (crustal thinning) calculated from the PBO data set (UNAVCO, 2008) coincide with the loci of earthquake seismicity and the general distribution of Holocene faults as mapped by U.S. Geological Survey, et al. (2006), both of which are concentrated in the topographic lows on either side of the province. The topographically higher, central part of the province is less active tectonically and has nearly neutral, or even slightly negative, 2D dilatation rate. This region correlates spatially with the newly defined mantle ''drip'' (West et al., 2009).
The fit between geodetic and permanent deformation in the Basin and Range is less good when examined at the scale of individual faults. Friedrich et al. (2003) argue for ''Wallace-type'' behavior for the Wasatch fault, which may have experienced a Holocene cluster of activity atypical of longer term behavior. A more glaring apparent discord between geodetic and permanent deformation can be found in north-central Nevada. Across the Crescent normal fault, an east-west instantaneous shortening rate between two stations has been explained as a transient effect due to fault unloading Wernicke et al., 2000), although more recent work indicates that the geodetic signal is debatably anthropogenic and due to mining and groundwater withdrawal (Gourmelen et al., 2007;Wernicke et al., 2008). Alternatively, when analyzed in two dimensions, the entire central Nevada region is one with slightly negative 2D dilatation rate indicating some more regional effect than a simple transient or local anthropogenic effect (e.g., West et al., 2009). In any case, these examples and others (Niemi et al., 2004;Oskin and Iriondo, 2004;Oskin et al., 2007) illustrate the potential dangers of assuming a priori that short-and long-term strain accumulation on individual faults are similar.
In some areas, interseismic deformation measured geodetically is conspicuously at odds with regional long-term deformation as, for example, in the Chilean forearc. Interseismic loading while the Nazca-South American plate boundary is locked produces large, margin-perpendicular shortening rate on the scale of the orogen. However, most of the structures in the Chilean Coastal Cordillera manifest margin-normal extension (Allmendinger and Gonzá lez, in press; Gonzá lez et al., 2003;Loveless et al., 2005;Niemeyer et al., 1996;von Huene and Ranero, 2003), which would appear to be more compatible with coseismic strain (Fig. 6). Though it is known that some permanent structures in the Coastal Cordillera form coseismically (Loveless et al., 2009), simple elastic dislocation models suggest that interseismic strain accumulation on the subduction interface can also produce a narrow belt of extension (Loveless, 2007). The network of GPS stations in northern Chile is, unfortunately, too sparse to test this hypothesis and evaluation of interseismic deformation from InSAR is complicated by the atmospheric effect resulting from persistent coastal fog (Loveless and Pritchard, 2008). Farther south in Chile, post-seismic transient deformation still dominates nearly 50 years after the great Valdivia earthquake of 1960 (Khazaradze et al., 2002;Wang et al., 2007).

Complementary nature of block modeling and strain analysis
Despite very different assumptions, inversion of velocity fields for strain and rotation rate, and the velocities predicted by elastic block modeling are complementary. Calculating the gradient of a geodetically constrained velocity field using a reasonably large distance weighting parameter yields a picture of the nominally interseismic, large-scale deformation at a plate boundary zone, while elastic block model analysis permits estimation of activity on discrete structures that accommodate the relative block motion. Fig. 7 demonstrates this relationship in Japan, where densely spaced geodetic observations show a northeast trending high strain rate zone coincident with the Niigata-Kobe Tectonic Zone (NKTZ, e.g., Sagiya et al., 2000). The structures comprising the NKTZ in a block model (Loveless and Meade, in press) are estimated to experience up to 15 mm/yr of dextral and reverse slip, consistent with the large magnitude shear and contraction rates given by the velocity gradient analysis.
The velocity gradient calculation makes no assumptions regarding rheology nor the existence and behavior of geologic structures, while block models inherently describe the geodetic velocity field as arising from processes occurring on a fixed geometry of interacting faults. At present, block models simulate the accumulation of interseismic strain on faults using geometric extensions of idealized earthquake cycle models (Savage and Burford, 1973;Savage, 1983) while additional lithospheric behaviors, such as visco-elastic relaxation (e.g., Savage and Prescott, 1978;Hetland et al., 2008;Hilley et al., 2009) are only beginning to be implemented into the theory.
The velocity gradient calculation provides information about vertical axis rotation, horizontal distortion, and, assuming constant volume deformation, 2D dilatation, which gives insight into the vertical tectonics of a region (England and Molnar, 2005;Allmendinger et al., 2007). The nature of the calculation yields maps of these parameters that vary smoothly throughout a region, with concentrations of strain occurring around prominent crustal structures. Block models, on the other hand, predict velocity fields that can vary smoothly or be characterized by discontinuities, depending on the balance of block rotation versus elastic strain accumulation effects. An abrupt change in the velocity field across a structure indicates a low degree of interseismic locking between the bounding crustal blocks, whereas a smoother velocity field gradient is the signature of a fault actively accumulating elastic strain throughout the seismogenic zone. Allmendinger et al. (2007, their Fig. 11) show that, given GPS stations spanning a fault and spaced similarly to typical modern networks, the strain field associated with interseismic locking on the fault is nearly indistinguishable from that resulting from a static offset, demonstrating consistency between interseismic and longterm stylesdand potentially magnitudesdof deformation. That is, a calculation of strain from GNSS data is similar to the predicted long-term strain field, at least across a single structure so long as the fault slip rate has not changed in time. Block models, on the other hand, are capable of distinguishing elastic deformation effects from the offset signal, thereby providing information about the locking depth of the fault or, more specifically, the spatial distribution of elastic strain accumulation on it. By differentiating the velocity field in the gradient calculation, using weighting parameters appropriate for evaluating regional tectonics, we potentially lose information about elastic strain accumulation on individual structures, although constraining the strain accumulation requires making the block model assumptions regarding fault geometry, processes, and rheology.
Finally, calculation of velocity gradients might be considered to correspond to the continuum class of models, but this is not necessarily the case. Where the continental blocks are large relative to the spacing of the GNSS network (and the weighting factor used in the case of the distance weighted algorithm), stable, internally little deformed blocks emerge naturally from the analysis. For example, the Tarim Basin clearly emerges as a low strain rate block relative to Tibet and the Tien Shan, as does the Sichuan Basin on the eastern edge of Tibet (Allmendinger et al., 2007). Unfortunately, Fig. 6. Coseismic deformation captured during the 1995 Mw 8.1 Antofagasta, Chile earthquake by GPS data (Klotz et al., 1999). Red line segments show the orientations of the local principal horizontal extension rate axes and color ramp depicts magnitudes. Heavy black lines show the locations of young forearc normal faults (tick marks on the down-thrown plate) and right lateral strike-slip faults. Although new surface cracks formed during this event, none of these faults are known to have moved during the event (Delouis et al., 1998;Gonzá lez et al., 2003). Units are years À1 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) within deforming regions, the spacing of current GPS networks remains too coarse to allow the rigid blocks to emerge naturally as the output of a velocity gradient analysis. Individual station spacing would need to be significantly less than the spacing of the geodetically active fault network.

Coseismic deformation most related to fault zone processes
Earthquakes can cause large ground displacements (up to [10s] of m in a few seconds for Mw 8-9 earthquakes) and so it is not surprising that many earthquakes have been observed geodetically (e.g., Yeats et al., 1997). The coseismic ground displacements involve a static (permanent) offset of the ground surface that has been well imaged by InSAR and GPS (e.g., Bü rgmann et al., 2000;Segall and Davis, 1997;Klotz et al., 1999). Additionally, a dynamic component (during the rupture propagation) can now be observed with high-rate GNSS observations (e.g., Larson et al., 2003), providing complementary and in some cases, superior information to seismometers, which can record near-field coseismic displacements on scale. Geodetic measurements have recorded significant post-seismic deformation (in some cases, exceeding the coseismic ground motion, e.g., Heki et al., 1997) starting immediately after the earthquake and lasting from several years to decades in the case of extremely large earthquakes.
The pattern of ground deformation from nearly all earthquakes has been successfully replicated by elastic dislocation sources representing finite slip on a fault plane. These models even seem to work for earthquakes from depths of over 100 km (Peyrat et al., 2006) because the materials still behave elastically during the short duration (seconds) of the seismic rupture. In one case, it was argued that the asymmetric deformation from a 1997 nearly-vertical strike-slip earthquake in Tibet required non-linear elasticity (Peltzer et al., 1999), however the asymmetry can also be fit with a linear elastic model that has a reversal of fault dip (e.g., Funning et al., 2007). At least a few earthquakes, however, show ground deformation related to triggered or induced motion on secondary structures that are weak and compliant Fielding et al., 2004;Hearn and Fialko, 2009). A major topic of research is to use precise geodetic observations to understand the complex spatio-temporal slip during the earthquake event, and how this slip is related to secondary motions, other earthquakes (and aftershocks), and the location and nature of post-seismic deformation.
The geodetic observations are complementary to seismic studies of earthquakesdby combining both, we get a more complete picture of the earthquake location, mechanism, and complexity. For example, Fig. 8a shows that by using InSAR data, it was possible to determine that the global earthquake catalogs mislocate a small earthquake (Mw 5.3) in southern Iran by up to 70 km, probably because of inaccuracies in the estimates in the seismic velocity structure. The precise earthquake location from the InSAR allows the earthquake to be placed in geologic context (Fig. 8b), and helps to assess the sub-surface dip of the fault. In a larger sense, earthquake parameters based on geodetic (and seismic) parameters facilitate comparison with geologic structures and allows assessment of any causal relations (e.g., Lohman et al., 2002;Loveless et al., in press).

Are differences between geodetic and geologic deformation rates real or an artifact of one or both types of measurements?
Several studies have shown that current geodetic (decadal) rates and patterns of deformation are different from deformation over the long-term (thousands to millions of years) measured by geologic mapping, paleoseismology, and seismic imaging (Friedrich et al., 2003Hsu et al., 2003;Jackson, 1999;Oskin and Iriondo, 2004;Oskin et al., 2007;Peltzer et al., 2001;Wright et al., 2004). In particular, these studies show that, in some areas, current deformation cannot explain recent geologic features like mountain ranges and rift zones, or even the direction of faulting. These studies also find that spatial variations in the rates of deformation are importantdfor example, only a subset of faults within a region might be active at any given time if fault motion alternates or migrates within a fault system (e.g., Jackson, 1999;Niemi et al., 2004). Furthermore, temporal variations in fault slip are also important. Fault activity can be clustered (e.g., Wallace, 1977) such that the observed geodetic rate may not be representative of a longterm average depending on whether measurements are made during a ''cluster'' or during a quiescent time (e.g., Oskin et al., 2008). In fact, earthquake clustering may be common, and is expected on theoretical groundsdeven in simple systems with a single earthquake on a single strike-slip fault can affect the frequency of earthquakes on nearby faults via visco-elastic processes (Kenner and Simons, 2005;Meade and Hager, 2005).
While the differences between geologic and geodetic deformation rates might reflect real variations in deformation over space and time, several research groups are also re-investigating whether apparent discrepancies might be resolved by considering earthquake cycle models with temporarily complex behavior. For example, fault slip models often assume that fault slip is occurring in an elastic homogeneous half-space, but the inferred fault slip rate can change if the there are variations in elastic plate thickness (e.g., Ché ry, 2008) or visco-elastic rheology (Johnson et al., 2007;Segall, 2002). In addition, the discrepancy between geodetic and geologic results is motivating new field investigations that find more complex fault interactions than originally suspected (Oskin et al., 2007).
6. Future directions 6.1. Technological advances enabling future directions Several trends in space-based geodesy promise a resolution of some of the ambiguities that plague current interpretations. Continued monitoring and increased densification of existing GNSS networks will provide more uniformly spaced observations which will help us to better evaluate whether a strain rate anomaly is real or simply reflects station clustering in areas of high strain. The on-going switch from campaign to continuous GNSS data will provide higher accuracy velocities with better resolved vertical components that will allow for three-dimensional strain analysis. However, higher errors in the calculation of vertical strain components will remain because the stations only occur on the surface of the earth and are thus very nearly coplanar except in areas of very high relief. The switch to continuous data will have another, more profound result: it will provide time series of deformation that will allow the capture of transients, as described below.
At least four technical innovations will increase the use of InSAR for tectonic studies. A primary limitation of the widespread use of InSAR in the last decades has been the lack of sufficient observations, but this is likely to change within the next decade. Currently, seven different SAR satellites (or constellations of satellites) are in earth orbit: European Space Agency's ERS-2 and Envisat; Japanese Space Agency's ALOS; Canadian Space Agency's Radarsat-1 and -2; German Space Agency's TERRASAR-X; and the Italian Space Agency's COSMOSky-MED. Several more missions are planned in the next few years, with an exciting prospect being NASA's DESDynI satellite which would be the first dedicated to InSAR, and with a likely launch date in 2017. The increased number of satellites will increase the spatial and temporal frequency of observations. While individual SAR scenes from the different satellites cannot (generally) be combined to form interferograms, temporally overlapping interferograms from different satellites (or orbital tracks of a single satellite) can be combined to allow multiple components of the three-dimensional deformation (or velocity) field to be recovered (e.g., Fialko et al., 2005) and provide constraints on the temporal evolution of deformation (e.g., Pritchard and Simons, 2006). With multiple components of the deformation field, it would be possible to do a strain analysis and do other types of modeling in combination with available GPS data. Another limitation for past InSAR studies has been that the observations were restricted to arid and urban regions because vegetation changes masked tectonic changes. The use of longer radar wavelengths (L-band or 23 cm like on the ALOS and DESDynI satellites) allows InSAR measurements in vegetated areas. The development of new computer algorithms like persistent-scatterer InSAR (Ferretti et al., 2001;Hooper et al., 2004) has also permitted InSAR measurements to be made in vegetated areas and the further development and application of these techniques along with increased frequency of SAR observations will further expand regions where InSAR measurements of tectonic deformation can be made.
A particularly frustrating aspect of current analyses of GNSS data is the difficulty of combining data sets collected by different groups of investigators. This difficulty arises from the fact that different research groups use different reference frames. While different universal reference frames exist, many older campaign networks  (Lohman and Simons, 2005) from southern Iran spanning 21 April to 26 May 1999, showing ground deformation from a Mw 5.3 earthquake at 3.2 km depth (red region), as well as the location and focal mechanism from the Global CMT catalog (mislocated in depth by 42 km and horizontally by 67 km) and the location in the ISC earthquake catalog (mislocated by 32 km in depth and 7 km horizontally). (b) Geologic cross section (McQuarrie, 2004) showing location and mechanism for precisely located earthquake by Lohman and Simons (2005) from (a). Black unfilled circles show earthquake located with synthetic seismograms (see references in McQuarrie, 2004), while white circles show microseismicity of Tatar et al. (2004). persist in publishing processed data in their own unique reference frames. Thus overlapping data from different surveys cannot be combined and analyzed together. In South America, for example, the data from the CAP network (Brooks et al., 2003;Kendrick et al., 2001) cannot be combined with the almost completely overlapping GFZ network (Khazaradze and Klotz, 2003;Klotz et al., 2001) and IPGP (Chlieh et al., 2004) networks. While there are mathematical approaches to combining data from different networks (Dong et al., 1998), increased collaboration between groups and open access to the raw observations (for example in the UNAVCO data archive) would ensure that disparate networks can be combined into a common reference frame.

Transients and their relations to processes
For a structural geologist, short-term geodetic measurement can reveal processes that might be valuable when interpreting the geologic record. Many different types of ''transient'' deformation signals involve the cycling of fluids or other dynamic processes that might not be apparent when studying fossilized systems. In addition to the inter-and coseismic periods, spatially and temporally dense geodetic data have revealed a broad spectrum of transient or intermediary behavior that occurs during the seismic cycle, including post-seismic and/or pre-seismic deformation caused by ductile, poro-elastic, frictional process, and fault zone healing (e.g., Bü rgmann and Dresen, 2008;Heki et al., 1997;Peltzer et al., 1996;Fielding et al., 2009), slow earthquakes or silent slip events (e.g., Schwartz and Rokosky, 2007), and partial interseismic coupling on seismogenic faults. Additionally, transitory deformation from glacier mass balance changes (e.g., Hooper and Pedersen, 2007;Pinel et al., 2007) and fluid related processes such as magmatic injection (e.g., Smith et al., 2004;Chang et al., 2007) have been documented. Comparing these behaviors to long-term deformation begs the question of whether transient phenomena leave a distinct permanent mark or are lost among signals of coseismic effects in the geologic record.
Most tectonic transient phenomena induce surface displacements in a direction similar to that of coseismic deformation. For example, the well documented periodic slow slip events on the Cascadia subduction zone interface produce a reversal of the interseismic trend recorded at forearc GNSS stations (Dragert et al., 2001(Dragert et al., , 2004. While earthquakes produce sudden offsets that are documented in stratigraphy, aseismic processes may act over such long time scales and be characterized by such small magnitude deformation that they do not produce a distinct geologic signal in the fault zone itself. Along the Parkfield segment of the San Andreas Fault, Toké et al. (2006 and references therein) note that while aseismic fault creep and earthquakes can produce similar styles of stratigraphic deformation, only earthquakes can result in fissuring, liquefaction, and colluvial wedges, thus introducing some distinction between the two types of slip behavior in the paleoseismic record. Practically, Toké et al. (2006) found only evidence of moderate (Mw6) earthquakes, considering aseismic creep as an unquantified factor in constructing sag ponds and other geomorphic evidence of fault slip.
Post-seismic deformation has been clearly separated from coseismic offsets in the geodetic record by its time dependent nature. Whereas the signal of an earthquake in a geodetic time series appears as a step function, post-seismic displacements often follow a non-linear pattern for some time period after the earthquake. Because the time scale of the decay of post-seismic deformation is typically on the order of years to decades, even for the greatest of earthquakes, it is difficult to separate it from the coseismic signal in the geologic record. Instead, the combined effects of co-and post-seismic deformation may be interpreted as resulting from an earthquake with an overestimated magnitude. For example, post-seismic GPS displacements following the 1994 Mw ¼ 7.4 Sanriku-oki earthquake offshore northern Japan have been interpreted as the result of after-slip on the subduction interface, the equivalent magnitude of which exceeds that of the earthquake (Heki et al., 1997). An idealized geologic record, recording 100% of the earthquake and after-slip, would be interpreted as a single event exceeding the assigned earthquake magnitude.
Transient processes have been discovered by using dense time series at a given location but have also been discovered when looking at InSAR snapshots of the spatial extent of deformation across a region. For example, Fialko et al. (2002) discovered motion on a series of faults 10 s of kilometers away from the main fault ruptures (indicating the characteristics of fault damage zones).

Conclusions
Space-based geodetic data have considerable potential to revolutionize our understanding of active structures. Many of the questions surrounding the interpretation of these datadis most of the deformation concentrated on networks of large faults separated by relatively stiff elastic blocks? How does spatial sample distribution affect our perception of deformation? etc.dhave direct parallels for structural geologists in the context of finite strain: how many thrust faults have to be included to capture the strain in a thrust belt? How do you measure strain in the middle of a heterogeneous simple shear zone? In one sense, however, there is a profound difference: geodetic data truly capture a geologically instantaneous snapshot in time. Because geodetic data span much less time than an individual earthquake cycle, we must use data from different locations as a proxy for time. That is, we are trying to understand the complete fault behavior by piecing together interseismic deformation from one fault with coseismic and postseismic behaviors from other faults in other areas. In geological study of deformed rocks, we can occasionally eke out the progression of events in a single area. That geologic record, however, may never fully record the transient events that are emerging from geodetic data. Likewise, the geodetic record is still too young to capture faults that may have had an important role in Cenozoic mountain building, but have ceased their activity, or are temporarily quiescent.