A new model for the growth of normal faults developed above pre-existing structures

Two end-member models that explain strain distribution during the evolution of fault systems are (1) the isolated model, whereby faults increase in length with displacement accumulation (the ratio of maximum vertical displacement to length, D_max/L, remains approximately constant) and fault tips propagate and interact incrementally; and (2) the constant-length model, whereby faults reach their final lengths near instantaneously, and fault tips interact by linkage with increasing displacement accumulation (e.g., Walsh et al., 2002; Rotevatn et al., 2019). These end-member models and the associated global compilation of D_max-L data (e.g., Walsh and Watterson, 1988; Peacock and Sanderson, 1991; Cowie and Scholz, 1992; Schultz and Fossen, 2002; Walsh et al., 2002; Rotevatn et al., 2019) have been widely used to evaluate fault populations in different settings, including mid-ocean ridges, magmatic rift systems, and extensional sedimentary basins (Gupta et al., 1998; Bohnenstiehl and Kleinrock, 2000; Polun et al., 2018).

Testing these models is problematic because (1) global compilations combine multiple data sets from noncomparable settings, span several orders of magnitude of fault throw, and have significant scatter, therefore the models do not have distinct distributions within D_max-L plots; (2) the style of early-stage fault growth requires data from syn-kinematic packages that are hard to observe or reconstruct (e.g., Walsh et al., 2002); and (3) the models need not be mutually exclusive (Rotevatn et al., 2019).

In this study, we investigate fault growth in syn-kinematic deposits using data from the northern Krafla fissure swarm (KFS), Iceland, a region where an established fault system was near-instantaneously magmatically resurfaced by the Storaviti lava flow at ca. 10 ka (Sæmundsson, 1991; Jóhannesson and Sæmundsson, 1998) and fractured by ∼20 subsequent rifting episodes (Björnsson et al., 1979; Buck et al., 2006). The data quantity and quality and the presence of syn-kinematic deposits provide an unparalleled opportunity to consider a single fault system with D_max-L spanning two orders of magnitude.

The KFS lies within the Northern Volcanic Zone (NVZ) (Fig. 1A), one of four rift zones across Iceland (Sigmundsson, 2006). The NVZ extends north from the Vatnajökull icecap in central Iceland to the north coast where it meets the Tjörnes Fracture Zone (Einarsson, 1991), and encompasses a series of north-northeast–striking, left-stepping en echelon fissure swarms following the plate boundary (Einarsson, 2008). The NVZ includes five main volcanic systems: Askja, Kverkfjöll, Fremrinámar, Krafla, and Theistareykir (Fig. 1B). The KFS covers a 5–10-km-wide and 100-km-long region and transects the 200 ka active Krafla central volcano and caldera (Sæmundsson, 1991).

The KFS contains fractures that form a continuum from pure extensional fissures through to well-developed normal faults (Hjartardóttir et al., 2012). The main deformation region is largely confined to a central zone of fissures, commonly flanked by graben-bounding normal faults (Gudmundsson, 1984; Opheim and Gudmundsson, 1989). The resulting KFS forms a set of graben structures, with the central graben extending from south of Hverfjall through the central caldera toward the coast (Angelier et al., 1997) (Fig. 1) and the density of fault- and fissure-like fractures decreasing with distance from the central caldera (Hjartardóttir et al., 2012).

We used a series of airborne lidar surveys acquired on 7 August 2007 and 5 September 2008 by the UK Natural Environment Research Council (NERC) Airborne Research and Survey Facility (ARSF) Dornier aircraft; these surveys have an optimal resolution of ∼0.5 m on the x- and y-axes and ∼0.2 m on the z-axis, resulting in a high-resolution (0.5 m) digital elevation model (DEM) of the KFS. Following initial post-processing of the raw lidar point cloud data to x, y, and z coordinates by the Unit for Landscape Modelling (ULM), University of Cambridge (UK), we used a convergent interpolation algorithm (Haecker, 1992) to build the DEM. The DEM provides an unprecedented view of the topographic surface (Figs. 1C and 1D). Along-fault measurement intervals of 2–6 m (average ∼4.8 m, decreasing to 2 m around fault tips and complex regions) are possible by using abrupt changes in surface expression to map hanging-wall and footwall cutoffs.

Given the subvertical nature of the faults in the KFS (Opheim and Gudmundsson, 1989), our compilation of vertical displacement profiles (the throw) allows us to compare faults across the study area, which includes fissure-like faults (here defined as primarily extensional fractures with significant segments of their vertical displacement profile at <1 m). We measured 775 faults with maximum displacement ranging from ∼0.5 to 37 m and surface lengths from ∼0.015 to 2 km. We extracted D_max/L for each measured fault and plotted them against the global population (Fig. 2A, bottom right inset), which reveals that there is a marked spread of data, particularly at lower displacements.

To investigate this wide distribution, we developed a quantitative approach for categorizing the fault profiles (detailed in Table 1) based on the percentage of fault profiles that have a vertical displacement >1 m. They range from category 1 fractures, which show no vertical displacement, to category 5 fractures, which show vertical displacements along their entire length (Fig. 2A, top inset). Category 5 fractures are subdivided into three subgroups (Table 1; Fig. 2B): fully displaced faults (5a, classical), exhibiting a single bell-shaped profile; faults with a flat-topped profile (5b, flat-top); and faults with local displacement minima (5c, linked; Fig. 2). We examined where within the range of the measured D_max/L data each category lies (Fig. 2A, main) and tested whether there is a relationship between fault geometry and a fault’s location within the distribution.

Category 1 and 2 fractures, representing fractures that are fissure-like (displacement <1 m) along 65% or more of their length, are preferentially distributed in the lower portion of the D_max/L plot (Fig. 2A). Noticeably, the distribution of fissure-like fractures has a well-defined upper limit in surface fault length, with only three fissure-like faults having lengths >400 m. Furthermore, with only four exceptions, all fractures with D_max <2 m are fissure-like (category 1 or 2), and these span lengths from ∼15 to 400 m. D_max/L for fissure-like fractures spans nearly the full global range (∼10⁻³ to 10⁻¹). Fractures that are predominantly fault-like (categories 4 and 5) plot with a much narrower range of D_max/L (∼10⁻² to 10⁻¹) and displacements >2 m for all but six faults. Category 3 fractures are both fissure- and fault-like, each for ∼50% of their surface length, and form a transition in the D_max/L distribution between category 1 and 2 fissure-like and category 4 and 5 fault-like fractures.

This change in the characteristic shape of fracture displacement profiles supports a model of fault growth in which fractures evolve from being “fissure-like” to “fault-like” as they increase their vertical displacement, without changing their surface length. We refer to this region of the D_max/L plot as the fissure to fault growth zone (FFGZ; Fig. 2A). The correlation between length and displacement profile for “fault-like” fractures suggests fractures can only increase their length once displaced along >75% of their surface length. The absence of category 1–3 fractures with longer lengths implies that further displacement would be required before they could grow in length. We propose that this represents an evolution from the fissure-dominated FFGZ to the fault-dominated population referred to in the D_max/L plot as the fault growth zone (FGZ; Fig. 2A).

Fault-like fractures with >1 m displacement along >90% of their length have a range of shapes (5a, classic; 5b, flat-top; and 5c, linked; Fig. 2B). Category 5a faults are more likely to plot toward the uppermost limit of the ratio D_max/L in the global distribution, with category 5b and 5c faults showing lower D_max values for their lengths. Category 5c faults have D-L profiles that are consistent with them having formed by linkage of shorter segments. We have divided a subset of category 5c faults into their possible pre-linkage fault segment lengths using displacement minima (example shown in Fig. 2C, inset), and their D_max/L values measured. We observe that all the constituent faults lie within the category 5 region in the overall D_max/L distribution (Fig. 2C). This could imply that faults are more likely to grow by linkage when they are fault-like fractures. Alternatively, this relationship could be attributed to faults with larger vertical displacement slipping and propagating faster, and therefore linking more rapidly than smaller faults.

By assuming that our analysis of a population of fissure- and fault-like fractures in the KFS represents a snapshot of D-L data of different stages of growth, we interpret that, in this location, faults form initially as fissures. The fissures then gradually increase their vertical displacement (FFGZ in Fig. 2) before entering a phase where both length and displacement increase (FGZ in Fig. 2). We propose that this pattern of fault growth is a result of the faults forming above pre-existing faults following resurfacing by the Storaviti lava flow (Fig. 3).

Collectively, our observations suggest four stages of fault growth. Initially there is a resurfacing stage (Fig. 3A), in which there is a geologically instantaneous episode of deposition of a homogenous package (in this case, the Storaviti lava flow) that buries the preexisting underlying faults. This is followed by the flexure and initial fracture stage (Fig. 3B), where slip on the underlying fault causes tensile bending stresses at the surface, resulting in horizontal opening of isolated fissures with minimal vertical displacements (<1 m). Subsequently, the fissures enter a constant length stage (Figs. 3C to 3D), with the faults propagating downward toward the underlying fault. The strike length of each fault that has displacement increases, as does the amount of total displacement, until the faults reach D-L profiles with displacements along >75% of the surface length (category 4 and 5 fractures). Faults then enter a mature growth stage (Figs. 3D to 3E) where they can grow by a combination of increasing both their displacement and length in isolation and by increasing their length by strike linkage and dip linkage (e.g., Rykkelid and Fossen, 2002) to increase their displacement accumulation. In this final stage, data distribution tends toward the higher boundary of the ratio D_max/L in the published global distribution. The process would be repeated following the next phase of resurfacing. We suggest that fault segments can remain unconnected from one another and from the underlying structure for most of their evolution, but can form dip linkage with the underlying fault in the final stages of the growth model. An important consideration is the rheology of basalt, where columnar jointing lends itself to accommodating flexure through fissuring. In contrast, weaker or unlithified sedimentary successions could be more likely to accommodate flexure through folding (e.g., Childs et al., 2017) in the flexure and initial fracture stage.

The wide distribution of global D_max/L data in previous compilations has been attributed to differences in tectonic setting, lithology, or maturity of the individual systems (Cartwright et al., 1995). Our results demonstrate that D_max/L data, even in a very well-constrained setting, have a wide natural variability. The combination of data resolution and syn-kinematic evolution involving both radial tip propagation and the linkage of initially isolated segments leads to this wide distribution.

Although our findings are most directly applicable to magmatic divergent plate boundaries, including mid-ocean ridges and subaerial magmatic systems such as those in Iceland or Ethiopia, they have wider implications for how D-L relationships are interpreted. D-L relationships have been widely used to estimate regional strain and to discriminate strain zones (Soliva and Schulz, 2008; Polun et al., 2018). Our results suggest that care must be taken when applying surface fault D-L relationships from topographic data to predict deeper fault-controlled mechanisms, e.g., when comparing surface fault relationships in slow- and fast-spreading mid-ocean ridges (Shaw and Lin, 1996; Bohnenstiehl and Kleinrock, 2000) where rates of resurfacing relative to tectonic strain may not be comparable. Application of our model could be critical in informing subsurface geological models that are commonly reliant on surface fault data, e.g., hydrothermal and geothermal systems in which underlying faults play a key role in controlling fluid flow (Hayman and Karson, 2007).

Our observations are applicable beyond divergent magmatic plate boundaries. Seismic hazard mapping in magmatic and sedimentary systems commonly relies upon surface fault expression to predict deeper fault patterns (Scholz and Gupta, 2000), yet we show that this may be a poor predictor of underlying faults where resurfacing occurs. On longer time scales, hybrid models have been applied to explain cumulative slip in sedimentary basin evolution (Rotevatn et al., 2019). Our model suggests such processes are present over the much shorter time scale of the seismic slip–syn-kinematic deposition cycle. Furthermore, our new model for fault growth can be applied to seismic reflection data, particularly aspects of fault growth such as early-stage evolution of fault zones that are not readily detected in conventional data sets (Hayman and Karson, 2007; Baudon and Cartwright, 2008; Welch et al., 2009).

We thank the UK Natural Environment Research Council (NERC) for supporting this research (grant NE/E007414/1 and studentship for Bramham), the NERC Airborne Research and Survey Facility (ARSF) for the lidar acquisition, Schlumberger (Gatwick, UK) for providing Petrel software, Freysteinn Sigmundsson and Asta Rut Hjartardóttir for support in Iceland, and John Walsh, Páll Einarsson, Chris Jackson, and Giovanni Camanni for their constructive reviews. The NERC Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET) is a partnership between UK universities and the British Geological Survey.