The concept and mathematics behind using a Monte Carlo approach to calculating dip slips were developed by Thompson et al. (2002). The motivation behind using this approach is that it accounts for uncertainty in key parameters required for calculating slip from structural data and topographic profiles (Fig. 1). Regression statistics are calculated for lines fit to the hanging wall, scarp, and footwall. Each line has a mean value of slope and y intercept, with 95 % confidence intervals for each component, and can thus be modeled as a normal distribution. The fault dip can be modeled as any type of distribution depending on how well constrained the values are (Fig. 1). For example, if fault dips can be measured in adjacent outcrops, it may be most appropriate to model the value as a normal distribution with a mean and 95 % confidence interval; if dip is constrained only by regional fault geometry data, the user may choose a uniform or trapezoidal distribution to reflect the additional epistemic uncertainty (Fig. 1). In the simplest case, the only other parameter required to calculate dip slip is the location of the fault projected up to the scarp.

Figure 1Slip statistics inputs. A schematic of inputs into the Monte Carlo Slip Statistics code, which include uncertainty in intercept, slope, fault dip, fault position on scarp, and the age of the geologic surface in question. For surface age, a uniform distribution is shown; however, the user could define a normal or trapezoidal distribution as well. Symbols used in the figure include slope (m), intercept (b), hanging wall (hw), footwall (fw), and position (x).

Download

This approach has been widely used on dip slip faults around the world (Thompson et al., 2002; Amos et al., 2010; Rood et al., 2011; Stahl et al., 2016; Stahl and Niemi, 2017). Thompson et al. (2002) first used their method to quantify slip rate variability in a transect across the Tian Shan of Kyrgyzstan. Amos et al. (2010) showed that it could be used to identify previously unrecognized offsets along the Kern Canyon fault in eastern central California, highlighting its role in accommodating internal deformation of the southern Sierra Nevada. Rood et al. (2011) analyzed a suite of geomorphic markers along the transition from the Sierra Nevada to the Eastern California Shear Zone–Walker Lane Belt to reveal interactions among multiple faults. Stahl et al. (2016) used manually surveyed GPS profiles along the Fox Peak fault in New Zealand to demarcate segment boundaries based on slip rate. Stahl and Niemi (2017) compared manually surveyed fault scarps to geodetic estimates of strain in the Sevier Desert, Utah, to discern a local magmatic vs. far-field tectonic extension regime along the Wasatch Front. The broad interest in calculating slip statistics from fault scarps around the world, and the increasing access to high-resolution DEMs and powerful coding environments, motivates our development of a semiautomated toolkit for analyzing bulk fault datasets.

The primary goals of creating the MCSST are to (1) lower the barrier to entry for scientists by limiting the required knowledge of coding languages or general programming techniques; (2) allow the user to check for consistency and accuracy in their interpretations of fault components (e.g., fault scarp, hanging wall, and footwall) before proceeding with the slip statistics calculations; (3) optimize the functionality of GIS environments by combining previously created plugins when possible; (4) introduce a methodology for efficiently analyzing numerous fault scarps and delineating distributions of spatial and temporal patterns for slip statistics; and (5) provide open-source software for those who do not have access to commercial products. We reference hanging wall and footwall throughout because this code was developed for use in a rift setting; however, upthrown and downthrown are equally viable terms considering the underlying math should be consistent for reverse faults as well.

MCSST was designed to leverage the power and broad code base of the open-source Quantum GIS (QGIS) plugin environment. It also utilizes the key data manipulation and visualization tools developed for comma-separated value (CSV) file formats analyzed in the Jupyter Notebook Python coding environment. Free, open-source applications, such as QGIS and Jupyter Notebook, provide increased access and functionality for scientific computing and displaying and analyzing spatial datasets. We also provide a version compatible with ArcGIS.

The data calculated in the previous step can be added to the original shapefile layer to allow for spatial analysis. As outlined in the user manual, the QGIS software allows for seamless “joining” of the CSV file with the original fault scarp profile shapefile layer. When using the ArcGIS software, the ArcGIS tool developed in this study, “attachStats2Profiles”, and its Python stand-alone version can also complete this task.

The Taupo Volcanic Zone (TVZ), in the central North Island of New Zealand, is a northeast–southwest-trending, 250 km long, 14 to 40 km wide active rift (Wilson et al., 1995; Wallace et al., 2004; Seebeck et al., 2014). Extension across the TVZ over the last 1.6 Myr is expressed primarily as regional volcanism, extensional fractures, and northeast–southwest-trending normal faults (Villamor et al., 2017).

The TVZ represents the ideal case study area to test our methodology because the Quaternary geology of the region has been extensively studied through detailed mapping of geological units and active fault locations. The regional kinematics and near-surface geometry of faults are known and contribute to one of the best paleoseismic datasets in the world, the results from which we can compare our findings (Villamor and Berryman, 2001, 2006). Many fault scarps are clearly imaged to offset a relatively planar geologic horizon which is easily correlated and has been mapped across the relatively narrow tectonic province. Additionally, many high-resolution datasets exist (e.g., the QMAP Geological Map of New Zealand Project, lidar-based DEM of 1 m horizontal resolution, and regional active fault map) (GNS Science, Leonard et al., 2010, Villamor et al., 2010) (Fig. 4).

Figure 4Taupo Volcanic Zone case study. Base images of (a, b, and c) are lidar hillshade maps of the North Island and Taupo Volcanic Zone, New Zealand. (a) Tectonic features of New Zealand where the green box represents the enlarged image in (b). The red fault traces are from the New Zealand Active Faults Database (NZAFD). Active faults have deformed the ground surface of New Zealand within the last 125 000 years. The blue arrow and rate represent geologic extension rates from geodetic and Quaternary fault data (Villamor et al., 2017). The yellow line represents the case study transect. (b) Fault scarp profiles used in this case study are shown in blue and are distributed along the case study transect shown in yellow. Mapped faults in the region are shown in red (NZAFD). For simplicity, only the Earthquake Flat Formation of the Okataina Group outcrops is shown from the Geologic Map of New Zealand Project (GMAP). This is the formation that was analyzed in this study. The light blue box represents the map extent of part (c). (c) Lidar hillshade map of the location of the profile shown in Fig. 2 and analyzed in Fig. 3.

We analyzed fault surface displacements on the Earthquake Flat Formation ignimbrite of the Okataina Group (Nairn, 2002), a volcanic formation of pyroclastic material, lapilli, and ash, which has been dated to 60–62 ka (Wilson et al., 2007). The constructional surface of the ignimbrite is well-preserved and represents a previously continuous, sub-horizontal, planar feature, which has displaced by normal faults in the TVZ.

Here, we applied the methodology for the QGIS-based tool to a transect across the central TVZ. By choosing to analyze fault displacements of a single age surface we can determine spatial patterns of the deformation rates across the region occupied by the surface for the time period represented by the surface chosen. In this case the Earthquake Flat ignimbrite covers practically the entire width of the active rift (only a couple of moderately high scarps on the western margin of the rift do not displace the ignimbrite). We obtained data from as close to the transect line as possible because scarp heights appear to vary along strike of the fault traces due to the complex nature of faulting in the TVZ.

We analyzed 33 faults in a 25 km transect from northwest to southeast in the direction of regional extension. We chose faults with significant offset (>5 m) because smaller faults contribute insignificantly to the overall extension of the region and often represent splays of the larger faults. Thus, because of the filtering small scarps and two missing moderately high-relief scarps on the west of the transect in this study, we report a minimum extension rate for the central TVZ below.

For each fault, we identified the key fault components (e.g., fault scarp, hanging wall, and footwall) and determined slip estimates using MCSST. For our fault geometry distribution, we chose a trapezoidal distribution, centered around 70–80^∘ in dip. This is consistent with data obtained from geological observations in trenches and superficial exposures within 40 m from the ground surface (Grindley 1959; Villamor and Berryman, 2001; Lamarche et al., 2006; Villamor et al., 2010, 2017).

The dip slip values obtained for each fault were in line with reasonable estimates based on extensive experience conducting field campaigns in the region (e.g., altimeter measurements along transects and detailed logging of faults in exploratory trenches) and visual inspection of the digital elevation model (Villamor and Berryman, 2001; Villamor et al., 2010; Villamor et al., 2017). Further, the values plotted in Fig. 3 are consistent with visual inspection of the distance vs. elevation profile (this study). We used recently derived age constraints on the Earthquake Flat Formation to convert these values into dip slip rates.

MCSST automatically converts dip slip rate estimates into horizontal extension rate values by using a trigonometric relationship defined by the fault angle distribution. This resulted in a minimum cumulative near-surface horizontal extension rate of 2.77 mm yr⁻¹ (±0.69 mm yr⁻¹) summed across the mean net slip rate value for all faults along the transect. The minimum cumulative near-surface horizontal extension rate when summed across the median net slip rate values for all faults along the transect is 2.73 mm yr⁻¹ with a 95 % confidence interval of 0.56–5.26 mm yr⁻¹.

This is similar to the values provided by the current best estimate of Quaternary extension rates derived from fault data of 2.4 mm yr⁻¹ (±0.4 mm yr⁻¹) for the central TVZ (Villamor and Berryman, 2001) and only 8 km southwest of this study's transect; 2.9 mm yr⁻¹ (±0.74 mm yr⁻¹) for the northern TVZ immediately offshore (Lamarche et al., 2006); and 2.7 mm yr⁻¹ (±0.3 mm yr⁻¹) for the southern TVZ at the Tongariro Volcano latitude (Gomez-Vasconcelos, 2017).

Note that near-surface fault-derived extension rates along the TVZ have been converted to higher (more realistic) total extension rates (from 4 to 15 mm yr⁻¹ from south to north) in the various studies by applying correction factors that include shallower fault dips and larger fault slips at depth and the contribution from small to moderate earthquakes (small faults that do not rupture the surface) (see Villamor and Berryman, 2001 for method). The differences in fault-derived total extension rate from south to north are mainly dependent on the average “deep” fault dip value chosen (60^∘ for the southern and central TVZ and 45^∘ for the northern TVZ), which remains one of the largest uncertainties for the TVZ faults. Geodetically derived extension rates show a clear extension rate increase from north to south (Wallace et al., 2004).

The utility of MCSST is best demonstrated by the efficiency with which this study was conducted and the agreement with the current geological understanding of the central TVZ. Once our methodology was defined and the workflow implemented, the analysis of the transect (33 faults across 25 km) was completed in 1 d at a work station. We note that installing and running the Python code packages, navigating the GIS software, and inputting parameters and gaining comfort with the interface for the Jupyter Notebook pose the greatest challenges to employing the MCSST toolkit, so we provide a detailed user manual to accompany the paper to lower the barrier of entry.

Scarp heights appear to vary along strike of the fault traces due to the complex nature of faulting in many geologic settings. This is one limitation of the MCSST. The user must define unique profiles along the fault scarp for each measurement and thus may not choose the position of maximum displacement, a location that can resolve the full 3-D slip vector, or reveal the complex nature of displacement on the fault (Mackenzie and Elliot, 2017). Additionally, the MCSST approach should be used with caution or not used when the fault is located at the base of a concave or complex slope, as it requires markers that were originally parallel and ideally demonstrated to be the same age. However, the approach allows the geologist to incorporate best judgement in selecting fault components and may select far-field, linear features if erosion or diffusion is present at the scarp.

All codes, as well as the user manual along with a copyright statement and disclaimer, can be found at this Github Repository: https://github.com/wolfefranklin/MCSST_2019 (Wolfe, 2020).

Datasets used in this study include the New Zealand Active Fault database, which can be found here: http://data.gns.cri.nz/af/, GNS Science, 2020a; the Geologic Map of New Zealand, which can be found here: https://www.gns.cri.nz/Home/Our-Science/Land-and-Marine-Geoscience/Regional-Geology/Geological-Maps/1-250-000-Geological-Map-of-New-Zealand-QMAP/Digital-Data-and-Downloads, GNS Science, 2020b; and elevation data, which can be found here: https://data.linz.govt.nz/data/category/elevation/ (Land Information New Zealand Data Service, 2020).

The supplement related to this article is available online at: https://doi.org/10.5194/esurf-8-211-2020-supplement.

FDW, TAS, and BL contributed to code development and testing. PV contributed to data analysis and case study development. All authors contributed to drafting the paper. BL and PV contributed datasets. PV and TAS contributed background material and context.

The authors declare that they have no conflict of interest.

This paper was edited by Wolfgang Schwanghart and reviewed by Christoph Grützner and Michael Hodge.