Untangling tectonic slip from the potentially misleading effects of landform geometry Open Access

Fault displacement estimates from landform offsets are used in a number of applications, e.g., estimating fault slip rates in combination with Quaternary dating (e.g., McCalpin, 2009), determining earthquake repeatability or fault slip history (e.g., Lindvall et al., 1989; Kondo et al., 2010; Klinger et al., 2011), and estimating maximum earthquake magnitudes through application of scaling laws (e.g., Scholz et al., 1993; Scholz, 1994; Wells and Coppersmith, 1994; Manighetti et al., 2007; Leonard, 2010; Stirling et al., 2013). In addition, repeat surveys have been used to study temporal variation in fault morphology, e.g., scarp degradation (Elliott et al., 2012a) and postseismic tectonic deformation (Wilkinson et al., 2010; Lienkaemper et al., 2016; Mackenzie et al., 2016).

In the past few years there has been a major increase in the availability of high-resolution topographic data [submeter sampling; e.g., structure from motion (a photogrammetric range imaging technique), terrestrial lidar (light detection and ranging), airborne lidar, and satellite stereophotogrammetry], providing access to the full three-dimensional (3D) topography. Such high-resolution data sets show great potential for the reconstruction of landform offsets (e.g., Hudnut et al., 2002; Salisbury et al., 2012; Zielke et al., 2015; Zhou et al., 2015a; Haddon et al., 2016). However, analysis of landform displacements using these data is usually still performed using profiles, reducing the data to two dimensions to measure offsets. This approach has arisen from conventional profile-based triangulation survey techniques, e.g., by topographic profiles (x-z) across a scarp to estimate fault throw, or planform projection (x-y) to estimate lateral slip. In each case a straight line is fit on either side of the fault and projected to the fault to measure an offset. This generally requires strong assumptions that any effects arising from the geometric configuration are negligible (often not stated explicitly, e.g., Zhou et al., 2015a; discussed in the following).

Efforts have been made to account for epistemic and aleatoric uncertainties in piercing line offsets (e.g., Lienkaemper, 2001; Salisbury et al., 2015), but even these fail to account for (and capitalize on) the error introduced by the confounding geometries of slip, slope, and profile azimuth.

There are standardized tools for using topography data sets to reconstruct lateral offsets (e.g., Zielke and Arrowsmith, 2012; Haddon et al., 2016) or dip-slip vertical offsets (e.g., Thompson et al., 2002; McCalpin, 2009; Burbank and Anderson, 2011), but these measurements generally lack a rigorous 3D framework, and the uncertainty introduced by unaccounted-for geometry is usually poorly constrained. There are no standard methods for characterizing the 3D displacements associated with oblique slip and/or imperfect geometry, which necessitate solving for multiple components of the slip vector simultaneously.

Scalar offset measurements are commonly used to build along-strike slip distributions (e.g., Kondo et al., 2010; Klinger et al., 2011) on the principle that each individual measurement can capture the throw or lateral slip. While these measurements are sufficient if the geometry is very simple (e.g., a vertically displaced horizontal surface, or nearly horizontal channel offset laterally), in the case of rough or steep topography, the geometry of the fault-topography intersection can give rise to an apparent slip much larger or smaller than the actual slip (Fig. 1). This is often mitigated by very careful choice of measurement location, but this filtering limits the amount of information gathered and potentially introduces a sampling bias. Furthermore, it does not make use of the information provided by the full topography field. The resulting slip distributions are almost always variable and noisy and it is often unclear whether this reflects real slip variation or artifacts of measurement noise, erosion, and slip geometry (e.g., Haeussler et al., 2004; Klinger et al., 2006; Elliott et al., 2012b; Schwartz et al., 2012; Rockwell and Klinger, 2013; Billant et al., 2016).

We aim to lay out the assumptions associated with present techniques for calculating fault slip distributions and slip vectors using profiles whether measured directly by differential global positioning system (GPS) or total station, or from a high-resolution digital elevation model (DEM). We derive an analytical expression linking the vertical separation a profile measures to fault slip under a generalized surface geometry. We use it to test the validity of the geometric assumptions and quantify the errors associated with breaking them. Based on the analytical expressions for fault offset, we suggest a Monte Carlo approach to uncertainty analysis of fault offset measurements similar to that of Thompson et al. (2002), but taking into account the 3D topography-fault intersection geometry.

Given the variation of apparent landform offset caused by fault-topography intersection geometry, we derive a new, inherently 3D method to solve for the full slip vector. Utilizing the variation of apparent offset with topographic attitude and slip geometry allows us to fully exploit high-resolution topographic data sets. We test the method on the terrestrial lidar data from the 2010 El Mayor–Cucapah (Baja California, Mexico) earthquake ruptures (Gold et al., 2013) and the satellite-photogrammetry elevation model of the scarps on the Hoshab fault (Pakistan; Zhou et al., 2015a).

Conventional methods will not always provide an accurate measure of fault slip. In general, measured separations are sensitive to the orientation of the displaced feature, so without a priori information or assumptions about the slip orientation, an individual displaced landform will often not be a good proxy for total fault slip (examples in the following).

Fault offsets are measured using geomorphic markers, i.e., features of known or estimatable prior geometry, displaced by surface faulting. These can generally be divided into two classes; quasi-linear features such as ridge crests or stream channels, and correlative surfaces such as alluvial fan surfaces or river terraces. Linear features provide the tighter constraint on fault slip because there are fewer degrees of freedom. They are only ambiguous along the line of the feature; e.g., a horizontal line perpendicular to a dip-slip fault cannot constrain the horizontal shortening or extension across the fault. Correlated planar surfaces are ambiguous along two directions within the plane and therefore give a slightly weaker constraint on fault slip; e.g., a planar slope with an aspect perpendicular to a pure strike-slip fault will not record fault offset.

In both cases, the separation vector in the direction orthogonal to the marker is the only measurement that can be uniquely made with no additional information or geometric assumptions. This is conventionally estimated with profiles in the case of marker surfaces, and with more sophisticated 3D line fitting for line markers. In order to relate this separation to the slip vector, the geometry of the fault plane relative to the separation vector must be known. In the case of nonparallel markers (lines or surfaces), the marker separation is measured at the estimated location of the fault surface trace, between the projected markers (e.g., Thompson et al., 2002).

Markers forming a simple geometry relative to the fault can provide a good measure of fault displacement; e.g., a line feature perpendicular to a fault trace will provide an accurate measurement of lateral slip without being affected by any dip-slip motion. However, for less perfect geometries the orientation of the fault relative to the topography affects the apparent offset measurement, disguising the slip vector.

The effect of marker geometry on the apparent offset for a given slip vector can be illustrated by several simple cases. Figure 1A shows an example map of two perpendicular fences intersecting an oblique right-lateral thrust fault, each at an orientation of 45° to the fault. The fences give apparent lateral offsets that depart from the true value of lateral slip (S_L) substantially as fault-perpendicular motion increases (in this case shortening). Similarly, Figure 1B shows a profile view example of a pure dip-slip thrust fault intersecting two different topographic slopes. The vertical separation (Δz) measured with a profile is the sum of two components: the fault throw (s_z) and a component contributed by the horizontal advection of the topography. Unless the second component is accounted for, two different vertical throw estimates will be made for an identical slip vector. The slope angle and the fault dip must therefore be known in order to accurately estimate the throw from a profile across a pure dip-slip scarp (e.g., Baljinnyam et al., 1993; Thompson et al., 2002).

Figure 1C shows a perspective diagram of a left-normal fault scarp cutting a slope with fissuring at the base of the scarp, demonstrating the effects of slope angle, fault dip, and oblique slip on the slip measurement. A profile across this scarp would measure a vertical separation of Δz = s_z – ω, underestimating the throw because the horizontal component of slip contributes its own apparent vertical separation by the advection of topography. Projecting the channel thalweg along its axis to the scarp free-face results in an overestimate of the lateral component of slip (S_L + ε). Projecting the end of the channel to the free face in the fault-perpendicular direction gives the correct lateral measurement, but relies on a precise identification of the end of the channel, a measure that is often difficult to identify in practice and very sensitive to erosion and scarp degradation. Particularly with dipping faults, the strike-perpendicular component of slip will distort the apparent offsets of fault-oblique landforms, which should therefore be avoided unless further constraints are available. This example shows a slope with an aspect perpendicular to the fault strike; nonperpendicular aspects lead to further complication.

In the general case, a marker can only constrain slip in the direction orthogonal to the marker; slip parallel to the marker is not recorded. Thus an individual marker separation is only a proxy for total fault slip in ideal or geometrically constrained cases.

In using profiles to measure lateral slip or throw, several important assumptions are inherently made, although not always stated and often overlooked. Slope and fault geometry is usually assumed to be simple, with fault strike perpendicular to slope aspect. Slip is usually taken to be pure dip slip or strike slip unless oblique motion is evident. The azimuth of the profile is chosen by the operator and assumed to be either perpendicular to the fault strike, or parallel to the slope aspect. In addition, the geometry of the fault (strike and dip) is usually taken as moderately well known. Not all of these assumptions are made every time, but usually at least one is. Plots of separation against profile azimuth and/or local fault strike are sometimes employed to check for geometric variation (e.g., Middleton et al., 2016). We form analytical expressions to test the effect of breaking the methodological assumptions (see discussion of Analytical Fault Throw). However, in addition to the methodological assumptions, it is also generally assumed that the surfaces measured are well preserved and linear or planar; markers are usually carefully chosen to avoid areas of erosion and deposition around the scarp, adding a degree of subjectivity and potential operator bias (e.g., Scharer et al., 2014; Salisbury et al., 2015). We present an alternative to measuring profiles that may reduce the sensitivity to scarp degradation using a more 3D approach (see the discussion of Slip Geometry: 3D Analysis).

In this paper we use a number of terms related to slip geometry that have different meanings and interpretations between communities. Our usage is based on Ragan (2009), extended to incorporate further mathematical terms and described here to avoid ambiguity. The fault slip vector is the slip in the fault plane, composed of three components: lateral slip (S_L, horizontal, along strike), heave (s_y, horizontal, strike perpendicular), and throw (s_z, vertical). The vector sum of heave and throw is the dip-slip component of slip, in the updip direction on the fault plane (Fig. 2).

A planar marker displaced by a fault can be described by its attitude: the dip (θ) and aspect (ψ) of the surface. The apparent offset of the marker in a specified direction is termed the separation and can be measured in a number of directions (Hill, 1959): strike separation and dip separation in fault-plane space (P and Q in Fig. 2B) and so-called stratigraphic separation when measured orthogonal to bedding planes. In our geomorphic focus on offset fluvial and alluvial surfaces we primarily refer to the latter, which for generality we term the orthogonal separation, the orthogonal distance, d, between the displaced planes (Figs. 2A, 2B). Similarly, a displaced line marker is described by its attitude and the orthogonal separation. Where the fault geometry is well known, linear and planar markers are commonly projected to the fault plane intersection, defining piercing points and piercing lines, respectively (e.g., Gold et al., 2013).

These two planes, r_FW and r_HW, are a simplified representation of a natural topographic slope offset by discrete fault slip (Figs. 2A, 2B).

exactly as described by Thompson et al. (2002) for the case of thrust scarps. Thus, assuming pure dip slip, a known slope angle (θ), and a known fault dip (δ), we can calculate the fault slip (e.g., Baljinnyam et al., 1993, p. 29). Similarly, in the more general case, if the slope aspect is known the fault slip can be calculated without further assumptions. These equations require us to know the entire geometry of the system to calculate the fault slip. Instead, it is often assumed that the geometry is close enough to simplest case with fault strike perpendicular to slope aspect (ψ ≈ 0). It is also commonly assumed that the fault dip is well known. In the next section we test the effects of breaking these assumptions.

Thus the slope angle will be systematically underestimated by a profile, except in the ideal case where the profile is in the direction of the slope aspect.

Figure 4 shows plots of the ratio s_{z, profile}/s_{z, true} for several different cases. Under these assumptions, with fault-perpendicular aspects (ψ = 0), a profile can perfectly recover the throw, provided that the fault dip is known and the profile gives a good estimate of slope angle. However, the measurement error increases with aspect, ψ, as it moves away from fault perpendicular, becoming larger with steeper slope angles and shallower fault dip. The dashed lines demonstrate that even a small degree of lateral slip significantly exacerbates this discrepancy. If we are willing to accept 5% errors in throw, then for steeply dipping (∼60°) dip-slip faults, a slope aspect varying by as much as ∼40°, even with reasonably steep slope angles, does not significantly affect the throw estimate. However, for shallower faults (e.g., ∼30°) this is reduced to 20°–30°, or even less if there is a small component of lateral slip (e.g., Zhou et al., 2015a; see following).

The plots of Figure 4 show several specific cases, but it is worth considering the full parameter space of Equation 11. Figure 5 shows a series of density plots for different cases of fault dip (30°, 45°, 60°) and different degrees of oblique slip (S_L/S_D = 0, 0.5, 1). For steeply dipping dip-slip faults, the effect of the slope attitude is relatively small; with δ = 60° with slope angles to θ ∼ 20°, varying the slope aspect relative to global fault strike by ∼40° does not affect the measured value at the 5% level.

However, as we might expect from the conceptual models, introducing a small component of lateral slip significantly reduces this range, an effect that is enhanced by shallower fault dips. Shallower fault dips and increased obliqueness of slip result in very narrow geometric tolerance; e.g., with S_L/S_D = 1 on a 30° dipping fault, and slope angles of >5°, a slope aspect deviating by as little as ∼10° yields a >5% systematic error. Therefore, with gentle slopes, measurements will be quite tolerant of breaking the assumption of simple geometry (ψ = 0), but even a small component of lateral slip can cause significant errors in the throw estimate. Importantly, these errors are systematic, and so cannot be characterized by a simple Gaussian uncertainty.

In cases where slopes are close to flat (θ ≈ 0) and/or fault-perpendicular (ψ ≈ 0), profiles can still provide a good constraint on fault slip, but we need a rigorous way to characterize the uncertainty in the measurements due to the fault, topography, and slip geometries. The geometric formalism laid out here (Equations 5–9) can provide this framework.

Thompson et al. (2002), Davis et al. (2005), Amos et al. (2010), and Gold et al. (2012) used Monte Carlo simulations on scarp profiles to model the errors associated with estimating fault offsets in 2D; assigning probability distribution functions (PDFs) to each of the important 2D geometric parameters, they sampled the PDFs 100,000 times to estimate the probability distribution for the fault throw. Using the geometry defined herein, we extend this method to a full 3D approach, providing a means to quantify the effect of imperfectly known 3D geometry on estimates of landform offsets.

In order to estimate the PDF for the fault throw (or total slip), each component of the fault scarp geometry must be attributed a PDF (Fig. 6). In the case of surfaces fitted with a least squares regression (e.g., slope angle and intercept), a normal distribution with standard error can generally be used. Similarly, parameters known from direct measurement (e.g., dip from fault exposure) are generally also assigned a measurement PDF. Parameters with asymmetric or otherwise constrained distributions (e.g., fault dip if unknown) are generally given uniform, triangular, or trapezoidal distributions depending on the individual case, although Scharer et al. (2014) noted that nonuniform PDFs must be carefully justified. Generalizing to 3D, we incorporate two further parameters (compared to Thompson et al., 2002): the slope aspect (ψ) and the strike-slip to dip-slip ratio (S_L/S_D).

We can relate vertical separations (Δz) to throw (s_z) with Equation 6; we could also write this in terms of total fault slip, or lateral slip in the case of a strike-slip fault. Thus, provided that we can constrain the PDFs of all the geometric parameters, we can estimate the PDF of the fault throw. Figure 6 shows a hypothetical case for profile measurements across a reverse fault scarp. In this case the fault dip is unknown but a possible range is estimated, which we represent with a uniform PDF. The slip is taken to be dip slip, but the rake is allowed to vary, introducing as much as 10% lateral slip, as we might expect such small amounts of lateral slip to be missed on inspection. We take the slope aspect and the strike-slip to dip-slip ratio to be normally distributed with zero mean (though other distributions could be used). The resulting throw PDF is asymmetric as a result of the trigonometric relations between the geometric parameters. The median and 95% confidence intervals are shown as dashed lines in Figure 6. The MATLAB (https://www.mathworks.com/) script for generating the throw PDF is provided as Supplemental Item 1¹. We suggest that this full error analysis should be carried out for any slip estimate from profiles in order to accurately characterize the geometric error.

An example of the errors associated with varied geometry can be seen from the fault scarps on the Hoshab fault in Pakistan, the source of the 2013 M_W 7.7 Balochistan earthquake (e.g., Avouac et al., 2014; Jolivet et al., 2014). The 2013 earthquake generated ∼200 km of surface ruptures, with as much as 11–14 m lateral slip and 1–4 m vertical slip (Zinke et al., 2014; Gold et al., 2015; Zhou et al., 2015a). Zhou et al. (2015a) used postearthquake stereo imagery to generate a high-resolution (1 m) DEM of the entire fault zone; they used planform projection to measure lateral offsets and profiles to measure vertical offsets, both for the 2013 event and cumulative scarps formed over multiple earthquake cycles.

We reassessed the fault displacement at the site described in Zhou et al. (2015a, Fig. 2 therein), focusing on a set of cumulative scarps in an old alluvial fan (Fig. 7A). Using the elevation point cloud (subset provided as Supplemental Item 2², full data set to be made available through OpenTopography; www.opentopography.org), we took swath profiles along the profile lines P1 and P2, and fit parallel lines above and below the scarp (Fig. 7A), estimating vertical separations of 12.1 ± 1.5 m and 7.3 ± 0.3 m, respectively. Our estimate for P1 is much larger than that of Zhou et al. (2015a), highlighting both the subjectivity of profile measurements and the tendency toward confirmation bias in expecting neighboring profiles to have similar apparent offsets.

Fitting a plane to the hanging-wall slope in each profile (yellow polygons, Fig. 7), we estimated the slope attitude for the surface represented by each. At P1, the hanging-wall surface has an aspect close to fault perpendicular (ψ = –3°) with a slope angle of ∼7°, so we expect measured separations to be relatively insensitive to lateral motion. The hanging-wall surface at P2 has an aspect of ∼32° and a slope angle of ∼7°, so we must consider the contribution of lateral slip to the vertical separation. Using these slope attitudes and the profile separations, we apply the Monte Carlo method to investigate the site with a rigorous treatment of the errors, comparing them to estimates made under simple geometry assumptions.

In Figure 7 we plot estimated PDFs for the fault throw under three different cases. We represent the measurement uncertainties on the vertical profile separations as Gaussian distributions with 1σ errors. Figure 7B shows the PDF for the approach of Zhou et al. (2015a), where it is implicitly assumed that the geometry is simple, i.e., θ ≈ 0, ψ ≈ 0, so the throw is exactly the vertical separation measurement. The two estimates of throw disagree at the 95% confidence level.

Figure 7C shows the PDF for the case of assuming ψ ≈ 0, but accounting for the slope angle θ. The fault dip is assumed to be 50°–70°, covering the values estimated by Barnhart et al. (2014) and Zhou et al. (2016) for this reach of the fault; lacking any further information, we assume a uniform probability distribution across this range. The estimated throw is systematically smaller than the measured separation (in this case by ∼0.5 m), and the resultant PDFs still disagree at the 95% confidence level between the two profiles.

Figure 7D shows the PDFs for the fault throw in both profiles taking account of the full orientation of the slopes measured. We assume a value of S_L/S_D = 4.5 ± 0.5 based on the estimate of Barnhart et al. (2014) for this stretch of fault in the most recent event, which is consistent within error to the average ratio based on measurements of Gold et al. (2015) and Zhou et al. (2015a). The orientation of the topographic surface at P2 is such that lateral slip reduces the vertical separation, so a larger (10.3 m) throw would be consistent with the measured vertical separation, while P1 is largely unaffected by lateral slip (11.2 m), but the measurement error is larger. These estimates for the throw agree well within the 1σ error bound, confirming that the change in vertical separation is primarily due to the change in orientation of the surface topography.

We conclude that in locations where slope aspect is approximately perpendicular to fault strike, profiles are a reasonable tool to estimate fault throw, with a mild (<10%) overestimate if the fault dip is unaccounted for. However, where slope aspects are not fault perpendicular, profiles do not constrain the throw; further geometric information is necessary. We describe an alternative 3D method to solve for the components of slip using the two slopes described and an offset river channel (see discussion: 3D Example: Hoshab Fault).

In general, for shallow slope angles and steep fault dips, variation in the expected vertical separation is relatively small with varying slope geometry (equivalent to moving along strike over different landforms). Therefore profiles can recover fault throw well and Monte Carlo simulations can give a reasonable estimate of the associated errors. However, with steeper slope angles and shallower faults, there is significant variation in vertical separation with change in slope aspect for a given fault slip, so individual 2D-projected profiles are insufficient to constrain fault motion. Similarly, adding even a small degree of oblique slip further increases the discrepancy. While this analysis was carried out for the case of throw estimates from profiles, a similar analysis can be performed for estimating lateral slip from planform projections.

To resolve statistically significant slip variation along strike, a rigorous assessment of these errors should be completed. Monte Carlo simulations based on our geometric framework provide a much needed method to rigorously characterize the uncertainty arising from poorly constrained 3D geometry in fault offset estimates from profiles or planform projection.

Instead of taking profiles, we can exploit the availability of high-resolution DEMs by approaching the analysis of surface faulting in an inherently 3D fashion. Elliott et al. (2012b) combined pairs of nearby planform markers across the Darfield fault to calculate the magnitude and orientation of the horizontal component of the slip vector, assuming that both markers were subjected to the same slip vector. Billant et al. (2016) amassed suites of apparent topographic offsets to define both the horizontal and vertical components of a single slip vector represented by adjacent contour separations. We extend these concepts to exploit the change in apparent landform offset with varying marker geometry. We build an inherently 3D framework to constrain the full slip vector by combining multiple along-strike measurements starting from the base assumption of a constant slip vector over the region of interest. This assumption is unlikely to be valid at the entire fault scale, but may be applicable at the scale of individual and/or neighboring geomorphic features (e.g., alluvial fans, river channels, and moraines). This will not always be applicable, as there are good examples of short wavelength variations in fault offsets that are not due to geometric effects (e.g., Rockwell et al., 2002; Rockwell and Klinger, 2013). However, even where there is minor variation over the region of interest, an average slip vector is still a useful constraint, e.g., in estimating paleo-earthquake magnitudes. The absence of precisely defined piercing lines often warrants reconstructions of broader features to obtain a representative slip estimate.

Instead of profiles, we use correlative surfaces on either side of the fault because these provide a natural estimate of the topographic slope attitude. The simplest case is to approximate the surfaces on either side of the scarp as two parallel planar surfaces. Nonplanar (e.g., conical fan surfaces) and nonparallel surfaces could similarly be used, but they require more free parameters, so we illustrate the concept with the simplest case. The incorporation of line markers in this model is also discussed in the following.

Each row of Equation 13 represents a single observation, independent of the measurement location, a hanging-wall plane displaced by d from the footwall. A point from the footwall plane translated by the slip vector (s) must be in that hanging-wall plane. Introducing a second observation (d_ii from a common origin), the end of the slip vector must be along the line of intersection between the two hanging-wall planes; with a third observation, the slip vector is uniquely constrained.

Thus, the problem is equivalent to solving for the intersection between all of the hanging-wall planes superimposed, where each plane is displaced from the origin by (s ⋅ _i)_i (Fig. 8A). For Equation 13 to have a unique solution, there must be at least three independent rows in the matrix N. There must be significant geometric variation between the measurements, i.e., at least three different slope aspects and/or dips. If all of the observation planes are parallel, the problem is completely unconstrained because the planes will not intersect. Similarly, if all the observed planes share an axis (or equivalently their normals are all within a single plane), then the slip vector is only constrained to be on that axis.

We showed that changing various geometric parameters results in different apparent vertical separations, confounding slip estimates; this same variation is necessary for solving for the slip vector in 3D. In areas with very shallow slope angles, slope aspects entirely perpendicular to fault strike, or steep fault dips and no lateral slip, the problem will likely be underdetermined. In these cases, with additional assumptions, or geometric knowledge (e.g., direct measurement of fault dip), we may still be able to solve a constrained problem, as we demonstrate herein for the Hoshab fault, Pakistan.

Alongside planar surfaces, line features such as ridge crests, moraine crests, or terrace risers are often displaced by faulting. These can provide a tighter constraint on the slip vector than a planar surface; they provide information equivalent to two intersecting planar surfaces (Fig. 8B). A line feature (fully captured in 3D) can be incorporated into the formalism of Equation 13 as two rows of N, i.e., two planes intersecting along the line of the marker.

For a simple straight line feature, the choice of intersecting plane pairs is arbitrary. Two parallel lines must be within a plane; in the following case studies, we choose the first pair of planes such that the line feature in both the hanging wall and footwall is within the same plane, i.e., two planes with a separation of zero. The second pair of planes is then perpendicular to the first, with the maximum orthogonal separation. Many anthropogenic or geomorphic line markers may often be better defined by the intersection of two planes on either side than as a line directly, because this is less sensitive to erosion and/or modification (e.g., ridge crests or curbstone crests).

Because a 3D line marker provides a tighter constraint than a planar marker, only one more offset feature is required to fully constrain the slip vector, provided that it is not a plane in which the line marker lies. However, it is important that a single line feature does not constrain any component of slip along the line. Knowing the precise orientation of the fault plane can provide this third constraint as the slip vector must be in the fault plane (s ⋅ = 0, for fault plane unit normal ⁠).

The scarps on the Borrego fault from the 4 April 2010 El Mayor–Cucapah earthquake provide an excellent test case for the methods described here and verify the method as a means of assessing paleo-earthquake scarps. Gold et al. (2013) studied the Borrego fault using terrestrial laser scans of two ∼200 m sections of the fault cutting a variety of different terrain geometries, and made very precise slip measurements throughout. We used the same terrestrial lidar data available through OpenTopography (http://dx.doi.org/10.5069/G9B85622, http://dx.doi.org/10.5069/G96H4FBB) for two of their sites, applying our method directly to the bare-earth point cloud. We attempt to measure an average slip vector on this stretch of the fault, with the assumption that it remains approximately constant. This is an inherently different approach from that of Gold et al. (2013), who looked for the change in slip across the sites, but we expect to measure similar slip to their average.

Using the KeckCAVES (W.M. Keck Center for Active Visualization in the Earth Sciences, University of California, Davis) open-source virtual reality software, LidarViewer (Kreylos et al., 2013), we selected sets of points on the raw point cloud in 3D, from correlated, approximately planar geomorphic surfaces in the footwall and hanging wall (Fig. 9A). We then fit parallel planes to each correlative pair of point sets, and determine the orthogonal distance between them. We discard pairs for which a parallel model does not give a good fit, i.e., apparently correlative surfaces with different attitude in the hanging wall and footwall. The varied slope angle and aspect along strike provide sufficient independent rows of N, so we solved for the slip vector s with a simple least squares fit. The MATLAB script used to fit parallel planes and solve for the slip vector is provided as Supplemental Item 3³. We find a vertical slip (throw) of 1.56 ± 0.02 m, lateral slip of 1.9 ± 0.3 m, and fault-perpendicular horizontal component (heave) of 0.8 ± 0.2 m. The errors quoted are the statistical error on the mean; the individual measurements vary by much more, but this represents the error in mean estimated from all the observations together. This captures the aleatoric uncertainties and highlights those components well or poorly constrained, but does not account for any epistemic uncertainty, such as the suitability of the landscape for reconstruction using planar features.

Figure 9B shows a plot of the surface-orthogonal separations measured and those predicted for the best-fit slip vector applied to the surface orientation at each location. Although the data are noisy, the model provides a reasonable fit, with a root mean square (rms) misfit of 20 cm. The black line shows the expected separations based on the average slip vector for each section measured by Gold et al. (2013) (site 1: vertical 1.4 m, lateral 2.5 m, site 2: vertical 1.2 m, lateral 2.1 m), giving a larger rms misfit of 25 cm.

Examining the variance of our three slip vector components, we find the vertical component very well constrained and in agreement with the previous estimate at the 2σ level. The lateral and extensional components are more poorly constrained but differ from those estimated by Gold et al. (2013) at the 3σ level, with less lateral motion but a larger heave.

In particular, the main difference in fit between the prior slip vector and our estimate occurs at site 1, where there was significant scarp collapse and secondary faulting in the hanging wall. The average slip vector of Gold et al. (2013) suggests a consistently smaller separation than we measure, likely arising from the wavelength of the problem. Projection of channel thalwegs to the fault plane measures the offsets at tens of centimeters away from the fault, whereas our estimate measures the offset of surfaces 1–10 m from the fault plane, capturing any distributed slip in the region closest to the fault.

A further reason for the discrepancy may be the absence of a measurement of the extensional component (heave) by Gold et al. (2013), who projected linear features on either side of the fault to the exposed nearly vertical footwall fault plane. This approach assumes a slip vector in the near-surface fault plane. Along the rupture, fissures as wide as ∼30 cm were observed, so significant extension did occur and the free face at the surface does not necessarily reflect the fault plane on which slip occurred. Therefore, in projecting non–fault-orthogonal linear features (e.g., channel thalwegs) across the fissure to the scarp-face plane, some of the heave can contaminate the apparent lateral and dip-slip measurements (Fig. 1C).

For the site on the Hoshab fault described herein (Fig. 7A), we can carry out a similar analysis, demonstrating the solution to the constrained problem with a known fault dip. Previous studies estimated a fault dip of 50°–70° on this part of the fault (Barnhart et al., 2014; Zhou et al., 2016), so we assume a dip of 60° to demonstrate the method.

The geometric interpretation of the known dip case is very similar to that of the full 3D case; the fault plane forms an extra plane in which the slip vector must be. We can therefore project each observation plane to the fault plane. Each hanging-wall plane is displaced from the footwall plane by d, and intersects with the fault plane along a line (Fig. 10A). Figure 10B shows the intersection lines of each hanging-wall and footwall surface in fault geometry and in their correct spatial locations. Shifting each pair to a common origin (Fig. 10C), the slip vector is measured from the origin to the intersection point, in analogy to the 3D case in Figure 8A. Solving for the best-fit slip vector required to give the observed separations is therefore equivalent to finding the intersection between the projected lines. We solve for the intersection (Equation 13) using weighted least squares, requiring the slip vector s to be in the fault plane (⁠⁠), by setting ⋅ s = 0 with a large weight (though a dip uncertainty could be used to weight the constraint relative to the observations).

At the site on the Hoshab fault (see Balochistan ruptures case study; point cloud provided as Supplemental Item 2 [^{footnote 2}]), we consider three markers; two pairs of parallel planes measuring the separation of the slopes at P1 and P2 and the displaced channel, which forms a line feature (Fig. 7A). At profiles P1 and P2, we select sets of points from the raw point cloud in the hanging wall and footwall of the correlated alluvial fan surface (yellow and orange polygons, respectively, Fig. 7A), fitting parallel planes through the points. We fit parallel lines to points picked in the river channel from which the lateral offset was estimated by Zhou et al. (2015a). The river channel has incised since the abandonment of the P1-P2 fan surface, so the vertical offset in the channel only represents slip in the most recent event. However, the position provides a measure of horizontal slip as the river channel location was locked in after the abandonment of the fan surface. We therefore represent the channel thalweg as a pair of vertical planes in the hanging wall and footwall (assuming the river has not migrated laterally since the abandonment of the fan surface). The river therefore constrains a single component of offset (one row of N). If each surface or line for each observation pair is fitted independently, the slope normals change by <2°, justifying the use of parallel planar markers in all three cases.

We project each marker to the fault plane (Fig. 11), and solve for the slip vector using a weighted least squares fit (weighted by plane separation uncertainties), estimating a lateral slip of 39 ± 4 m and 14.2 ± 2.0 m of updip slip (12.3 ± 1.7 m throw). This throw matches our earlier best estimates from the two profiles at the 1σ level. Incorporating the assumed fault dip, this gives a ratio of lateral to vertical slip of 3.1 ± 0.6 (or equivalently 2.7 ± 0.5 strike slip to dip slip), significantly lower than that estimated by Zhou et al. (2015a) and that of the most recent event (Barnhart et al., 2014; Gold et al., 2015; Zhou et al., 2015a).

The extra constraint of the fault geometry can resolve the slip vector with fewer observations necessary compared to the unconstrained case. However, given that we only have three observations, the uncertainty is not well constrained. We also note that the slip estimate is strongly dependent on the fault dip; we have chosen a fixed dip of 60° to demonstrate the dip-constrained method, but allowing a range of dips would further quantify the uncertainty.

Salisbury et al. (2015) discussed the influence of operator decisions on the measurement of offsets on strike-slip faults, and concluded that discrepancies generally arise from misinterpretation of features and their pre-event geometries. They also found that standardizing remote measurement methods is crucial in obtaining useful and repeatable offset measurements. In this paper we have provided solutions to standardize two elements of fault offset measurement: the rigorous characterization of errors in 2D profiles, and a method to characterize fault slip in the case of nonoptimal slope geometry and/or oblique slip when the assumptions are violated for conventional methods. Both approaches still require the user to identify correlated surfaces on either side of the fault scarp and address the associated epistemic uncertainty in feature correlation, but we hope that they provide a framework under which measurements can be standardized.

In cases where profiles do not perform well, slip measurements have generally been very limited. In some cases, with fresh ruptures, distinct point to point vector displacements can be measured (e.g., Treiman, 2002), or where the entire geometry is well constrained, features may be projected appropriately to estimate offsets (e.g., Gold et al., 2013). In an increasing number of cases, the full topography field is available both for pre-event and post-event, allowing the measurement of the full surface displacement field (Borsa and Minster, 2012; Nissen et al., 2012, 2014; Glennie et al., 2014; Zhou et al., 2015b).

We suggest an alternative to profiles or planform projection for measuring discrete fault offsets at the surface. Combining multiple measurements along strike (without reference to their individual locations), we can solve for the best-fit slip vector to give these separations, based on an assumption of constant slip over the region of interest. Such a method can be used over regions in which it is reasonable to assume an approximately constant slip vector, such as individual gullies, ridge crests, or anthropogenic buildings, roads, and curbs. Using multiple planes to represent offset linear markers permits constraint by more sample points and will often be a better representation of the marker; e.g., discrepancies among measurements of an offset curbline in the 2014 South Napa (San Francisco Bay area, California, USA) earthquake (Morelan, 2015) could be mitigated by capturing the line as the intersection between the planes on either side.

The El Mayor–Cucapah earthquake scarps on the Borrego fault provide a case study demonstrating the utility of combining multiple offset planes to estimate an average slip vector in contrast to using individual line markers. However, even in this case of relatively varied and steep topography, we note large errors on two components of the slip, suggesting that the unconstrained problem may often be poorly resolved. It may therefore often be more appropriate to use a constrained form to characterize fault slip with some prior geometric constraints, as demonstrated by the example of the Balochistan ruptures. The formalism described here is a mathematically rigorous framework under which the constrained problem is a special case.

A detailed understanding of the 3D surface geometry highlights a number of common pitfalls associated with landform offset measurement. We have shown that in the case of relatively simple geometries, the common assumptions hold well. For dip-slip faults of known dip and surfaces that have fault-perpendicular aspect, 2D topographic profiles provide an easy and accurate method by which to estimate fault throw. Similarly, for pure strike-slip faults, planform projection of linear markers can provide an excellent measure of fault slip. Numerous previous studies have used these methods to good effect, using careful selection parameters to choose sites with good geomorphic preservation and simple geometries (Scharer et al., 2014). Where the assumptions hold, profiles and planform projection can give a good estimate of fault displacement, justifying their continued use provided that the assumptions and limitations are recognized.

In general, however, profile surveys do not characterize the errors associated with these assumptions. We provide a best-practice recommendation for rigorously characterizing fault slip and associated aleatoric uncertainties using a Monte Carlo simulation method. Using analytical geometric relations, a PDF can be generated for the component of fault slip based on the uncertainties in all the geometric parameters. Previous studies (e.g., Baljinnyam et al., 1993; Thompson et al., 2002) performed this analysis in an inherently 2D framework, but our sensitivity tests show that incorporating uncertainty in the full 3D geometry is also important. The resultant fault displacement PDFs will generally be asymmetric, so uncertainties should generally be reported either as a PDF or a range covering the possible values (or 95% confidence bounds). Use of this Monte Carlo method will naturally highlight cases where a profile or planform projection based approach is inappropriate, presenting as large uncertainty (broad PDF), e.g., in cases of nonzero slope aspects (|ψ| >> 0) and significant oblique slip (|S_L/S_D| > 0).

Many studies compile offset measurements along strike to estimate a fault slip distribution (e.g., Haeussler et al., 2004; Klinger et al., 2011; Haddon et al., 2016). This treats measured marker offsets as independent proxies to the local fault slip at the surface. However, any marker only constrains the component of slip in directions perpendicular to the marker (e.g., Figs. 1A–1C). The method we present combines multiple markers to constrain more components of the slip vector but still requires operator judgement in feature selection.

Scharer et al. (2014) and Salisbury et al. (2015) discussed the effects of operator biases and experience in the quality of the measurements made; they found that the quality of the geomorphic reconstruction is often the largest source of uncertainty and introduced a quality measure whereby each observation is assigned a quality rating (1–3, 1 being highest quality). Salisbury et al. (2015) noted that in general, the more experienced researchers can make more consistent estimates from a given feature, but the preservation and shape of the geomorphic marker often limit the reconstruction quality. Thus, both the quality of the reconstruction and the measurement error must be taken into account when comparing landform offset measurements. The methods we describe reduce or characterize the aleatoric uncertainties associated with the measurement geometry but do not account for epistemic uncertainty arising from poor-quality geomorphic markers. In practice, the operator must still identify suitably correlated surfaces on either side of the fault, with the uncertainty that introduces.

Scharer et al. (2014) also noted that a significant trade-off exists in building a slip distribution from multiple observations along strike, between making as many observations as possible and the quality of those observations. Our methods allow the inclusion of sites that would traditionally be excluded due to poor marker geometry, enabling more measurements along strike. In this framework, markers that are not fault perpendicular or slopes that are cut obliquely are incorporated in a mathematically rigorous manner to provide additional constraints on fault slip.

With the availability of high-resolution topographic data with point uncertainties of <5–10 cm (e.g., Gold et al., 2013), the formal measurement uncertainties are typically very low (often <10 cm), excluding any errors due to the geometric assumptions. However, in slip distributions compiled along strike, we typically see much greater variation over short wavelengths (hundreds of meters), and it is unclear whether this reflects slip variation or reconstruction quality, measurement errors, or geometric errors. Our sensitivity tests show that significant noise can arise from the breaking of the geometric assumptions and therefore be removed by correctly accounting for the geometry. However, noise arising from the marker reconstruction must still be evaluated on a site by site basis (Scharer et al., 2014). Using our methods with careful identification of high-quality markers should allow slip variation to be resolved over a variety of length scales.

The recent explosion in availability of high-resolution topography data presents great potential for the measurement of fault slip at the surface. In general, for shallow slope angles and steep fault dips, conventional analysis tools can capture fault offsets accurately. However, with steeper terrain and shallower fault dips, we get significant change in apparent offset as the slope attitude changes; this is worsened by even a small (e.g., 10%) component of oblique slip. The rigorous 3D framework described here can be used to characterize fault slip uncertainties through Monte Carlo methods, incorporating geometric uncertainty into the resultant slip estimate. A case study for this method on the obliquely slipping Hoshab fault demonstrates that profile measurements can recover the fault throw well when the slope aspect is fault perpendicular, but give biased estimates when this assumption is violated.

Where the assumptions of conventional profile based methods are violated, we propose to exploit new high-resolution topographic data sets, combining multiple measurements to solve for a 3D slip vector. With 3D marker surfaces (here planar) we use the variation in apparent offset with terrain and measurement geometry to determine the slip vector under a simple linear formulation. We have demonstrated the method in estimating a slip vector from the El Mayor–Cucapah rupture scarps on the Borrego fault, where we estimate a throw similar to previous estimates, but less lateral slip and greater (though poorly constrained) heave. A further example from the Hoshab fault demonstrates the constrained case where the fault dip is relatively well known. The geomorphic marker surfaces can be projected to the fault plane and an analogous 2D problem solved.

In some regions of significant topography or oblique slip, 2D-projected measurements will generally give biased estimates of fault slip. In contrast, our inherently three-dimensional approach can resolve the true slip vector and is enhanced by oblique slip and topographic variation.

The data acquisition and data processing and classification of the Borrego fault data set at University of California Davis KeckCAVES (W.M. Keck Center for Active Visualization in the Earth Sciences) was supported by the Southern California Earthquake Center, which is supported by the United States Geological Survey (USGS) and the National Science Foundation (NSF) (http://dx.doi.org/10.5069/G96H4FBB, http://dx.doi.org/10.5069/G9B85622). This work and the data acquisition and processing of the Hoshab fault photogrammetry derived data set at the University of Oxford were supported by the UK Natural Environment Research Council (NERC) through the Looking Inside the Continents from Space (LiCS) project (NE/K011006/1), and the Centre for the Observation and Modeling of Earthquakes, Volcanoes and Tectonics (COMET, COME30001, http://comet.nerc.ac.uk). We thank Yu Zhou for providing us with the Hoshab fault data set (to be made available via OpenTopography presently [www.opentopography.org]; subset available as Supplemental Item 2; see ^{footnote 2}) and his aid in examining the Hoshab fault case study. This material is based on data services provided by the OpenTopography Facility with support from NSF Awards 1226353 and 1225810. Discussions with many colleagues have helped in developing the methods presented here; we especially thank Barry Parsons for many helpful comments and suggestions. We thank Ramon Arrowsmith and two anonymous reviewers for giving helpful feedback that improved this manuscript.