Thin-skinned foreland fold-and-thrust belts are characteristically found at the margins of contractional orogenic belts created by oceanic subduction or continent–continent collision (see review in McClay and Whitehouse, 2004; Figure 1A). The structural style of these fold-and-thrust belts is described as thin-skinned as it lacks any significant involvement of crystalline basement. Orogenic belts such as the Pyrenees, the Alps or the Appalachians are typically asymmetric with marginal fold-and thrust belts, often called “doubly vergent orogenic wedges”, the formation of which has been the focus of much research (e.g., Willett, 1999; Burbidge and Braun, 2002; McClay and Whitehouse, 2004; Naylor et al., 2005; Butler et al., 2019; Butler, 2020). Thrust faults and fault-related folds are now recognised as being an important part of the surficial morphology and structure of many such mountain belts. The manner in which individual thrust faults initiate, grow, and interact with other thrusts faults to build a wedge is complex and the subject of much previous research (e.g., Hardy et al., 1998; Selzer et al., 2007; Hardy et al., 2009).

Analogue models of thrust wedge formation typically provide a side-wall and map-view picture of the wedge initiation and evolution but give limited insight into the dynamic evolution of such systems (e.g., Figure 1B), although they highlight the influence of syn-tectonic sedimentation and erosion (see e.g., Storti and McClay, 1995; Mugnier et al., 1997; Wu and McClay, 2011; Graveleau et al., 2012). However, the syn-tectonic “sediments” in many of these studies are homogeneous, lacking mechanical stratigraphy and inter-layer slip-surfaces, quite unlike natural sedimentary sequences (e.g., Figure 1C). Numerical modelling studies, both kinematic and mechanical, of thrust wedge formation are varied and numerous and have produced some quite interesting results whereby the simulated thrust wedges in both weak, frictional and “strong” materials are generally quite similar to those observed on Earth and in analogue modelling studies (Willet et al., 1993; Hardy et al., 1998; Burbidge and Braun, 2002; Naylor et al., 2005; Selzer et al., 2007; Hardy et al., 2009; Simpson, 2010; Dean et al., 2013; Fillon et al., 2013). However, many of these numerical studies have not included syn-tectonic sediments or “growth strata”. Those that do either treat them as passive markers of fault/fold growth (e.g., Hardy and Ford, 1997), have a very coarse resolution (100 s of m) such that detailed, local interactions of faulting and sedimentation are not discernible (Selzer et al., 2007; Fillon et al., 2013), or are accretionary wedge-type finite element simulations where only hinterland sediment supply and diffusion are simulated (e.g., Selzer et al., 2007; Simpson, 2010).

Conceptually, experimentally and observationally, syn-tectonic sedimentation (and erosion) is thought to have a significant influence on local thrust tectonics and wedge development in general (Nichols, 1987a; Nichols, 1987b; Storti and McClay, 1995; Storti et al., 1997; Mugnier et al., 1997; Benesh et al., 2007; Bigi et al., 2010; Fillon et al., 2013; Dean et al., 2015; Butler et al., 2019; Butler, 2020; Noda et al., 2020). This syn-tectonic sedimentation may be the result of local erosion, transport and sedimentation, sediment supply from the hinterland, or more regional basinal processes such that thrusts propagate into a growing (aggrading) foreland sedimentary sequence. Whilst previous numerical and analogue studies have addressed the interaction of sedimentation and tectonics, their resolution has typically not allowed the examination of detailed, single structure, tectono-stratigraphic relationships. Here, we present a high-resolution (average 9.7 m element radius) discrete element model of thrust wedge development (Figure 1D). We use frictional-cohesive assemblies of discrete elements, grouped in multiple, pre-defined layers separated by weak interfaces, to represent probable, mechanically variable, cover rheologies. This also allows for the development of folds and fault-related folds within the thrust wedge. Homogenous frictional-cohesive models produce results that are very similar but are dominated by shear zones rather than fault-related folds. The syn-tectonic strata added here are also frictional-cohesive and we call them “mechanical growth strata” as they interact with the growing thrust fault array. The model allows us to examine the detailed evolution of thrusts and thrust wedges as they propagate into a growing (aggrading) basinal sedimentary sequence in the foreland. Implications of the results for the interpretation of coupled tectono-sedimentary systems are discussed.

Here we use a 2D discrete element numerical model to simulate the growth of a thrust wedge in a frictional-cohesive cover overlying a weak, horizontal, decollement. Modelling of deformation to high strain is an ideal candidate for the discrete element technique as it is well-suited to studying problems in which discontinuities (shear-zones, faults, fractures, voids, etc.) are important. It also allows deformation involving unlimited relative motions of individual elements and complex, abrupt and changing boundary conditions (Cundall and Strack, 1979; Finch et al., 2004:; Egholm et al., 2007; Hardy, 2008, Hardy, 2011, Hardy, 2019a; Thompson et al., 2010), particularly germane to the simulation/formation of thrust faults and wedges. The discrete element code used here is “cdem2D” which has been used previously to undertake a wide variety of structural geology modelling in both 2D and 3D (as cdem3D). The basic numerical methodology and its application to a variety of problems (fault-propagation folding, orogenic wedge growth, caldera collapse, dike intrusion on Mars, viscous flow, Gilbert deltas, etc.) has been published previously (e.g., Botter et al., 2014, Botter et al., 2016; Hardy, 2008; Hardy et al., 2009; Hardy, 2016; Hardy, 2018; Hardy, 2019a; Hardy, 2019b). The details of this code applied to frictional and frictional-cohesive materials can be found elsewhere (e.g., Hardy et al., 2009; Hardy 2011, Hardy, 2015, Hardy, 2019b, Hardy, 2021), however an overview is given below. Here, we will illustrate the use of this modelling scheme to simulate the evolution of thrust faults/wedges and ongoing syn-tectonic, “growth”, sedimentation.

The discrete element technique, in common with other numerical techniques, has both advantages and disadvantages when considering its application to any specific geological problem. As discussed before, modelling of large deformations and the possibility of applying complex boundary conditions are critical for modelling thrust wedges. However, one disadvantage of the technique lies in the necessary, but time-consuming, calibration of micro-particle parameters to bulk physical properties (cf. Egholm et al., 2007; Botter et al., 2014). The interaction of many tens of thousands of particles, both locally and globally, also leads to situations wherein our ability to explain why a particular fault or fracture grew at the expense of a neighbouring one is limited. Such issues also exist in analogue modelling where repeated experiments under the same boundary conditions are reasonably reproducible but not in the finer details of the fault and fracture systems (e.g., van Gent et al., 2010; Schreurs et al., 2016). Computational limitations on element size and/or model resolution, whilst important previously, are no longer a particular concern due to rapid advances in computational power and the parallelisation of discrete element codes.

The discrete-element model described here is a variant/development of the lattice solid model (LSM) (Mora and Place, 1993, Mora and Place, 1994; Place et al., 2002). Here, a sedimentary rock is treated as an assembly of circular elements, which interact as if connected by breakable elastic springs (bonds) and that undergo motion relative to one another. The behaviour of the elements assumes that the particles interact through a “repulsive–attractive” force (Mora and Place, 1993) in which the resultant (normal) force, F _n , is given by: F n = K n ( r − R ) r < r 0 i n t a c t   b o n d F n = K n ( r − R ) r < R b r o k e n   b o n d F n = 0 r ≥ R b r o k e n   b o n d (1) Here, K _n is the elastic constant (normal spring stiffness) of the bond, R is the equilibrium distance between the particles (sum of particle radii), r ₀ is a breaking threshold, and r is the current distance between the particle pair. Particles are bonded until the distance (r) between them exceeds the defined breaking threshold (r ₀ ) at which time the bond breaks. The force acting on a bond at this threshold represents the force necessary for a bond to fail or yield, or, alternatively, can be cast as the stress acting on a particle’s bond at failure. After this breaking threshold, the particle pair experiences no further attractive force and the bond is irreversibly broken. However, if the two particles return to a compressive contact (r < R), a repulsive force still acts between them. This is the case with the frictional-cohesive models discussed herein where all elastic bonds are initially broken.

For a frictional-cohesive material, in addition to treating the normal force (F _n ) between particles, we also calculate the tangential (shear) force, F _s , as a result of displacement (X _s ) perpendicular to the vector connecting the particles’ centroids. This frictional force acts in a direction opposite to that of the relative tangential velocity and is modelled as a threshold-limited elastic spring with a cohesive force term (C ₀ ) in parallel with that used to calculate the normal forces (cf. Cundall and Strack, 1979; Mora and Place, 1994). The magnitude of this force is limited to be less than or equal to the shear force allowed by Coulomb friction: F s = K s X s + C 0 F s m a x = μ F n F s = F s m a x i f F s > F s m a x (2) where K _s is the elastic constant (shear spring stiffness) of the contact, F _smax is the maximum (limiting) shear (frictional) force, F _n is the normal force at a contact, and µ is the inter-particle coefficient of friction. If a contact is ‘‘lost’’ between two touching elements (i.e., they separate), then all the forces between the elements are set to zero. The total elastic force, F _i,α exerted on a particle is thus obtained by summing the normal and tangential forces on each contact/bond that links a specific particle to its neighbours, calculated by: F i , α = ∑ j = 1 , α f i , j (3) where f _i,j is the elastic force (normal and shear) experienced by particle i from its neighbouring particle j. However, we include a viscous damping term (proportional to the particle velocity) that acts to dampen reflected waves from the rigid edges/boundaries of the model, preventing a build-up of kinetic energy within the closed system, a standard technique to ensure numerical stability (cf. Mora and Place, 1994; Place et al., 2002).

Finally, gravitational forces, F _g , acting on each element are calculated in the vertical direction, increasing the vertical stress with depth. Therefore, the total force (F) on any particle is given by: F = F i , α − υ x · + F g (4) where υ represents the dynamic viscosity and ẋ is the velocity of the particle. At each time step, the particles are advanced to their new positions by integrating their equations of motion using Newtonian physics and a velocity-Verlet-based scheme (Allen and Tildesley, 1987). In discrete element models such as the one used here, parameters such as strength, coefficient of friction etc of the assembly are emergent properties and do not relate directly to the element properties. This relation must be assessed through the use of angle of repose and/or unconfined/confined biaxial numerical tests, a process known as calibration (cf. Oger et al., 1998; Finch et al., 2004; Holohan et al., 2011; Botter et al., 2014).

The key additional feature presented here is that, unlike simple, homogeneous, discrete element assemblies, we can choose to have element-element contact properties (such as cohesion, friction or bonding) only activated within a set of (predefined) layers and not between these layers. In this paper, we incorporate this heterogeneous layering into the models in order to facilitate flexural slip deformation along layer interfaces and thus better approximate the typically anisotropic mechanical stratigraphy of sedimentary cover rocks, be this a layered sandstone-shale sequence or salt layers etc (cf. Figure 1C). These inter-layer contacts have no cohesive bonds or friction, in order to make these interfaces as weak as possible. These interfaces cannot be frictionless in the true sense, due to the “roughness” of circular elements and thus there is always some residual, intrinsic friction at the interface. However, it remains much weaker than a frictional-cohesive interface. All other numerical techniques and approaches used on the elements are identical to those discussed in previous papers to which the interested reader is referred (e.g., Hardy et al., 2009; Hardy, 2011, Hardy, 2015).

A second feature presented here is the use of the discrete element technique not only to model deformation of a pre-existing sedimentary sequence but also to model the deposition (and subsequent deformation) of a “growth” sedimentary sequence deposited during active tectonics. After syn-tectonic sediments are added to the evolving model, these elements must be allowed to equilibrate/adjust to a final, initial deposited state. This is done by pausing any external boundary movement after discrete element addition for a predefined (small) number of time-steps until the added elements are in a quasi-static location and effectively stable. This is a crucial computational step. The experiment then continues as before but with a new sediment load added to the system. Compaction per se is not considered, but the greater load, as a result of sedimentation, changes forces and element packing at depth. The specific heuristic used here is to add discrete elements of variable size (identical to the pre-growth and with an average radius of c. 9.7 m), filling up any available submarine accommodation space beneath a defined baselevel. These elements are chosen randomly and then added sequentially with each additional element attempting to minimise the remaining accommodation space until it is filled. They are carefully placed to ensure that they maximise the filling of available accommodation space while just touching the preexisting deposited elements. These growth elements can either be frictionless to represent weak, unconsolidated (un-lithified) sediments or can be frictional-cohesive to represent rapidly lithified (via compaction or cementation, and thus much stronger) growth strata. Here, we consider the growth elements to be frictional-cohesive and as strong as the pre-growth elements. In addition, they can also experience flexural-slip between their defined layers. We do not consider the transport of these elements or their erosion in this study. The spatial resolution chosen here is arbitrary but it is appropriate to simulate the thrust wedges and to achieve realistic runtimes. However, with current rapid advances in computing power, model resolution is becoming less of an issue.

In this paper we have simulated, using the discrete element method, the growth of thrust wedges under increasing amounts/rates of syn-tectonic sedimentation. These experiments have allowed us to assess the influence that such sedimentation has on tectonics rather than the usual perspective of “tectonics and sedimentation” where the growth strata are some kind of passive recorder of fault or fold activity. We can perhaps think of them as “mechanical growth strata” as they interact with, and influence, the growing thrust wedge. While similar analogue and numerical studies have been undertaken previously, here we have examined the detailed interaction of sedimentation and tectonics at a high-resolution, allowing structural and stratigraphic relationships to be observed at the local, almost outcrop, scale.

The first-order result of these experiments is that they confirm that, indeed, syn-tectonic sedimentation can have a profound effect upon thrust wedge development (cf. Wu and McClay, 2011; Graveleau et al., 2012; Butler et al., 2019; Butler, 2020). In the experiment of thrust wedge development without syn-tectonic sedimentation (experiment 1), a forward-breaking sequence of thrust activity is seen: producing a typical thrust-wedge structure, consistent with many analogue, numerical and conceptual models (Figure 2; cf. Storti and McClay, 1995; Burbridge and Braun, 2002; Wu and McClay, 2011). In general, the inclusion of syn-tectonic sedimentation produces thrust wedges composed of fewer major forward-vergent thrusts and with only minor thrust activity in the foreland (experiments 2 to 5; Figures 3–6). In most models the sequence of thrust activity is complex and not simply forward-breaking. Very high syn-tectonic sedimentation (experiment 5) results in the formation of a single basin-bounding thrust fault and no thrust-wedge per se (Figure 6). Thus, with increasing syn-tectonic sedimentation, the frontal thrust has greater offset/displacement and overrides a much thicker package of earlier syn-tectonic sediments. Clearly, tectonics remains the driving force behind contraction and the formation of a thrust wedge (in nature and in the models). However, sedimentation clearly has a more local effect upon the nature of the wedge itself and particularly on structural style.

In a sense, these results for the dynamic development of thrust wedges are predictable/understandable from a minimum work perspective (e.g., Hardy et al., 2009; McBeck et al., 2020). In straightforward terms, minimum work is based on the principle that there is a continual tradeoff/balance between frictional and gravitational work during wedge growth, and that the “path” of minimum work or resistance is favoured during deformation (e.g., Masek and Duncan, 1998; Del Castello and Cooke, 2007; Cooke and Madden, 2014; McBeck et al., 2020). Without syn-tectonic sedimentation (or erosion), the minimum work model predicts propagation of thrust fault activity towards the foreland through time as a result of increasing topographic (gravitational) loads and frictional work on deformed hinterland faults. As regards our models, stated simply, with increasing syn-tectonic sedimentation, there is a “penalty” for propagating forward along the decollement to create outboard structures rather than simply uplifting or imbricating close to the moving end-wall. The physical basis for this “penalty” is that more (frictional and gravitational) work is required to displace a thicker pile of sediment, which thus favours reactivation of, or continued activity on, pre-existing faults that are overlain by less overburden (i.e., a thinner pile of sediment). Consequently, these experimental results support the idea that minimum work optimization techniques can provide critical insights into crustal tectonic evolution (e.g., Hardy et al., 1998; Masek and Duncan, 1998; Del Castello and Cooke, 2007; Cooke and Madden, 2014; McBeck et al., 2020).

A simple chronology, or sequence, of fault activity in our models is not observed. Out-of-sequence thrusting and folding occurs during overall foreland propagation of deformation. The detailed spatial and temporal activity of faults is complex and reflects the interaction between the development of topography, the contrast between basal (décollement) and internal coefficients of friction, and the loading effect of sedimentation. The key factor is the relative rates of sedimentation and tectonic deformation. As with all complex systems, these two rates cannot be considered independently. In particular, and the specific focus of this study, rates of sedimentation are shown to be important controls on fault activity both spatially and temporally. Sedimentation acts to increase gravitational and frictional work on a given fault, and therefore has the potential to “blanket” faults and render them inactive. However, sedimentation clearly does not “shut down” all faults, as many of the frontal thrusts here show continued activity while overriding the sedimented footwalls. Model results emphasise the complex, and subtle, feedbacks that can exist between tectonic and surficial mass transport processes.

At a more local, outcrop, scale the effects of contemporaneous sedimentation and tectonics are quite marked in our experiments. Figure 10 shows detailed zooms of the “thrust-front” in each of the models with syn-tectonic sedimentation. What is obvious is that a wide range of distinct geometries and relationships can be produced due solely to changes in the magnitude of syn-tectonic sedimentation. The first, obvious, effect is the increase in displacement/offset on one or two major faults with increasing syn-tectonic sedimentation. Secondly, and as a consequence, is the greater area and thickness of syn-tectonic sediments that are over-ridden at the thrust front. Thirdly, the main, frontal, thrust in each experiment also becomes steeper with increasing syn-tectonic sedimentation (compare Figures 10A–D).

What is also clear from our results is that the manner in which sediments are introduced to analogue or numerical models is important. Many studies have included only hinterland sediment supply and/or diffusive erosion/sedimentation, or sedimentation that thins and is wedge-like towards the foreland (e.g., Selzer et al., 2007; Wu and McClay, 2011). This approach is entirely reasonable, but in a sense produces very straightforward, and perhaps biased, results. This is because by adding sediments to the wedge in this manner, deformation is being driven towards the foreland. On the other hand, basinal, banking aggradational sedimentation as implemented herein is a very different set-up/approach that produces a buffering effect on thrust propagation rather than a driving effect. Natural systems, most likely, lie somewhere between these extremes and we must consider both possibilities when examining exhumed, ancient basins. The objective of this study was not to compare our results directly to specific structures and, in addition, the total amount of displacement in our experiments (2.5 km) is low compared to natural orogenic belts. However, some general comparisons can be made. At a local level, Nichols (1987a, b) noted that in the Western External Sierras of the Pyrenees (N. Spain) complex, high-displacement frontal thrust structures were associated with localised thrust-front fanglomerates appearing to “buffer” forward propagation. At a much larger scale, Butler (2020) highlighted that in the northern Apennines thrust systems are ramp-dominated where syn-tectonic sedimentation was high. In contrast, in the southern Apennines, the Lagronegro allochthon is a high displacement structure (>65 km) as the foredeep basin received little sedimentation. Our experiments are consistent with these observations.

The excess-area analysis (Epard and Groshong, 1993) of our experimental results has tested its applicability in analysing thrust-fault arrays where we know the boundary conditions, depth to detachment and displacement (i.e., the answer). This serves to illustrate that its application to typical (incomplete) surface and subsurface geological data can provide robust estimates of several key parameters. Finally, these model results are only a first step in unravelling the complex interaction between sedimentation and tectonics. There are many other internal and external facets which can be important in conditioning the growth of a thrust wedge - climate change, erosion, isostasy, more complex mechanical stratigraphy, basal friction etc. These and other topics are subjects of ongoing research.