Broadhaven revisited: a new look at models of fault–fold interaction

Assessment of the relationships between faults, folds and fractures in deformed rock is important for understanding deformation processes and the kinematic evolution of structure (e.g. Ramsay 1967; Ramsay & Huber 1987). Conceptual understanding of the relationships between thrust faulting and folding can be traced back to early work based on outcrop (e.g. Buxtorf 1916) and subsurface observations (e.g. Rich 1934). Early concepts were later developed into kinematic models of fold-thrust evolution (e.g. Suppe & Medwedeff 1990; Erslev 1991; Allmendinger 1998; Shaw et al. 2005). These kinematic models – fault-bend fold, fault-propagation fold, detachment fold and trishear (Fig. 1) as outlined in Shaw et al. (2005) – provide a classification structure for natural examples by which they can be understood. They are useful in ensuring that fold-thrust interpretations are geometrically valid (e.g. Suppe 1983), and have been used to relate deformation history to fold- and fault-related fractures (e.g. Bellahsen et al. 2006; Ghosh & Mitra 2009; Watkins et al. 2015).

These kinematic models of fold-thrust evolution rely on a series of assumptions: a basal detachment or thrust flat (that requires displacement to increase with depth); folding that occurs as a consequence of faulting; and concentric folding. Further, these models emphasize a rigid or undeformed footwall and do not account for the influence of mechanically heterogeneous stratigraphy in defining fold shapes. Fold shape is determined by fault geometry in fault-bend folds, by fault geometry and fault slip/propagation ratio in fault-propagation folds, and fault displacement and décollement layer thickness in detachment folds (e.g. Jamison 1987; Mitra 1990). Trishear does not obey the area or line length balance of conventional kinematic modelling – folding ahead of a propagating thrust tip is controlled by a complex set of parameters including the dimensions of the ‘trishear triangle’, fault geometry and the slip/propagation ratio of the fault. Despite these differences, the overarching theme is that folding occurs as a consequence of faulting. Similarly, therefore, any fractures associated with folding are also linked to faulting.

Recognition that displacement along thrusts is related to their map length (Elliott 1976) and that maximum displacement occurs at the centre of faults (Elliott 1976; Muraoka & Kamata 1983) led to conceptual descriptions of thrusts that propagate radially in 3D. These models of thrust-fault development by Elliott (1976), Williams & Chapman (1983), Chapman & Williams (1984) and Pfiffner (1985) differ from the later standard kinematic approaches in that fault planes propagate radially from a point source, and do not require basal detachments or fault propagation in a single direction. In these models a ductile bead deforms the volume of rock ahead of the radially propagating thrust tip (Fig. 2). It is this concept that underpins the bow and arrow rule of Elliott (1976), and that Williams & Chapman (1983) developed in their explanation of the geometries they observed in fault-related folds.

The ductile bead model of Elliott (1976), developed by Williams & Chapman (1983) for Broadhaven and similar sites, relies on internal strain in the form of stratigraphic thickening ahead of a propagating fault tip; for compatibility, there is thinning along the displacing fault. Williams & Chapman (1983) recognized the importance of the multilayered stratigraphy in developing asymmetrical folds ahead of the fault tip, verging in the direction of thrusting. Indeed, compression itself will not automatically result in folding in a homogenous medium; some heterogeneity is required for folding to initiate e.g. (Biot 1957). Once the fault has propagated through the instantaneous ductile bead, the strain is ‘frozen’ and any hanging-wall deformation is simply carried in the over-riding thrust sheet unless the hanging-wall material is translated over a bend in the thrust plane (e.g. a ramp–flat).

Williams & Chapman (1983) argued that the slope of a displacement–distance plot is related to the value of internal strain, and they developed this concept in Chapman & Williams (1984) to show how different fold thrusts may be categorized into different groups based on the slope of the best-fit line to displacement–distance data from a single thrust (Fig. 3). Uniform displacement of units along a thrust results in a displacement–distance plot with a gradient of zero, and consequently no requirement for ductile deformation in the rock volume around the fault (relative stretch = 1). Conversely, folded strata that show no evidence of faulting will yield a relative stretch of 0 and fall on the vertical axis of the distance–displacement diagram (Fig. 3). While this classification allows displacement gradients to be related to groups of fold-thrust structures, it is limited in that it assumes no change in displacement–distance relationships and does not allow for hybrid models of structural development (e.g. initial folding followed by fault nucleation and propagation).

Implicit in any kinematic model of fold-thrust development is a specific deformation sequence. McConnell et al. (1997) showed that fault displacement gradients and, thus, slip/propagation ratios of faults were also implied by these models. If fault slip/propagation ratios impact the deformation imparted to the ductile bead, it follows that kinematic models should imply fault-related ductile deformation (with the exception of fault-bend folds, which are essentially force folds formed due to translation across a bend in the fault plane). Retention of bed length and thickness, as in the constant-thickness fault-propagation fold of Suppe & Medwedeff (1990), does not, however, allow for this ductile behaviour (Fig. 1b). Kinematic models that do not allow for fault-related thickness changes are thus incompatible with the ductile bead. Where models accommodate folding through internal strain in front of the propagating fault tip (Fig. 1c, d) – for example, fixed-axis fault-propagation folding (Suppe & Medwedeff 1990) or trishear (Erslev 1991) – these are compatible with the ductile bead. The form of the ductile bead in these examples, however, is fairly arbitrary: it is generally straight edged and triangular in shape, an artefact of the geometrical and numerical modelling that underpins these kinematic fold-thrust models (Fig. 1c, d).

Direct evidence for the relationship between faulting and folding, and the nature of the ductile bead, relies on observations of structural geometries and, where possible, structural measurements. The use of quantitative geometrical structural data from field outcrops to inform models of structural evolution and strain partitioning was key to Ramsay's approach, and was used to refine understanding of how fault, fold and fracture processes link (Ramsay 1967 and others: e.g. Sorby 1908; Peacock & Sanderson 1991). These studies show that precise quantification of geometries is important to accurately constrain structural development. Supplementary measurements, such as bed thickness (e.g. Deng et al. 2013; Cawood & Bond 2018) and fracture attributes (e.g. Stearns 1964; Price 1966; Fischer et al. 1992; Ghosh & Mitra 2009), can improve understanding of fault–fold interaction and structural development, particularly where faults are not exposed. Here we employ the virtual outcrop method to digitally extract fault-displacement, bed-thickness and fracture-attribute measurements to develop the original work of Williams & Chapman (1983), and refine the existing concept of the ductile bead.

Williams & Chapman (1983) used this well-known outcrop as one of three examples with which to examine the relationship between faulting and folding (Fig. 4). Mapping of the outcrop structures, presumably through a combination of field sketches and photographs – as most are exposed on physically inaccessible cliff sections – allowed the authors to create fault-displacement profiles. Using these displacement–distance plots, they predicted slip/propagation ratios and distance to the tips of the measured thrust faults. In combination with observations of fold-thrust geometry, these displacement–distance data formed the basis for predicting deformation ahead of an advancing fault tip. These relationships formed the basis for a formal description of the ductile bead, which relates fault-displacement gradients to the amount of internal deformation (including balanced contraction and extension) in the hanging wall and footwall to the thrust (see the earlier subsection ‘The ductile bead model’). Subsequent work focused on defining groups of fold-thrust-structure-based displacement–distance patterns at Broadhaven and other outcrop sites (Chapman & Williams 1984).

In a reassessment of the Den's Door outcrop (fig. 9 of Williams & Chapman 1983), we have created a virtual outcrop model of the fold–thrust pair. The data acquisition, processing, analysis and interpretation are summarized in the following subsections. Like Williams & Chapman (1983), we created a displacement–distance profile of the main thrust that cuts the outcrop. We have augmented our assessment of the displacement–distance profile with quantified structural observations from the virtual outcrop model, including bed-thickness characterization and fracture analysis.

The classic fold-thrust outcrop at Den's Door (51°47′11″ N, 5°6′14″ W), 1 km north of Broadhaven, lies on the western edge of the Pembroke peninsula, SW Wales (Fig. 5a). This and other local structures record complex Variscan deformation of the South Wales Lower Coal Measures Formation (Hancock et al. 1982). The outcrop exposes and repeats 25–30 m of fluvio-deltaic to marginal-marine units of the late Westphalian (313–312 Ma) age (George 2008), interpreted to have been deposited in the foreland basin of the advancing Variscan front (Powell 1989). Locally, folds are generally northwards-facing, asymmetrical and bound or decapitated by northwards-propagating thrusts which commonly display curved geometries (Williams & Chapman 1983). Relay zones are common (Nicol et al. 2002) and add to structural complexity. Within fold-thrust structures, sandstone units are deformed disharmonically, with irregular spaces between these units accommodated by duplexing and thickening of shaly, fine-grained material (Hancock et al. 1982; Smallwood 1985). Locally, shortening is estimated at 40% over a 1.5 km section, a somewhat higher value than for regional estimates of 25–30% (Smallwood 1985).

The outcrop is dominated by a thrust fault and fold pair which deforms the multilayer sequence of sandstones, siltstones and shales (Fig. 5b). The main sandstone unit in the middle of the cliff displays a prominent fracture pattern and apparent thickening in both the footwall syncline and hanging-wall anticline. Disharmonic folding of interbedded siltstone and shale units is apparent in both the core and outer arc of the fold. Several thrusts cut through the succession, with minor splays branching off the main thrust. A number of minor, isolated thrust faults are present within the exposed structure, some of which serve to form relay zones between these and larger, more persistent, faults (Fig. 5b).

A total of 248 terrestrially acquired digital images were collected at the study site for subsequent processing and digital reconstruction of the outcrop. Images were acquired with a 24 megapixel Nikon D5300 DSLR fitted with a 50 mm f/1.8 fixed focal length objective. The collection of close-range (3–15 m) imagery of the outcrop was achieved either by handheld photography or through the use of a 2 m-long camera pole, to increase the proximity of the camera positions to the outcrop. An automatic shooting mode with variable shutter speeds and apertures to account for changes in lighting was combined with a fixed ISO value of 100 during acquisition. During the in-field acquisition of imagery, consideration was given to sufficient image overlap (60–80%) for digital photogrammetry and orthogonality of camera stations to outcrop surfaces, following established Structure-from-Motion protocols (James & Robson 2012).

Digital reconstruction of the study site was achieved by alignment and processing of acquired images using digital photogrammetry or Structure-from-Motion (see James & Robson 2012). This step was performed using the Structure-from-Motion (SfM) software package Agisoft Photoscan Professional 1.3.5. Image alignment and data processing (see Bemis et al. 2014 for the detailed processing steps) yielded a final photorealistic 3D mesh of 1.5 million faces (Fig. 6a). Final mesh resolution of the dataset yielded an 8 mm point spacing in 3D, with a photorealistic texture resolution of 20 mm ground pixel resolution.

Georeferencing and scaling of the virtual outcrop was achieved by the inclusion of ground control points in the survey. Control point locations were recorded by differential GPS, which allowed a 3D accuracy of 0.012–0.023 m to be achieved for this survey. Further calibration of model scaling, orientation and location was performed by comparing this dataset from high-resolution LiDAR scans of the study site. In addition, data was ground-truthed with in-field orientation and scale measurements following established workflows (e.g. Cawood et al. 2017). Calculations from applied corrections provide an estimate of >1° rotation and a scaling ratio of 0.998 of the virtual outcrop relative to control measurements.

Subsequent to processing and georeferencing of the photorealistic 3D mesh, the virtual outcrop was transferred to Midland Valley's Move software for interpretation and interrogation in a digital 3D environment. Initial interpretation of the structure was performed by polyline interpretation along key horizons and thrust segments. This interpretation of key horizons through zones of structural complexity and across thrust traces was helped by the use of high-resolution images used for photogrammetric reconstruction, recorded field observations and summary sedimentary logs recorded while in the field. In addition to key horizons and thrusts, fracture traces within the prominent sandstone beds in the outcrop were digitized to allow examination of changes in fracture attributes through the structure (Fig. 6b).

Following interpretation, 3D polyline interpretations of horizons, faults and fractures were projected onto a north–south cross-section plane (Fig. 7a), to allow analysis in 2D. A vector of azimuth 261° and plunge 10° was chosen for the projection of 3D polyline interpretations to a north–south-orientated cross-section. This vector was defined by the orientation of the π-axis, taken as normal to the best-fit π-girdle through poles to bedding and fault-surface orientations (Fig. 7b). These measurements were collected by a combination of fieldwork and digital data extraction from the virtual outcrop, and projected to a 2D section following established protocols (e.g. Martín et al. 2013; Cawood et al. 2017). Projected interpretations of horizons, fault traces and fractures (from a 3D to a 2D section) were used for subsequent measurements of bed thicknesses, fault displacements and 2D fracture orientations. These were made in 2D to avoid measurement bias on non-planar surfaces in 3D.

Digital interpretation of the study site using the virtual outcrop (Fig. 6b) provided nine interpreted units, of which units h4 and h5 are the prominent sandstone layers (Fig. 7; see Fig. 5 for a field photograph). When projected onto a north–south-orientated cross-section (Fig. 7), the interpretation outlines the hanging-wall anticline and footwall syncline geometry of the outcrop, with prominent thrust faults that truncate horizons. There is a degree of structural complexity and disharmonic strain accommodation through the deformed units, particularly in the hanging-wall anticline. In general, we record an overall structural geometry (Fig. 7) that is similar to the original interpretation (Fig. 4a) of Williams & Chapman (1983). In spite of similarities between interpretations, there are a number of readily apparent differences that should be noted:

• Within equivalent areas, this study records the presence of 18 thrust segments (Fig. 7) compared to 11 in the original study (Fig. 4a). This difference may be attributed to improved viewpoints and virtual access to zones on the structure not accessible on foot, such as those high on the cliff (see Fig. 5b for a field image) using the virtual outcrop method. Equally, however, low-offset thrust segments may have been omitted from the original study because they were seen to be insignificant or to improve the clarity of the interpretation.

• At the base of the outcrop, Williams & Chapman (1983) interpreted a low-angle thrust that truncates the lower part of the footwall syncline (Fig. 4a). When interpreting the virtual outcrop (Fig. 6b), we could not find any evidence of this lower thrust, or of the disharmonically folded units directly above it, as shown in the original interpretation (Fig. 4a). Deformation in the footwall to the dominant thrust is thus recorded in this study as being less complex than in the original interpretation: harmonic folding dominates deformation in the footwall syncline.

The technique of using displacement–distance profiles to determine slip/propagation ratios and distance to fault tips for thrust faults was established by Williams & Chapman (1983). Here we construct a displacement–distance profile from digital measurements by plotting unit displacements against distance of units, in the hanging wall of the thrust, from a reference point (see Williams & Chapman 1983 for displacement–distance profile construction). To compare our dataset with that from the original study, a scaled version of the Williams & Chapman (1983) interpretation was re-measured to produce a displacement–distance profile in which fault-displacement gradients could be directly compared (Fig. 8a). 2D measurements on projected interpretations yielded a virtual-outcrop-derived distance–displacement profile that, in agreement with the data of Williams & Chapman (1983), records an increase in fault displacement downwards (Fig. 8a). This suggests that maximum fault displacement occurs lower down in the section, at an unexposed position.

Several differences in displacement patterns are apparent, however, when these datasets are overlain: our profile records a lower overall displacement gradient than the original study, and a greater variability in displacement along the thrust (Fig. 8a). Best-fit lines to the displacement–distance data (Fig. 8b, c) represent the average fault-displacement gradient, which relates to the internal strain ahead of the propagating fault tip (see Figs 2 & 3). Both datasets record some deviation of data away from the best-fit line (R² = 0.91 for the Williams & Chapman 1983 profile and an R² = 0.64 for the virtual-outcrop-generated profile in this study). A lower R² value for the virtual-outcrop-derived data records greater variability in displacement gradient along the thrust, and thus slip/propagation ratio, than shown in the original study. This deviation of virtual-outcrop-derived data from the best-fit line (Fig. 8b) is particularly pronounced within the sand-rich units h4 and h5 (see Fig. 7 for the interpretation), where a roughly constant displacement (or flat displacement gradient) is recorded.

In order to assess the interaction between folding, thrusting and bed-thickness changes, projected polyline interpretations were used to measure bed thicknesses in each mapped unit, across the fold-thrust structure (Fig. 9). Unit thicknesses were measured at the following four structural positions: syncline long-limb; footwall syncline hinge; hanging-wall anticline hinge; and anticline back-limb. Data record a trend for increased unit thicknesses with proximity to fold hinges: units appear to be thickest where they are crossed by mapped axial traces of the syncline or anticline (Fig. 9). It should be noted that while the measurements in Figure 9 record a general trend for thickening of units towards fold hinges, the mechanism for unit thickening is unknown, although is interpreted to be as a result of folding strain.

A systematic assessment of thickness variability according to structural position was performed by measuring the thickness of each unit at regularly spaced 1 m intervals, along the base of that unit, through the cross-section. Thickness measurements were made orthogonal to the bedding surface of each unit and parallel to the north–south-orientated cross-section through the length of each deformed horizon. Thickness measurements show general agreement with widely spaced data (Fig. 9) in that they record a general trend for thickening of units through fold hinges (Fig. 10). In the footwall of the main thrust, units retain relatively constant thicknesses through the long limb of the syncline (Fig. 9) and thickening coincident with the syncline hinge zone (Fig. 10). The majority of units in the footwall record a return to average thickness (see Fig. 10) with immediate proximity to the thrust. A return to average thickness of units adjacent to the thrust fault appears to correspond with increased distance of the syncline axis to the thrust; we thus interpret unit thickening in the footwall as fold-related deformation rather than fault-related thickening along the main thrust.

Units in the hanging-wall anticline to the thrust similarly record a general trend for thickening in the hinge zones of the fold, although with greater variability (Fig. 10) than in the footwall counterpart. As with the footwall, units retain relatively constant thicknesses in the long limb of the fold followed by a general increase in thickness through the hinge zone (Fig. 10). Patterns of thickness change, however, are more complex in the hanging wall, with several units recording alternations between thickening and thinning along bed lengths through the anticline. Interpretation of the structure reveals two hinge zones (Fig. 9) within the hanging-wall anticline, which broadly correlate with increased unit thicknesses in the deformed succession (Fig. 10). As several low-offset, minor thrusts are present through the fold-thrust structure, measured thickening of units was separated into fault-related and fold-related (ductile) components (Fig. 10). This approach allowed: (a) a direct comparison of fault-related and fold-related thickening; and (b) an area balancing correction to be applied to the restored cross-section (see the ‘Cross-section balancing’ subsection later in this section).

Interpretations of undeformed thicknesses were derived by fitting straight lines to thickness measurements where there is no clear pattern of tectonic thickening (see Fig. 10 for interpretations). It should be noted that this approach is limited in that it does not account for the influence of sedimentary processes on unit thicknesses. Unit h2 records a gradual decrease in thickness from north to south (Fig. 10) that does not appear to correspond to directly observed tectonic features. This trend may record thinning in more ductile units in the long limb of the fold or variability in the original depositional thickness. Erosive horizons and channel sandstones of variable thickness have been recorded at this locality (George 2008) and elsewhere in units of equivalent age and setting in Pembrokeshire (Cawood & Bond 2018). Estimates of original depositional thicknesses and the amount of subsequent shortening by tectonic thickening (Fig. 10) may thus be biased by disregarding the influence of stratigraphic geometries and the possible lateral thinning or thickening of sedimentary bodies by depositional processes. Similarly, this approach does not account for penetrative strain undergone by units before the onset of folding and faulting, and where measured thicknesses fall below the undeformed average (e.g. unit h9: Fig. 10) it is difficult to differentiate sedimentary from tectonic thinning, or a combination of both.

Our virtual-outcrop-derived thickness measurements show several differences to the interpretation of the structure by Williams & Chapman (1983) (Fig. 4). Of particular note are differences in interpreted bed thicknesses: the graphic interpretation of Williams & Chapman (1983) (Fig. 4a) suggests that significant thickening of units occurs right up to the thrust surface in the footwall syncline, with no apparent influence of axial-trace positions. Our data suggest this is not the case (Fig. 10); a difference in interpretation due to virtual access provided by the virtual outcrop method. Apparently, greater footwall unit thicknesses (Fig. 4a) in Williams & Chapman (1983) than in this study (Fig. 7) are likely because of a foreshortening effect when the structure is viewed from the ground. Williams & Chapman (1983) also noted thinning of beds in the hanging wall of the structure sub-perpendicular to the thrust, and suggest that this stretch balances the internal strain and thickening of strata in the fold hinges, but they do not attempt to quantify either element.

The prominent sandstone units in the outcrop, units h4 and h5 (Fig. 9), show distinct fracturing which fans around the folded strata (Fig. 5b); to better understand deformation within this sandstone unit we have calculated the fracture intensity (e.g. Price 1966). To avoid the effects of relative unit competences and thicknesses on estimated fracture intensity (e.g. Harris et al. 1960), we used circular scanlines along the interface between units h4 and h5 only. Scan circles of 1 m radius were spaced at 1 m intervals along this interface. Fracture intersections with scan circles were counted (Mauldon et al. 2001) at each station (Fig. 12a) and plotted as a function of distance along the scan transect (Fig. 12b). As all scan circles are of equal radius, we do not take into account the effects of scan-circle size on the estimated intensity (Mauldon et al. 2001), but present intersected fracture count (n) as a function of distance along the chosen interface.

Circular scanlines along the h4–h5 interface demonstrate an increase in fracture intensity with proximity to the core of the fold-thrust structure (Fig. 12a). The highest numbers of intersected fractures appear to coincide with fold hinges and zones where minor faulting occurs, with the lowest numbers in the relatively undeformed parts of the section. Scan circles at the southern end of the cross-section record the lowest fracture counts. Scan-circle fracture counts record general agreement with unit thicknesses: unit thickening in the core of the structure (Fig. 10) is broadly matched by a trend for increased fracture intensity through the same part of the structure (Fig. 12). Based on this agreement between datasets, we suggest that unit thickness variability is primarily related to tectonic processes at the study site, rather than thicknesses changes within the stratigraphic template.

An investigation of fracture orientation was performed by assigning colours to fractures based on their 2D orientation in the north–south-orientated cross-section (Fig. 13). Deviation in fracture orientation from a vertical reference datum was performed in FracPaQ 2.0 (Healy et al. 2017) by modification of the workflow for fractures mapped on a horizontal surface. Fracture orientations are observed to relate to the orientation of bedding surfaces through the fold-thrust structure, with the majority of fractures near-perpendicular to bedding (Fig. 9), but with some refraction in the footwall syncline hinge (Fig. 13, inset a). Abrupt changes in fracture orientations coincide with mapped thrusts (Fig. 13, inset a; see Fig. 9 for the interpretation). Several areas of the outcrop were identified as containing fractures that provided direct evidence for slip surfaces, not directly observed, along bed interfaces (Fig. 13, insets b–d).

This two-stage approach yields a final restored cross-section length of 52–56 m, and a calculated shortening of c. 36%. Although not stated in the original study, we derive a shortening value from the restoration of the Williams & Chapman (1983) cross-section of c. 33%. The higher estimate of shortening amount provided by our data is likely to be due to the excess area correction: our line-length restored section provides a shortening value of 34%, similar to the amount calculated from the Williams & Chapman (1983) interpretations. The overall geometry of our balanced cross-section (Fig. 14a) broadly agrees with that presented by Williams & Chapman (1983) (Fig. 14b) in that two prominent, shallowly dipping (30°) north-propagating thrusts are present in both restored cross-sections, as are a number of subsidiary, low-offset thrusts, back thrusts and layer-parallel slip horizons.

The similarities in shortening values and restored geometries between this and the original study appear to support the compatibility of our digital approach. Digital interpretation and high-resolution data do not, however, allow us to circumvent some of the well-known limitations of cross-section restoration (see Hossack 1979 for a review). Adjusting restored line lengths by combining undeformed thicknesses with excess unit areas (Fig. 10) is limited in that the method assumes plane strain, area conservation and layer-cake stratigraphy. Potential stratigraphic thickness variations at the study site (e.g. unit h2: Fig. 10), evidence for ductile thickening (Fig. 10), the possibility of ductile thinning and the potential for penetrative strain during progressive deformation (e.g. Koyi et al. 2004) are just some of the elements that may challenge the line-length, area-balanced approach to cross-section restoration. Given these limitations, cross-section restoration was performed as a means of comparing our data to that of Williams & Chapman (1983), rather than as an attempt to definitively capture the exact nature of the pre-deformational template.

Displacement–distance data (Fig. 8) record relative stretch values of 0.93 (this study) and 0.83 for the re-measured Williams & Chapman (1983) interpretation (Fig. 4a), which yield slip/propagation ratios of 0.07 (this study) and 0.13 (original study). As noted by Williams & Chapman (1983), this low slip/propagation ratio suggests rapid thrust propagation through the structure, with limited deformation of the ductile bead. Based on the work of Chapman & Williams (1984), calculated relative stretch values classify this structure as a Group 2 (Fig. 3) fold thrust: hanging-wall and footwall cut-offs cannot be exactly matched across the fault (due to a sloped displacement gradient) and thus folds are inferred to have grown in advance of the propagating thrust. Thrust-perpendicular limbs of the fold pair (Fig. 9) and significant thickening of units (Fig. 10), however, are suggestive of a significant amount of ductile deformation. We record a total undeformed thickness of 7.55 m for mapped units and a maximum deformed thickness of 12.5 m for the equivalent stratigraphy in the hanging-wall anticline. Based on the ductile bead of Williams & Chapman (1983), if this thickening is to be accounted for by fault propagation alone, the gradient of the displacement–distance profile would need to coincide with an internal stretch value of c. 0.6 (Fig. 3), which is not supported by our data (Fig. 8).

Our data support the proposal by Williams & Chapman (1983) that significant folding in the hanging wall and footwall coupled with a low slip/propagation ratio for the thrust is indicative of early folding before the onset of thrusting (Fig. 15a). As noted by others (e.g. McConnell et al. 1997), low slip/propagation ratios of thrusts through folds may indicate fault segments that propagated with little associated folding or thrusts that post-date the early folding of strata. Our data show that displacement is not uniform along the thrust, suggesting that some of the observed deformation is related to thrust propagation. The fold-thrust structure at Den's Door thus probably represents a hybrid between an asymmetrical buckle fold, a fault-propagation fold (e.g. Suppe & Medwedeff 1990) and a break-thrust fold, as defined by Willis (1894). It is important to note that the use of displacement gradients to classify fold-thrust structures (Chapman & Williams 1984) is limited as straight-line displacement–distance relationships cannot exist along an entire fault. Any fault with a finite geometry must record a bow and arrow (Elliott 1976) geometry: displacement must increase from zero at the fault tips to a maximum displacement near the centre of the fault. The overall shape of the displacement–distance profile (i.e. how much variation in displacement gradient exists) will thus determine how useful information from a straight-line fit on displacement–distance data is.

Generalized displacement–distance relationships and their relationship to fold-thrust classes or groups (Fig. 3) do not take into account the effect of mechanical anisotropy of the deformed strata – they implicitly assume homogenous media. In contrast to the original study, our data show an irregular pattern in displacement along the thrust (Fig. 8b) with three distinct segments on the displacement–distance plot: a relatively steep segment in mixed lithology units h6–h9 (relative stretch = 0.86), to a zero gradient in sandstone units h4 and h5 (relative stretch = 1), followed by an increase in gradient in mudstone units h1–h3 (relative stretch = 0.75). Detailed measurements of normal-fault displacements through multilayers (e.g. Muraoka & Kamata 1983; Wilkins & Gross 2002; Roche et al. 2012) record similar stepped patterns in displacement profiles. Muraoka & Kamata (1983) classified these as ‘M’- or mesa-type profiles, where the displacement–distance plot records steep gradients in less competent material and zones of constant displacement in rigid or competent beds (see Fig. 15c). Our displacement–distance plot resembles part of an ‘M’-type profile, suggesting that fault propagation through the sandstone layers h4 and h5 was faster, with little or no associated ductile deformation, than through the finer-grained units above and below.

‘M’-type displacement profiles (Muraoka & Kamata 1983) record a ductile–competent–ductile transition and, as a consequence, little or no fault-related thickening or thinning of the competent unit (Fig. 15b). Correspondingly, using an area-balanced ductile bead model, there is a requirement for increased ductile deformation, either by thickening or thinning, in adjacent incompetent units (Fig. 15b). This suggests that displacement profiles through deformed multilayers should record steepest displacement gradients and corresponding increases in thickness variability through incompetent material (see Fig. 15b, c). Our data appear to support this hypothesis: we record a decrease in the coefficient of thickness variation (Fig. 11) through sand-rich units h4 and h5, where fault-displacement gradients are also lowest. Thickness measurements (Fig. 10) record a degree of complexity, with greater variability in the hanging wall to the thrust (Fig. 11), particularly higher in the stratigraphy. This hanging-wall thickness variation is likely to imply that subsequent to initial folding, deformation by fault propagation was restricted to the hanging wall of the structure.

The presence of extension and contraction in a ‘quadrantal distribution’ (Williams & Chapman 1983) around the fault-nucleation point (Fig. 2) assumes deformation of a homogenous rock volume. Where a mechanically heterogeneous volume makes up the ductile bead, thickness variability in less competent strata is required (Fig. 15b) for compatibility of variable displacement gradients (Muraoka & Kamata 1983). This is manifest as extension and contraction of the ductile bead at the tips of the ‘M’-type fault segments, which will modify the geometry of the ductile bead and its internal mechanical behaviour. Our data record an array of bedding-parallel slip horizons (Fig. 13), brittle deformation and outer arc extension within competent sand-rich units (Fig. 12), isolated thrust and back-thrust segments (Fig. 9), and disharmonic folding, particularly in the hanging wall to the thrust (Fig. 9). Heterogeneous distribution of strain and variable mechanisms of deformation within the ductile bead highlight the complexity of fault-related deformation in a multilayer system, particularly when a structure evolves through multiple phases of deformation.

We interpret the structure at Den's Door as having developed by initial folding followed by fault propagation and finally by thrust breakthrough (Fig. 15a). The structure has therefore developed on the displacement–distance classification plot of Chapman & Williams (1984), redrafted in Figure 3, through different fold-thrust classification groups. At t1 and t2 (Fig. 15a) the structure falls into the Group 0 classification (Fig. 3), followed by a transition to Group 2 at t3 and Group 3 at t4. Progressive deformation highlights the problem with assigning a single fold-thrust grouping to a structure based on final displacement–distance characteristics: the technique is not applicable for fold thrusts that developed through hybrid mechanisms. Categorizing a fold-thrust structure into a single group may lead to erroneous predictions of structural evolution and understanding, particularly when limited exposure only allows displacement–distance plots for parts of faults to be produced.

The original mechanical stratigraphy, and any deformation that occurs before fault nucleation and propagation, will influence the initial shape, extent and behaviour of the ductile bead (Fig. 15a). As such, a simplified box-like pre-deformational template (Fig. 2) is not applicable for multilayered stratigraphy and/or where multiphase deformation has taken place. Thickening of strata in fold hinges will modify the geometry of the ductile bead, as it propagates ahead of the thrust, making the ‘quadrantal distribution’ (Williams & Chapman 1983) of extension and contraction difficult to identify. Thus, even when a volume appears to be homogenously deformed, reliably relating deformed thicknesses to a generalized ductile bead will be difficult if significant deformation occurred before fault nucleation. The nature of the ductile bead is thus influenced internally by mechanical heterogeneity, arising from multilayered strata and pre-thrust deformation.

Improved measurement resolution afforded by remote data acquisition to build a virtual outcrop model has allowed a detailed examination of the Den's Door fold-thrust structure, resulting in a refinement of previous work undertaken on the outcrop by Williams & Chapman (1983). Our virtual-outcrop-derived measurements record small-scale variations in fault displacement, bed thickness and fracture attributes – this has allowed detailed quantification of the interaction between fault-displacement gradients ahead of a propagating thrust tip. Our results record a degree of structural variability and complexity not recorded in the original study (Williams & Chapman 1983). We attribute this to improved measurement data provided by remote acquisition, particularly as the majority of the outcrop is inaccessible on foot and views of the structure are foreshortened from below. Virtual access to geological data and improved measurement availability are key advantages of the virtual outcrop method (e.g. Buckley et al. 2008), which allowed us to quantify, in detail, small-scale structural relationships across the whole structure, and significantly develop the work from the original study.

The availability of high-resolution (e.g. Vollgger & Cruden 2016), accurate (e.g. Cawood et al. 2017) measurements is a key advantage of remotely acquired geological datasets. An accurate projection of data and the removal of viewpoint bias is an advantage of the virtual outcrop method (e.g. Tavani et al. 2016) that can significantly impact final results. From an examination of the original outcrop interpretation (Fig. 4) by Williams & Chapman (1983), which was likely to have been constructed from field observations and images, we infer that viewpoint bias and a lack of direct access to the outcrop in the original study led to erroneous interpretations of bed thicknesses that impacted the collection of high-resolution fault-displacement data. This led to the significant differences shown, in fault-displacement data, between our dataset and that from the original study (Fig. 8), and consequently to differences between structural interpretations. As noted by Pfiffner (1985), displacement–distance measurements and the derived interpretations of process are reliant on both accurate data projection and cross-sections that are precisely orientated parallel to the transport direction. Here we have circumvented viewpoint bias in the original study (Fig. 4) by both remote data acquisition and geometrical projection of data (Fig. 7) along a precisely defined vector.

Small changes to measured displacement gradients may significantly impact structural interpretations and inferences about the mechanical behaviour of deformed strata (e.g. Muraoka & Kamata 1983). Our data record thrust displacements that differ from the original study (Williams & Chapman 1983) by only a few tens of centimetres (Fig. 8). But the impact on the resultant interpretative models is significant. Recorded non-linearity in fault displacement (Fig. 8b) highlights, in our interpretation, the role of mechanical heterogeneity in defining heterogeneous slip/propagation ratios along faults. Measured changes in bed thicknesses across the structure (Fig. 10), along with fracture analysis (Figs 12 & 13), allowed us to constrain both the mechanical behaviour of the multilayer stratigraphy and the progressive deformation of these units. This detail was crucial for understanding the processes associated with structural development; measured changes in structural attributes across the structure provided important information about the influence of multilayer heterogeneity on structural processes. Ramsay (1967) demonstrated that this type of detailed, careful approach to data collection at outcrop is of critical importance in the understanding of structural processes.

The non-linearity in fault displacement that we have documented in this study (Fig. 8b) has been recorded in deformed multilayers in both extensional (e.g. Muraoka & Kamata 1983; Ferrill & Morris 2008; Roche et al. 2012) and contractional (e.g. McConnell et al. 1997; Deng et al. 2013) settings. Changes in the displacement gradient through mechanically heterogeneous strata are predicted by ductile bead theory (Fig. 15b, c; after Muraoka & Kamata 1983; Williams & Chapman 1983). Variability in slip/propagation ratios and associated changes in the geometry and behaviour of the ductile bead must thus be the rule, rather than the exception, in faulted multilayers. Where mechanical contrasts between layers are sufficient and fault displacement data are of high enough resolution, we predict these trends to be recorded at other locations. The influence of multilayer properties on fault displacement gradients and the ductile bead are essentially manifestations of the mechanical stratigraphy on deformation partitioning and structural development: multilayer heterogeneity imparts a structural style onto deformation patterns.

Many authors have noted that variability in structural style is defined by the stratigraphic template at a range of scales, from outcrop (e.g. Wilkins & Gross 2002; Roche et al. 2012; Cawood & Bond 2018) to mountain-scale (e.g. Butler 1989; Pfiffner 1993; Dominic & McConnell 1994; Cooley et al. 2011). The stratigraphic template has been shown to determine structural style throughout Variscan Pembrokeshire, from outcrop (e.g. Nicol et al. 2002) to regional scale (e.g. Smallwood 1985; Powell 1989); any model of fault-related fold development in Pembrokeshire or analogous foreland settings should thus consider the impact of mechanical stratigraphy. We have shown that the shape, size and internal structure of a multilayer ductile bead is likely to be complex: its properties will vary according to both the stratigraphic template and the deformation history before the onset of thrust faulting.

Our results show that it is not necessary for faulting to precede folding in fold-thrust systems. Physical experiments (e.g. Dixon & Tirrul 1991), conceptual models (e.g. Ramsay 1992; Mitra 2002) and field evidence (e.g. Fischer et al. 1992; Willis 1894) for folding that precedes faulting has been documented from a number of settings, including locally in Pembrokeshire (e.g. Cawood & Bond 2018), and in units of similar age and deformation history (e.g. Lloyd & Chinnery 2002). Similarly, a large volume of work exists on the formation of buckle folds, either symmetrical or asymmetrical, that form without the requirement for faulting (e.g. Currie et al. 1962; Ramsay 1974; Pfiffner 1981). It has been shown by a number of authors that hybrid models of process are more likely to capture the complexity observed at outcrop (e.g. Mitra 1990; Dominic & McConnell 1994; Erslev & Mayborn 1997; Cooley et al. 2011; Hughes & Shaw 2014) than single end-member kinematic models. Our results and interpretations agree with this conclusion: we suggest that multiple phases of deformation acted on the fold thrust at Den's Door to achieve the present-day structural geometry (Fig. 15a).

Temporal variation in the shape, size and mechanical behaviour of the deforming rock volume results in the present-day geometries observed at Den's Door. Given the spatial and temporal complexity observed here, elsewhere in Pembrokeshire and in deformed multilayers more generally, it is to be expected that single end-member kinematic models are thus not appropriate for this type of deformation. A combination of detailed measurements, hybrid kinematic models and non-kinematic approaches to structural data and observations are more likely to realistically capture deformed multilayer geometries in a meaningful way.

Using digital photogrammetry to construct and interrogate a high-resolution virtual outcrop, we have developed earlier work by Williams & Chapman (1983), Chapman & Williams (1984) and Pfiffner (1985) to reassess the link between folding and faulting at a classic fold-thrust outcrop. Based on detailed measurements from the digitally reconstructed Den's Door outcrop, we make the following conclusions:

• Faulting does not have to be a precursor to folding in fold-thrust systems. Kinematic models of fault–fold interaction emphasize the importance of faults in defining fold geometries; detailed structural measurements provide a quantified alternative model of development at this site.

• We suggest that early folding occurred before fault nucleation and growth at the study site. This was followed by a final stage of fault breakthrough, which translated the hanging wall of the thrust. Given our model of structural development, hybrid models are likely to be more appropriate for classifying natural structures, rather than single end-member kinematic models.

• We record variable thrust displacement gradients through this deformed multilayer. Variability in displacement gradient is likely to be the norm in faulted multilayers, as will changes in the size, shape and internal mechanical behaviour of any ductile bead ahead of a propagating thrust-fault tip in mechanically layered stratigraphy.

• Whether approaches to better understand deformation are kinematic or non-kinematic, the availability, accuracy and precision of measurements used for initial interpretation are of critical importance.

Thanks to Yuki Totake for field assistance and for providing images used in the virtual outcrop construction. Midland Valley Exploration is thanked for academic use of the Move 2016 software. We thank the two anonymous reviewers whose careful work and constructive comments improved the manuscript.

This study was carried out as part of a University of Aberdeen PhD supported by the NERC Centre for Doctoral Training in Oil & Gas (grant No. NE/M00578X/1 awarded to A.J. Cawood).