Observations of fault geometry and cumulative slip distribution serve as critical constraints on fault behavior over temporal scales ranging from a single earthquake to a fault’s complete history. The increasing availability of high-resolution topography (at least one observation per square meter) from air- and spaceborne platforms facilitates measuring geometric properties along faults over a range of spatial scales. However, manually mapping faults and measuring slip or scarp height is time-intensive, limiting the use of rich topography datasets. To substantially decrease the time required to analyze fault systems, we developed a novel approach for systematically mapping dip-slip faults and measuring scarp height. Our MATLAB algorithm detects fault scarps from topography by identifying regions of steep relief given length and slope parameters calibrated from a manually drawn fault map. We applied our algorithm to well-preserved normal faults in the Volcanic Tablelands of eastern California using four datasets: (1) structure-from-motion topography from a small uncrewed aerial system (sUAS; 20 cm resolution), (2) airborne laser scanning (25 cm), (3) Pléiades stereosatellite imagery (50 cm), and SRTM (30 m) topography. The algorithm and manually mapped fault trace architectures are consistent for primary faults, although can differ for secondary faults. On average, the scarp height profiles are asymmetric, suggesting fault lateral propagation and along-strike variations in the fault’s mechanical properties. We applied our algorithm to Arizona and Utah with a specific focus on the normal Hurricane fault where the algorithm mapped faults and other prominent topographic features well. This analysis demonstrates that the algorithm can be applied in a variety of geomorphic and tectonic settings.

Scientists use a tectonic fault’s three-dimensional geometry and the surrounding damage zone to probe research questions about fault mechanics and hazard [16]. Testing relationships about fault mechanics including fault growth and linking this behavior to seismic hazard [711] requires observations of faults over many length scales and with varying properties including slip mode and structural maturity [12, 13]. Until recently, making even a few observations required to probe these relationships was very resource intensive, often requiring a large amount of time collecting field data and/or making measurements at a computer. Over the past two decades, a proliferation of global and high-resolution topography (HRT; more than one observation per square meter) has greatly increased the number of observations along fault zones [1417]. However, fault mapping and measuring 3D geometry remains very time intensive.

Given the recent increase in topographic data availability and quality, there have been recent efforts to automatize portions of dip-slip fault mapping and scarp height calculations on the assumption that scarp height indicates the cumulative vertical fault slip [10]. Automatic approaches for mapping strike- and dip-slip faults are often applied to optical imagery and include edge-detection methods (e.g., [18]), “ant-tracking” methods [19], and an increasing variety of deep-learning approaches (e.g., [2022]). Methods for estimating fault height or surface slip often use topography data. For example, Zielke et al. [23] and Stewart et al. [24] each developed algorithms that estimate lateral and vertical fault slips from offset markers. Hodge et al.’s [25] Scarp Parameter Algorithm (SPARTA) maps normal fault scarps from calibrated slope and curvature parameters. Sare et al.’s [26] algorithm uses a curvature template and a hillslope diffusion model to map dip-slip fault scarps and calculate height and morphologic age. Howe et al.’s [27] approach constrains slip by mapping paleolake shorelines from topographic inflection points. Despite these advancements, many algorithms either map faults or measure slip and are often applied to a single geomorphic setting, which limited their broad usage.

Here, we present a new computational algorithm that maps normal faults and estimates scarp height from high-resolution and global topography datasets. Our algorithm significantly reduces the time required to build large fault geometry and scarp height datasets. Our approach identifies fault scarps based on length and slope parameters calibrated from a manually produced fault map. We applied the algorithm to the Volcanic Tablelands of eastern California, where faults are well-exposed and their length varies over several orders of magnitude. We validated the results at two test zones and compared results to those generated by a second algorithm. We derived a height versus length scaling relationship for 152 faults which is consistent with the relationship derived from ~850 faults from the literature in different settings and over a variety of spatial scales. This demonstrates the robustness of our measurements and hence of our algorithm. We applied our algorithm to the normal Hurricane fault of western Arizona and Utah and to the full states of Arizona and Utah to demonstrate its broader utility. We anticipate that our algorithm could successfully map normal fault systems in a variety of tectonic contexts.

2.1. Geologic Setting of the Volcanic Tablelands

The Volcanic Tablelands near Bishop, California, are located in the western portion of the Basin and Range province of western North America (Figure 1). The province is in the Walker Lane Belt of eastern California that includes a mix of Miocene to present normal and strike-slip faulting (e.g., [28]). Owens Valley forms a graben that is bounded to the west by the Sierra Nevada range front normal fault system and to the east by the White Mountain normal fault, each with 3-4 km of vertical throw [29, 30]. The Long Valley volcanic system, located 70 km to the NNW of the Volcanic Tablelands, erupted and produced the Bishop tuff in 758.9±1.8ka [31, 32]. In the Tablelands, the Bishop tuff is 150 m thick and consists of moderately welded ignimbrites with polygonal cooling joints located stratigraphically above poorly welded ignimbrite [33, 34]. The tuff buried the preexisting topography with densely spaced normal faults. Likely, many of these faults were reactivated posteruption. These faults and their scarps are our focus.

The modern Tableland topography includes normal fault scarps, ancient fumaroles, and fluvial channels. Faults typically strike 335°W to 025°E, dip steeply in both east and west directions, have scarp heights that can exceed 100 m and lengths of meters to ten kilometers. The ancient fumaroles form meters to tens of meters long pillow-shaped structures. As shown in Figure 1(b), cooling joints in the tuff are exposed along the scarps’ steep faces [34]. The eroding scarps are buried by debris flow deposits and wash-slope sediments [35], as shown in Figures 1(b) and 1(c). As the faults grow, they likely reactivate the cooling joints [33].

Owens Valley is still seismically active: the 1892 oblique-right lateral M7.5-7.9 Owens Valley (Figure 1(a)) earthquake ruptured 120 km of the strike-slip and normal Owens Valley fault with offsets still preserved a few tens of kilometers south of the Volcanic Tablelands [36]. The 1986 M6.2 Chalfant earthquake (Figure 1(d)) ruptured the White Mountain Fault with a primarily right-lateral focal mechanism [37] and produced dilatational fractures with 1-2 cm of opening distributed over 10 km in the Volcanic Tablelands [38].

2.2. Topographic Fault Scarps and Fault Mechanics

Fault growth research has focused on constraining the scaling between cumulative maximum displacement (dmax) and fault length (L) with the power-law relationship,

While there has been an extensive literature discussion on the value of “n” and therefore how fault growth depends on the fault’s current stage (e.g., [57]), researchers generally agree on a linear scaling [2, 39]. We solve for Equation (1) with our new measurements of scarp height and length in the Volcanic Tablelands. We use the terminology maximum displacement and fault slip interchangeably to describe general fault behavior, often when referring to literature studies. We use the term scarp height to refer to our measurements; surface slip can be calculated from scarp height with knowledge of fault dip.

We applied our algorithm to four topographic datasets in the Volcanic Tablelands with spatial resolution varying over three orders of magnitude. “High-resolution topography” includes small uncrewed aerial system (sUAS) photogrammetric data, airborne lidar, and Pléiades satellite topography (Figure 1(d)). We use the 30 m global Shuttle Radar Topography Mission (SRTM) topography dataset in the Volcanic Tablelands (Figure 1(a)) and for the Arizona and Utah mapping.

3.1. Small Uncrewed Aerial System (sUAS)

We collected 2334 images over a single fault system in the southern test zone (Figure 1(d)) on April 12, 2019, with a sUAS that flew at ~110 m above ground level. We produced a topographic point cloud using structure-from-motion techniques [15, 16, 40] and Agisoft Metashape software. We used external georeferencing with differential Global Navigation Satellite System (GNSS) and ground control points. The average 3D root-mean-square error relative to the GNSS positions is 2 cm for the six ground control points. The point cloud covers 1.6 km2 and has a point sampling of 278 points/m2. The algorithm input is the point cloud resampled to 28 pts/m2, which reduces the computational burden for the mapping. Scott et al. [41] published the dataset which is available from OpenTopography (https://www.opentopography.org/).

3.2. Airborne Light Detection and Ranging (Lidar)

The airborne lidar data were acquired in 2014 [42] along predominantly ESE-oriented flight swaths with an ~375 m spacing. This dataset spans the southern half of the Tablelands and has 15.5 pts/m2. Like many lidar datasets, absolute errors are unreported, although surveys typically have 5-15 cm vertical errors and five times larger horizontal errors [4345]. The error is commonly observed as an apparent topographic offset along the flight line overlap. The algorithm input is the full resolution point cloud. The dataset is available from OpenTopography.

3.3. Stereosatellite Topography

The stereosatellite imagery was acquired on August 17, 2017, by the Pléiades satellites. To produce digital surface models (DSMs), we follow the same method as Mattéo et al. [22]. Each dataset has a one panchromatic image at a 50 cm resolution and a multispectral image with four bands at a 2 m resolution. Using the panchromatic tristereo images, we produced a 50 cm DSM with a 1 m vertical resolution (MicMac software, free open-source solution for photogrammetry; [46, 47]). The algorithm input is the DSM with points stored as individual x, y, and z coordinates, like in a point cloud. Figure 2 shows an example topographic profile from the Pléiades data.

3.4. Shuttle Radar Topography Mission (SRTM)

We use the 30 m global SRTM data [48] accessed from OpenTopography. The topographic dataset has a vertical precision of ~10 m. Like the Pléiades data, algorithm input is the DEM with points stored as individual x, y, and z coordinates.

4.1. Manual Mapping

A modest-sized manually produced fault map determines the algorithm’s mapping scale based on a calibration step that optimizes the slope and length mapping parameters based on this manual map. Manual maps must be completed at the scale(s) of interest for the fault mapping (i.e., coarse or fine) and must cover an area large enough to include a variety of features at the relevant scale. A single area can have more than one manual map if the mapper would like to explore the impact of scale on the results, as discussed in Section 6.1. Our maps for the Volcanic Tablelands are several kilometers along-strike and ~1 km long perpendicular to the strike.

We provided three input fault maps: (1) a coarse- and a (2) fine-scale map at the Fish Slough fault test zone based on the Pléiades topography and (3) a fine-scale map at the southern test zone based on the lidar topography, as shown in Figure 1(d). The coarse-scale mapping includes only primary faults, while the fine-scale mapping includes primary and secondary faults and ancient fumaroles. The fine-scale mapping ensures sensitivity to faults with a similar geometric dimension to the fumaroles.

The scarps can be mapped consistently at their top, steepest (approximate middle), or bottom location. Generally, a fault is best mapped at its bottom. In the Volcanic Tablelands, however, the scarp bottom is often obscured with a debris and/or wash flow. An alternative approach is to map the trace at the intersection between the scarp’s free face and the debris or wash flow. We mapped the scarp at its bottom (coarse-scale map at the Fish Slough fault) and top (fine-scale map at both test zones). While the scarp top is unusual to map as it does not represent the physical fault trace, this location is more apparent in the Volcanic Tablelands and hence easy and useful to track. We resampled the manual fault maps to a 1 m scale.

4.2. Fish Slough Fault Test Area

The Fish Slough fault is the largest Tableland fault and has a scarp height over 100 m. The fault has large breached relay ramps (Figures 1(d) and 3(a)), and the majority of the hanging wall forms a low relief surface covered by welded tuff [33]. The test zone spans 5 km along-strike and contains primary and secondary faults. The coarse- and fine-scale manual maps are shown in Figures 3(b) and 3(c), respectively. Using the Pléiades-derived topography, we mapped the faults and compared the algorithm and manual scarp heights derived from the coarse- and fine-scale calibration maps (Figures 4 and 5).

4.3. Southern Test Area

The southern test zone includes a segmented west-dipping primary fault, secondary faults, and ancient fumaroles, which are all included in the fine-scale mapping. We applied the algorithm to all four topography datasets (SRTM, Pléiades, lidar, and SfM) with the same manual mapping, as shown in Figures 6 and 7.

Our MATLAB algorithm maps scarps based on fault-perpendicular topographic gradients. In a rocky landscape such as the Volcanic Tablelands, individual boulders often cause more local topographic relief and roughness than the major faults. The scale of our manual mapping indicates the algorithm’s sensitivity to features at the fault scarp-scale and not boulders.

Conceptually, the algorithm maps scarps as regions where a high topographic slope is maintained for a specified distance or longer. These slope and length criteria based on the manually produced map ensure that the scarps are prominent features in the topography. The slope and length parameter values are determined based on a wide-range grid search. For each slope and length parameter set, a fault map is produced by the algorithm and compared to the manual calibration map through a scoring metric described in Section 5.5. The algorithm then calculates the scarp height from the offset of the hanging wall and footwall flats at the scarp location. The resulting maps depend on the manual map scale: a coarse manual map results in only major faults while a fine scale resolves major and secondary faults. The mapping parameter set that has the highest score over the calibration zone is then applied to produce a final fault map and calculate scarp height over the full area of interest. Table 1 presents the mapping algorithm steps, which are discussed as follows.

5.1. Topographic Swath Extraction

The algorithm extracts topographic transects perpendicular to the main fault strike. The swath width is 5 m and 90 m for the HRT and SRTM topography, respectively. The several-pixel swath width dampens noise, although we did no testing of different widths. Figure 2 shows an example topographic profile along the Fish Slough fault from Pléiades topography.

5.2. Topographic Slope

For each transect, the algorithm calculates the absolute value of the topographic slope over the full 5 m (HRT) and 90 m (SRTM) width every 1 m with a 2 m sliding window (HRT) and every 30 m over a 90 m sliding window (SRTM). The transect width and moving window dampen the noise. See Figure 2(a) for the slope differences between scarps and flats.

5.3. Identify Topographic Scarps and Flats

Based on a set of slope θslopeand length θlength calibration parameters which are later optimized (see Section 5.5), the algorithm identifies topographic scarps and flats from a parameter grid search. For the HRT, θslope ranges from 0.01 to 0.75 and θlength from 1 to 72 m. For SRTM topography, θslope ranges from 0.02 to 0.4 and θlength from 30 to 540 m. Conducting the grid search on a narrow set of parameters can result in local minima.

Individual 1 m segments along the transect are part of a scarp based on the following conditions: for each meter (xi) along the transect, the slope parameter φslope is the mean slope over the distance θlength,
The square brackets indicate extracting the slope over the closed interval surrounding xi. xi belongs to a scarp when
There are many gentle benches or flats located along more prominent fault scarps. The calibration map detail determines when a bench that divides a scarp is grouped with the larger scarp or remains a bench. A potential bench or flat with length of lengthflat remains a flat when

If lengthflat<θlength, the potential flat does not become a flat, but instead, the potential flat and two bordering scarps are grouped into a single scarp.

5.4. Calculate Scarp Height

The algorithm fits lines to the topographic flats bordering each fault (Figure 2(c)) using least-squares. The hanging wall and footwall slopes bordering a fault can differ. The best-fit lines are projected to the scarp’s steepest location, and their vertical offset is the scarp height. The algorithm records the top, steepest, and bottom scarp positions.

5.5. Slope and Length Mapping Parameter Calibration

Along each topographic transect, we conducted a grid search over the mapping parameters (θlengthand θslope) based on the manual maps, as shown in Figures 4 and 5. The ideal mapping parameters depend on dataset resolution and noise, mapping scale, and the overall fault size. We evaluate the mapping parameters with two criteria that penalize for manually mapped faults missed by the algorithm and for algorithm mapped faults not mapped manually. A higher score indicates better performance. For the HRT, we used a 10 m buffer for the offset of the algorithm and manually mapped faults, which allows for some variability in fault position.

The correct fault score (CorrFault) is the fraction of the algorithm’s fault points with a corresponding manually mapped fault within the distance buffer:

This score penalizes for false positives, in other words, for faults mapped by the algorithm but not manually. With many false positives, much of the algorithm fault will lack a corresponding manual fault, and CorrFault will indicate low performance.

The complete fault score (ComFault) is the fraction of manually mapped fault points with a corresponding algorithm fault point within the distance buffer:
A low ComFault indicates low algorithm performance and that many manually mapped faults lack a corresponding algorithm fault. Total score (TS) sums CorrFault and ComFault:

The best mapping parameter sets have the greatest TS.

5.6. Scarp Height Error

There are three main error sources in the scarp heights: (1) the choice of mapping parameters, (2) topographic dataset error, and (3) nonfault features identified by the algorithm. (1) The parameter calibration depends on the scale and fault map quality, which determine the landforms identified as scarps. The Fish Slough mapping results at the two scales indicate the large impact of the mapping parameters (Figure 3). However, this error is difficult to quantify. Due to (2), we interpret some heights with caution, for example, the EW “scarps” in the lidar-derived heights (particularly visible in the southern half of Figure 6(e)). (3) Non-fault feature like ancient fumaroles and river channels are approximately symmetric so their net height is often close to zero. Asymmetric features can be removed by a geologist (Section 5.8).

To calculate height measurement uncertainty, we identify the mapping parameters with

The scarp height error is the 16th and 84th percentiles (i.e., 1σ) of the height ensemble that satisfies Equation (8), which are along the same transect and are within 10 m (HRT) or 225 m (SRTM) along the transect path from the TSbest measurement of focus. These distances allow for variability in scarp location with different mapping parameters.

5.7. Fault Mapping over the Entire Area of Interest

With calibrated parameters, the mapping can be completed over the full area of interest. To calculate the scarp height error, the mapping must be done with the ensemble of parameters that satisfy Equation (8).

The algorithm mapped scarps over ~290 km2 in the Volcanic Tablelands from the Pléiades topography with θslope=20m and θlength=0.17, as shown in Figure 8. This parameter set is an approximate average of the test zone calibration parameters. With one parameter set, any biases introduced by the parameter choice are consistent throughout the mapping.

5.8. Analysis of the Faults and Their Height Derived from the Fault-Mapping Algorithm

Following the scarp mapping over ~290 km2 in the Volcanic Tablelands, we manually extracted scarp height profiles. Our algorithm does not group or label individual height measurements into faults, so this process must be done manually.

To do so, we extracted the height and length profiles for 152 individual and composite fault systems. For each of the 152 faults, we manually encircled the measurements that compose each fault scarp and no other scarp or geomorphic feature. We manually mapped a single trace which denotes the location of the fault or the main fault in the case of a composite scarp. We used the length of this manually drawn strike line to calculate the fault length because the algorithm can miss the fault tips and thus underestimate length (Figure 4). This process of labelling and simple mapping of faults is greatly expedited by the algorithm’s results. Along composite arrays, we summed the height profiles (as in Figure 7). The 152 fault maps and their profiles are shown in the supplemental material (available here). We summarized this analysis with normalized cumulative scarp heights (Figure 9) and fault length to height scaling relationships (Figure 10).

5.9. Validation: Manual Estimation of Scarp Height

To make manual scarp height measurements, we used the MATLAB code published by Scott et al. [63] for Bello et al. [64]. We extracted 2 m wide transects of the Pléiades topography separated by 50 m and perpendicular to the average fault strike. We manually selected two points along each of the footwall and hanging wall flats and one point at the scarp middle. Identical to our algorithm, scarp height is the offset of the best-fit hanging wall and footwall lines projected to the fault location. Scarp height error is calculated from a propagation of uncertainty based on fault height and location error as shown in Equation (4) of Bello et al. [64]:

mfoot and mhanging are the best-fit slopes to the footwall and hanging wall lines, respectively, and bfoot and bhanging are the y-intercepts. ΔFx is the fault position error, equal to 25% of the scarp height.

6.1. Fish Slough Fault Test Zone

For validation, we compared the algorithm maps and scarp heights to each other and to our manual results. In the X-X’ transect of Figure 3, the two faults are measured similarly by the coarse- and fine-scale mapping parameters, although their errors do not quite overlap. In contrast for the Y-Y’ transect, the coarse-scale mapping parameters resolved a single 101 m tall scarp, while the fine-scale parameters resolved three scarps with a total 70 m height. This lower height reflected topographic warping between scarps that is not counted as height.

The mapping scale impacts the resulting fault map. The coarse-scale map has few false positives (Figure 4(g)) while the fine-scale map has false positives along secondary faults (Figure 4(j)). For both calibrations, the algorithm missed ~100 m of fault length at the tips (Figures 4(h) and 4(k)). Overall, the coarse-scale map captures the bulk fault geometry and therefore may often be preferable. The fine-scale results illustrate that the algorithm can capture additional detail, although the scarp heights miss warping.

The coarse-scale scarp height profile (Figure 5(a)) shows a southern and several overlapping northern scarps. The errors are small except in areas of fault overlap. The algorithm and manual heights agree well (Figure 5(b)). For the fine-scale map, the main fault and secondary faults have a larger error with a high upper bound (Figure 5(d)), because the mapping parameters with 95% TS did not all divide the large scarp into the smaller scarps. For the northern faults (Figure 5(e)), the manual and algorithm measurements generally agree because the two faults are distinct. Along the southern fault, the algorithm captures smaller-scale fault steps that were too small to be measured manually.

6.2. Southern Test Zone

As shown in Figures 6 and 7, the algorithm identified the major NNE striking and west dipping faults in each dataset. With SRTM data, the algorithm mapped the central portion of the two major west-dipping faults but missed many secondary faults. With the Pléiades and lidar data, the algorithm mapped the primary faults, many secondary faults, and ancient fumaroles. The lidar flight line offset produced east-west-oriented strips of “scarps” separated by 300-400 m.

The optimal mapping parameters (shown in the lower right corner of each panel in Figure 6) have a similar θslope for all datasets and small variations reflect resolution and noise. θlength varies more: SRTM’s θlength of 30 m equals the data resolution. For the lidar and SfM datasets, θlength is 5 m. SfM’s lack of extra sensitivity may reflect that the manual mapping is done on the lidar data, the algorithm ignores small-scale details, and/or the SfM data was resampled by a factor of 10.

The main fault has four segments whose height decreases northward. Using all four datasets, the algorithm resolved the ~30 m maximum scarp height and the second largest scarp whose apparent height decreases with dataset resolution. The algorithm identified the two northern segments from the HRT and a single broader segment from the SRTM topography.

6.3. Comparison to the Brigham2021 Algorithm

We compared our algorithm to Brigham and Crider’s [65] algorithm, abbreviated here as Brigham2021. Brigham2021 automatically maps fault scarps and estimates scarp height by measuring the elevation difference between the fault scarp’s lower and upper extremes (the toe and crest, respectively) that bound the free face and the talus slope, as shown in Figure 11(a). The algorithm inputs are scarp profiles that are roughly perpendicular to the fault strike. To determine the local curvature minima and maxima, each topographic profile is fit with a sum of B-spline functions using the Python package NURBS-Python [66]. To map the scarp crest position, the algorithm considers the angle a between each curvature maximum and its proximate points (situated m=20 profile points away from the maximum); the crest is the maximum with the smallest a. The toe is determined in a similar fashion, as the point with the smallest angle b between the crest, the curvature minimum, and the profile starting point. The scarp height is the difference in elevation between the crest and toe.

The Brigham2021 (Figures 11(b) and 11(c)) and our algorithm both captured similar major fault traces and height profiles. However, there are some differences: The along-strike continuity of the fault is less pronounced in Brigham2021 where there are 100 m wide gaps between fault segments and some secondary faults obscure the primary fault. Both algorithms extracted similar height profiles for the two largest peaks and the maximum summed heights are within 2-5 m. The northern area differs as follows: Brigham2021 mapped a several-hundred-meter-long scarp with summed height peak of 20 m, while our algorithm measured an ~1 km broad area with several local 20 m height maximums. This difference may reflect the fact that Brigham2021 characterizes composite scarps formed by successive earthquakes for morphologic dating and may join several local peaks into a larger peak.

7.1. Hurricane Fault

The normal Hurricane fault [67, 68] crosses the western Arizona and Utah state border in the physiographic region between the Basin and Range and Colorado Plateau, as shown in Figure 12. The high-angle fault is 250 km long with west-side down synthetic scarps and east-side down antithetic scarps. We chose this as a second location for our mapping algorithm because of the prominent fault and other features such as the Grand Canyon. The Hurricane fault is larger than the Tableland faults, so this analysis demonstrates that our algorithm can be successfully applied to faults of different scales. We used SRTM topography because this dataset is ubiquitous below 60° latitude [48] and therefore covers the Hurricane fault as well as other locations where our algorithm could be applied in future studies.

For the calibration, we manually mapped the synthetic and antithetic scarps surrounding the Shivwits section (Figure 12(b)). With the optimized mapping parameters of θlength= 360 m and θslope= 0.21, the algorithm detected the synthetic trace yet missed many of the antithetic segments. Other mapping parameters would have mapped more antithetic scarps as well as synthetic fault steps and the cinder cones that were not in the manual map.

We applied the algorithm to the entire Hurricane fault (Figures 12(d) and 12(e)). The resulting trace follows the trace visible in the topographic hillshade. The algorithm detects nontectonic features such as the Grand Canyon’s walls and tributaries, as expected due to its sensitivity to topographic scarps independent of their direct relationship to faulting.

7.2. Arizona and Utah

We applied the algorithm to SRTM topography in Arizona and Utah using θlength=360m and θslope=0.21, as optimized for the Hurricane fault and shown in Figure 13. The algorithm detected many faults in the US Geological Survey’s Quaternary fault database [69], for example, the Hurricane and Wasatch faults. It mapped some Quaternary and older faults in southern Arizona as well as the WNW-oriented transition between the Basin and Range province and the Colorado Plateau in central Arizona. The Grand Canyon’s tributaries highlight the algorithm’s sensitivity to both tectonic and nontectonic features. In Utah, the algorithm distinguished between the western Basin and Range landscape and the Canyonlands and Uinta Mountains to the east.

8.1. Optimal Fault Mapping Parameters

The ideal fault mapping parameters depend on the dataset resolution, noise, and feature size. Without surprise, at the Fish Slough fault, the fine-scale mapping (Figure 3) requires a lower θlength than the coarse-scale mapping. At the southern test zone, the optimized θslope varied from 0.09 to 0.13 m/m with the highest value for the Pléiades data, likely due to the high pixel to pixel noise. The SRTM data had the largest θlength equal to 30 m, the pixel size. The lidar and SfM datasets had identical mapping parameters, likely because the pixel size is much smaller than the dimension of the faults and boulders. The Hurricane fault mapping parameters exceeded those at the southern test zone in the Volcanic Tablelands with SRTM topography due to the Hurricane fault’s large scale.

8.2. Fault Cumulative Height Profiles

As shown in Figure 9, the envelope shape of the cumulative height profile for the average Volcanic Tablelands fault is asymmetric and not elliptical or bell-shaped. When the height profile is normalized and plotted with the maximum height between 0 and 0.5, the maximum height occurs from 0.3 to 0.45. The profile is approximately straight prior to the maximum height and curved afterwards. These observations are consistent with other faults worldwide, whose cumulative slip profiles are roughly triangular and systematically asymmetric independent of fault length [10, 70].

The nonelliptical shape of the many individual and average fault height profiles in the Volcanic Tablelands indicates that these faults, like most other faults worldwide, do not behave as elastic cracks in a homogeneous elastic medium [10, 71]. While different mechanical factors impact a fault’s slip distribution [1], tectonic damage surrounding major faults is likely responsible for the asymmetric subtriangular profiles [10, 11, 72]. Greater slip occurs on the fault section embedded in more compliant rock that is typically the most mature and therefore has had the longest slip activity. The fault profile asymmetry attests to its history, in particular the long-term lateral fault propagation and lengthening that produces a maturity gradient as the fault zone becomes damaged [8, 11, 72].

8.3. Fault Slip-Length Scaling Relationships

We compared the fault length to height scaling relationships for the 152 Bishop faults that we measured to 847 faults from other studies in the Volcanic Tablelands and globally whose lengths span ten orders of magnitude, as shown in Figure 10. For the length to maximum displacement scaling law (Equation (1)), we solve for a best-fit n of 0.62 for the 152 Volcanic Tableland faults. There is considerable scatter, and the n=1 line passes through the observations, particularly for the larger faults. A recognized challenge is that there is often a low number of larger faults, and therefore, the few large faults have a disproportionate impact on the n-value. We conclude that our dataset is consistent with a large fault compilation from the literature which shows linear scaling.

In detail, the small faults (length<~300400m) have a higher Dmax/L ratio. The faults that were active prior to the 758.9±1.8ka Long Valley eruption which deposited a 150 m thick layer of volcanic tuff will now be relatively long relative to younger faults because we can only measure height accumulated post-eruption. Dawers et al. [73] made a similar observation on Tablelands faults. The variability may also indicate that faults grow by an initial lengthening followed by a period of a dominant increase in slip [74]. Still, we conclude an overall Dmax to length scaling from the Tablelands faults that is consistent with the literature compilation. Therefore, while fault growth includes phases of irregular slip accumulation, overall, when considered over multiple seismic cycles, fault growth occurs through generic, length-dependent slip accumulation. And semiautomatic processing of topography which serves as an offset marker is a powerful tool for examining the geometric and offset properties of faults.

We developed a new algorithm that maps dip-slip faults and measures scarp height from well-preserved faults using the growing archive of global and high-resolution topography datasets. The algorithm detects fault scarps from the spatial extent of the scarp’s high slope area based on relatively small calibration maps produced by an expert geologist. We focused the application of our algorithm on normal faults in eastern California’s Volcanic Tablelands. We used topography datasets with a decimeter to tens-of-meters resolution. After applying our algorithm to the entire Volcanic Tablelands, we solved for a fault slip to length scaling relationship that is broadly consistent with literature studies and further demonstrates a fairly linear scaling between maximum slip and length. The cumulative fault height profiles are typically asymmetric, consistent with observations of other faults worldwide. This asymmetry is likely related to the fault’s mechanical properties, the surrounding damage zone, and growth by lateral propagation. We successfully applied the algorithm to the Hurricane fault in western Arizona and Utah and to the entire states of Arizona and Utah. We anticipate that our algorithm will serve as a powerful tool for mapping and measuring 3D fault geometry along normal faults using a variety of high-resolution and global topography datasets.

The sUAS (https://doi.org/10.5069/G97S7KXT), airborne lidar (https://portal.opentopography.org/lidarDataset?opentopoID=OTLAS.042016.26911.1), and SRTM (https://doi.org/10.5069/G9445JDF) topography datasets are available from OpenTopography. The Pléiades satellite orthorectified image IDs used to calculate the DSM are DS_PHR1A_201709151853278_FR1_PX_W119N37_0712_01909, DS_PHR1A_201709151853179_FR1_PX_W119N37_0712_01933, and DS_PHR1A_201709151852520_FR1_PX_W119N37_0712_01913. The Pléiades images can be accessed from https://www.intelligence-airbusds.com/optical-and-radar-data/.

The authors declare that they have no conflicts of interest.

C. Scott was funded by the US National Science Foundation Postdoctoral Fellowship 1625221, the School of Earth and Space Exploration at Arizona State University and by the Université Côte d’Azur, France, where C. Scott spent one month as a visiting scientist. The Pléiades satellite image acquisition was supported by public funds received in the framework of GEOSUD, a project (ANR-10-EQPX-20) of the program “Investissements d’Avenir” managed by the French National Research Agency. A portion of the field work and Giampietro, Leclerc, Manighetti, and Mattéo’s contribution were funded by the French National Research Agency (ANR Grant FAULTS_R_GEMS # ANR-17-CE31-0008). D.A. Laó-Dávila was supported by Oklahoma State University with a sabbatical leave at Arizona State University. In addition, we thank Stephane Dominquez, Jacques Malavieille, Tyler Scott, Simone Bello, Federica Ferrarini, Chris Milliner, and Andrea Donnellan for joining us in the field at the Volcanic Tablelands to collect the sUAS imagery. We thank Dione Perkins at the Bishop Bureau of Land Management office for helping us to gain access to the Volcanic Tablelands. The scarp mapping algorithm is available from GitHub (https://github.com/cpscottasu/Fault_mapping_height). This GitHub site provides the mapping algorithm, shape files for the fault maps for the southern text zone in Bishop, and the Hurricane fault as well as links to process and download the lidar and SRTM topography datasets for these areas from OpenTopography (http://www.opentopography.org/).