Elucidating fault-related fold mechanics: a 2D finite element analysis of bending, slip, and buckling mechanisms

1 Introduction

Fault-related folds, evident in both compressional and extensional regimes (Mitra, 1993; McClay, 2011; Morley et al., 2011), are of significant geological interest, particularly due to their potential as structural traps for hydrocarbon resources (McClay, 1995; Kent and Dasgupta, 2004). Owing to their economic importance, numerous researchers have investigated fault-related folds through various approaches such as field studies, seismic research, sandbox experiments, and numerical modelling, as illustrated in works by McClay (1995), Allmendinger (1998), Hughes and Shaw (2015) and Ghanbarian et al. (2021). Various numerical modeling techniques have been employed in the study of geological structures. These include finite element analysis (Smart et al., 2012; Yang et al., 2017; Li et al., 2019; Khalifeh-Soltani et al., 2021a; Khalifeh-Soltani et al., 2021b; Khalifeh-Soltani et al., 2021c; Khalifeh-Soltani et al., 2023a; Khalifeh-Soltani et al., 2023b), boundary element analysis (Cooke and Pollard, 1997; Desai et al., 2011; Maerten et al., 2014), and discrete element analysis (Finch et al., 2003; Hardy and Finch, 2006). Moreover, numerical models that utilize velocity vectors have been developed for fault-related folds (e.g., Johnson and Berger, 1989; Hardy and Poblet, 1995; Hardy et al., 1996; Allmendinger, 1998; Cardozo, 2005).

Fault-related folds consist of three end-members: fault-bend folds (e.g., Berger and Johnson, 1980; Suppe, 1983; Brandes and Tanner, 2014), fault-propagation folds (e.g., Mitra, 1990; Suppe and Medwedeff, 1990; Brandes and Tanner, 2014), and detachment folds (e.g., Epard and Groshong, 1995; Homza and Wallace, 1995; Poblet and McClay, 1996; Homza and Wallace, 1997). The fault-bend fold developed above a flat-ramp-flat that is called a single-step fault-bend fold. Fault-propagation folds are caused by the gradual reduction of slip on a propagating thrust fault (Suppe and Medwedeff, 1990). These folds are characterized by steep or overturned forelimbs (Suppe and Medwedeff, 1990; Hughes et al., 2014; Hughes and Shaw, 2015). Detachment folds may form above or below a detachment or a decollement (e.g., Poblet and Hardy, 1995; Poblet et al., 1997; Peacock et al., 2000). The detachment fold grows by moving the detachment layer’s particles toward the centre of the fold (Homza and Wallace, 1995; Poblet and McClay, 1996).

Since fault-related folds can be hydrocarbon traps, a correct understanding of their kinematics can help us determine the temporal relationship of structure growth with maturity and migration of hydrocarbon resources. On the other hand, the strain distribution in the reservoir can control its permeability (Hughes and Shaw, 2015). Therefore, there is a need for a kinematic model that can clarify the mechanisms governing these folds and provide its stress and strain history, which, unfortunately, most of the kinematic-geometric models that have been presented so far, despite their advantages, have some limitations: 1) The primary assumption of these models is based on some geometric simplifications; for example, they consider the length and thickness of the layers to be constant during folding, which is not the case in real examples, and hence they cannot accurately represent the geometry of real examples, 2) The role of effective parameters on folding, such as the mechanical properties of materials, variation in stress and strain, and the interface characteristics of layers, are not considered (Allmendinger, 1998; Logan, 2010; Smart et al., 2012), and 3) Most presented kinematic models are based on the classical kink-band migration mechanism. The kink band migration model also has limitations: i) The conditions for forming kink folds are also ideal. According to Ramsay and Huber (1987), when the contrast in thickness and viscosity of the layers is high, kink folds are created. In natural examples, except for the viscosity contrast of the decollement and other layers, the viscosity and thickness of the sedimentary layers forming the fault-related folds are usually similar, ii) The simple shear parallel to the layering creates the kink folds (e.g., Ramsay and Huber, 1987). While buckling mechanism, which plays a role in the formation of fault-related folds, is neglected in these models, especially in detachment and fault-propagation folds (Barani, 2012; Li et al., 2018; Butler et al., 2020; Ghanbarian and Derakhshani, 2022; Rashidi et al., 2023).

Unlike kinematic-geometric models, geomechanical models do not have these limitations: 1) These models have no geometric simplifications, the length and thickness of the layers change during folding, so they create folds closer to actual examples (Smart et al., 2012), 2) They consider the influence of the mechanical properties of materials, the characteristics of the interface of layers and faults, 3) they can provide the evolution of stress and strain at any moment and any location of the model, and 4) Their setup does not limit the participation of interlayer sliding, bending and buckling mechanisms. Therefore, geomechanical models provide a more appropriate insight into the mechanism of fault-related folds than kinematic-geometric models. This study examines the mechanics of fault-related folds through a set of geomechanical models for the first time.

Since inter-layering slip, buckling, and bending mechanisms influence the development of folds (Twiss and Moores, 1992). Here, we investigate the contribution of each of them in the development of fault-related folds through five series of 2D finite element models. For this purpose, we first simulate three end members of fault-related folds (reference models) using the 2D finite element method (Series R). Then, we fold single-layers and multi-layers by bending and buckling mechanisms with two different elastic and elastic-plastic rheologies (Series A and B). Also, by removing interlayer slip, we examine its effect on buckling and bending (Series C and D). Finally, the strain ellipses of these models (Series A- D) are compared with reference models (Series R).

2 Modeling

In this study, five series of 2D elastic and elastic-plastic finite element models to determine the contribution of each of the mechanisms of bending, buckling and inter layering slip in creating fault-related folds are presented (R, A, B, C and D). All the models of this study ran in ABAQUS™ software (version 2017).

The first series (R) consists of three elastic-plastic R₁, R_2, and R₃ models, which simulate the three end members of fault-related folds, i.e., detachment, fault-propagation, and fault-bend folds. In fact, these are reference models to which the results of other models are compared. Series A is a set of elastic and elastic-plastic models that simulate bending at single-layer and multi-layer (with interlayer slip). The elastic and elastic-plastic models are listed in Table 1. Series B simulates buckling at single-layer and multi-layer (with interlayer slip) by a set of elastic and elastic-plastic models. As well as series C simulates non-slip buckling in single-layer and multi-layer modes. Series D investigates the role of non-slip bending and buckling. The difference in the buckling models of series C and D is in the equations used for modeling, which will be described in the following sections.

TABLE 1. Rheology and folding mechanism of the models.

2.1 Modeling equations

In all the elastic-plastic models in series R, A, B and D (Table 1), Hooke’s law for the portion of elastic behaviour of the models is used, and the standard Mohr-Coulomb criterion for the part of elastic behaviour of the models is used. Only Hooke’s law for the elastic models in series A, B, C, and D is used (Table 1). In series C, the buckling equation is used instead of the standard Mohr-Coulomb criterion to simulate buckling folding. It is described below.

Abaqus software uses an incremental loading pattern, Q^N, to calculate the eigenvalue buckling:

Where K₀^NM is the stiffness matrix corresponding to the base state, considering that the models presented in this study include one step (buckling step), their base state is the same as the pre-deformation model, so the preload effects (P^N) are not present in them. K_Δ^NM is the initial differential stress and load stiffness matrix due to the incremental loading pattern (Q^N). λ_i is the eigenvalues. υ_i^M is the buckling mode shapes (eigenvectors). M and N refer to degrees of freedom of the whole model, and i refers to the ith buckling mode.

Here, we examine three outputs, elastic and plastic strain and Mises stress of elastic-plastic models, and two outputs of elastic stress and Mises stress of elastic models. Elastic and plastic strains represent reversible and permanent deformation of the models, respectively. Mises stress is a measure of the combined magnitude of all stress components (tensile, compression, and shear) at any point. Therefore, it is extremely useful for predicting the failure modes in structures.

2.2 Geometry and mechanical properties of models

The dimensions (length and thickness) and mechanical properties of all the models are similar (Figure 1; Table 2). The material properties of the layers (i.e., density, Young’s modulus, Poisson’s ratio, cohesion, internal friction angle and dilation angle) are listed in Table 2. Like our previous works, these data were derived from the results of rock mechanics experiments carried out by the National Iranian South Oilfields Company (Khalifeh-Soltani et al., 2021a; Khalifeh-Soltani et al., 2021b; Khalifeh-Soltani et al., 2023a).

FIGURE 1. Boundary conditions of series R models; (A) detachment fold model (R₁), (B) fault-propagation fold model (R₂), and (C) fault-bend fold model (R₃).

TABLE 2. Material properties are used in all finite element models (Khalifeh-Soltani et al., 2021b; Khalifeh-Soltani et al., 2021c).

Series R contains five layers. The ramp dip in the models of the fault-propagation (R₂) and fault-bend folds (R₃) is 45° and 25°, respectively (Figure 1). While series A, B, C and D include four layers (Figure 2). Because layer 1) does not play a role in folding in reference models (R), it has been omitted.

FIGURE 2. Boundary conditions of the models; (A) single-layer bending at three points (A₁ to A₈ models), (B) multi-layer bending at three points (A₉, A₁₀, D₁, and D₂ models), (C) multi-layer bending at four points (A₁₁, A₁₂ and D₃ models), (D) single-layer buckling (B₁, B₂, C₁ to C₄ models, (E) multi-layer buckling (B₃ to B₅, C₅ and D₄ models).

It is necessary to mention here that, since fault-related folds are mainly created in sedimentary rocks, and except for detachment levels, the difference in mechanical properties in layers is not significant, the results of this study can be valid even for other regions. However, to generalize the results to other areas, the lithology influence should not be ignored.

2.3 Friction coefficient

For the interaction of fault and layers in series R, friction coefficients of 0.01 and 0.25 are considered, respectively. A friction coefficient of 0.25 is considered for the interaction of layering surfaces in multi-layer models of series A and B, while multi-layer models of series C and D have no interlayer slip.

2.4 Setup and boundary conditions

The models presented here include one step, and displacement is applied to simulate folding in that step. The boundary conditions of the reference models (R) are shown in Figure 1. As shown in all three models, the model’s left side is not allowed to move horizontally during folding (X = 0; Where X, Y, and Z are displacements along the X, Y, and Z-axes, respectively). At the same time, displacement is applied to the right side of the model (Figure 1). The displacement in the detachment, fault-propagation, and fault-bend fold models is 100, 75, and 150 m, respectively. Footwalls in the fault-propagation and fault-bend fold models and also layer (1) (the basement) in the detachment fold model are fixed (X=Y=0).

Figures 2A–C shows the boundary conditions of the models that simulate single-layer and multi-layer bending with slip (series A) and non-slip (D₁ to D₃ models). Figures 2A, B show the boundary conditions of single-layer and multi-layer bending at three points. Figure 2C shows the boundary conditions of multi-layer bending at four points. In all three cases, the right and left sides of the model are not allowed to move horizontally and vertically (X = Y = 0). Also, in all three cases, two rigid circular bodies with a diameter of 40 m are fixed on the model top at a distance of 120 m from the model’s edge (X = Y =0). In the case of bending models at three points, a rigid circular body is located under the centre of the models, and the bending at three points is simulated by its upward movement (100 m) (Figures 2A, B). In the case of bending models at four points, two rigid circular bodies are located under the model centre at a distance of 200 m from each other; the bending at four points is simulated by their simultaneous upward movement (100 m) (Figure 2C). Figure 2D shows the boundary conditions of buckling in single-layer models of series B and C. In the same way, Figure 2E shows the boundary conditions of buckling in multi-layer models for both sliding (series B) and non-sliding conditions (series C). In all models shown in Figures 2D, E, the model’s left side cannot move in the horizontal direction (X = 0), while the displacement is applied to the model’s left side to simulate buckling (Figures 2D, E). Of course, the amount of displacement is different in these models (Actually, models with larger displacements have been added to examine the results more precisely).

All models’ meshes are square, with 10 × 10 m dimensions. We use the deformation of these squares to estimate the final strain ellipse of the models. Of course, in the basement and footwalls of the reference models (series R), the dimensions of meshes are large and irregular because these parts are fixed during folding.

3 Modeling results

Here, we present five series of 2D elastic and elastic-plastic finite element models to determine the contribution of interlayer slip buckling, and bending mechanisms to develop fault-related folds. The results are as follows.

3.1 Reference models (series R)

Series R includes three reference models, which simulate detachment (R₁), fault-propagation (R₂), and fault-bend folds (R₃). Then, the geometric evolution, distribution of elastic and plastic strains, and Mises stress in these models are investigated.

The detachment fold model creates a symmetric fold in which its amplitude and wavelength increase with displacement increases (Figures 3A–C). The elastic strain is distributed throughout the model and is mainly concentrated at the crest and fold limbs (Figures 3A–C). Mises stress is distributed mainly in the third and fourth layers (especially the limbs) (Figure 3D). The plastic strain is concentrated at the detachment fault tip and distributed towards the axial surfaces (Figure 3E). The maximum fault slip occurs at the model’s boundary but decreases towards the detachment fault tip. On the contrary, the amount of interlayer slip is zero at the borders of the model and increases towards the model’s centre (Figures 3A–C). In addition, the fold crest moves forward during folding (Figures 3A–C).

FIGURE 3. Distribution of elastic strain and geometry of detachment fold model (R₁) in; (A) 33 m of displacement, (B) 66 m of displacement, (C) 100 m of displacement, (D) distribution of Mises stress in 100 m of displacement, (E) plastic strain distribution in 100 m of displacement. The maximum and minimum slip points during folding are shown in turquoise colour (A–C). Layers one to five are shown in (A)

In the fault-propagation fold model, the amplitude and wavelength of the fold increase while the interlimb angle decreases during the folding, and the fold becomes tighter (Figures 4A–C). The elastic strain is distributed throughout the model and is mainly concentrated in the crest and limbs at the fourth and fifth layers (Figures 4A–C). Mises stress is more concentrated in the fourth and fifth layers (Figure 4D). Plastic strain is localized on the fault surface, the fault tip, and the forelimb (Figure 4E). Similar to the detachment fold model, the maximum fault slip occurs at the model’s boundary (75 m), decreasing towards the model’s centre (38 m) and reaching zero at the fault tip (Figures 4A–C). On the contrary, at the model border, interlayer slip is zero and increases towards the centre.

FIGURE 4. Distribution of elastic strain and geometry of fault-propagation fold model (R₂) in; (A) 25 m of displacement, (B) 50 m of displacement, (C) 75 m of displacement, (D) distribution of Mises stress in 75 m of displacement, (E) plastic strain distribution in 75 m of displacement. The maximum and minimum interlayer slip points during folding are shown in turquoise colour (A–C). Slip changes on the flat-ramp are shown in white, and blak shows the ramp tip. Layers one to five are shown in (A)

In the fault-bend model, the fold amplitude and wavelength increase with increasing shortening (Figures 5A–C). In the early stages of folding, the elastic strain is distributed throughout the fourth and fifth layers (Figure 5A). But with the progress of folding, it is concentrated in the fold limbs (Figure 5C). Mises stress is focused on the second to fourth layers (Figure 5D). The plastic strain is concentrated on the fault surface, fault tip and forelimb of the fold (Figure 5E). The maximum fault slip is observed at the model boundary (150 m), decreasing towards the ramp (80 m) and reaching zero in the upper flat (Figures 5A–C). While the interlayer slip is zero at the model’s boundary and maximized at the forelimb (Figures 5A–C).

FIGURE 5. Distribution of elastic strain and geometry of fault-bend fold model (R₃) in; (A) 50 m of displacement, (B) 100 m of displacement, (C) 150 m of displacement, (D) distribution of Mises stress in 150 m of displacement, (E) plastic strain distribution in 150 m of displacement. The maximum and minimum interlayer slip points during folding are shown in turquoise colour (A–C) (A–C). Slip changes on the flat-ramp-falt are shown in white (layer 2) and yellow (layer 3). Layers one to five are shown in (A)

3.2 Bending-sliding folding in multi-layer and bending in single-layer (series A)

Series A consists of twelve models, and the first eight models (A₁ to A₈) simulate bending in the single-layer with both elastic and elastic-plastic behaviours. The odd models are elastic (i.e., A₁, A₃, A₅, and A₇), while even models are elastic-plastic (i.e., A₂, A₄, A₆, and A₈) (Table 1). Figure 2A shows the boundary conditions of these eight models. The results show that all the elastic models are entirely run until the end of the applied displacement (Figures 6A, C, E, G). While the elastic-plastic models stop after some displacement (Figures 6B, D, F, H). In elastic models, the Mises stress and elastic strain are concentrated in the model’s centre and the hinge zone (Figures 6A, C, E, G). In elastic-plastic models, the Mises stress and elastic strain are localized in the central part of the model and at the contact surface of the rigid circular bodies with the layers (Figures 6B, D, F, H). While in elastic models, the Mises stress and elastic strain are concentrated on the opposite side of the rigid circular bodies (Figures 6A, C, E, G). In elastic-plastic models, plastic strain is also localized in the centre of the model (Figures 6B, D, F, H).

FIGURE 6. Distribution of elastic, plastic strains, and Mises stress in single-layer elastic (A, C, E, G) and elastic-plastic (B, D, F, H) bending models. The odd models A₁, A₃, A₅ and A₇ are elastic, the even models A₂, A₄, A₆ and A₈ are elastic-plastic.

The next four models (A₉ to A₁₂) create bending in a multi-layer (Figures 7, 8). The rheology in these models is different; the odd and even models are elastic and elastic-plastic, respectively (Table. 1; Figures 7, 8). The A₉ and A₁₀ models simulate bending at three-point with interlayer slip. Figure 2B shows their boundary conditions. The elastic strain is distributed throughout the A₉ model (Figure 7A). The Mises stress is mainly distributed in layers (2) to (4) (Figure 7A). The elastic strain and Mises stress are increased in the A₁₀ model compared to the A₉ model (Figure 7). The plastic strain is concentrated around the rigid circular bodies (Figure 7B). The geometry of the elastic and elastic-plastic models is almost the same; only in the elastic-plastic models the rigid circular bodies are indented at the contact surface with the layers, indicating the model’s plastic rheology (Table. 1; Figure 7). Both elastic (A₉) and elastic-plastic (A₁₀) models are run until the end of the applied displacement. In contrast, the single-layer elastic-elastic models stopped after a few displacements. Maybe this is due to the effect of interlayer slip in multi-layer models; the results of the next series will reveal this problem (series D).

FIGURE 7. Distribution of elastic strain and Mises stress in multi-layer bending model at the three-point for; (A) elastic (A₉), and (B) elastic-plastic models (A₁₀). Layers one to five are shown in (A)

FIGURE 8. Distribution of elastic strain and Mises stress in multi-layer bending model at the four-point for; (A) elastic (A₁₁), and (B) elastic-plastic models (A₁₂). Layers one to five are shown in (A)

In the same way, A₁₁ and A₁₂ models simulate bending at four points with interlayer slip (Figure 8), whose boundary conditions are shown in Figure 2C The elastic strain and Mises stress are distributed throughout the A₁₁ model and are localized at the contact surface of the rigid circular bodies with the layers (Figure 8A). The Mises stress and elastic strain are increased in the A₁₂ model compared to the A₁₁ model (Figure 8). The plastic strain is concentrated around the rigid circular bodies (Figure 8B). Also, like the A₁₀ model, in the A₁₂ model, the rigid circular bodies are indented at the contact surface with the layers, indicating the model’s plastic rheology (Figures 7B, 8B).

Investigation slip variations in multi-layer models showed that in the hinge surface and the borders of the A₉ and A₁₀ multi-layer models, the slip is zero, and it towards limbs becomes maximum (Figure 7). While in models A₁₁ and A₁₂, the slip is zero only at the borders of the models, and it without a specific pattern is distributed throughout the model (Figure 8).

3.3 Buckling-sliding folding in multi-layer and buckling in single-layer (series B)

B series consists of five models that simulate single-layer and multi-layer buckling (with slip). The B₁ and B₂ single-layer models are elastic and elastic-plastic, respectively (Table. 1; Figure 9A). To summarize, only layer 2) is modelled for single-layer bending. The boundary conditions of these models are shown in Figure 2D. The B₁ elastic model is run until the end of applied displacement (Figure 9A). Nevertheless, any fold is not created, and only the thickness of the layer increases as pure shear (Figure 9A). Mises stress, and elastic strain are uniformly distributed throughout the fold and are localized at the model’s left border (Figure 9A). The B₂ elastic-plastic model stops after a few displacements and does not form any fold (Figure 9A). Mises stress, elastic strain, and plastic strain are concentrated at the model’s left border as in the B₂ model (Figure 9A).

FIGURE 9. (A) Distribution of elastic and plastic strains and Mises stress in single-layer buckling models with elastic (B₁) and elastic-plastic rheologies (B₂). Distribution of elastic and plastic strains and Mises stress in multi-layer buckling models for; (B) elastic-plastic rheology (B₃), (C, D) elastic rheology (B₄ and B₅). Layers one to five are shown in (B). The blue arrow shows the half-wavelength from which the strain ellipse presented in Figure 13 was extracted.

The B₃, B₄ and B₅ models simulate multi-layer buckling with interlayer slip (Figures 9B–D). The boundary conditions of these models are shown in Figure 2E. The rheology in the B₃ model and B₄ and B₅ models are elastic-plastic and elastic, respectively (Figures 9B–D). The difference between the B₄ and B₅ models is in the amount of displacement; their displacements are 100 and 150 m, respectively (Figures 9C, D). The B₃ multi-layer model, like the B₁ and B₂ single-layer models, stops after a few displacements and does not create any fold; Just the thickness of the layers increases as pure shear (Figure 9B). The elastic strain is distributed throughout the model (Figure 9B). The Mises stress is the maximum in the third layer. The plastic strain is localized at the base of the model’s left side (Figure 9B). In the B₄ model, a fold with a full wavelength is created. The Mises stress and elastic strain are concentrated in the hinge of the anticlines and synclines (Figure 9C). They are minimal at the fold limbs (Figure 9C). The distribution of elastic strain and Mises stress in the B₅ model is the same as in the B₄ model (Figures 9C, D). Instead of a full wavelength, only 3/4 of a wavelength is created. Its amplitude is also increased compared to the B₄ model (Figures 9C, D).

Investigation slip variations in multi-layer models showed no interlayer slip in the B₃ elastic-plastic model. In the B₄ and B₅ models, like the A₉ and A₁₀ multi-layer models, slip is zero in the hinge surface, and the models’ borders, and it towards limbs becomes maximum (Figures 7, 9B–D–D).

3.4 Pure buckling folding in single-layer and multi-layer (series C)

Series C consists of five models. The C₁ to C₄ models simulate the buckling of layers 2) to 5), respectively (Figure 10). The C₅ model also simulates multi-layer buckling non-slip (Figure 11). In addition, the solution method in this series is different; here, the buckling equation is used instead of the standard Mohr-Coulomb criterion (Eq. 1). The boundary conditions of series D models are the same as series B models (Figures 2D, E).

FIGURE 10. (A) distribution of elastic strain in the six buckling modes of the second layer, (B) distribution of Mises stress in the six buckling modes of the second layer, (C) distribution of elastic strain in the sixth mode of buckling for third layers (C₂), fourth (C₃) and fifth (C₄), (D) distribution of Mises stress in the sixth mode of buckling for third layers (C₂), fourth (C₃) and fifth (C₄).

FIGURE 11. (A) distribution of elastic strain in the six buckling modes in a multi-layer (C₅), (B) distribution of Mises stress in the six buckling modes in a multi-layer (C₅). Layers one to five are shown in (A). The red arrows show the divergence of the boundary of the deformed squares in the two adjacent anticlines is equal. The pink arrow shows the half-wavelength from which the strain ellipse presented in Figure 13 was extracted.

Six buckling modes are shown for model C₁ (Figures 10A, B), and only the sixth mode of buckling is shown for models C₂ to C₄ (Figures 10C, D). Each buckling mode produces 1/4 of the fold wavelength, so the fourth mode creates a full fold wavelength, and the sixth mode produces 1.5 of the fold wavelength. In all modes of all models, the elastic strain is maximum in the outer arch of the hinges zone in both anticlines and synclines (Figures 10A, B). Inversely, it decreases towards the centre of the limbs (Figures 10A, C). Mises stress is also maximum in the hinges zone (both anticline and synclines) and is minimum in limbs (Figures 10B, D). The created folds are symmetrical and uniform. However, the layers have been thinned in the hinge zones and thickened in the limb zones; these changes in thickness from the first to the sixth mode of buckling intensify (Figure 10). The C₅ model shows the six modes of pure multi-layer buckling (non-slip). The elastic strain is localized in the outer arch of the hinges zone (Figure 11A). Mises stress is concentrated in the outer arc and the mid-layer inflection zone (Figure 11B).

As expected, non-slip is observed at the boundary of the layers of this model (Figure 11). In other words, the primary squares do not slide at the boundary of the layers, and folding is created by stretching the squares in the outer arc and their compression in the inner arc. As a result, the boundary of the squares creates an outward-diverging and inward-converging pattern. As shown in Figure 11A, the degree of their divergence in the outer arc of the adjacent anticlines is the same (Figure 11A).

3.5 Buckling and bending folding non-sliding (series D)

Series D consists of four models. In fact, they are multi-layer models that investigate the role of buckling and bending non-slip. The D₁ and D₂ models simulate bending at three-point and are elastic and elastic-plastic, respectively (Table 1; Figures 12A, C). Their boundary conditions are similar to Figure 2B. The F₃ model simulates bending at four-point and has elastic rheology (Figure 12B). The D₄ model also models multi-layer buckling non-slip (Figure 12D); its boundary conditions are shown in Figure 2E The difference between the D₄ and C₅ models, which simulate multi-layer buckling non-slip, is in the used equation. The C₅ model uses the buckling relationship, and the D₄ model uses the Hook relationship.

FIGURE 12. Distribution of elastic strain and Mises stress in non-slip multi-layer bending elastic models at three points (A) and four points (B, C) distribution of Mises stress and elastic and plastic strains in non-slip multi-layer bending elastic-plastic model at three points, and (D) distribution of elastic strain and Mises stress in non-slip multi-layer elastic buckling model.

A half-wavelength fold is created in D₁ and D₃ models, as in models A₉ and A₁₁ (Figures 7, 8, 12A, B). The elastic strain is localized in the hinge of the anticline (Figures 12A, B). The Mises stress is also localized in the hinge of the fourth layer (Figures 12A, B). No interlayer slip is observed throughout the two models. Bending occurs by extension of squares in the outer arc and compression in the inner arc. The squares’ boundary shows a weak divergence towards the outer arc (Figures 12A, B). The D₂ elastic-plastic model stops after a few meters of displacement (Figure 12C); its stress and strain distribution are the same as the early stages of deformation of the D₁ model (Figures 12A, C).

In the D₄ model, no fold is created; only the thickness of the layers increases as in pure shear (Figure 12D). The elastic strain and the Mises stress are concentrated in the second and third layers, respectively (Figure 12D). In this model, no interlayer slip is observed. Only the primary squares are compressed in the displacement direction (Figure 12D).

4 Discussion

This research investigates the roles of sliding parallel to the layering (or along faults), buckling, and bending mechanisms in the formation of fault-related folds through the use of five series of 2D finite element models, encompassing a total of 29 models. Below, we compare our findings ith those of prior studies.

In series A, which simulates single-layer and multi-layer bending folding, the single-layer elastic-plastic models stop after applying a few meters of displacement (Figures 6B, D, F, H). In contrast, the multi-layer sliding model (A₁₀) is run until the end of applied displacement (Figure 7B). This indicates the effect of sliding parallel to the layers in the elastic-plastic model of bending folding. The results are confirmed by stopping the multi-layer elastic-plastic no-slip model D₂ that is presented to evaluate the results (Figure 12C). On the other hand, single-layer, and multi-layer elastic models are run until the end of the displacement (Figures 6A, C, E, H, 7A). Since the multi-layer elastic model has the interlayer slip (Figure 7A). We ran the D₁ model to investigate the role of slip in elastic multi-layer bending (Figure 12A). This model is run until the end of the applied displacement (Figure 12A). Therefore, multi-layer bending elastic models are implemented in non-sliding and sliding conditions (Figures 7A, 12A). This issue shows the influence of rheology in bending. Therefore, the results of series A show that two parameters of rheology and sliding parallel to the layers are effective in bending folding.

Series B simulates buckling in single-layer and multi-layer. Any fold is not created in its single-layer models (elastic and elastic-plastic) and multi-layer elastic-elastic models (Figures 9A, B). Only its multi-layer elastic models can simulate buckling folding (Figures 9C, D). This issue shows the influence of rheology in buckling. To investigate the role of sliding parallel to the layer to buckling, the D₄ model is run, a multi-layer elastic non-slipping model. As shown in Figure 12D, no fold is developed in this model. Hence, interlayer slip and rheology are controlling parameters buckling.

Earlier, Ramberg. (1964) investigated the effect of rheology characteristics on buckling. He found that competent materials responded better to buckling than incompetent materials. Buckling in competent materials creates a fold train, while incompetent materials only shorten and do not produce a fold train. This is consistent with the results of this research’s elastic and elastic-plastic models. This study’s multi-layer elastic buckling models (B₄ and B₅) generated a train of folds, while the multi-layer elastic-plastic buckling model (B₃) only shortened and did not create any fold (Figures 9B–D). Therefore, sliding parallel to the layer and materials’ rheology characteristics play an essential role in the development of buckling folds.

4.1 Evaluation of strain ellipse of models

To compare the deformation pattern in this research’s models with each other and the models presented in previous works (Ramsay, 1967; Ramsay and Huber, 1987; Twiss and Moores, 1992), we manually extracted the strain ellipse of the multi-layer models of this research (Figure 13), based on the deformation of the primary squares, as shown in Figure 14 for the detachment fold model. The strain ellipse in the multi-layer bending-sliding models of the A series shows an all-over stretching in the entire model, which is increased in the contact surface of the rigid circular bodes with layers, even in the inner arc of the fold, the strain ellipse stretching is observed, and this tension in the bending model at four points is larger than the bending model at three points (Figure 13). While in the multi-layer non-slip bending models of series D, the strain ellipse is stretched in the outer arc and compressed in the inner arc (Figure 13). Of course, an amount of turbulence can be seen in the contact surface of the rigid circular bodes with layers (Figure 13).

FIGURE 13. The strain ellipses of this study’s models are estimated based on the deformation of the initial squares. The half wavelengths from which these strain ellipses were extracted are marked with arrows in the previous figures.

FIGURE 14. Overview of strain ellipsoid extraction from the deformation of square meshes in a detachment fold model, as follows: (A) Initial square mesh setup, (B) Progressive deformation illustrating strain ellipsoids, and (C) Enlarged detail of a fully developed strain ellipsoid within the fold apex.

In the multi-layer buckling models of series B (B₄ and B₅), the strain ellipse is stretched in the outer arc and compressed in the inner arc. In the non-slip buckling model of series D (D₄), no folding occures, and like pure shear, the strain ellipse is compressed throughout the model (Figure 13).

In the buckling model of series C (C₅) that use the buckling relationship, the strain ellipse is strongly stretched in the outer arc and strongly compressed in the inner arc (Figure 13). In other words, except for the A series in all multi-layer models, the strain ellipse is stretched in the outer arc and compressed in the inner arc. However, this issue is more significant in theC₅ model (Figure 13). The strain ellipse of this model, is the same as the buckling multi-layer folding model presented by Ramsay. (1967) (Figures 13, 15A).

FIGURE 15. (A) Strain ellipse in the multi-layer buckling fold model (Ramsay, 1967), (B) slip changes in the flexural-slip fold model, and the size of the arrows shows the amount of slip (Twiss and Moores, 1992).

The strain ellipse in the detachment fold model is similar to (D₅) buckling model (Figure 13). Except for the second layer of the detachment fold, which is deformed by the movement of material towards its centre and has not experienced buckling, The strain ellipse of the third to fifth layers is the same as the buckling model of series C (C₅), whose strain ellipse is stretched in the outer arc and compressed in the inner arc. Perhaps it is reasonable to expect a bending mechanism from a detachment fold created by the upward pressure of its underlying layer, similar to bending sedimentary layers by the upward pressure of a magmatic mass. Contrary to expectation, the strain ellipse pattern in the upper layers of the detachment fold model is the same as the strain ellipse pattern in buckling (Figure 13). Of course, this is consistent with Butler et al.'s opinion (2020) about detachment folds. They believe that in addition to slipping parallel to the layer, the buckling mechanism also plays a critical role in developing detachment folds. The kinematic models presented for this folding ignore this issue.

In the fault-propagation and fault-bend fold models, the strain ellipses compress on the fault surface and forelimb while stretched on the back-limb (Figure 13). In fact, it seems that the presence of the fault has changed the place of compression from the centre of the fold to the fault tip and forelimb and the place of stretching from the crest to the back-limb (Figure 13). Therefore, in the geomechanical model of the fault-propagation and fault-bend folds, the effects of buckling can be observed in limbs of the folds. Of course, Butler et al. (2020) did not consider buckling influential in developing fault-bend folds. Still, the modelling results of this study showed that the buckling mechanism is also effective in the development of fault-propagation and fault-bend folds. However, its role is not as bold and significant as detachment folds.

In the multi-layer bending models at three-point of series A, and multi-layer buckling models of series B, the slip is zero at the model terminals and hinge surfaces, and it is maximized in the inflexion point (Figures 7, 9C, D). This is similar to the pattern of slip changes in the flexure slip folding mechanism presented by Twiss and Moores. (1992) (Figure 15B). While in the three fault-related fold models, the interlayer slip increases towards the model’s centre (hinge zone), inversely, the fault slip is reduced towards the hinge zone. In fact, in the model’s centre, the reduction of fault slip is compensated by the increase of interlayer slip, thickening of the layers, and buckling (Figures 3–5). Therefore, the main factor of deformation in the fault-related folds is the presence of the fault. The reduction of fault slip causes buckling in the hangingwall of the fault-propagation and fault-bend folds and the upper layers of the detachment layer in the detachment fold. Thus, variations in fault slip, slip parallel to the layer, and buckling are the controlling parameters of the development of fault-related folds.

On the other hand, the results of geomechanical models show that all three types of fault-related folds are created due to the reduction of fault, and their difference is in the place where the slip reaches zero. The slip reaches zero in the lower flat, or the detachment fault tip in the detachment folds. While in the fault-propagation and fault-bend folds at the ramp tip and upper flat, respectively, the slip reaches zero. Slip reduction in these places is compensated by increasing interlayer slip, increasing the thickness of layers, and buckling the upper layers.

5 Conclusion

In the presented research, the role of sliding mechanisms parallel to layering, buckling, and bending in the formation of fault-related folds were comprehensively investigated using five distinct 2D elastic and elastic-plastic finite element models. The key findings from our analysis are summarized as follows.

1) The rheological properties of the material and sliding parallel to layering are the dominant factors influencing bending in both single-layer and multi-layer models. Single-layer elastic-plastic models are not capable of folding via the bending mechanism. However, multi-layer elastic-plastic models that allow slip parallel to the layering do exhibit folding by bending. Conversely, where interlayer slip is absent (non-slip), these models do not exhibit any folding. Both single-layer and multi-layer elastic models fold by bending, regardless of whether sliding parallel to the layers is present or not.

2) Buckling is dependent on the rheological characteristics of the material and sliding parallel to the layer. The single-layer and multi-layer elastic-plastic models do not produce any fold by the buckling mechanism, Upside-down, only multi-layer elastic models are folded by the buckling mechanism at sliding conditions, and no folds are formed at no-slip conditions.

3) Variations in interlayer slip in bending and buckling models are the same as the pattern of slip changes in the flexure slip folding, where slip is absent at the hinge surfaces, then increases towards the limbs, and becomes maximum at the inflection points.

4) The strain ellipse and slip variation pattern in fault-related fold models show that sliding parallel to the layer and fault and buckling are two key mechanisms of fault-related folding. The reduction in fault slip creates all three types of fault-related folds; their difference is where this reduction occurs. The slip diminishes to zero at the detachment fault tip, the ramp tip, and the upper flat in the detachment, fault-propagation, and fault-bend folds, respectively.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

AK-S: Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Software, Visualization, Writing–original draft, Writing–review and editing. SA: Conceptualization, Methodology, Project administration, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing. MRG: Conceptualization, Methodology, Project administration, Supervision, Validation, Visualization, Writing–review and editing. MG: Methodology, Project administration, Validation, Visualization, Writing–review and editing. RD: Funding acquisition, Investigation, Resources, Writing–review and editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References