## Abstract

The growth of normal faults in mechanically layered sequences is numerically modelled using three-dimensional Distinct Element Method (DEM) models, in which rock comprises an assemblage of bonded spherical particles. Faulting is induced by movement on a pre-defined normal fault at the model base whilst a constant confining pressure is maintained by applying forces to particles lying at the model top. The structure of the modelled fault zones and its dependency on confining pressure, sequence (net:gross) and fault obliquity are assessed using various new techniques that allow (a) visualization of faulted horizons, (b) quantification of throw partitioning and (c) determination of the fault zone throw beyond which theoretical juxtaposition sealing occurs along the entire zone length. The results indicate that fault zones become better localized with increasing throw and confinement. The mechanical stratigraphy has a profound impact on fault zone structure and localization: both low and high net:gross sequences lead to wide and relatively poorly localized faults. Fault strands developing above oblique-slip normal faults form, on average, normal to the greatest infinitesimal stretching direction in transtensional zones. The model results are consistent with field observations and results from physical experiments.

Faults are generally mapped and modelled as single surfaces across which there is a discrete displacement. In reality faults are complex zones comprising multiple slip surfaces that contain variably deformed rock volumes, ranging from intact fault-bound lenses to fault rock (Childs *et al.* 2009). This complexity in fault zone structure can be of significance in a variety of application areas, for example, by providing across-fault connectivity of flow units. Fault zone complexity and architecture will vary depending on the nature of the faulted sequence and the prevailing deformation conditions (e.g. Patton *et al.* 1998; Micarelli *et al.* 2005; Ferrill *et al.* 2007; Ferrill & Morris 2008; Kettermann & Urai 2015). Empirical constraints on the 2D and, in particular, 3D geometry and content of fault zones are, however, relatively sparse, in the sense that insufficient data are available on fault zones developed within multi-layered sequences, with different stacking patterns and rheological properties, and under different deformation conditions.

The published literature contains many accounts of conceptual models for the evolution of fault zone structure based mainly on outcrop studies of fault geometry (see review articles by Crider & Peacock 2004; Wibberley *et al.* 2008; Faulkner *et al.* 2010). These are complemented by a range of numerical models that attempt to replicate the processes invoked in these conceptual models (e.g. Schöpfer *et al.* 2016). Numerical modelling of fault zone growth, which requires simulating the formation of fractures, followed by detachment, rotation and comminution of fracture-bound volumes, is however cumbersome or impossible in many of the standard modelling approaches used in geomechanics (e.g. finite elements). The Distinct Element Method (DEM) is ideally suited to the modelling of fault zones but, as detailed in this paper, requires significant effort to achieve the realistic rheological properties that are routinely defined by constitutive equations in other modelling approaches.

The purpose of the present study is to investigate the impact of a selection of factors, such as confining pressure and layer stacking (net:gross), on the structure of normal faults and associated sequence juxtaposition using 3D DEM modelling. In earlier studies we have already shown the usefulness of this numerical technique to model fault growth in layered sequences in two dimensions (Schöpfer *et al.* 2006, 2007*a*, *b*, 2009*a*). These earlier studies are somewhat limited because fault propagation and fault zone evolution are inherently 3D processes and, for example, 2D modelling in the fault slip direction fails to capture features such as lateral fault propagation and relay ramp development. Therefore, while the DEM models presented in this study are in some respects a 3D extension of earlier models, they represent a significant advance that incorporates a wider range of processes allowing consideration of the lateral variability in structure along fault zones.

In this paper we first provide an overview of the DEM and existing models of normal faulting. Then, for the 3D normal faulting models presented here, we describe in detail the model materials and their properties, the model generation and boundary conditions. These models are used to investigate the importance of confining pressure, mechanical stratigraphy (net:gross) and slip direction on fault zone structure and localization. We conclude that many of the model results are consistent with field observations and results from physical experiments.

## Distinct Element Method

The DEM comprises a broad class of methods that models finite displacements and rotations of discrete bodies (Cundall & Hart 1992). The DEM is capable of modelling the growth of discontinuities, such as faults and joints, without the limitations of continuum descriptions (Cundall 2001). The elements interact with each other via a force displacement law and can be of arbitrary shape, circular discs and spheres being the most common ones in 2D and 3D, respectively.

### Two-dimensional DEM models of normal faulting

To our knowledge, Saltzer & Pollard (1992) were the first to apply the DEM to model normal faulting in sedimentary sequences. Although these early models were run using 900 non-bonded particles only, they already highlighted the potential of the DEM to model faulting. Finch *et al.* (2004) investigated the effects of cover strength and fault dip on the style of fault-propagation folding above a normal fault using layers comprising bonded particles. Schöpfer *et al.* (2006) simulated the propagation of a normal fault through a vertically confined, mechanically layered sequence. In later studies, Schöpfer *et al.* (2007*a*, *b*, 2009*a*) explored the impact of sequence strength, confining pressure and net:gross (fraction of strong layers in the sequence) on fault zone structure. Egholm *et al.* (2007) reproduced sandbox experiments of basement-induced normal faulting using a modified DEM approach (Egholm 2007), which was also used for modelling normal faulting in sand–clay sequences (Egholm *et al.* 2008). Abe *et al.* (2011) and Hardy (2011, 2013) modelled the propagation of normal faults in vertically unconfined bonded particle models. Botter *et al.* (2014) simulated normal fault growth using non-bonded particle models and generated, on the basis of volumetric strain and empirical relations, synthetic seismic images.

### Three-dimensional DEM models of normal faulting

Owing to the high computational demand of 3D particle simulations few studies on normal faulting exist in the literature. Walsh *et al.* (2001) sheared a cube of bonded particles and studied how displacement is accommodated in a volume containing a predefined fault system. Imber *et al.* (2004) sheared a cube containing two overlapping predefined faults and investigated relay ramp growth and fault linkage. Carmona *et al.* (2010) combined the DEM with a process-based model of sedimentation to model syntectonic sedimentation. Recently, Schöpfer *et al.* (2016) presented 3D DEM models of normal faults in layered sequences that underpin a geometric model for the evolution of the internal structure of fault zones (Childs *et al.* 2009).

In summary, numerous studies on the growth of normal faults exist in 2D, but 3D DEM model studies of the growth of normal faults are rare. An overview of the DEM, including definition of particle/bond properties, the calculation cycle and model generation, is provided in the remainder of this section.

### Definition of particle and bond properties

We use the commercially available Particle Flow Code in Three-Dimensions (*PFC3D*; Itasca Consulting Group 2008), which implements the DEM. Rock is represented by an assemblage of randomly packed spherical particles that are bonded together (Fig. 1a). These bonds can break if either their normal or shear strength is exceeded, which corresponds to fracturing of the material. Table 1 summarizes the microparameters that need to be assigned for a (parallel) bonded particle model. A more detailed definition of these parameters can be found in Potyondy & Cundall (2004).

Elasticity is captured by allowing overlap at particle–particle contacts. A linear contact law is used so that the normal contact force is proportional to the amount of overlap; the proportionality constant is the normal stiffness *K*_{n} (Fig. 1b), which is inherited from the two contacting particles assuming that both particles’ stiffness act in series, i.e. with and being the particle normal stiffnesses. By relating the behaviour of a contact between two particles to that of an elastic beam with its ends at the particle centres, an expression can be obtained for the particle normal stiffness as a function of the so-called ‘contact Young's modulus’ *E*_{c}, so that *k*_{n}=4*rE*_{c}, where *r* is the particle radius. This modulus-stiffness scaling ensures that the bulk elasticity becomes particle size independent (Potyondy & Cundall 2004). The shear force is computed in an incremental fashion (Fig. 1b). Each subsequent relative shear-displacement increment produces an increment of elastic shear force that is added to the current shear force (which is, when a contact is formed, initially zero). The constant that relates shear displacement and shear force is the contact shear stiffness *K*_{s}, which is also inherited from the two contacting particles assuming that both particles’ stiffness and act in series. The shear stiffness *k*_{s} is assigned to each particle via the particle stiffness ratio, *k*_{n}/*k*_{s}. The contact stiffness ratio (*K*_{n}/*K*_{s}) is an important factor that influences, for example, Poisson's ratio of the model (*K*_{n}/*K*_{s} should be >1, otherwise a negative Poisson's ratio may be obtained; Bathurst & Rothenburg 1988; Rothenburg *et al.* 1991; Wang & Mora 2008). The absolute shear force at a contact is limited by Coulomb friction so slip occurs at a contact when the ratio of contact shear to normal force is equal to μ_{c} (Fig. 1b). If the two contacting particles have differently assigned friction coefficients, the smaller value is used as a contact friction coefficient.

Particles can be bonded with massless, linear elastic ‘cement’, which acts in parallel with the above-described linear contact law, hence the bond model is referred to as the parallel bond model (Potyondy & Cundall 2004). A parallel bond can be envisioned as a set of elastic springs uniformly distributed over a circular cross-section lying on the contact plane. In the present study these disc-shaped bonds have a radius equal to the smaller of the two bonded particles (Fig. 1a). Relative motion at parallel bonded contacts causes axial and shear forces as well as bending and twisting moments to develop (Fig. 1c). The maximum normal and shear stress acting on the bond periphery can be computed via beam theory (see equation 16 in Potyondy & Cundall 2004): the bond normal stress increases with increasing normal force and bending moment, whereas the bond shear stress increases with increasing shear force and twisting moment (Fig. 1c). A consistent means of setting the elastic bond microparameters is to define the bond Young's modulus and the ratio of bond normal to shear stiffness ratio (). The modulus-stiffness relation for the bonds is , where *r _{A}* and

*r*are the radii of the two contacting particles. Once either the maximum normal or shear stress carried by the bond exceeds the normal (tensile) or shear strength, and respectively, the bond breaks (Fig. 1c) and is immediately removed from the system along with its accompanying force, moment and stiffnesses. The fact that the bond and the slip model act in parallel has two important consequences: (a) slip does occur at bonded particle–particle contacts if the maximum shear force is exceeded; and (b) in compression load is carried by both the particles and the bond, whereas in tension load is carried by the bond only. These two features lead to pronounced non-linear elastic bulk behaviour (Schöpfer & Childs 2013).

_{B}### Calculation cycle and damping

The calculation cycle in *PFC* is a time-stepping algorithm that requires the repeated application of the second law of motion to each particle (force=mass times acceleration), a force displacement law to each contact (e.g. Fig. 1b), and a constant updating of particle and wall (rigid boundaries) positions. The equations of motion are integrated in *PFC* using a centred finite-difference scheme, which requires a time-step. *PFC* automatically calculates the critical time-step (which depends on mass and stiffness of the system) that is necessary for ensuring that the solution produced is stable and uses an actual time-step which is a fraction of the critical time-step for safety reasons.

*PFC* provides a true dynamic solution. However, deformation under quasi-static conditions requires dissipation of energy via damping. In the present study, so-called local damping is used, where a damping-force term is added to the equations of motion (damping constant is 0.7; Potyondy & Cundall 2004). For compact assemblies, this form of damping is the most appropriate form to establish equilibrium and to conduct quasi-static deformation simulations.

## Numerical rock laboratory and properties of model materials

As described above the properties of a DEM model material are defined by a suite of particle and bond properties (Table 1 and Fig. 1). The bulk mechanical behaviour emerges from the interaction of particles via bonds and their contacts. Bulk mechanical properties have to be determined using numerical rock laboratory experiments, such as unconfined and confined compression tests. These numerical tests on DEM model materials are typically performed using rigid, frictionless loading and confinement platens (Potyondy & Cundall 2004), which can lead to significant boundary effects. An improved loading and confinement method, recently developed by Schöpfer & Childs (2013), is adopted for determining the properties of the materials used in the present study. In this section, first the particle/bond properties used are given, followed by a description of the numerical rock laboratory boundary conditions and the bulk mechanical properties of the model materials.

### Particle and bond properties of model materials

In the present study, relatively simple DEM materials are used which are entirely defined using the set of parameters listed in Table 1. A drawback of using simple model materials is that, although the form of the failure envelope of natural rock can be reproduced (Schöpfer *et al.* 2013), the ratio of unconfined compressive to tensile strength and the friction angle is lower in the DEM material than in most rocks. The principal reason for the different behaviours of rock and DEM materials is that the latter comprises perfect spheres that offer no resistance to rolling. There are various techniques that have been implemented to improve the bulk properties of parallel-bonded particle models (Schöpfer *et al.* 2009*b*), but these often lead to a significant increase in model run time and/or lead to otherwise unrealistic behaviour.

In the normal faulting models two different materials are used to represent weak and strong rocks. The different particle and bond properties of these two materials are given in Table 1. To enable normalization of model parameters for comparison between models and the application of model results to natural conditions, the two materials used have identical stiffness/strength ratios and identical bond strength ranges (expressed as the coefficient of variation, which is standard deviation normalized by mean value).

### Numerical rock laboratory boundary conditions

A limitation of existing loading procedures that are used for determining bulk mechanical properties of DEM model materials is the use of rigid platens to apply both loading and confinement. Loading under unconfined (and confined) uniaxial compression or extension tests can however be achieved using distortional periodic space (Schöpfer & Childs 2013). For this type of testing a bonded particle model is first generated within a cubic periodic domain (Fig. 2a) using the dynamic model generation procedure described in Potyondy & Cundall (2004). Then a core (cylindrical sample) is carved so that the sides of the cylinder do not straddle the periodic boundaries (Fig. 2b). The top and bottom of the cylinder are therefore periodic, i.e. the cylindrical sample is effectively infinitely long. The core diameter and resolution are identical to the thickness and resolution of the layers used in the normal faulting models (ratio of layer thickness to average particle diameter is 5.3).

The core sample can either be shortened (uniaxial compression test) or lengthened (uniaxial extension) by distorting the periodic space. Since loading platens are absent, the applied stress cannot be calculated by force/area, so a so-called ‘measurement sphere’ is located within the sample in which the average stress tensor is computed using particle contact forces and volumes (Potyondy & Cundall 2004). A sensitivity study revealed that at strain rates of 0.01 s^{−1} or less no significant changes in the bulk strength occur, suggesting quasi-static conditions. These unconfined tests provide the sample Young's modulus (*E*), Poisson's ratio (ν), the unconfined compressive strength (UCS) and the tensile strength (*T*).

Confinement of the samples without the use of rigid platens is achieved by applying inwards-directed forces to the lateral particles that comprise the surface of the cylinder, thus simulating a jacket (Fig. 2c). The lateral surface particles and the direction and magnitude of the applied forces are found using the so-called ‘shining-lamp’ algorithm, which is explained in detail in Appendix A. These confined tests provide the strength under confined compression or extension loading conditions and are important for defining the failure envelope of the model materials.

### Bulk mechanical properties of model materials

Stress–strain curves for the two model materials (‘strong’ and ‘weak’) tested under both triaxial compression and triaxial extension test conditions are shown in Figure 2d. Since both particle/bond modulus and bond strength in the strong material are 10 times higher than in the weak material (see Table 1), the stress–strain curves are much steeper and reach higher stress differences in the strong material. The secant Young's modulus is extracted from these curves at an axial strain of 0.125 and 0.0625% for the triaxial compression and extension tests, respectively, and plotted v. confining pressure in Figure 2e. These tests show that under unconfined conditions Young's modulus in extension is about half the modulus in compression because load in compression is carried by both particles and bonds, whereas in tension load is only carried by bonds. With increasing confining pressure the difference in Young's modulus in extension and compression decreases, because the proportion of contacts that carry compressive load increases. Poisson's ratio also exhibits a confining pressure dependency (not shown here; see Schöpfer & Childs 2013) and is, for the conditions used in the normal faulting models, *c.* 0.18.

Another striking feature is that, for the range of confining pressures used, the strong material exhibits pronounced post-peak stress weakening, whereas the weak material exhibits little or no weakening (Fig. 2d). In fact, it appears that the strong material is brittle and the weak material is ductile. In the following we shall therefore use the definition that brittle behaviour exhibits a significant stress drop after peak stress, while ductile behaviour shows no such strain softening. The reason for this different behaviour is that the majority of tests on the weak material are conducted at confining pressures that are greater than the pressure at which the brittle–ductile transition, as defined by Byerlee (1968), occurs. The emergence of this type of brittle–ductile transition, or more precisely the brittle–cataclastic-flow transition, was recently explored by Schöpfer *et al.* (2013) in 3D stress space using the DEM. At low confining pressure peak-stress data fall onto a non-linear envelope that intersects (at the transitional pressure) a linear envelope (Fig. 2f). Under triaxial extension stress conditions this kind of brittle–ductile transition occurs in the weak material at a confining pressure of *c.* 21 MPa (Fig. 2f). The shape of the failure envelope and the location of the brittle–ductile transition characteristic of the numerical model materials used in this study can be better appreciated by plotting peak stress data normalized by the bond strength (Fig. 2g). This type of normalization can, in principle, also be applied to the normal faulting models presented in later sections. Identical results are expected if all intrinsic (modulus, strength) and extrinsic (confining pressure) parameters in stress units are scaled by the same factor. For clarity, however, all model results are presented in terms of their actual, rather than dimensionless, mechanical properties and applied confining pressures.

### Discussion of material properties

Before presenting the normal faulting results it is worthwhile discussing the results of the numerical rock laboratory experiments, summarized in Figure 2. The purpose of determining the bulk mechanical properties of the DEM model materials is not to calibrate the materials against a certain rock. We rather try to ensure that some properties (Young's modulus, unconfined strength) are realistic, i.e. comparable to rock. The range and variability of rock mechanical data is however enormous, e.g. sandstones can have a uniaxial strength ranging from almost zero up to 300 MPa (Lockner 1995). A useful dimensionless parameter is the modulus ratio, defined as the ratio of Young's modulus to UCS, for which values between 200 and 500 are considered as ‘average’ (Deere & Miller 1965). The modulus ratio of the two model materials is *c.* 400 and hence typical for rock. Moreover, in the classification scheme proposed by Deere & Miller (1965), the UCS values of the strong and weak material fall within ‘high strength’ and ‘very low strength’ fields, respectively. However, the ratio of unconfined compressive strength to tensile strength (UCS/*T*) is too low in our model materials (*c.* 3.8) when compared with natural rocks, which have UCS/*T*>10 (Schöpfer *et al.* 2009*b*). Since it is not possible using the present approach to match UCS/*T*, the question arises which one of the two strengths is believed to be more important under normal faulting conditions. In fact both *T* and UCS are important. The tensile strength controls the stress at which localization occurs (in the form of opening mode fractures under the conditions used in the present study, indicated by arrows in Fig. 2f). The compressive strength controls to some extent the comminution of fault-bound volumes. Practically speaking, if we were to match *T* then the very early stage of faulting would be realistic, but later fault zone evolution would be characterized by pervasive fracture and rapid comminution of fault-bound volumes. To generate interesting (and realistic) looking fault zone structure, we chose to match the compressive strength with the knowledge that the equivalent tensile strength is too high by a factor of >*c.* 2.5.

## Geometry and boundary conditions of normal faulting models

The normal faulting modelling is basically a three-step process; (a) model generation; (b) model confinement; and (c) application of the normal faulting boundary conditions.

### Model generation

In mesh-based continuum mechanics methods, such as the Finite Element Method (FEM), the geometry of the mesh is of utmost importance. In the DEM, when applied to the simulation of rock fracture, an initial dense particle packing is crucial. Moreover, in the FEM smaller elements are often used in volumes of anticipated large strain. A similar technique, in which smaller particles are used in the volume in which higher resolution is required, is used in the current DEM approach (Fig. 3a).

In each normal faulting model three different materials are used (Fig. 3b): the periodically layered sequence in which faulting is analysed comprises strong and weak beds. This sequence is overlain by a layer of non-bonded particles that have the same particle properties as the weak material and act as a buffer layer between the layered sequence and the particles that apply the confining pressure at the top of the model.

The model is generated by first filling the model volume, bound by planar, frictionless walls, with a dense packing of frictionless spherical particles with radii chosen randomly from a uniform size distribution. Particles located in the volume above the future fault are refined by means of nested regions in which the refinement level *N _{R}* increases inwards. Each particle in a refinement region is replaced by two smaller particles that have the same total volume as the original particle, and this replacement process continues until the specified refinement level

*N*is reached, so that the average radius decreases by

_{R}*r*=

*r*

_{0}exp[−ln(2)

*N*/3], where

_{R}*r*

_{0}is the average radius of the base material. The lateral boundaries of these refinement regions are parallel to the future pre-defined basal fault (Fig. 3a). Particle stiffnesses are assigned via the aforementioned (particle size-dependent) modulus–stiffness scaling using contact Young's modulus values assigned to each layer. Then the system is cycled to equilibrium, a low isotropic stress is installed and the number of ‘floating’ particles (with fewer than three contacts) is reduced (see Potyondy & Cundall 2004, appendices A5 and A6). Finally, particle friction coefficients are assigned. Bonds are installed after confinement, as described below.

### Model confinement

Once a model is generated, it is confined to the desired stress level that will be held constant throughout the model run. This confining pressure simulates the load of an overlying rock volume. We assume that this hypothetical overlying rock volume is in an isotropic stress state and therefore exerts no shear tractions to the upper model boundary, so that only normal tractions act on the upper boundary. Gravity is absent in the models (stress differences arising from the vertical stress gradient are negligible at the scale of our models) so that the magnitude of the confining pressure is depth independent. Despite this rather simple definition of boundary conditions, the difficultly lies in identifying particles that are at the top of the model and computing both the direction and magnitude of the applied force in such a way that a constant confinement is maintained while the geometry of the top surfaces changes during normal faulting.

The ‘shining-lamp’ algorithm (described in detail in Appendix A) is used to identify and apply appropriate forces to the uppermost particles so that the prescribed confining pressure is achieved in a robust way (Fig. 3c–e). Model confinement is applied in a step-wise manner, so that the applied force increases gradually (in order to prevent shock to the system). During initial model confinement the stress acting on the bottom wall is monitored and only when the basal stress level is stable (within a pre-defined tolerance level) is the model considered to be under static stress conditions. Since the model is confined by applying forces to its top surface only, the initial (pre-faulting) stress is not isotropic, meaning that the lateral stress (σ* _{h}*) is a fraction of the vertically applied stress (σ

*=*

_{v}*P*

_{conf}); since the model is loaded uniaxially the horizontal principal stresses have the same magnitude. The magnitude of this lateral pre-stress is given by the relation σ

*=*

_{h}*K*

_{0}σ

*, where*

_{v}*K*

_{0}is the so-called Earth pressure coefficient at rest, which is

*c.*0.42 in these DEM models. The magnitude of

*K*

_{0}is a function of the friction angle (ϕ) of the particle assemblage and can be estimated by the Járky expression,

*K*

_{0}=1−sinϕ, which gives a friction angle of 35°.

Bonds are installed after model confinement (except within the top layer) so that the pre-faulting bond stresses are zero, which, from a geological point of view, implies post-consolidation cementation. Post-confinement ‘cementation’ was necessary as application of confining pressure to bonded particle models can give rise to bond breakages, which would have led to different bulk mechanical properties of materials under different confining pressures (Schöpfer *et al.* 2009*b*).

### Normal faulting boundary conditions

Once a model is confined to the desired stress level the bottom platen is replaced by frictionless footwall and hanging-wall platens and a frictionless predefined ‘fault’ is generated in the ‘basement’ (Fig. 3e). The dip of this pre-defined fault is always 70°, because at lower fault dips antithetic faulting frequently develops in the overlying sequence (Horsfield 1977; Nollet *et al.* 2012). The lateral platens are replaced by frictionless infinite walls (which ensure that no particles ‘escape’ and also increase computing speed). Since the normal stiffness of the newly generated walls is identical to that of the initial walls, no disturbance of the force network occurs. During a model run the hanging-wall platen moves down with constant velocity to a final throw of between one and three times the strong layer thickness (final throw varied depending on the specific questions to be addressed with a certain set of models). Model states are saved in regular intervals of times the strong layer thickness for later analysis.

A recurring question and problem in DEM modelling is how fast a model can be deformed so that dynamic effects are negligible (so called quasi-static conditions). The optimal velocity, which ensures both quasi-static conditions and fast model runs, depends on model dimensions, properties and boundary conditions. For the normal faulting boundary conditions the most reliable diagnostic criterion is the stress acting on the hanging-wall platen (other measures, such as unbalanced forces or the bond breakage evolution are also useful). The basal stress evolution, for example, is illustrated in Figure 3f for four models run at different confining pressures. These curves illustrate that initially the normal stress acting on the footwall platen increases, whereas the stress on the hanging-wall platen decreases. This differential stress is obviously a result of the bulk strength of the overlying sequence. Once the strong beds are entirely offset (which depends on the nature of the sequence and confining pressure as described later) the footwall and hanging-wall platen normal stresses converge back to their initial (pre-faulting) stress level. By systematically varying the displacement rate of the pre-defined fault, we found the optimal hanging-wall platen velocity to be 0.1 m s^{−1}. This final velocity is preceded by an acceleration phase at the beginning of normal faulting, to avoid an abrupt change of boundary conditions that may lead to dynamically induced damage in DEM models.

### Model runs

A list of all models run for the present study is provided in Table 2. A constant strong layer thickness of 0.67 m is used in all models. Four models comprising *c.* 380 000 particles were run at confining pressures of 10, 20, 30 and 40 MPa (Fig. 3). The periodically layered sequence in which faulting is analysed comprises six strong and six weak layers of equal thickness. If we assume that the strong layers are reservoir rocks, then these models have a net:gross of 1:2. Various techniques for extracting fault zone geometry and quantitative characteristics and their confining pressure dependence are illustrated using these four models (Figs 4, 5, 6). Later in this paper we explore, using smaller models comprising *c.* 200 000 particles, the impact of net:gross and fault obliquity on fault zone structure (Figs 7, 8, 9, 10, 11).

## Fault zone visualization and quantification

### Passive markers and fault throw profiles

DEM model results are often represented by simply plotting particles, which are sometimes coloured according to their displacement. Plotting the locations of broken bonds is also common practice. For ease of comparison of our model fault zones with geological data we extract deformed passive markers from our models, effectively faulted bedding planes. Particle centres in the undeformed (pre-faulting) stage are triangulated using a Delaunay triangulation dividing the model volume into tetrahedra. A series of passive horizontal marker planes are defined through the centres of the strong layers and their intersection points with the edges of the tetrahedra are found. The passive markers are consequently discretised into three- and four-sided polygons. The equivalent polygons in the deformed state define the faulted passive markers and the geometrical distortions of the tetrahedra define the strain distribution over the marker. From these maps, faults are defined on the basis of a maximum shear strain cut-off and throughout this paper a tetrahedron is considered to be ‘faulted’ if it has a maximum shear strain (*s*_{1} and *s*_{3} are the maximum and minimum principal stretch, respectively). Polygons within ‘faulted’ tetrahedra are displayed in black, while unfaulted polygons are plotted according to their dip (e.g. Fig. 4).

On each passive marker the boundaries between faulted and unfaulted polygons define the fault/horizon intersection polygon or cut-off polygon. The difference in elevation between the upthrown and downthrown polygon cut-off defines the throw profile along the marker horizon. A detailed description of how fault throw profiles can be extracted from these fault polygons is provided in Appendix B. It is important to note that in the DEM models the individual fault strands are volumes of pervasive fracture and pronounced particle rearrangement (loosely speaking fault rock). The fault ‘heave’ polygons extracted with our method hence not only reflect the heave of the individual fault strands, but also include the width of the fault strands. Extraction of fault throws on which our analysis is based is not affected by the width of the fault strands.

### Quantitative measures of throw partitioning

For extraction of quantitative measures of throw partitioning fault throw profiles along marker horizons are resampled along regularly spaced sections (a spacing of 0.25 m is used throughout this paper, resulting in 39 sections per layer). Fault throw partitioning along each section is quantified using three measures (see also Ferrill *et al.* 2011): (a) the throw on the largest fault strand; (b) the cumulative throw on all smaller faults; and (c) the cumulative throw accommodated by drag and rotation. The sum of (a) and (b) is the discontinuous deformation in the model and the sum of all three is equal to the model throw. When normalized by the total fault zone (=model) throw, these three measures can be displayed in a ternary diagram (Fig. 5). For clarity, each point plotted in these ternary diagrams is the average of 39 readings for each layer. Data falling outside of the triangular plot are due to the presence of antithetic faults, which sometimes develop in the models as a result of the pre-defined fault at the base of the sequence and can give rise to displacements that are actually larger than that on the basement fault.

The normalized values for the three throw components are referred to as *fractional throws* and enable meaningful comparisons of fault zone structure at different stages in model. The *fractional cumulative fault throw*, for example, is the throw accommodated by all faults along a section divided by the total zone throw and is a useful measure of fault localization. In the discussions below we also refer to the *normalized model throw* (*T*_{norm}) which is the throw on the pre-defined fault divided by the thickness of the strong layers (0.67 m in all models shown in the present study).

The measures of throw partitioning in the following sections are presented as average values. While average values serve to illustrate particular trends in, and controls on, fault throw partitioning and localization, they are not necessarily the best measures for important application areas. One particular application of throw partitioning is in considering the across-fault connectivity of flow units for which minimum throw values are often of utmost importance; some of the conclusions of our analysis are recast in terms relevant to across-fault connectivity in the discussion section.

## Model results

### Variation in fault zone structure with displacement

The primary control on throw partitioning in the models is fault displacement. At low displacements, with normalized model throws of 0.25, bending of beds is the main contributor to the total throw accommodated across the fault zone and average fractional throws accommodated by drag and lens rotation are *c.* 0.8 (Fig. 6c). Fractional throws on the largest fault progressively increase and are generally >0.8 when the zone throw is three times the thickness of the strong layers (Fig. 6a). The average contribution of smaller faults to the overall throw is generally low and the highest average contribution is *c.* 15% at low confining pressures and normalized model throws of ≤1.0 (Fig. 6b). Once normalized fault throws are >1, the fault throw component is strongly concentrated onto the largest strand in the fault zone (Fig. 6). Therefore, while new minor faults can initiate to accommodate the imposed model throw at any time during fault development, all models record a progressive localization of throw onto a single through-going fault. The precise relationship between fault throw and throw partitioning varies between models depending on the range of parameters described in the following sections.

It is important to note that, while the average values shown in Figures 5 and 6, for example, give the impression that smaller synthetic faults are of little significance, these average values hide the fact that at particular locations over the fault surface they may be significant. For example, Figure 4 shows that at a normalized throw of 1.25 the fractional throw on the largest fault strand is locally, owing to the presence of a fault-bound lens (at the right model edge), <0.5.

### Influence of confining pressure on fault zone structure

The impact of confining pressure on change in fault throw partitioning during fault growth is illustrated in Figures 5 and 6. In general the models show that localization is most efficient under higher confining pressures with fractional throw values on the largest fault strand of >0.8 achieved when model throws are smaller than the thickness of the strong layers in the models with confining pressures ≥20 MPa. At the lowest confining pressure considered (10 MPa), significant localization only occurs when the throw approaches the thickness of the strong layers, i.e. at a normalized throw of 1. The spatial variation in throw partitioning, reflected by the standard deviation, is also strongly dependent on confining pressure. At higher confining pressures the variability in throw partitioning decreases rapidly at normalized model throws of 1 while at the lower confining pressure of 10 MPa the variability in throw partitioning does not decrease with increasing throw (Fig. 6). This feature of the data reflects the fact that, at lower confining pressure, initially formed fault-bounded blocks are not comminuted with increasing displacement, a feature discussed again later.

### Influence of net:gross on fault zone structure

The thicknesses of the strong and weak layers in the DEM models presented so far are identical. In this section the influence of net:gross on fault zone structure is discussed. The periodically layered sequence in these models is always 4 m thick. The thickness of the strong beds is always 0.67 m, so that by varying their number the following four different net:gross sequences are obtained (Fig. 7a; Table 2): 1:3, 1:2, 2:3 and 5:6. Confining pressure in all models is 20 MPa. The geometry of the DEM model fault zones is quite variable under these pressure conditions (Fig. 6) so that four realizations are run for each net:gross value. Models are saved in 0.25 normalized model throw intervals and final model throw is two strong layer thicknesses.

Cross-sections through a selection of these models in which particles are coloured by their vertical displacement illustrate that the best-localized faults appear to develop in the net:gross=2:3 model (Fig. 7a). This conclusion is also illustrated by the passive markers and associated throw profiles extracted from these four models (Fig. 8): In the net:gross=1:3 model a wide fault zone containing multiple fault-bound lenses develops. In outcrop such a fault zone may well be categorized as a ‘forced fold’ (e.g. Cosgrove & Ameen 1999). In fact, at a normalized model throw of 0.5 (as shown in Fig. 8), *c.* 50% of the total throw is accommodated by drag and lens rotation in the net:gross=1:3 model. However, with increasing net:gross the model fault zones become better localized and the proportion of throw accommodated by lens rotation decreases (Fig. 8). However the highest net:gross model (5:6) exhibits again pronounced rotations and a geometry similar to a kink-band.

The trends in the different net:gross models are quantitatively illustrated in fractional cumulative fault throw graphs in Figure 9, where each data point is the average of 39 sections per layer and four realizations. The curve that most rapidly converges towards a value of 1 (meaning that the entire model throw is accommodated by faulting) is for the bottom layer in the net:gross=5:6 models. The main reason for this trend is that the weak layers in the net:gross=5:6 models are very thin (see Fig. 7a), so that the bottommost strong layer is offset via a discrete fault at a very early stage. In that respect, the fractional throw curve for the bottom layer in the net:gross=5:6 is not representative of the general trend seen in the different net:gross models. Indeed, all other fractional throw curves support the aforementioned observation that the best localized fault zones develop in the net:gross=2:3 models.

### Oblique-slip faulting models

In this section, the impact of fault obliquity on fault zone structure is investigated using net:gross=1:2 models, each of which comprises three strong and three weak layers of equal thickness (Fig. 7b; Table 2). All other dimensions are identical to those used in the large models (Fig. 3).

Oblique-slip faulting is modelled by changing the trend of the pre-defined fault over the entire model width (Fig. 7b) so that obliquity ranges from 0° (=dip-slip) to 40° with 10° increments. The obliquity (ω) is here defined as the angle between the dip- and slip-direction of the pre-defined fault (see insets in Fig. 7b). The apparent dip, measured in a section containing the slip-vector, of the pre-defined fault is 70° in all models. The refinement regions (model volumes that contain smaller particles for better resolution) are adjusted accordingly (Fig. 7b). Confining pressure ranges from 10 to 40 MPa with 10 MPa increments. Models are saved in 0.25 normalized model throw intervals and final model throw is one layer thickness. Fault analysis is again made on the basis of passive markers (=horizons) through the centres of the strong layers (Fig. 10).

The fault zone dependencies (e.g. fractional throws) of the oblique models on confining pressure are, as expected, similar to those observed in the previously described, larger models (compare Figs 11a and 6a). Perhaps the most interesting feature of the different obliquity models is, however, the dependence of fault zone structure on obliquity, particularly the fault strand orientations (Fig. 10). The trends of individual fault strands are estimated by fitting straight lines to the fault median lines of strands with a length >0.5 m (see Fig. B1c). The map-view orientations of fault strands are measured as angles relative to the slip-normal trend; positive values are anticlockwise. These fault strand trend data are shown for the topmost strong layer and for the full range of confining pressures in Figure 11b. Despite the scatter, the data suggest that the mean (and median) fault strand trends increase with increasing fault obliquity by a factor of 0.5. In other words, a fault obliquity of ω leads to an average fault strand trend of *c.* 0.5ω. This relation can be rationalized by infinitesimal strain theory for zones of transtension in which faults nucleate normal to the greatest instantaneous stretching axis (e.g. McCoss 1986). A graphical construction by McCoss (1986), which is based on the model by Sanderson & Marchini (1984), illustrates the theoretical fault strand orientations as a function of obliquity (Fig. 10). The fact that a kinematic model for transtension offers a rationale for the observed fault strand orientations suggests that, in the DEM models, fault segments nucleate within a pre-cursory monocline, as already demonstrated in earlier 2D DEM models (Schöpfer *et al.* 2006).

## Summary and discussion of model fault zones

### Throw and confining pressure dependencies

The DEM model fault zones become better localized with increasing throw and confinement, trends which are consistent with earlier DEM normal faulting models (Schöpfer *et al.* 2007*a*, *b*, 2016) and physical experiments comprising cohesive powder embedded in cohesionless sand (Kettermann & Urai 2015). The DEM model fault zones initially comprise segmented fault arrays (Fig. 4). Model rock volumes enclosed by faults and fault plane irregularities (asperities) become progressively comminuted with increasing throw, giving rise to more planar faults. All models become better localized with increasing throw (Figs 5, 6 & 9) with strongly localized faults developing when fault throw is greater than one or two bed thicknesses. The modelled fault geometries are limited in scale from a few particle diameters to the model size, and therefore segmentation and bifurcation occur within a limited scale range. It is likely that faults in larger models would have greater scale geometrical complexity and would therefore require increased throws to achieve localization.

We consider two plausible end-member explanations for increased fault localization with increasing confining pressure: (a) fault zones are initially (i.e. at low throw) more segmented at lower confining pressures than at high confining pressures so that localization does not occur until higher throws are attained; and (b) initial fault segmentation is independent of confining pressure, but existing fault segments are more readily abandoned and/or fault-bound volumes are more easily comminuted at higher confining pressure. The DEM model results favour scenario (b), because various fractional throw measures at low model throws appear to be confining pressure independent (Figs 5 & 6).

It appears therefore that the DEM models support a fault zone evolution model in which faults are segmented over a wide range of confining pressures and that faster localization with increasing confining pressure is due to abandonment of earlier segments at lower throws and to more efficient communition of fault-bound volumes. The question is why at low confining pressure, fault segments which are perhaps not optimally orientated to accumulate large displacements are not abandoned and the rock volumes they bound remain to a greater extent intact. There are two analogies that may offer a mechanical explanation for this behaviour. (a) It is a well-known fact that in granular shear the deformation mechanism is shear zone normal stress dependent, so that at low normal stress grains remain largely intact and exhibit only surface abrasion, whereas at higher normal stress grain splitting occurs (e.g. Mair & Abe 2011). Therefore we would expect that fault-bound volumes with a given size will break up more readily at higher confining pressures. (b) Mechanical analysis (based on the Coulomb criterion) of slip surfaces with a saw-tooth profile suggests that at a certain normal stress the shear stress required for reactivation is greater than for failure of the intact material, resulting in a bi-linear criterion for the shear strength of rock fractures (Patton 1966). We would therefore expect that fault surface irregularities (in the fault profile plane) are more readily removed with increasing confining pressure. We suggest, by analogy, that the low fractional throws accommodated by faulting persist at low confining pressures (Fig. 6) because the strength of the unfaulted volumes is significantly greater than the differential stress required to reactivate existing fault surfaces, even if they are poorly oriented for reactivation. At high confining pressures existing fault surfaces are subject to high normal stresses so that it is more likely that they will become by-passed. The two aforementioned mechanisms probably both occur in the DEM models, where scenario (a) occurs after both tip-line and asperity bifurcation and scenario (b) occurs during asperity bifurcation (Childs *et al.* 1996).

### Impact of fault obliquity

The fault strand orientations extracted from DEM models deformed by an oblique-slip fault at the model base suggest that faults nucleate normal to the maximum instantaneous stretching axis (Figs 10 & 11b). This observation is consistent with results from physical experiments conducted with sand and silicone (Richard 1991) and clay (Schlische *et al.* 2002). Both of these studies are worthwhile discussing in light of the DEM model presented here.

Richard (1991) used an experimental apparatus for oblique-slip faulting consisting of a large table divided in two halves. In the oblique-slip normal faulting experiments one half was moved along a 45° dipping basement fault with a strike-slip to dip-slip ratio of 1. The pitch of the slip vector upon the fault plane is therefore 45° and the obliquity ω, as defined in Figures 7b & 10, is atan(2^{0.5})=55°. If faults in the overburden were to develop normal to the maximum instantaneous stretching direction, then the expected angle between the basement fault trend and the faults is 27°. In the experiments with a silicone layer at the base of the sand pack, fault zones develop in which the maximum fault strand trend of normal faults is indeed about 27° (fault strands with a smaller angle to the basement fault have a significant strike-slip component). An important result is also that wide fault zones (width increases with increasing silicone layer thickness) comprising en échelon normal fault segments only develop if a weak substratum is present (Richard 1991, fig. 4b–d). An experiment conducted without a weak substratum (Richard 1991, fig. 4a) led to the formation of a graben, bound by normal faults parallel to the basement fault, with a strike-slip fault in its centre. Effectively, oblique motion is partitioned into slip on faults with different senses of motion (Bowman *et al.* 2003). Weak layers (which are present in our DEM models) therefore promote the formation of wide transtensional zones.

Schlische *et al.* (2002) performed oblique-slip normal faulting experiments (obliquity ω=45°) on thin and thick clay layers and quantified the fault trend frequencies in the ‘overburden’. Both in the thin and thick clay layers the most frequent fault strand trends are halfway between the so-called master-fault and displacement-normal trend (Schlische *et al.* 2002, fig. 4). Again, this result is consistent with infinitesimal strain theory for transtensional faults. A notable difference between these clay models and our DEM models, however, is relay ramps, which are frequent and temporarily persistent in the former and rare in the latter (see Fig. 10). There are several reasons for this difference; of the two most likely, the first is that the along-strike dimension of the models is not large enough so that well-developed relay ramps can develop (i.e. the probability of a large relay to form somewhere along strike is too low). This explanation could in theory be tested by running larger models, which, however, would require impractical model run times. The second reason is partly one of definition. The DEM models do display rotation of ramps but, because large permanent strains can only form by fracture (i.e. bond breakage), these ramps are entirely fault bounded at low throws and are therefore fault-bound lenses rather than relay ramps. If the maximum shear strain cut-off applied to define a fault (see Fig. B1a) is increased, then many structures will display relay geometries. Accommodation of large permanent strains without bond breakage would require a rate-dependent bond model that mimics micromechanical processes, such as precipitation-dissolution creep. This line of research would be, in our opinion, demanding but potentially extremely useful.

### Influence of net:gross

The DEM models illustrate that, for the mechanical properties and confining pressure used in the models presented, the most rapidly localizing faults develop in the net:gross=2:3 models (Figs 8 & 9). Our preferred mechanical explanation for the observed net:gross dependence is that at low net:gross fracture-bound blocks can easily rotate within a ‘ductile’ matrix, which promotes continuous throw accommodation (drag and lens rotation). As net:gross increases, the adjacent layers begin to impinge on one another so that they are no longer bounded by a thick weak ‘matrix’ and rotation becomes more subdued, resulting in faster localization. At the highest net:gross, however, layer-parallel shear across the thin weak layers becomes suppressed and brittle layers impinge on each other early during fault zone growth to the extent that the assemblage of closely spaced beds operates as a single, thick mechanical unit which rotates as a whole, therefore suppressing localization. The ‘optimal’ net:gross value at which the best localized faults develop (which is 2:3 in our models) will most likely depend on other factors such as confining pressure and strength contrast. However, we would expect that the general trend (pronounced rotations at low and high net:gross) would also emerge under different confining pressure and/or strength conditions because fault zone structure in layered sequences depends mainly on relative properties (e.g. layer strength relative to confining pressure; Schöpfer *et al.* 2007*b*).

The DEM models are consistent with several field studies. For example, Ferrill & Morris (2008) demonstrate that the structure of normal faults in Cretaceous carbonates in the Balcones fault system (Texas) greatly depends on the mechanical stratigraphy. As in our DEM models, a high proportion of ‘incompetent’ material promotes fault-related folding. On the other hand, fault zones in massive (poorly bedded) sequences appear to exhibit complex internal structure, which may reflect that rotation of fault-bound volumes is, in the absence of bedding-parallel slip, inhibited (e.g. Bonson *et al.* 2007).

### Implications for across-fault reservoir connectivity

In this paper we have primarily been concerned with average values of throw-partitioning parameters as measures of fault zone geometrical evolution (Figs 6, 9 & 11a). However, one of the main potential applications of this modelling is in evaluating the likelihood of connectivity of flow units across fault zones, which requires consideration of the extreme rather than average values. The key measure is the minimum throw on the largest continuous fault strand. If this minimum throw is larger than the flow unit thickness, then it is disconnected across the modelled fault zone length (if across-fault flow is considered in 2D only, then the throw on the largest fault strand relative to the thickness of the flow unit determines juxtaposition sealing; Fig. 12a). The minimum throw on the largest continuous fault strand is determined in an iterative fashion; a sample result is illustrated by means of a strike-projection and a throw profile in Figure 12b. For the models presented here we consider the strong layers to be the flow units and express our results in terms of the model throw at which individual flow units are offset along the length of the fault divided by the flow unit thickness, i.e. the model throw at which the minimum throw exceeds the thickness of the strong layers; we refer to this measure as the *normalized critical model throw*. A normalized critical model throw of 2 means that the model throw must be twice the thickness of the flow unit before it is entirely offset along the length of the fault and a value of 1 means that all throw is concentrated onto a single continuous fault strand (Fig. 12c). The majority of natural faults and also our models exhibit, to a greater or lesser degree, displacement partitioning onto multiple strands and bed rotations. Consequently the throw at which juxtaposition sealing occurs is typically greater than the thickness of the flow unit so that the normalized critical model throw is greater than one (Fig. 12c).

Plots of normalized critical model throw v. both confining pressure and net:gross (Fig. 13) show similar relationships to those established from average values of fractional throw (Figs 6 & 9). The normalized critical model throw decreases with increasing confining pressure and approaches a value of 1 (the throw at which juxtaposition occurs is only slightly greater than the thickness of the strong layers at 40 MPa confinement; Fig. 13a). At a confining pressure of 10 MPa the normalized critical model throw is, on average, >2; in this case the throw must be greater than twice the flow unit thickness so that no juxtapositions occurs over the entire fault zone width, reflecting that a significant proportion of the total throw is accommodated by smaller faults and bed rotations. Similar to the conclusion derived from fractional throw measures (Figs 6 & 11a), the variation in normalized critical model throw decreases as the confining pressure increases so that the model fault is not only more localized but also more uniform.

The normalized critical model throws show a general decrease with increase in the proportion of strong layers within the model and values of nearly 1 are achieved at net:gross of 5:6 (Fig. 13b). It appears that the variability in normalized critical model throw does not change with net:gross and there is significant variability between different horizons within the same model. The average values of normalized critical model throw are lowest at a net:gross value of 2:3 in line with the qualitative impression of fault zone width and localization provided by Figure 8 and the fractional cumulative throw curves shown in Figure 9.

The present analysis is based on detecting the minimum throw on the largest continuous fault strand that spans the entire model width (Fig. 12b). For complex fault zones we expect that the minimum throw on the largest continuous fault will decrease with increasing fault zone length, so that the total throw at which juxtaposition sealing occurs over the entire fault zone length will increase with increasing model, or sampling, length. In other words, the probability of juxtaposition sealing along the entire fault zone decreases with increasing sampling length because the probability of encountering a structure that promotes across-fault fluid flow (e.g. a relay ramp) increases with increasing sampling length. A study on the scale-dependency of juxtaposition sealing, perhaps based on the algorithms presented here, would be a useful line of future research.

## Conclusions

The model fault zones suggest that the fraction of strong material in a layered sequence (net:gross) has a profound impact on throw partitioning and localization:

models with net:gross of 1:3, 1:2, 2:3 and 5:6 deformed at a confining pressure of 20 MPa suggest that the most rapidly localizing fault zones develop in the net:gross=2:3 models;

in low net:gross sequences significant proportions of the total zone throw are accommodated by rotations, not unlike forced folding;

in high net:gross sequences bedding parallel slip is suppressed and strong layers impinge on each other early during fault zone growth, leading to wide and poorly localized fault zones.

Model fault zones developing above pre-defined oblique-slip normal faults suggest the following:

the average fault strand orientations are consistent with predictions based on infinitesimal strain theory for transtensional zones;

the fact that a kinematic model for transtension offers a rationale for the observed fault strand orientations suggests that fault segments nucleate within a pre-cursory monocline.

Analysis of DEM models deformed at four different confining pressures and juxtaposition sealing analysis support the following conclusions:

fault zones become better localized with increasing throw and confinement;

at low confining pressure the fault zone throw must exceed, on average, twice the bed thickness so that across-fault flow is not possible (assuming juxtaposition sealing only);

the critical throws beyond which no across fault flow is possible decrease with increasing confining pressure;

the ranges of critical throws are considerable, because critical throws are controlled by the local, rather than the average, fault zone structure.

Many of the aforementioned conclusions concerning the growth of fault zones are consistent with field observations and results from physical experiments. The DEM may therefore provide a basis for evaluating the likely complexity of fault zone structure and associated sequence juxtapositions in a given sequence faulted under a certain confining pressure.

## Acknowledgments

The DEM modelling presented here was performed during the course of the multi-company project ‘QUAFF’, brokered by ITF (Aberdeen) and funded by Anadarko, BG International, BP Exploration, ConocoPhillips (U.K.), Eni, ExxonMobil, Marathon Oil Corporation, Statoil, Total E&P U.K. and Woodside Energy. C. Childs is funded by Tullow Oil. Some methods were developed during the course of the earlier multi-company project ‘FaultZone’, funded by ExxonMobil and Shell. We also thank the staff of Itasca Consulting Group for their support, in particular Dave Potyondy for fruitful discussions and for providing his ‘shining-lamp’ algorithm, and Sacha Emam and Peter Cundall for clarifying the concept and usage of distortional periodic space. Thorough reviews by Michael Kettermann and Kevin Smart are gratefully acknowledged.

## Appendix A: Pressure boundary condition via the ‘shining-lamp’ algorithm

A recurring problem in 3D DEM models comprising spherical particles is the application of a constant pressure boundary condition without the use of rigid, servo-controlled platens. In the DEM, forces can be applied to particles, but a robust pressure boundary condition that mimics, for example, the overburden pressure exerted on a rock volume at depth requires two important features: (a) all particles that are at the surface of the model must be detected, both at initiation of the model and as the surface deforms and becomes non-planar; and (b) both the magnitude and direction of the force applied to each surface particle must be determined so that the force vectors are normal to the surface. The so-called ‘shining-lamp’ algorithm, developed by D. Potyondy (Itasca Consulting Group, unpublished memorandum ICG7233-L, 2012) provides these features.

The pressure-application procedure is referred to as the shining-lamp algorithm because it defines the surface particles as those that receive direct illumination from light sources that shine toward the body. Each surface particle is assigned an externally applied force (**F**) that is related to the specified pressure (*P*_{conf}) and the amount of light that it receives (which is related to its illuminated surface area).

The principles of the algorithm are illustrated by means of a simple assemblage of discs that receive light laterally from plus and minus infinity (with inward-directed unit-normal vectors **p** and **p′**, respectively; Fig. A1a). For each particle the applied force (acting at its centre) is the sum of the force contributions from each light source **F**=*P*_{conf} *A***p**, where *A* is the projected illuminated area. Since the sum of these projected areas is identical on both sides of the assemblage, the applied forces balance.

The preceding description assumes that the illuminated surfaces can be computed exactly. The shining-lamp algorithm discretizes the light sources into rays emanating from regular grids (square grids are used in the present 3D models) and thereby provides an approximation of each illuminated surface (Fig. A1b); the approximation converges to the exact surfaces in the limit as the grid size approaches zero. The algorithm is very efficient and only requires one pass through the list of particles (in the normal faulting models only particles comprising the top layer are checked; Fig. 3). For each particle the light-ray intersections are computed (sphere-line intersections in 3D) and the positions of these intersections are stored at each corresponding grid-point. If an intersection is closer to the light source than an earlier computed intersection, then the current particle becomes a surface particle (the addresses to surface particles are also stored at each corresponding grid-point). A consequence of this discretization is that the applied particle forces are proportional to the number of light–ray particle intersections (Fig. A1b), so that **F**=*P*_{conf} *n _{p}a*

**p**, where

*n*is the number of intersections and

_{p}*a*the grid-spacing in 2D and the grid-spacing squared in 3D (in the case of square grids). Note that the discrete area (

*a*) is illustrated as line segments, centred on the intersection points, in 2D (Fig. A1b) and as squares in 3D (black squares in Fig. 2c; cyan and magenta squares in Fig. 3c). The computation of ray–particle intersections and hence applied particle forces is performed in regular intervals (e.g. every 100 time-steps) to account for changes of the model's surface geometry.

In the present study the shining-lamp algorithm is implemented in two model environments, (a) the triaxial extension/compression tests (Fig. 2c) and (b) the normal faulting models (Fig. 3b–e). In the triaxial models light rays are emanating from square grids with inward-directed normal vectors parallel to the global *x-* and *z*-coordinates (four grids in total; grid spacing is 0.5*r*_{min}). Depending on its position, a particle may therefore receive light from more than one light source; hence forces are applied both in the *x-* and *z*-directions, so that the force vectors are normal to the surface of the cylindrical sample, as they should be.

In the normal faulting models a constant overburden pressure is applied via the shining lamp algorithm. The efficiency of this procedure obviously depends on how well the model surface is illuminated by light, hence two sources (i.e. discrete grids with a spacing of 0.2*r*_{max}) with rays inclined at 45° and 135° (counter-clockwise from the positive *x*-direction) are used (cyan and magenta, respectively; Fig. 3b). Computationally, however, it is more efficient to use light rays that are parallel to the global coordinate system. Hence a coordinate transformation is performed when the light–ray particle intersections for the top layer are computed: each particle position is rotated clockwise by 45° around the origin (located in the model centre) and then intersected with light rays emanating from the positive *x*- and *y*-directions (the latter points upwards). Obviously, the applied forces computed in the transformed coordinate system have to be adjusted to the global coordinate system. Since only particles lying at the model surface should receive light, care has to be taken in the positioning of the grids so that the lateral surfaces of the model are not illuminated. Consequently, the grid from which rays are emanating at 45° (cyan in Fig. 3b, c) must be translated according to the throw of the moving hanging-wall platen.

### Appendix B: Extraction of fault throw profiles from passive markers

As described in the main text, horizon maps through the centres of each strong layer are extracted by intersecting passive marker planes in the pre-faulting stage with the edges of a 3D Delaunay mesh obtained from triangulated particle centres. Horizon maps are hence represented by polygons (triangles and four-sided polygons) for which the 3D coordinates of the vertices and the strain tensor are known. Measures derived from the strain tensor, such as maximum shear strain, can then be used to visualize faults (e.g. Fig. 4). In this Appendix a new technique that allows extraction of displacements from faulted horizon maps is described (Fig. B1).

The first step is to generate fault outlines, commonly referred to as fault polygons. This step is achieved by calculating the horizontal polygon centres and by drawing a contour using one contour level (which is equal to the aforementioned maximum shear strain cut-off). These contours provide the map view coordinates of the fault polygons (Fig. B1a). The vertical coordinates are readily computed by interpolation. This first step provides 3D outlines of the faults, viz. the footwall and hanging-wall cut-offs. The original contour may be noisy and hence small isolated contours and small peaks parallel to the trend of the pre-defined fault are removed.

The next step is to classify the contours. Three types can be distinguished (Fig. B1b): (a) isolated contours; (b) contours within contours (which must enclose fault-bound volumes); and (c) contours that straddle the model boundaries. Then local maxima, parallel to the trend of the pre-defined fault, are identified and classified. Three types of local maxima are distinguished (Fig. B1b): (a) fault tips; (b) branch points; and (c) model edges. Tips and branches are identified on the basis of simple geometric rules that assume that the trends of individual fault strands do not vary too dramatically (which is typically the case in our model faults). For example, a local contour maximum pointing to the positive trend-parallel direction (abscissa in Fig. B1) must be a fault tip if a point shifted slightly from this position to the positive direction falls outside the fault polygon. In contrast, if this shifted point would fall inside the fault polygon then the local contour maximum must be a branch point. Note that ‘contours within contours’, which correspond to fault-bound volumes, may also have both tips and branches.

The contours are then intersected by vertical planes, or sections, whose normals are parallel to the trend of the pre-defined fault. First, sections that pass through the previously identified tips and branches are defined (a selection of these traces is shown in Fig. B1b). This is an important step, because the number of fault strands is constant in-between tips/branches (Fig. B1d). Second, additional regularly spaced sections are defined to enhance the resolution of the final displacement profiles. Third, the section positions are checked and additional sections are defined, if necessary, to ensure that there is at least one section between each tip/branch location (otherwise no throw profile can be constructed if, for example, there are no sections in-between two tips).

Then the intersection points (3D coordinates) of the contours with the previously defined sections are computed and the median points (halfway between the footwall and hanging-wall cut-offs) are calculated (Fig. B1c). Fault strand throws and widths (vertical and horizontal differences between the footwall and hanging-wall cut-offs, respectively) are also computed. Fault strand throws are positive if the fault is synthetic to the pre-defined fault, whereas antithetic faults have a negative throw.

The final step is to extract fault segments (Fig. B1c). As stated earlier, in map view a fault zone can be ‘sliced’ by lines normal to its ‘average’ strike into sub-areas of constant fault strand number (Fig. B1d); the slicing locations are either tips or branches. For each of these sub-areas fault segments are extracted. The two lateral ends of each segment can be classified as follows: (a) fault tip; (b) branch point (merging or splitting branch points are distinguished); (c) model edge; and (d) no fault end. For each fault segment the 3D coordinates of the median line, the associated throws and widths, and segment end identifiers are used for later plotting and analysis.

- © 2017 The Author(s). Published by The Geological Society of London. All rights reserved