cdem: A macOS program for discrete element modeling of tectonic structures

Since the pioneering work of Cundall and Strack (1979), the discrete element method (DEM) has been used to model small- and large-scale tectonic structures. Some examples include shear zones (Morgan and Boettcher, 1999; Mair and Abe, 2008), fractures (Schöpfer et al., 2011; Virgo et al., 2014), fault zones (Egholm et al., 2008), fault-related folds (Finch et al., 2003, 2004; Benesh et al., 2007; Hughes et al., 2014), normal faults (Schöpfer et al., 2007; Finch and Gawthorpe, 2017), fold-and-thrust belts (FTBs; Hardy et al., 2009; Dean et al., 2013), strike-slip faults (Liu and Konietzky, 2018), and salt structures (Hardy, 2018; Pichel et al., 2019).

The DEM simulates the rock as an assembly of elements, which interact under gravity and contact forces. This mesh-free technique is ideal for modeling large deformation and discontinuities (e.g., faults), which are naturally emerging features in the model (Gray et al., 2014). Computational limitations on element size and model resolution and small time steps required to maintain system stability, while important, are no longer a particular concern due to recent advances in computer power and parallelization. The main challenge in DEM is perhaps the need to constrain independently the bulk properties of the assembly, given that these cannot be defined a priori from the element properties. Several studies have documented best practices for assembly calibration (e.g., Schöpfer et al., 2009, 2013), and some have tried to overcome this limitation (Egholm, 2007).

Although the DEM is a well-established technique within the modeling community, it is not broadly used by geologists in academia and industry. There can be several reasons for this, including the numerical complexities of the technique; a high-entry computational level; poorly designed interfaces that don’t integrate pre-processing, processing, and visualization; engineering codes that don’t yield realistic geological results; and prohibitive costs.

In this paper, we introduce the macOS program cdem to model tectonic structures in two dimensions (2-D) using the DEM. cdem is based on the DEM code cdem2D by Stuart Hardy, which has been applied to a variety of problems including fault-propagation folding (Hardy and Finch, 2007), orogenic wedge growth (Hardy et al., 2009), caldera collapse (Hardy, 2008), dike intrusion on Mars (Hardy, 2016), salt flow (Hardy, 2018), and Gilbert deltas (Hardy, 2019a). cdem integrates pre-processing, processing, and visualization in a friendly document-based architecture. Several documents corresponding to different simulations can be run and compared simultaneously. cdem is fast because it is designed to target all available CPU cores. Most importantly, cdem was made by geologists for geologists, and it produces realistic tectonic structures.

We illustrate the application of the program to assembly calibration (biaxial and collapse tests), extensional and contractional structures, and syn-tectonic (growth) strata. cdem is less powerful than its precursor cdem2D, but it can import and visualize cdem2D simulations, making these two programs a robust suite for modeling complex structures at fine resolution. We illustrate this functionality as well. Overall, cdem provides geologists with an intuitive and powerful framework for mechanical modeling of tectonic structures, which can help validate geological interpretations.

The DEM in cdem is a variant of the lattice solid model (LSM) of Mora and Place (1993, 1994) and Place et al. (2002). In 2-D, a rock is treated as an assembly of circular elements, which interact as if connected by breakable elastic springs (bonds) and undergo motion relative to one another. The elements interact through a “repulsive-attractive” force (Mora and Place, 1993) in which the resultant (normal) force, F_n, is given by:

Here, K_n is the elastic constant (normal spring stiffness) of the bond, r is the current distance between the elements, R is the equilibrium distance between the elements (sum of element radii), and r₀ is the breaking strain. Elements are bonded until the separation between them (r − R) exceeds the breaking strain times the equilibrium distance (r₀ × R), at which time the bond breaks. The force acting on a bond at this stage represents the force necessary for a bond to fail or yield or, alternatively, can be cast as the stress acting on an element’s bond at failure. After this breaking threshold, the element pair experiences no further attractive force, and the bond is irreversibly broken. However, if the two elements return to a compressive contact (r < R), a repulsive force still acts between them. This initially bonded material is the simplest model in cdem.

A second, more realistic model in cdem is the frictional-cohesive model. In this model, all elastic bonds are initially broken, and in addition to treating the normal force (F_n) between overlapping elements (Equation 1), we also calculate the tangential (shear) force, F_s, because of displacement (X_s) perpendicular to the vector connecting the element centroids. This frictional force acts in a direction opposite to that of the relative tangential velocity and is modeled as a threshold-limited elastic spring with a cohesive force term (C₀) in parallel with that used to calculate the normal force (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:

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-element 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 an element is thus obtained by summing the normal and tangential forces on each contact or bond that links a specific element to its neighbors, calculated by:

where f_{i, j} is the elastic force (normal and shear) experienced by element i from its neighboring element j. However, we include a viscous damping term (proportional to element velocity) that acts to dampen reflected waves from the rigid 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 element is given by:

where υ represents the dynamic viscosity (damping constant) and x is the velocity of the element. At each discrete time step, the elements 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). Table 1 lists the default values of the spring constants, time step, and damping constant in cdem. Note that the normal and shear spring stiffnesses (K_n and K_s) have the same value, the time step (dt) is proportional to the square root of the minimum mass of the elements, and the damping constant (υ) is proportional to dt. Also, because the element mass is proportional to its cross-sectional area, dt and υ are proportional to the minimum radius of the elements.

The core of cdem is a stripped-down version of cdem2D. We translated the original cdem2D code from C to Swift, designed the interface following macOS guidelines, and parallelized the program via Grand Central Dispatch (GCD) to use all available CPU cores. cdem is a document-based application. Several documents, each one encapsulating a simulation, can be run and compared at the same time.

A document consists of three views: the Set-up view, the Summary view, and the Results view (Fig. 1; see Supplemental Material¹). The Set-up view is the pre-processing module of the program. In the first part of the Set-up view, the user can enter the assembly width and height (1 in Fig. 1A), the unit scaling (2), the total runtime (3), the display interval (4), the equilibration time (5), gravity (6), and density (7). The assembly in cdem is rectangular, and it can be between 10 and 100 units in width and between 2.5 and 7.5 units in height. The elements have four different radii of 0.05, 0.075, 0.1, and 0.125 units. The unit scaling defines the length scale of a model unit, and it can be between 100 and 1000 m. For the default width, height, and unit scaling (Fig. 1A), the assembly size is 10 × 5 km and the element radii are 50, 75, 100 and 125 m. Note that every run includes an equilibration phase (5 in Fig. 1A) before any boundary displacement begins. This phase is necessary because “a priori” the assembly does not know what its properties will be, and therefore it needs to “feel” the new conditions and equilibrate to a static state that is appropriate to the run.

The second part of the Set-up view deals with the material model. Here the user can choose between the granular frictional (Equations 1 and 2) or the initially bonded material (Equation 1) and assign parameter values (8 in Fig. 1A). For the granular frictional material, there is one checkbox that, when checked, indicates the cohesion is lost after the first slip event at any element pair contact; when unchecked, the cohesion is lost after the two elements physically separate. In general, first slip (the checkbox on) works well in contraction, while first separation (the checkbox off) works better in extension.

The third part of the Set-up view deals with the heterogeneity of the assembly. Here the user can set the layer contacts to be frictionless (if granular material) or unbonded (if bonded material) to simulate flexural slip (9 in Fig. 1A). The next option (10) makes the friction (if granular material) or breaking strain (if bonded material) of a lower or upper portion (0%–100%) of the assembly equal to a fraction (0%–100%) of the specified values in the material model (8). Entering 0% results in a frictionless or unbonded interval. In addition, it is possible in this option (10) to set the density of this interval to be a fraction (50%−100%) of the density specified in the first part of the Set-up view (7). The “weaker” and/or “lighter” interval is highlighted as semi-transparent in the model. The number of pre-growth layers in cdem is 24. The lowermost option (11) allows setting the friction coefficient of layers 5–8, 13–16, and 21–24 to the value entered in the related text field. This is useful to simulate a sequence of competent (e.g., friction coefficient = 0.25) and incompetent (e.g., friction coefficient = 0.1) layers of equal thickness. Note that this option and the previous one (11 and 10, respectively) are mutually exclusive, and the lowermost option (11) is not available for the bonded material.

The fourth part of the Set-up view deals with the boundary conditions and faulting modes. Here the user can choose between a planar (or non-planar by ticking the Listric CH or CD checkbox) normal or reverse fault, two vertical faults in extension (caldera) or contraction (piston), a detachment fault under contraction, a repose-collapse test, and a biaxial test (12 in Fig. 1A). Note that for the planar (or non-planar) fault, it is possible to include tectonic inversion at a specified stage of the run.

Finally, the fifth part of the Set-up view allows the user to introduce overburden pressure (13 in Fig. 1A), erosion of elements above a defined base level (14), or syn-tectonic sedimentation of elements below a static or rising sea level (15).

Once the model parameters and boundary conditions are defined, the user can initialize the model by pressing the Initialize button (Fig. 1A). This writes the model parameters and boundary conditions to the Summary view (Fig. 1B). The text in this view is important. Here, key information—such as the number of elements, the model dimensions, element radii and masses, physical constants, calculation parameters and constraints, runtime details, boundary displacements, display increments, fault configuration, material properties, and growth sedimentation properties—are all reported. Reading this text is an effective way to check any possible inconsistencies before the run and provides a definite record of the simulation.

Pressing the Run button in the Summary view (Fig. 1B) runs the model in the background and displays the evolution of the assembly in the Results view (Fig. 1C). A counter in the lower left corner of the Results view shows the current display increment and the net displacement (for faulting), strain (for the biaxial test), or time (for the repose-collapse test). This is the processing part of the program, but with the additional advantage that one can visualize the output in real time and, if desired terminate the simulation by pressing the Stop button (Fig. 1C).

Visualizing and analyzing the results after the run is completed (post-processing) is easy. Any increment can be chosen using the slider at the base of the Results view or the text field next to it (Fig. 2A). Pressing the play-pause button next to the slider plays and pauses the simulation. Pressing the rewind button brings the simulation to its first increment. The assembly can be zoomed, and if desired, a scale bar can be displayed (Fig. 2A). Displacement vectors and strain and stress axes can also be displayed, and the elements can be colored by displacement, strain, or stress (Figs. 2B–2D). Incremental displacement or strain is calculated between the current increment and five increments before, while total displacement or strain is calculated between the current increment and the increment at the end of the equilibration phase. Strain is calculated at the center of each element using neighboring elements within a distance twice the maximum radius of the elements in the assembly (minimum three elements are required to compute strain; Cardozo and Allmendinger, 2009). Stress is calculated from the element forces (Equations 1–4) and areas.

The simulation can be saved as a document file with .cdem extension, or increments can be saved as image files or animated GIFs. Increments can also be exported as text files for opening in other programs such as SSPX (Cardozo and Allmendinger, 2009). Finally, biaxial (strain-stress) data can be saved as text files. cdem files store all display increments in the simulation, and thus they are large (the default simulation in the Demo program is 10 MB). Therefore, one must exercise care when saving these files. A sensible strategy is, if possible, to increase the display interval (4 in Fig. 1A), thus reducing the number of increments that are saved.

General settings that apply to any document can be determined via the Preferences panel (Fig. 3). Here, one can set the colors of the pre-growth and growth strata as well as highlight a group of layers (1 in Fig. 3). Note that the growth strata can be divided in two sequences whose thickness can be specified by entering the number of layers in the lower growth sequence 1. One can also set the quality of the animated GIFs and reduce their size by including only odd increments (2). Stress is calculated if the Calculate stress checkbox (3 in Fig. 3) is on (by default it is not). The next checkbox (4 in Fig. 3) controls whether the spring and damping constants are calculated as in Table 1 (checkbox on) or are modified by unit scaling (checkbox off; this is the default). If the checkbox is off, the spring and damping constants for unit scaling <1000 m are calculated as follows:

where R_avg is the average radius of the elements. This provides shorter and stable runtimes, but it results in overlap of the elements and a “weaker” material (see the Assembly Calibration and Simulation of Tectonic Structures sections).

The second-to-last setting in the Preferences panel increases the height of the lateral walls (5 in Fig. 3); this is convenient in contractional simulations (e.g., detachment simulations) where the thickness of the assembly can get high and elements may fall out of the lateral walls. The last checkbox allows reducing the element radii by half (6). Toggling on this checkbox makes the element radii 0.025, 0.0375, 0.05, and 0.0625 units, thus producing a denser assembly with four times as many elements.

Table 2 includes the spring and damping constants and element radii general settings for all the cdem simulations included in the paper.

In cdem, the element properties and the emergent bulk material properties can be assessed through biaxial and angle-of-repose tests. In the biaxial test, gravity is turned off and an assembly of 10 × 5 units is shortened while a user-defined confining pressure is applied to the elements on the free edges of the assembly (these are highlighted in green). During equilibration and before shortening, the confining pressure is raised in a ramp fashion, from zero at the middle of the equilibration period to its full value at the end of equilibration. Therefore, it is important to use a long equilibration time. We normally use twice the equilibration time for biaxial tests than for other types of simulations.

We start with the default unit scaling (1000 m, 10 × 5 km samples), frictional-cohesive element properties (μ = 0.25, C₀ = 5 × 10⁹ N, first slip on; Fig. 1A), and element radii. The assembly has 2121 elements and is equilibrated for 800 s, and the total runtime is 2220 s to reach an axial strain of 0.2. We run biaxial tests at 0, 5, 10, 15 and 20 MN/m (MPa with consideration of across-section area). Figure 4A shows the results of the tests as axial strain versus differential stress (left) and Mohr’s circle (right) diagrams. The assembly displays a realistic stress-strain behavior, and strain localization is evident in the sample geometry and total maximum shear strain (γ_max) (Fig. 4A, left). The bulk behavior of the assembly is characterized by cohesion (C) of ~4 MPa and friction angle of ~34°.

We now look at the effect of reducing the radii of the elements by half (turning on the half radii checkbox in the Preferences panel; option 6 in Fig. 3). This results in a denser assembly of 8342 elements, but because the time step (dt) is proportional to the minimum mass (and radius) of the elements (Table 1), dt and the shortening rate are halved. Thus, we use twice the previous equilibration and total runtimes, namely 1600 and 4440 s, respectively, to reach an axial strain of 0.2. Figure 4B shows the results of the biaxial tests from 0–20 MPa confining pressure. The stress-strain behavior of the denser assembly is similar, but the assembly is “stiffer” (lower axial strain at failure). Also, because the elements are smaller, faulting is more localized (geometry and γ_max in Fig. 4B, left). However, the bulk behavior of the assembly is about the same with C of ~4 MPa and friction angle of ~34°.

To explore the effect of unit scaling and the spring and damping constants, we run biaxial tests with a smaller unit scaling of 200 m (2 × 1 km samples) and again default element radii (2121 elements). However, because in this case the elements are smaller, we reduce the inter-element cohesion (C₀) to 1 × 10⁷ N. We first run the tests using the default spring and damping constants (turning on checkbox 4 of the Preferences panel, Fig. 3). Under this condition, dt and the shortening rate are five times smaller than at unit scaling 1000 m (because the elements are five times smaller; Table 1), but the same equilibration and total runtimes of 800 and 2220 s, respectively, are required to reach 0.2 axial strain (because the sample is five times shorter). Figure 5A shows the results of the biaxial tests from 0–20 MPa confining pressure. The assembly is weaker, less stiff, and less cohesive (C ~1 MPa) than the larger assembly at unit scaling 1000 m (Fig. 4A). However, the bulk behavior of the assembly is still frictional, with a similar friction angle (~34°) and shear band (fault) orientations as the larger assembly at unit scaling 1000 m (compare Figs. 5A and 4A).

We run the biaxial tests again but this time modifying the spring and damping constants by unit scaling (turning off checkbox 4 of the Preferences panel, Fig. 3). This results in a larger dt and shortening rate (about twice as much as the default values; Equation 5). Therefore, shorter equilibration and runtimes of 400 and 1050 s, respectively, are required to reach 0.2 axial strain. Figure 5B shows the results of the biaxial tests. In comparison to the previous case with default spring and damping constants (Fig. 5A), the modified spring and damping constants make the simulations run faster, but they result in more element overlap, and a weaker and more ductile assembly, with less strain localization, and a lower friction angle of 25° (Fig. 5B). This illustrates the importance of the spring and damping constants. If running time is not a major concern, we recommend not varying these constants by unit scaling (using Table 1), although varying the constants as in Equation 5 can result in some interesting behavior in terms of strain localization behavior (see Planar Faults in the Simulation of Tectonic Structures section).

A faster yet reasonable way to calibrate the assembly is by performing a repose-collapse test. In this test, the assembly is equilibrated within a box delimited by walls, but at the end of the equilibration stage, one or both lateral walls are removed, producing the collapse of the assembly. Figure 6 shows repose-collapse tests (right wall removed) for the assemblies of Figures 4 and 5. All tests are run with the same dimensions and parameters as the biaxial tests but with longer runtimes to ensure slope stability and very small incremental displacements at the end (Fig. 6). For the more brittle assemblies of Figures 4 and 5A, the repose angles from the collapse tests (Figs. 6A–6C) are close to the friction angles from the biaxial tests. For the more ductile assembly of Figure 5B, this equivalence breaks down, and the repose angle from the collapse test (Fig. 6D) is higher than the friction angle from the biaxial tests. Also in the more ductile assembly, the slope is unstable for a longer time (~6× the total runtime of the biaxial tests) than in the brittle assemblies (~3× the total runtime of the biaxial tests) (Fig. 6).

The displacement boundary conditions in cdem allow for movement along a planar normal or reverse fault, two vertical faults under extension (caldera) or contraction (piston), and a horizontal contractional detachment (12 in Fig. 1A). Figure 7 shows simulations of a low-angle 30°-dipping normal fault (Fig. 7A), a high-angle 50°- dipping reverse fault (Fig. 7B), a caldera (Fig. 7C), and a detachment fault with the lateral walls moving inwards (Fig. 7D). In all cases, we use unit scaling 1000 m; equilibration and total runtimes of 1600 and 8000 s, respectively; gravity, density, and frictional-cohesive element properties as in Figure 1A; and half element radii. The initial size of the assembly is 10 × 5 km (8342 elements) for all cases except the detachment model, where the assembly is 15 × 3 km (7600 elements). For the extensional cases (normal fault and caldera), we use first slip off, while for the contractional cases (reverse fault and detachment) we use first slip on (8 in Fig. 1A). In addition, for the detachment model, we set no friction or cohesion between the elements and the base wall to simulate a “weak” detachment.

cdem produces realistic geometries and finite strains (Fig. 7). From the applied basal (Figs. 7A–7C) or lateral (Fig. 7D) displacement boundary conditions, faults propagate through the assembly and deform and offset the layers. The low-angle normal fault generates a graben limited by both synthetic and antithetic faults (Fig. 7A). These results are consistent with sandbox and finite differences models by Exadaktylos et al. (2003). The reverse fault produces a monocline with a triangular zone of deformation that expands outwards (Fig. 7B). This resembles clay box models by Bonanno et al. (2017, their isotropic experiment) and finite element models by Albertz and Lingrey (2012). Caldera collapse is accomplished by two high-angle faults that propagate and splay outwards; displacement is normal on the caldera margins but is reverse in the central part of the caldera (Fig. 7C). These results are compatible with analogue models by Roche et al. (2000) and DEM models by Hardy (2008). Finally, shortening and sliding of the assembly along the detachment produces two thrust-related anticlines of opposite vergence (Fig. 7D). Similar structures are observed in FTBs with a weak basal detachment, e.g., the Niger delta outer FTBs (Corredor et al., 2005). In all four of these simulations, much can be learned about the evolution of the structures by displaying both incremental and total displacement and strain as well as stress and observing these variables through time in cdem.

To further illustrate the impact of the spring and damping constants, we run the caldera model (Fig. 7C) but with a smaller unit scaling of 200 m, reduced inter-element cohesion (C₀) of 1 × 107 N, and as in Figure 7C, first slip off. Figure 8A shows the geometry and finite strain of the simulation with the spring and damping constants as in Table 1. Here, the caldera-bounding fault zones are relatively wide and there is folding of the uppermost layers. Figure 8B shows the simulation with the spring and damping constants modified by unit scaling as in Equation 5. Here, the assembly is more brittle and displays narrower fault zones. The surface of the caldera is irregular and exhibits several fault scarps as opposed to the more continuous folded surface of Figure 8A. The first model (Fig. 8A) is perhaps more realistic (Roche et al., 2000), but by using in Figure 8B a larger time step and approximately the same damping value as in Figure 8A, it is possible to simulate more strain localization and a more brittle caldera.

In cdem, non-planar faults are introduced by “sculpting” the left wall into an irregular fault. This also removes the footwall elements, so we focus only on the hanging-wall deformation. The default geometry of the non-planar fault is listric with a lowermost fault dip that can be specified in the “Fault dip” field (12 in Fig. 1A) and an uppermost dip of 90°. However, it is also possible to introduce any irregular fault shape using coordinates of the fault vertices from a text file (Model → Irregular fault menu). In any case, there are two possible displacement boundary conditions. The first one is constant heave (CH checkbox in Fig. 1A; Verrall, 1981), where movement along the fault is simulated by displacing the fault elements horizontally and the hanging wall vertically to fit the displacement along the lowermost fault segment. In this case, it is important to set the element-fault friction coefficient and cohesion to zero (8 in Fig. 1A). The second one is constant displacement (CD checkbox; Williams and Vann, 1987), where movement along the fault is simulated by displacing the fault elements parallel to the fault. This is analogue to the mylar sheet condition typically used in sandbox experiments (Ellis and McClay, 1988). In this case, it is important to have non-zero element-fault friction and cohesion (8 in Fig. 1A).

Figure 9 shows simulations of a listric normal fault under constant heave (Fig. 9A) or constant displacement (Fig. 9B), and a normal fault with irregular geometry under constant displacement (Figs. 9C and 9D). In all three simulations, we use unit scaling 1000 m; equilibration and total runtimes of 1600 and 10,000 s, respectively; gravity, density, and element properties as in Figure 1A; first slip off; and half element radii. Also, before removal of the footwall elements, the assembly is 15 × 5 km. Movement along a listric normal fault results in a rollover structure (Figs. 9A and 9B). However, the choice of displacement boundary conditions is important. Under constant heave, the hanging wall collapses mainly by movement along antithetic normal faults, some of them displaying fault scarps (Fig. 9A). Constant displacement leads to more distributed deformation and a gentler fold surface, although there is more strain localization on a zone between the listric and planar fault domains (Fig. 9B). Movement along the irregular normal fault produces a hanging-wall anticline (Figs. 9C and 9D). However, the evolution of this structure is more complex than in the listric fault case. From the start to about halfway through the simulation, a hanging-wall syncline forms in the uppermost layers close to the fault (Fig. 9C). This syncline collapses as fault displacement continues (Fig. 9D). Hanging-wall deformation is accomplished by both synthetic and antithetic normal faults. Incremental shear strain (γ_max, not included in Fig. 9) indicates that antithetic faults nucleate from the concave-upwards fault bends, while synthetic faults nucleate from the convex-upwards fault bends. This is consistent with the kinematic model of Xiao and Suppe (1992).

So far, we have experimented with a homogeneous assembly. Now we look at the influence of assembly heterogeneity by making either the layer contacts frictionless (flexural slip; 9 in Fig. 1A) or the lower part of the assembly frictionless (10 in Fig. 1A). Figure 10 shows the low-angle normal fault simulation of Figure 7A but with no friction between the layers (Fig. 10A) or with the lower half of the assembly frictionless (Fig. 10B). In addition, in both simulations, we use no friction and/or cohesion between the body and the wall elements (8 in Fig. 1A).

No friction between the layers promotes layer-parallel slip. For the low-angle normal fault, this results in more strain localization (narrower fault zones) and a better-defined graben than in the case of the homogeneous assembly (compare Figs. 10A and 7A). A weak frictionless substrate results in more distributed deformation and no graben formation (compare Fig. 10B). Normal faults (60°–70° dip) offset the upper frictional part of the assembly but tip out in the lower frictionless part, which is mainly dragged along the basal master fault. The mechanical stratigraphy of this simulation resembles brittle sandstone over weak over-pressured shale, and the results are quite like actual structures and sandbox experiments of such fault systems (Gabrielsen et al., 2019).

Here, it could be interesting to try also layered friction (11 in Fig. 1A) by, for example, having competent layers with friction coefficient = 0.25 and incompetent layers with friction coefficient = 0.1. These results are not included here, but a similar DEM study on the role of mechanical stratigraphy on extensional fault-propagation folding was conducted by Hardy (2019b).

Syn-tectonic (growth) strata not only record the evolution of tectonic structures but also impact how the structures develop. In cdem, syn-tectonic sedimentation can be introduced by checking the Syngrowth checkbox in the Set-up view (15 in Fig. 1A). This results in the addition of new elements at a given sediment interval and the equilibration of these elements during a given settling time when there is no boundary displacement. The choice of sediment interval and settling time is not trivial. In cdem, they are defined as follows:

Thus, the sediment interval is proportional to the display interval (4 in Fig. 1A), and the settling time is proportional to the equilibration time (5 in Fig. 1A). Therefore, one must select carefully these two parameters such that elements are added not too often (not enough accommodation space) and not too seldom (space to fill is too big) and the elements have enough time to equilibrate without unnecessarily slowing down the simulation. Also, one may end up with the wrong scenario that the time for settling the sediments is larger than the time between sediment intervals. If this is the case, the program will report an error, which should be corrected by either increasing the display interval (to increase the time between sedimentation events) or decreasing the equilibration time (to reduce the settling time).

The upper boundary of the space available for new elements (the accommodation space) is sea level. Initially, sea level is equal to the maximum height of the assembly after equilibration (y_max). Sea level can be static or rise throughout the simulation at a low, mid, or high rate (15 in Fig. 1A). After equilibration, sea level is calculated as:

where u̇ is the displacement rate (displacement per time step), and factor is 0.0 for null, 0.25 for low, 0.5 for mid, and 0.75 for high sea-level rise.

Figure 11 shows the 30°-dipping normal fault models of Figures 7A and 10 but with syn-tectonic sedimentation and static sea level. Sediment growth leads to steeper normal faults and a narrower graben in both the homogeneous assembly (compare Figs. 11A and 7A) and the heterogeneous assembly with flexural slip (compare Figs. 11B and 10A). In the assembly with a weak frictionless substrate, sediment growth leads to more deformation of the substrate and more strain localization around the basal master fault (compare Figs. 11C and 10B). Sediment growth increases the overburden stress, which enhances strain localization and shear failure in all three models (Kettermann and Urai, 2015) and leads to more deformation of the weak substrate along the master fault in the third model of Figure 11C (Gabrielsen et al., 2019).

Figure 12 shows the final geometry and total maximum shear strain of FTB models without (Fig. 12A) and with (Figs. 12B–12C) syn-tectonic sedimentation. In all models, we use unit scaling of 1000 m; gravity, density, and element properties as in Figure 1A (except for no friction on the base wall); first slip on; flexural slip; and half element radii. The initial size of the assembly is 25 × 2.5 km (initial number of elements is 10,520), and the assembly is shortened by moving the right wall to the left a total amount of 7.5 km (30% shortening). In the simulation without syn-tectonic sedimentation, a forward-breaking sequence and a thrust-wedge geometry are observed (Fig. 12A). The introduction of syn-tectonic sedimentation with low (Fig. 12B) or mid (Fig. 12C) sea-level rise produces a wedge composed of fewer major forward-breaking thrusts, with almost no fault activity in the foreland. With increasing sea-level rise and sedimentation from Figure 12B (final number of elements = 13,944) to Figure 12C (final number of elements = 18,679), the frontal thrust has greater displacement and offsets a thicker package of growth strata. These findings are consistent with those of Hardy and Cardozo (2021), who documented a change of FTB geometry from a typical wedge in the case of no syn-tectonic sedimentation to a frontal steep ramp accumulating most of the FTB displacement in the case of high syn-tectonic sedimentation. These authors used cdem2D simulations with a lower unit scaling and element radii (about one-tenth of those in Fig. 12) and a denser assembly (~10× as many elements as in Fig. 12).

When the unit scaling is small (~100 m, meter-size elements) and the assembly is very dense (tens of thousands of elements), the simulation in cdem can take very long, and it can be better to use a program such as cdem2D. As mentioned earlier, cdem2D is a 2-D DEM code written in C that uses the OpenMP application programming interface for parallelization and effective use of all cores in a machine. Compilation of the code (using the Intel compiler) produces an executable. This executable requires a text file that defines the assembly (this file is also internally present in cdem) and a runtime text file that defines the model parameters (which are like those in cdem, Fig. 1A). Running the executable from a terminal or command window produces a model parameters text file with the details of the simulation (like Fig. 1B) and, throughout the simulation, increments or text files at a user-specified display interval. Each one of these files contains the center x and y coordinates, radii, and layer number of the elements at the increment. In addition, if stresses are calculated, they are also included in the text file.

This information is essentially the same as that defining the element arrays in cdem. Therefore, we added to cdem the capability to import the increments of a cdem2D simulation. This is achieved by choosing in cdem the File → Import cdem2D simulation menu option and selecting a folder containing the increments. The program imports and displays sequentially the increments in the Results view as well as the details of the simulation (from the model parameters file) in the Summary view. Once all the increments are imported, they can be examined as in any cdem simulation. Any increment can be chosen and visualized in terms of displacement, strain, or stress (if computed), and the simulation can be played forward and backward. Although it is possible to save the imported simulation as a single cdem file, we don’t recommend doing so because the file can be very large and take a long time to open. Importing cdem2D increments into cdem does not take long, so the user may keep the original simulation folder.

This functionality allows cdem to work as a post-processor for cdem2D. Figure 13 illustrates the possible interaction between these two programs. The first step is to prepare the input (runtime file) for cdem2D (gray circle, Fig. 13); the second step is to run cdem2D, which produces a series of increments. Then, either while cdem2D is running or when it has terminated, one can import and analyze the increments in cdem. Thus, cdem can be used to check and, if desired, stop ongoing cdem2D simulations or process completed simulations. This process can be repeated by modifying the runtime file for an alternative scenario or condition and going through the loop again.

Figure 14 shows a cdem2D simulation of positive tectonic inversion, imported into cdem. Unit scaling is 125 m, the initial assembly size is 6.25 × 2.8 km, and the initial number of elements is 46,335 (average element radius is ~10 m). The assembly is frictional-cohesive (μ = 0.25, C₀ = 6 × 10⁷ N), and there is syn-tectonic sedimentation and low sea-level rise. Three stages are shown: at the maximum normal displacement just before inversion (54,288 elements; Fig. 14A), after 50% inversion (55,665 elements; Fig. 14B), and after 100% inversion (55,855 elements; Fig. 14C). Note that this model is more refined than the cdem models we have included before, and because the elements are quite small (meter size), the faults, folds, and growth strata are better defined. Processing of the simulation reveals useful information. For example, coloring of the assembly by uplift (red) or subsidence (blue) displays well the transition of the structure from extension to compression and the movement of the compressional deformation (and fault null points) downwards as tectonic inversion increases (Williams et al., 1989).

Denser and more refined simulations are not the only advantage of cdem2D. The code offers other functionalities not present in cdem, such as viscous, time-dependent deformation (salt deformation) and more complex boundary conditions (e.g., keystone graben, regional rotation) and syn-tectonic sedimentation (e.g., deltas) (Hardy, 2018, 2019a). In addition, one has more control over the model parameters via the runtime file, or even through the code itself, which is more flexible because there is no user interface.

We have described the functionality of cdem, an intuitive and powerful macOS program for discrete element modeling of tectonic structures and their associated deformation. DEM modeling of geological structures is not trivial, but cdem facilitates this process by including a rectangular assembly that can be resized, scaled, and sculpted. In addition, critical computation parameters, such as the time step and damping constant, are calculated in terms of more intuitive parameters, such as unit scaling and element density (Table 1). This offers great flexibility, but it can also lead to unexpected, undesired results, so one should be careful. Element properties, either initially bonded or frictional-cohesive, result in realistic bulk behavior and the natural development of faults and folds. Calibrating the assembly is a difficult problem, but the program offers intuitive tools to do so via biaxial (Figs. 4 and 5) or repose-collapse (Fig. 6) tests. Complex structures can be simulated with relatively simple boundary conditions (12 in Fig. 1A); for example, a graben can be formed by displacement along a planar low-angle normal fault (Fig. 7a), or a FTB by shortening and sliding of the assembly over a weak detachment (Fig. 12). Overburden load, erosion, and syn-tectonic sedimentation are easily introduced by defining an overburden pressure (13 in Fig. 1A), a base level for erosion (14 in Fig. 1A), or a static or rising upper boundary (sea level) for sedimentation (15 in Fig. 1A). In any case, these three parameters impact the evolution of the structure (Figs. 11 and 12). Besides these core functionalities, perhaps the greatest advantage of cdem is its visual display. There is a lot of information encapsulated in the increments of a DEM model. Being able to display this information in terms of displacement, strain, or stress and examining the evolution of the structure by playing these variables through time is a game changer. The application of this functionality to imported cdem2D simulations has allowed us to establish a robust and friendly modeling environment for denser and more refined DEM simulations (Figs. 13 and 14).

Because cdem (and any DEM) is computationally intensive, it is not the right tool to precisely fit the geometry of a structure or to model the structure’s uncertainties (Cardozo and Oakley, 2019). Kinematic models and simpler mechanical models (e.g., elastic models) are better for those applications. The benefits of DEM for structural modeling are nicely summarized by Gray et al. (2014), and they boil down to: (1) testing the mechanical validity of geological interpretations and (2) constraining the properties of the rock materials that compose these structures, both at the scale of observation (e.g., seismic) and at smaller scales (e.g., sub-seismic). Understanding the mechanical evolution of a structure and the impact boundary conditions, rock properties, and mechanical stratigraphy have on the structure also helps in interpreting the structure. This is particularly important for seismic interpreters working in structurally complex areas where seismic imaging may be poor, and much can be gained by validating interpretations against mechanical models. Finally, DEM models contain large amounts of data, and these data could be used to derive mechanically consistent yet simpler and more efficient models using data-driven science and engineering (Karpatne et al., 2019; Brunton and Kutz, 2022). This is an exciting area of research.

cdem is continuously evolving. Future improvements include tracking the connectivity of elements along fault zones (fault sealing), modeling a viscous, time-dependent substrate (salt), and introducing lateral facies changes and unconformities, and modeling several faults that can act simultaneously or at different times. Better processors and parallelization technologies in the future will mean a faster program. cdem is written in Swift, and therefore it is relatively simple to implement it in mobile iOS devices. Finally, we are not limited to 2-D. Author Hardy has implemented a three-dimensional DEM code called cdem3D. We could use cdem to read 2-D slices from cdem3D, or alternatively we could implement cdem3D as an additional functionality in cdem. cdem is robust and flexible enough to explore these possibilities.

cdem has been our hobby during these past six years. We thank our families for their patience, and our employers for allowing us to collaborate informally at a distance. We thank the Validé-Plogen program in Stavanger (Norway) for seed funding, which made possible the implementation of the program. cdem is free for academic, non-profit use, and it can be obtained by contacting Nestor Cardozo at nestor.cardozo@uis.no. Likewise, for access to cdem2D, please contact Stuart Hardy at stuart.hardy@icrea.cat.