Characterizing Complex Surface Ruptures in the 2013 Mw 7.7 Balochistan Earthquake Using Three‐Dimensional Displacements

We use satellite‐derived high‐resolution topography and orthoimages, namely, preearthquake Advanced Land Observing Satellite World 3‐D data and postearthquake Pleiades data, to retrieve 3‐D displacements in the 2013 Balochistan, Pakistan, earthquake. Previous studies of this earthquake have revealed many complex rupture patterns, such as off‐fault deformation (near‐field strike‐slip displacements smaller than the far field) and inelastic surface deformation (horizontal shortening on the surface significantly larger than that from simple elastic model). In this paper, we reanalyze the complexities of surface ruptures in a systematic way using the newly derived 3‐D displacements. In doing so, we made use of the vertical component of the curl and the horizontal divergence of the displacement field. By comparing the off‐fault deformation and curl field along the rupture, we found that curl is a good measure of the width of the deformation zone. The curl width ratio (CWR)—the ratio of the basic resolution of the curl field and curl width—was used to quantify the degree of surface slip localization. When CWR ≥ 0.9, there is no or very little off‐fault deformation, whereas when CWR ≤ 0.6, surface deformation is almost 100% distributed. The distributed deformation observed is controlled by both the fault geometry and local material types. Despite the overall strike slip with some thrusting, the divergence shows localized extension or enhanced shortening in the near field due to fault geometric variations at a spatial scale of tens to hundreds of meters, consistent with the 3‐D displacements and geological interpretations.

InSAR can give 3-D displacements (Fialko et al., 2005;Wright et al., 2004) but often no near-fault (within hundreds of meters from the primary rupture trace) estimates in large earthquakes due to steep phase gradients and intense ground shaking. Optical image (e.g., Landsat-7/8 images) matching alone can measure near-fault displacements but is commonly restricted to the horizontal components of displacement. Differences are often seen between these two types of measurement. For example, Fialko et al. (2005) studied a few large strike-slip earthquakes using InSAR and observed that slip gradient decreases toward the surface (shallow slip deficit), whereas Xu et al. (2016) reanalyzed the same earthquakes using optical image matching but did not observe this phenomenon. Far-field (more than 1 km away from the primary rupture) deformation determined by InSAR is useful for resolving patterns of slip at depth, and near-field deformation help us constrain shallow slip distribution.
To improve our understanding of earthquake physics and hazards, we need to be able to get the 3-D displacements including both the far field and near field. Matching and differencing light detection and ranging (LiDAR) topography has been a valuable tool for determining 3-D displacements in earthquakes (e.g., Glennie et al., 2014;Nissen et al., 2012Nissen et al., , 2014Oskin et al., 2012;Scott et al., 2018), but its relatively high cost and low availability has limited its applicability, particularly as preearthquake LiDAR acquisitions are rare Figure 1. Regional map of the 24 September 2013 M w 7.7 Balochistan, Pakistan, earthquake. Active faults in the eastern Makran are mapped in black. USGS focal mechanism for the main shock is shown in blue. N-S, E-W, and vertical components of surface displacements are derived from correlating and differencing the 5-m Advanced Land Observing Satellite and Pleiades (downsampled) DEMs. The elongate red polygon denotes the extent of the preearthquake and postearthquake DEMs (5 km wide × 240 km long). Positive displacements represent northward, eastward, and upward motion. Blue squares along the rupture indicate the locations of the corresponding figures. USGS = U.S. Geological Survey; DEMs = digital elevation models. . New global topographic data sets that were released in 2016 have opened up new opportunities for mapping surface deformation in 3-D. These new global data sets will allow comprehensive analyses for all major continental earthquakes in suitable terrains from now onward.
In this paper, we use satellite-derived topography, that is, the preearthquake Advanced Land Observing Satellite (ALOS) World 3-D data set and a postearthquake Pleiades digital elevation model (DEM), to retrieve 3-D displacement fields in the 2013 M w 7.7 Balochistan, Pakistan, earthquake. Previous studies of this earthquake (Avouac et al., 2014;Gold et al., 2015;Jolivet et al., 2014;Vallage et al., 2015Vallage et al., , 2016Zinke et al., 2014;Zhou, Walker, Elliott, & Parsons, 2016) have analyzed the surface rupture in great detail using various sources of imagery, but most of them focused mainly on the horizontal motion, except  where they directly measured vertical offsets based on postearthquake Pleiades topography for the whole of the rupture. Here we combine measurements of horizontal and vertical displacements to quantify fault geometric variations at shallow depth and show how they can affect surface slip at the fault. Another aim of this paper is to demonstrate simple ways of using high-resolution displacement fields to characterize complex surface deformation in earthquakes.

Data
A limited number of global high-resolution topographic data sets are now becoming available recently, potentially enabling us to measure 3-D earthquake deformation worldwide. One of our goals in this paper is to assess the capability of these data sets and discuss how these new satellite-derived observations of 3-D coseismic displacements can help characterize earthquake ruptures.
We use the 2013 Balochistan, Pakistan, earthquake (Figure 1) as a case example. To retrieve 3-D surface deformation in earthquakes, both preevent and postevent topography are needed. The 1-arc sec Shuttle Radar Topography Mission (SRTM) DEM (Farr et al., 2007), the 30-m Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) GDEM (ASTER, 2009), the 30-m ALOS World 3-D DEM (referred to as the ALOS DEM hereafter), and the 5-m DEMs based on the same ALOS Panchromatic Remote sensing Instrument for Stereo Mapping (PRISM) images (Tadono et al., 2014) were all constructed from data acquired before the 2013 earthquake. For the purpose of comparison (see Figure 2), we derived 3-D displacement fields using each of these four data sets in turn as preevent topography and the Pleiades DEM from  as postearthquake topography. Note that all the 30-m DEMs are freely available, but the 5-m ALOS DEM must be purchased to use under license; we also purchased the 2.5-m ALOS orthoimagery.
The global (covering 56 ∘ S to 60 ∘ N) 1-arc sec SRTM DEM was released in 2014. Before then, only 3-arc sec data were available outside the United States. The vertical error of the 1-arc sec DEM is less than 16 m 10.1029/2018JB016043 They are artifacts in the ALOS DEM, caused by using integer values in transforming the data from ellipsoid heights to orthometric heights. The artifacts get lost in the hanging wall (west of the fault) due to noises in the vertical displacement field. (b) The jump is ∼1 m as shown by our swath profile (black polygon). ALOS = Advanced Land Observing Satellite; DEM = digital elevation model. (Farr et al., 2007). The current ASTER GDEM (covering 83 ∘ S to 83 ∘ N) was released by National Aeronautics and Space Administration and Japan's Ministry of Economy, Trade and Industry in 2011. The vertical accuracy of GDEM is claimed to be 20 m at 95% confidence (ASTER, 2009). The ALOS DEM (covering 82 ∘ S to 82 ∘ N) was generated by the Japan Aerospace Exploration Agency using the archived data (2.5-m resolution stereo images) of the PRISM that were acquired from 2006 to 2011 (Tadono et al., 2014). The nominal relative height accuracy of the ALOS DEM is 5 m (1 ; Tadono et al., 2014). The Pleiades DEM (1-m resolution) was produced from 13 stereo data sets of 0.5-m nominal resolution, panchromatic Pleiades images acquired by Airbus. The images were processed using the LPS module of the ERDAS Imagine 2013 software (version 13.00.00, Build 281). A pixel-by-pixel matching procedure was implemented with a window size of 5-by-5 pixels and a correlation coefficient of 0.3 to 0.7 . The relative vertical accuracy of Pleiades-derived DEMs was assessed to be 0.3 m .
When we differenced the ALOS (both the 30-and 5-m versions) and Pleiades DEMs, we found some linear artifacts in the ALOS DEM. As shown in Figure 3, we noticed some linear steps in the DEM difference map about every 20-25 km. These lines do not follow any geological structures, so they are unlikely to be associated with earthquake deformation. We estimated a uniform ∼1-m elevation change across the lines. After reexamining the data and discussions with NTT DATA (the suppliers of the ALOS data), the cause was identified as due to the choice of number representation used in constructing the ALOS DEMs. When the 30-m ALOS DEM is first derived from satellite data, the height estimated is the ellipsoidal height, the vertical distance from a reference ellipsoid (i.e., a mathematical model of the Earth). What is required is the orthometric height or elevation measured relative to the geoid in the direction of gravity at the Earth's surface. The DEMs must be converted from ellipsoidal height to elevation using the known distance between the geoid and the reference ellipsoid. The 30-m global DEM used integer values in making the conversion, in which case the geoid heights were stored in integer format; hence, the correction was rounded, resulting in the 1-m jumps as its magnitude varies. The problem can be fixed by specifying the data type as float when purchasing the 5-m product, but the freely available 30-m ALOS World 3-D data retain this artifact.

Three-Dimensional Coseismic Displacement Fields of the 2013 M w 7.7 Balochistan Earthquake
Previous studies of the 2013 Balochistan, Pakistan, earthquake have used various types of optical imagery, including Landsat-8 (Avouac et al., 2014;Jolivet et al., 2014;Zinke et al., 2014), SPOT-5 (Vallage et al., 2015, Pleiades , and WorldView (Barnhart et al., 2014;Gold et al., 2015), to study the 2-D coseismic deformation patterns in detail. Here we show how the complete 3-D deformation field, that is, adding the vertical component of displacement, can contribute to our understanding of earthquake deformation.

Overview of the Balochistan Earthquake
The 2013 M w 7.7 Balochistan, Pakistan, earthquake ruptured a 200-km section of the arcuate and moderately dipping (50-70 ∘ N) Hoshab fault in the eastern Makran accretionary wedge (Avouac et al., 2014;Jolivet et al., 2014; Figure 1). The rupture was bilateral but propagated much further to the south of the epicenter than to the north where fault creep (∼3.4 mm/year) has released a large fraction of the accumulated stress (Fattahi & Amelung, 2016). The southward propagation ended in a relay stepover (Zhou, Walker, Elliott, & Parsons, 2016). Previous studies have measured an average of ∼8-m left-lateral strike-slip (Avouac et al., 2014) and ∼1.6-m vertical  motion during the earthquake. Although the surface trace of the 2013 earthquake appears to be relatively continuous on a spatial scale of tens of kilometers (despite two were measured by fitting parallel lines to topographic profiles along rivers and alluvial fan surfaces . FF offsets (black line, 1 error) were measured from the vertical displacement field (see Figure 2e for an example). In general, the NF and FF vertical offsets are consistent. Nonetheless, we found a number of places where the NF and FF vertical offsets are in the opposite sense (e.g., Figure 5) and also places where the NF vertical offsets are larger than the FF offsets (e.g., Figure 6). Red boxes show the locations of corresponding figures. small restraining stepovers), high-resolution displacements along the surface ruptures have revealed some complex deformation patterns.
The first notable feature as pointed out by Zinke et al. (2014) and Gold et al. (2015) is the significant off-fault deformation (Dolan & Haravitch, 2014) in the fault-parallel, strike-slip motion. Off-fault deformation is a phenomenon where near-field (0-10 m from the main fault) strike-slip displacements are smaller than the far field (more than a few hundred of meters from the main fault). Zinke et al. (2014) measured offset geomorphic features, that is, near-field strike-slip offsets, and compared the offsets to the Landsat-8 matching results, that is, far-field strike-slip offsets. They found that on average, the near-field offsets were ∼45% smaller than the far field. Gold et al. (2015) reanalyzed the near-and far-field strike-slip offsets using WorldView images and estimated the discrepancy to be ∼28%, implying that earthquake deformation is not always concentrated on the fault. It has been suggested that whether or not earthquake deformation is concentrated on the fault is related to several factors, such as the immaturity of the fault (Dolan & Haravitch, 2014), and the type of material at the surface through which the fault cuts Zinke et al., 2014).
Another interesting feature is that at some places, horizontal shortening at the surface is much larger than that predicted from simple elastic models. Using the surface displacements derived from SPOT-5 imagery, Vallage et al. (2015) showed that the far-field fault-normal surface slip (1-2 m) may be significantly smaller than the near field (reaching 6-7 m). They described such feature as inelastic surface deformation, possibly caused by subsurface flattening of the fault (Vallage et al., 2015).
In the following sections, we use our 3-D displacement fields to characterize these complex surface deformation patterns in a systematic way.

Three-Dimensional Surface Deformation Derived From Preearthquake and Postearthquake DEMs
To derive the 3-D displacement field, we first downsampled the postearthquake Pleiades DEM from 1 m to match the spatial resolution of the preearthquake DEMs. We then correlated the preearthquake and postearthquake DEMs using the COSI-Corr software package (Ayoub et al., 2009;Leprince et al., 2007; to estimate E-W and N-S displacements. The correlation window size was set to be 32 × 32 pixels with a step of 8 pixels. The horizontal displacement field contains holes (Figures 1 and 2) due to filtering out noisy displacements. The holes were filled in using bilinear interpolation. The processed horizontal displacements were then used to correct for the horizontal advection of topography between the DEMs. We differenced the corrected DEMs to obtain unbiased measurements of vertical displacements (Figures 1 and 2).
Another way to determine 3-D displacements using preearthquake and postearthquake topography is iterative closest point (ICP). This method is designed to process dense and accurate LiDAR point clouds. Nissen et al. (2012Nissen et al. ( , 2014, Glennie et al. (2014), and Scott et al. (2018) obtained good 3-D deformation fields based on ICP. In this study, the ALOS DEM that we purchased is in raster format and has a much lower resolution and vertical accuracy than LiDAR data (5 m versus tens of centimeters). Hence, instead of converting the ALOS DEM to point clouds and using ICP, we use a simple cross-correlation method.
The resulting displacement fields from the ALOS DEMs delineate the 2013 rupture very clearly, but the 2013 rupture is much less clear in the displacement fields derived from the ASTER and SRTM DEMs ( Figure 2). Swath profiles were drawn to estimate the offsets across the fault (Figures 2c-2e). Results from the 30-m ALOS DEM Cartoon illustrating how the geometric variation controls the near-field deformation. We projected the horizontal displacements onto the fault-parallel (c) and fault-perpendicular (e) directions assuming a mean strike of 135 ∘ . (d) The fault-parallel offset along swath profile S1-S1 ′ (80 m wide), that is, the far-field strike-slip component, was estimated to be 7 m. We did not find any evident near-field strike-slip offsets, suggesting that the near-field offset is much smaller than the far field. (e) Fault-perpendicular displacement field shows an area of localized extension at the bend (marked by the black box). (f ) The displacement profile shows shortening (1.1 m) in the far field and localized extension (1.5 m) in the near field. (g) A map of vertical displacements at this site. (h) The FF vertical offset at the bend was estimated to be ∼1.7 m (west side up corresponding to thrusting). The NF vertical motion was measured to be 3 m (west side down corresponding to normal faulting; see Figure 4 for location). FF = far-field; NF = near-field.

10.1029/2018JB016043
are consistent with those from the 5-m ALOS DEM, though much smoother due to its coarser spatial resolution. The ASTER and SRTM DEMs yielded very noisy estimates, suggesting that their vertical accuracy is not as good as the 30-m ALOS DEM with comparable resolution.
Given that we also have both the ALOS and Pleiades orthoimages with higher resolution, we then calculated the horizontal displacement fields by matching the 2.5-m ALOS and Pleiades orthoimagery. We see stripes about every 60 m in these displacement fields (like those shown in Figures 5c and 5e); they are due to the misalignment of the charge-coupled device chips inside the ALOS PRISM sensor. The resulting periodic effect across the stripes is ≤1.2 m. We made no attempt to remove this noise by filtering in order to avoid smoothing the data, and we did not use them to correct for the horizontal advection of topography. Zinke et al. (2014) and Gold et al. (2015) have compared near-field (0-10 m from the primary surface traces) and far-field (>350 m from the primary surface traces) fault-parallel offsets across the fault and found some differences between them (near-field fault-parallel offsets smaller than the far field). Here we reanalyze the near-and far-field offsets by including our vertical measurements. The near-field vertical offsets were measured by fitting parallel lines to topographic profiles along rivers and alluvial fan surfaces within ±200 m from the fault . The far-field vertical offsets were measured by taking swath profiles (1.5 km wide and 4 km long) from the vertical displacement field (e.g., Figure 2e). We determine linear fits to the far-field data on either side of the fault and extrapolate to the fault, where the measurement is taken (following, e.g., Milliner et al., 2016). The uncertainties were calculated based on the fitting errors. The measurements are plotted in Figure 4. In general, the near-and far-field vertical offsets agree well, both showing an average vertical motion of 1.5 m (west side up due to thrusting). However, at some places, we do notice differences between the near-and far-field measurements.

Comparing Near-Field and Far-Field Displacements
One example is shown in Figure 5, at the location along the rupture marked in Figure 4, where the far-field vertical offset is west side up given the overall thrusting, but the near-field vertical offset is west side down.
Here the strike-slip motion has been turned into localized normal faulting by a bend in the fault, the jog to the left producing a mini basin (Figure 5b). To confirm this interpretation, we projected the E-W and N-S displacements onto the fault-parallel and fault-perpendicular directions (Figures 5c-5f ). The far-field strike-slip offset (7 m) is larger than the near field (we did not find any evident laterally offset features across the surface ruptures), indicating decreasing strike-slip motion at the fault. The strike-slip offset is distributed over about 200 m, a distance corresponding to the shadowing effect of the fault jog along strike. The fault-perpendicular motion reveals extension in the near field, despite the overall shortening. These observations are consistent, together showing transformation between strike-slip and vertical motion in the near field as a result of geometric variations in the fault trace. The west-side-down motion is localized within 500-600 m from the fault, implying that the deformation associated with this feature is confined to shallow depths.
Conversely, at some restraining bends, the near-field vertical motion is larger than the far field. This is because strike-slip motion has been turned into localized shortening, thus adding extra thrusting to the overall deformation at the fault. Figure 6 shows an example of such geometric variation. At the bend, as shown in Figure 6a, the far-field strike-slip offset is estimated to be 9 m. We measured ∼4 m offset at the fault (also in Zinke et al., 2014), suggesting that the near-field fault-parallel slip is smaller than that in the far field at the bend. The fault-perpendicular displacement profile (Figure 6f ) reveals 4.2-m shortening in the near field, twice as much as the far field (2.1 m), implying extra thrusting and hence more vertical motion at the bend. The vertical displacement profile (Figure 6h) shows a larger near-field vertical offset (3.5 m) than in the far field (1.8 m).
Pull-apart basins and elevated topography are well-known long-term expressions of, respectively, left-stepping and right-stepping jogs on a left-lateral strike-slip fault. The examples in Figures 5 and 6 show the processes giving rise to such features as they occur. Both examples clearly show interactions between the horizontal and vertical motion on the fault in response to shallow geometric variations, emphasizing the importance of analyzing surface deformation in 3-D.

Curl and Divergence
To analyze the displacement field (F) and deformation patterns along the fault trace, we computed the vertical component of the curl ( ) and the horizontal divergence ( ): where u and v are the displacements in the x (east) and y (north) directions, respectively. Positive corresponds to anticlockwise rotation viewed from above ( conforms to the right-hand rule); positive indicates an increase in area, that is, net extension, and negative indicates a decrease, that is, net shortening.
Since the horizontal displacement fields from matching the DEMs have gaps in flat areas with few topographic features, to calculate and , we used the E-W and N-S displacements derived from matching the 2.5-m ALOS and Pleiades orthoimagery. As can be seen in Figure 7, which shows and for the examples in Figures 5 and  6, delineates the fault trace clearly (a similar effect has been observed in the rotational component of ICP deformation fields, e.g., Nissen et al., 2014); left-lateral strike-slip results in anticlockwise rotation and hence positive at the fault. It provides a simple means of visualising and analyzing the surface ruptures. Note that the noise level in and is ∼0.02, mainly due to the stripes in the horizontal displacement field.
We also found some interesting deformation patterns at the fault from the divergence. It varies along the rupture, and, despite the overall shortening, does not remain negative, but at some locations becomes positive, suggesting that extension occurred at the fault locally. The change in (extension) in Figure 7d is consistent with the change in the near-field vertical motion (west side down) as shown in Figure 5, both of which are caused by fault geometric variations at subsurface. Figure 7j reveals extra shortening at the fault bend, which is also consistent with our geological interpretations and observations as presented in Figure 6.
There are differences between the distribution of deformation in the two cases. Where there is additional shortening, this can be easily taken up by thrusting at the fault; that is, the shortening remains concentrated at the fault (Figure 7j). In the example where extension occurs at the fault, because overall there is net shortening, that shortening must occur off the main fault. This can be seen in Figure 7e to the left of the location where the extension occurs.

Curl Width
The curl ( ) is determined by the surface motion at the fault. We know that the surface rupture in the 2013 earthquake varies from a single, linear trace (localized deformation) to a wide zone of distributed cracks (distributed deformation). It is of interest to explore whether can represent such variations. Figure 8 shows three typical scenarios of surface ruptures. For simplicity, we consider pure left-lateral strike-slip motion. In the first scenario ( Figure 8a) where surface slip is entirely localized on the fault, we would expect a single surface rupture. In other cases, the surface deformation may be partly distributed (Figure 8f ) or 100% distributed in the extreme case (i.e., no detectable displacements on faults within the error, as shown in Figure 8k). To see how the curl varies accordingly, we extracted profiles (swath profiles through the curl data with a width of 200 m) along the fault. Here we present one example for each of the scenarios.
In the first example (Figures 8b-8e), we see a clear linear trace from the orthoimage and the strike-slip displacement field. Localized strike slip will result in maximum rotation at the rupture (i.e., maximum ). Figure 8d shows anticlockwise rotation along a single, linear rupture. It increases from 0 to the maximum within a very narrow zone (Figure 8e). In the limit of deformation concentrated on the fault, the width would be 0, but because our displacement field only has a resolution of 20 m, represents the deformation within tens of meters. To quantify the width of the deformation zone (hereafter referred to as the curl width), we use the distance between the points where is 50% of maximum (see Figure 8j for illustration). The curl width is 43 m in this first case. Note that the curl width cannot reveal the level of localization below this limit that is set by the image resolution.
Figures 8g-8j present an example of the second scenario. This site is also shown in Gold et al. (2015) as an example of off-fault deformation. At this site (see Figure 6 in Gold et al., 2015), the far-field strike-slip offset is 7.2 m, whereas the near-field strike-slip offset is 4.7 m, 34% smaller. Because the deformation is only partly distributed, we can still see a peak in the curl profile (Figure 8j), corresponding to the main fault on which slip is localized. However, the curl width (90 m in this example) is larger than that with 100% localized slip, indicating partly distributed deformation.
In the extreme case (Figures 8l-8o), the deformation is accommodated in the form of extensional or compressional cracks at the surface without any detectable displacements on faults. At the site in Figure 8l, we can map many small surface cracks, but none of them shows significant offsets within the resolution (0.5 m) , the orthoimage, strike-slip displacement field (calculated using a mean strike of 160 ∘ ), and curl reveal a single, linear rupture at the surface. We plotted along swath profile C1-C1 ′ (200 m wide) in (e). The peaks at the fault; the width of the zone (i.e., curl width) defined by points where is 50% of maximum is 43 m. (f-j) Surface slip that is partly localized and partly distributed. The site in (g) is also shown in Gold et al. (2015) as an example of off-fault deformation. At this site, the far-field strike-slip offset is 7.2 m, whereas the near-field strike-slip offset is 4.7 m, 34% smaller. In the case of partly localized surface deformation, we would expect some anticlockwise rotation within a wider zone. As shown in (i) and (j), indicates deformation within a zone of 90 m, larger than that with 100% localized deformation. (k-o) The extreme case where surface deformation is 100% distributed; that is, there are no detectable offsets at faults within the image resolution. In this scenario, we would expect ≈ 0 within the measurement error. (l-o) An example of such scenario. From visual analysis of the 0.5-m Pleiades imagery, we can map many discontinuous surface cracks (red lines) but no significant offsets, suggesting distributed deformation. (m) The fault-parallel (strike-slip) displacement field also shows that the deformation occurs within a wide zone (a gradual change in slip across the fault). Unlike (e) and (j) where peaks at the main fault, is almost 0 within error along profile C3-C3 ′ . There appear to be a central feature (∼300 m wide) in the profile, but the peak values follow the stripes. We therefore treat this feature as noise. Even if it was real signal, the width would still be much larger than those in the first two scenarios. Green dashed lines show the noise level (profiles were extracted >2 km away from the fault). of the Pleiades imagery, suggesting almost 100% distributed deformation. In this scenario, we did not see a clear peak in the profile, but instead, we found ≈ 0 within the measurement error.
The three examples show that the curl width does correlate with the surface deformation patterns (i.e., localized and distributed surface slip). In section 4.2, we discuss the variation in curl width along the entire 2013 Balochistan rupture.

Comparing Surface Displacements Derived From Different Sources
The surface displacements in the 2013 Balochistan earthquake have been measured using various types of images in this and previous studies (Avouac et al., 2014;Barnhart et al., 2015;Jolivet et al., 2014;Vallage et al., 2015;Zinke et al., 2014), and the question of the consistency between the different results is a pertinent one. We plotted all the measurements in Figure 9 for a quantitative comparison. The measurements derived from SPOT-5, Landsat-8, ALOS, and Pleiades data are in good agreement in general, regardless of the data type (optical imagery and DEMs) and resolution (ranging from 0.5 to 15 m). Any inconsistencies could result from data quality (e.g., noise level and orthorectification error in the case of optical imagery) or different subpixel matching algorithms, for example, Zinke et al. (2014), and we used COSI-Corr (Leprince et al., 2007), Vallage et al. (2015) used MicMac (Rosu et al., 2015), and Barnhart et al. (2015) used ROI_PAC (Rosen et al., 2004).
The Landsat-8 result (Avouac et al., 2014) is smoother due to its coarser resolution (240 m). The WorldView-derived estimates appear to be more scattered compared to all the other estimates. This could result from its orthorectification since the 0.5-m WorldView images were orthorectified with a much coarser DEM (the 30-m SRTM DEM). Topographic errors can introduce large artifacts (a few to tens of meters) in apparent horizontal displacements, depending on the incidence angle, flight height, topography, and ground displacements (Barišin et al., 2015;Hollingsworth et al., 2012;Zhou, Walker, Hollingsworth, et al., 2016). As shown in Figure 10, the WorldView-derived horizontal displacements are strongly correlated with topography. For example, there should not be any displacement along profile A-A ′ in Figure 10, but the WorldView estimates show an apparent E-W offset of 3.5-4 m caused by the 4-m topographic relief across the alluvial fan that is not visible in the SRTM DEM (Figure 10c). Such apparent displacements may be even larger across the fault with greater topographic relief.
It is necessary to use orthorectified images for matching, particularly when using high-resolution imagery such as Pleiades and WorldView. However, many studies that use optical images to retrieve earthquake surface deformation ignore orthorectification errors, but, as shown in our example, this effect may be significant. Image matching based on DEMs is free from orthorectification errors, but it would be challenging in areas of flat topography with few topographic features (similar to image contrast in the case of optical imagery).

Quantifying the Degree of Surface Slip Localization
As we mentioned in section 3.1, studies of the 2013 Balochistan earthquake Zinke et al., 2014) have found significant off-fault deformation; that is, deformation is sometimes distributed within a wide zone (off fault) rather than localized on the fault, and such distributed deformation may be too small to be seen from imagery or in the field. This phenomenon has been observed in many other earthquakes as well, for example, Rockwell et al. (2002), Oskin et al. (2007), Dolan and Haravitch (2014), .
In section 3.5, we illustrated the use of curl width as a measure of distributed deformation (Figure 8). When the rupture is along a single, linear trace, the curl width is ∼40 m (ranging from 40 to 45 m) set by the basic resolution of the displacement fields. We define 40 m as the localization width, meaning that the surface slip is 100% localized when curl width is 40 m. In cases where there were multiple surface ruptures at a site, the curl width was much larger than 40 m (e.g., Figures 8g-8j). We computed the ratio of the localization width to the curl width, that is, curl width ratio, CWR, and compared it to the off-fault deformation from Zinke et al. (2014). We see a reasonably good correlation between CWR and the difference between the near-and far-field displacements ( Figure 11) with the exception of locations near the sites in Figure 8g where the CWR is ∼0.5, suggesting wide fault zone, but the near-field and far-field measurements of Zinke et al. (2014) agree. We have checked their near-field estimates and believe that this particular site may have been misreported or misinterpreted . Otherwise, when CWR ≥ 0.9, the near-and far-field strike-slip offsets are in good agreement, suggesting no or very little off-fault deformation. On the other hand, when CWR < 0.8, the near-field strike-slip offsets become increasingly smaller than the far-field measurements. We also noticed that in extreme cases where strike-slip offsets are almost invisible in the 0.5-m Pleiades and WorldView imagery, CWR is less than 0.6.
For comparison, we also defined a deformation width using fault-parallel displacement profiles along the length of the rupture such as those illustrated in Figures 5d and 6d. The deformation width was defined as the distance between the main breaks in slope along the profiles. This deformation width is compared in Figure 12a with the deformation width obtained by Gold et al. (2015, their Figure 3d), which was the width of the zone encompassing 2013 surface ruptures observable in their WorldView imagery. Over the distance range 30-120 km the two measurements are in reasonable agreement. The deformation width obtained from the displacement profiles here is limited by the 20-m resolution of the displacement fields. Where our deformation width is at this limit (e.g., 0-20 km), and that in Gold et al. (2015) is smaller, this suggests that faulting is highly localized in these places. However, at the eastern end of the rupture (170-200 km), the deformation width from the displacement fields is well above the resolution limit and is considerable larger than the width of surface ruptures from Gold et al. (2015). This suggests a lack of observable surface ruptures along this part of the fault due to the diffuse deformation and small fault-parallel displacements. The CWR at the site in Figure 8g seems to be inconsistent with the near-field estimates from Zinke et al. (2014). We have checked their near-field estimate and believe that this particular site may have been misreported or misinterpreted. FF = far-field; NF = near-field.
In Figure 12b, we compare the deformation width from the displacement profiles to the curl width. Both measures are subject to the same resolution limit of the displacement fields. It is clear that there are quite a few locations where the curl width is close to or equal to the resolution difference, but the deformation width from the displacement profiles is considerably larger. The reasons for these differences will become clear in the next two sections, where we take a closer look at the displacement fields for a few additional examples.

Geometric Variations Controlling Surface Slip: Areas of Local Extension
Although it has been widely observed that large-scale (tens to hundreds of kilometers) fault geometric variations can play an important role in accommodating fault motion, small (tens to hundreds of meters) features due to geometric complexities are often only well mapped by field geologists (e.g., mole tracks, pressure ridges, tension fissures, sags, and ponds). In addition, these geometric complexities are seldom compared with satellite measurements of surface slip. Our study of the 2013 Balochistan earthquake shows that local geometric variations can dominate surface slip at the fault, the vertical motion in particular. As demonstrated in Figures 5-7, high-resolution 3-D displacements provide a systematic means for analyzing the effects of complex geometric variations.
We calculated the divergence ( ) along the length of the 2013 rupture, in the same way as was done in Figures 7e and 7j, and plotted the peak value within the fault zone in Figure 13 as a function of distance.
We then compared to the near-field vertical motion. The 2013 earthquake is a strike-slip event with some thrust motion. As shown in Figure 13, the divergence mostly shows shortening at the fault, with the vertical motion west side up, conforming to the fault geometry (dipping to the northwest). Nonetheless, both the vertical offsets and divergence exhibit considerable variation along the rupture. The most notable feature is that where the divergence changes sign, that is, the fault-normal motion at the rupture changes from shortening to extension, the near-field vertical offset changes its sense of motion from west side up to east side up. As  Figure 5 showed, this is a result of part of the strike-slip motion being turned into localized extension in the near field in response to shallow geometric variations in the fault trace, that is, slip readjustment.
The example in Figure 5 has a CWR of about 0.5, and we explained the distributed deformation here as a shadowing effect of the extensional jog. However, we note that there is another area of local extension and subsidence at a distance of 90 km along the rupture (Figure 13), just to the east of the main restraining stepover, which apparently has a narrow zone of deformation with a CWR close to 1. To understand the difference between the two examples of local extension, we looked at the displacement fields and curl and divergence profiles of this new example ( Figure 14). From the fault-parallel displacement profile, the deformation zone is indeed about 350 m in total, but the strain is concentrated in two narrow zones, the main rupture and secondary faulting to its north, with an area of low strain in between. In this instance, the method of calculating curl width picks out the peak in curl on the main rupture, resulting in a low value of curl width. It would appear that the region of secondary faulting may take up part of the overall shortening that must occur, although the values of divergence here are for the most part not far from the noise level.
There are two differences between the examples in Figures 5 and 14. In the latter, the jog over seems to occur over a longer distance along strike and is located within or at the edge of consolidated material, whereas the former appears to be located within alluvial material. Off-fault deformation has been inferred to be controlled by factors including fault structural maturity (Dolan & Haravitch, 2014) and the type of surface material Zinke et al., 2014). The examples shown here suggest that off-fault deformation at a stepover may occur in different types of material but that material properties influence the style of distributed faulting.

Geometric Variations Controlling Surface Slip: Areas of Shortening
The inelastic surface deformation described by Vallage et al. (2015) is another example of slip readjustment. They found an area of geometric complexity on the fault where the displacement field rotated and near-field   fault-perpendicular displacements were significantly larger than predicted by simple elastic modeling. To illustrate this example, we present the 3-D displacement fields for this area derived from our ALOS and Pleiades data, together with the associated curl and divergence, in Figure 15. The surface rupture mapped from the Pleiades orthoimagery shows evident geometric variations. To estimate the strike-slip and shortening, we projected the E-W and N-S displacements onto the fault-parallel and fault-perpendicular directions using a mean strike of 139 ∘ . Although the far-field strike-slip offset is large (11 m along profile I-I ′ ), we did not find any offsets at the fault, suggesting significant off-fault deformation, which is also confirmed by the CWR (0.2 in Figure 15j). Figure 15f shows that the near-field shortening (4.7 m) is much larger than the far field (1 m), that is, inelastic surface deformation. As shown in Figure 15l, the large divergence also indicates extra shortening at this site. Extra shortening in the near field results in more uplift at the fault (4.5 m) than the far field (2.6 m; Figure 15h). Using the near-field shortening (4.7 m) and uplift (4.5 m), we can estimate the fault dip at shallow depths (arctan = 69 ∘ ). The dip estimate at greater depths agree with the inversion result (60-70 ∘ ) of global low-frequency W-phase waveforms (Jolivet et al., 2014). The dip estimates suggest that the fault flattens when it comes to the surface, supporting the hypothesis in Vallage et al. (2015). This change in fault geometry may be the cause of the distributed deformation zone at the surface.
A second example of diffuse shortening and irregular fault geometry is that shown in Figure 8l. The complete set of displacement fields and profiles, together with curl and divergence profiles for this location, is shown in Figure 16. In this case, too, shear strain in the fault-parallel component is greater at both edges of the zone of deformation, as shown in the curl profile (Figure 16j). Shortening occurs over a wider zone than the strike-slip component, continuing south of the main rupture (Figure 16f ). We show an enlargement of the area near the intersection of the profile and the fault rupture in Figure 16b. At the edge of the zone of shortening on the profile (point C) we identified a thrust fault running for a short distance across the alluvial fan. From the map of fault-normal displacements (Figure 16e), we see that the shortening south of the main rupture is quite limited in lateral extent, suggesting that the shortening south of the fault is confined to thrusting within the alluvial fan. The shortening here is quite diffuse, with divergence values mostly within the noise and no concentrated shortening at the fault. However, points A and C correspond to two places where the divergence is just outside the noise bound.

Distribution and Causes of Fault Deformation for the Balochistan Earthquake
In the example in Figure 15 discussed in the previous section, the change in dip must occur within 1 km of the main rupture, that is, at depths less than 1 km. This implies that distributed deformation occurs in the top layer at the front of the thrust from 0 to 1 km in depth, that is, with an average thickness of ∼500 m. We also observe a horizontal length scale of 400-500 m in the mini basin shown in Figure 5, again suggesting faulting or deformation confined to a similar depth range.
The examples in Figures 7f-7j and 15 show that part of the strike-slip offset is taken up on secondary faults indicated by a secondary maximum in the curl profiles. The location of the secondary faults is such as to try to straighten the fault, eliminating the deviations that have led to the distributed deformation in the first place. These examples also show that, behind the enhanced shortening occurring on the main rupture, there is a region of almost uniform extension. In Figure 15f, the extensional strain is ∼0.003 (3 m over 1 km). While this strain is below the noise level (±0.02) in the divergence ( ) plots, it has consequences elsewhere. Assuming the extension occurs incompressibly in a layer 500 m thick, one would expect subsidence of ∼1.5 m behind the uplift at the fault. This seems to be visible in the vertical displacement plot Figure 15h, although more than 500 m from the rupture the vertical displacements become very variable.
By looking at strike-slip displacements only, off-fault deformation in the 2013 Balochistan earthquake has been inferred to be mainly controlled by the type of surface material Zinke et al., 2014). However, as we have shown here, slip readjustment due to geometric variations can also cause the near-field strike-slip motion to decrease and to be smaller than the far field, leading to apparent off-fault deformation, in agreement with Dolan and Haravitch (2014) who pointed out the role of structural maturity. We do not know in general whether fault geometry, or the type of surface material, is the more important factor controlling the surface slip. Comparing the examples in Figures 5 and 14 suggests that both the nature of the fault geometry and of the local material may play a role in the type of distributed deformation observed. We note that much of the work on the width of the deformation zone has focused on the fault-parallel component of motion, a natural thing to do given that this component is larger than the fault-normal component for the 2013 earthquake. However, both examples in section 4.4 showed a fault-normal component distributed over a wider zone than the strike-slip component, suggesting in some places the two components may be taken up on different faults. The use of the fault-normal component led to the identification of thrust faulting away from the main fault rupture in the example in Figure 16.
The ability to derive detailed 3-D deformation fields from high-resolution stereo images allows a rich range of phenomena at the fault to be revealed. Some caveats are needed, however. It is not uncommon for surface fault traces to be used to constrain the fault geometry at depth in inverting geodetic, for example, InSAR, measurements. In this case, perhaps for faulting with a thrust component in general, doing so is clearly problematic. In addition, it is possible that the off-fault deformation is a result of deformation of a shallow layer of fractured material and only indirectly reflects slip at depth. Nonetheless, careful examination of the 3-D displacements fields enables the far-field offsets to be extracted; these can be used to constrain deeper slip on the fault.

Conclusions
We have shown that satellite-derived high-resolution DEMs provide new opportunities for producing 3-D earthquake deformation fields. By analyzing the newly derived 3-D coseismic displacements of the 2013 M w 7.7 Balochistan earthquake, we found that localized variations in fault geometry play an important role in controlling the surface slip at the fault. Several observations suggest that the distributed deformation occurs in a relatively shallow layer. Using the 3-D displacement field and its properties, we are able to identify subtle faulting that has previously been overlooked. Our measurements of fault motion are consistent with the geological interpretations, demonstrating that 3-D displacement fields, combined with the associated curl and divergence, are useful for characterizing complex surface deformation patterns in earthquakes. We also demonstrate that the 5-m ALOS DEM is a valuable global benchmark for studies of ground deformation post-2011.