## Abstract

We study high-amplitude folding in layered rocks with two-dimensional numerical simulations. We employ the finite-element method to model shortening of an incompressible multi-layer with power-law viscous rheology. The Lagrangian numerical mesh is deformed and re-meshed to accurately follow the layer interfaces. Three settings are considered: (i) pure shearing of a confined multi-layer, (ii) simple shearing of a multi-layer above a detachment, and (iii) slump folding owing to gravity sliding. In our pure shear simulations, finite-amplitude folds always develop despite confinement and thin weak interlayers. The fold shapes can be significantly irregular, resulting from initial geometrical heterogeneities that are perturbations of the layer interfaces and differences in layer thickness. The bulk normal viscosity of the multi-layer decreases significantly with progressive folding. This structural softening decreases the bulk normal viscosities by a factor of 2–20. For simple shear, the multi-layer does not develop asymmetric fold shapes significantly. Fold axial planes in the multi-layer are mostly curved and not parallel. For slump folding, fold shapes can be significantly asymmetric exhibiting strongly curved fold axial planes and overturned fold limbs. The rheology of the competent layers has a major impact on the fold shapes for gravity-driven multi-layer folding.

## 1. Introduction

Rocks are often mechanically layered and fold when compressed in a direction sub-parallel to the layering. Folds are common structures in deformed rock and range in scale from sub-millimetre to hundreds of kilometres [1–3]. Understanding the mechanics of folding is essential for a better comprehension of tectonic processes, for example mountain building. The mechanics of rock folding have been studied with mathematical methods for more than 100 years [4]. Mathematical models for folding are usually based on the concept of continuum mechanics [5,6], and among the existing mathematical models one may distinguish two general approaches.

The first approach is based on the theory of stability of structures or elastic stability [7–9]. In this approach, the simplified so-called thin-beam equations are used to study two-dimensional folding (or thin-plate equations to study three-dimensional folding). These equations are based on the assumption that the thickness of the folding layer is significantly smaller than its lateral dimension. The main advantage of the thin-beam equation is that the mathematical description of a two-dimensional plane strain folding process is reduced to a one-dimensional mathematical model (correspondingly, the thin-plate equation is a two-dimensional mathematical model describing a three-dimensional folding process). The main disadvantage is that the treatment of shear deformation within the layer is only basic, which yields inaccurate solutions especially when the viscosity ratio (or any other ratio of material properties describing the mechanical strength, e.g. Young's modulus) between the layer and the embedding matrix is small [10]. The first study applying the equation for a thin beam on a foundation to folding of rock is presumed to be by Smoluchowski [4]. In his analysis, the foundation term represents the gravity force acting against folding of the Earth's crust. Later, it was Biot [11–13], in particular, who applied the elastic stability approach to rock folding. Among other things, Biot provided thin-beam solutions for a half-space foundation term (representing the resistance of the medium embedding the layer, i.e. the matrix) and also solutions for beams and matrices exhibiting viscous and viscoelastic rheologies. Biot derived the dominant wavelength theory [11,12], which is one of the most important results for understanding rock folding. This theory can explain several fundamental features frequently observed in single-layer folds, such as the proportionality between fold arc length (i.e. distance between neighbouring fold hinges measured along the folded layer) and layer thickness. Many studies elaborate on the linear thin-beam solutions of Biot by introducing nonlinearities in the thin-beam equation. Two types of nonlinearities can be distinguished, namely geometric and material nonlinearities. Solutions considering geometric nonlinearities mainly focus on folding up to high amplitudes and limb dips [14–16]. These nonlinear finite-amplitude solutions show that, for folding with velocity boundary conditions, the fold amplification rate decreases with increasing fold amplitude, which means that the folding instability becomes weaker with progressive folding [15–17]. The fold amplification rate is constant in the linear Biot solution [12]. The second type of nonlinearities, material ones, has been mainly analysed to investigate non-periodic fold geometries and localization during folding. Analytical solutions show that material softening can cause localization [18] and several studies demonstrate that material nonlinearities, in particular due to material softening, can be a possible explanation for strongly localized and irregular fold geometries [19–21]. Furthermore, some studies applied so-called shear-deformable beam equations (e.g. Timoshenko beam) to study folding. These theories, sometimes referred to as thick-plate theory, provide more accurate (but still not exact) results for the shear deformation within the layer [14]. Thin-beam equations have also been applied to study folding in multi-layers [12,22,23] and thin-plate equations to study folding in three dimensions [24,25].

The second approach to mathematically study rock folding is based on the theory of fluid dynamics, and the first studies applying this theory to rock folding are presumed to be by Ramberg [26,27]. The fluid dynamics approach has been combined with the hydrodynamic stability theory and applied to folding by Fletcher [10,28] and Smith [29,30]. In the fluid dynamics approach, the two-dimensional flow of a fluid layer can be accurately calculated from a so-called stream-function equation, which is a fourth order partial differential equation (so-called biharmonic equation) [31,32]. The fluid dynamics solution is not based on any beam or plate approximation and predicts the shear deformation within the layer accurately even for small viscosity ratios. The fluid dynamics solution is sometimes also referred to as the thick-plate solution [33]. The fluid dynamics solution was applied to folding of fluids with linear and nonlinear (power-law) viscous rheology. Strictly speaking, the analytical fluid dynamics solutions are valid only for infinitesimal amplitudes because they are solutions of linearized equations. However, numerical simulations [34] as well as laboratory experiments [35,36] have shown that these solutions are sufficiently accurate until fold limb dips between 10° and 20°. This is important because the wavelength selection process during folding takes place for limb dips smaller than 10–20° [28,35,37]. For limb dips smaller than 10°, layer shortening is dominated by uniform layer thickening but the initial (infinitesimal amplitude) irregular waviness in the layer shape is simultaneously selectively amplified into a more or less regular train of low limb-dip folds [38]. Between limb dips of 10° and 20°, a transition between uniform layer thickening and folding at more or less constant thickness takes place [38]. Around 20°, the hinge positions in the amplifying layer are ‘locked in’ and the fold arc length varies little with additional shortening [32,35]. The wavelength selection acting while limb dips are small is strongly sensitive to rheological properties of the layer and matrix and has a strong impact on the final fold geometry. Therefore, the fluid dynamics solutions predict the finite-amplitude fold shapes quite well although they are linear and, strictly speaking, valid only for infinitesimal amplitudes. Because the fluid dynamics solution is also exact for nonlinear rheologies, it can be used to validate and calibrate numerical algorithms (i.e. determine the necessary numerical resolution, the magnitude of the time step and the convergence criteria for iterations to numerically solve for nonlinear rheologies) used to simulate folding of viscous layers [39] or viscous flow in general. To study high-amplitude folding accurately, the linearized fluid dynamics solution can be extended to higher order [32], which basically corresponds then to a numerical spectral method [40,41]. Fluid dynamics solutions have been derived for folding of multi-layers [32,42–45] and for three-dimensional folding of linear viscous [46] and power-law viscous layers [47]. The fluid dynamics approach does not impose any constraints on the layer's aspect ratio (i.e. the ratio of length to thickness) and is therefore widely applicable. It has been applied to study ductile rock deformation, which generates structures other than folds, such as mullions (i.e. folding of a single interface between two materials of different properties) [30], pinch-and-swell structures on the small scale [29,48] and on the large scale (considering gravity) [49,50] and gneiss domes [51]. Furthermore, the fluid dynamics approach is also applied in geodynamics to study diapirism and thermal convection [52].

Both the abovementioned approaches are well suited to investigate folding and the mathematical models are useful because they establish a link between the observed fold shapes and the tectonic history [33] and enable the assessment of: (i) the bulk shortening in the folded rock [53], (ii) the effective viscosity ratio (or ratio of elastic properties) between different layers, which is important for assessing the rock rheology and the strength distribution in stratified rock [54], (iii) the folding mechanism, i.e. active folding (or buckling) or passive folding (e.g. formation of a monocline or bending owing to the movement of rock layers along non-planar surfaces such as fault bend folding) [2], (iv) the type of active folding mechanism, e.g. folding of consolidated rock owing to tectonic compression or folding of unconsolidated sediments owing to gravity sliding (i.e. slump folding), and (v) the bulk deformation style, such as pure shear or simple shear.

In nature, multi-layer folds (figure 1) are more frequent than single-layer folds, but researched significantly less than their single-layer equivalents. Here, we study the mechanics of two-dimensional, finite-amplitude multi-layer folding with the finite-element method. We use the fluid dynamics equations for slow viscous flow [31]. The low-amplitude stages of two-dimensional multi-layer folding have been studied analytically [32,42,44,45,55], and the finite-amplitude stages have mainly been investigated with laboratory experiments [56,57] and with numerical simulations [40,41,58]. While these simulations have been performed for linear viscous materials, rheological arguments as well as single-layer folding research [28,30,59] indicate that a power-law viscous flow law, employed in this study, would be more applicable. The power-law flow law is used here to describe the effective rheology of the folded layers without specification of a particular microscopic rock deformation mechanism. In general, the power-law flow law can be applied to study rock deformation on many scales as it effectively describes a wide variety of deformation mechanisms: (i) diffusion creep for *n*=1 (*n* being the power-law stress exponent) [60], (ii) the pressure solution for *n*=1 [61], (iii) dislocation creep for 3<*n*<5 [60], (iv) the so-called exponential flow for high stress levels where flow laws are determined by fitting laboratory creep experiments usually for *n*>7 (a satisfactory theoretical microscale model is lacking; also, exponential flow laws have been suggested for the pressure solution) [62–64], (v) low-temperature plasticity based on the Peierls mechanism [63,65], (vi) von Mises plastic deformation for *n*≫10 [49], (vii) viscosity strain softening [66], (viii) micro-scale brittle deformation resulting in an effective macroscale ductile deformation for 5<*n*<15 [67], (ix) orthotropic linear viscous fluids [28,39], and (x) the effective rheology of large-scale brittle–ductile tectonic units such as the entire lithosphere [68].

We study multi-layer folding for three different settings: (i) bulk pure shear of a confined multi-layer, (ii) bulk simple shear of a multi-layer above a detachment, and (iii) gravity sliding of a multi-layer. Finite-amplitude multi-layer folding for these three settings has been little studied numerically; however, setting (i) has been investigated analytically for linear viscous rheology [55,69] and experimentally [42] and setting (iii) has been studied experimentally [70]. All of these settings occur frequently in nature. Applications of (i) and (ii) can be found in figure 1, while an example for (iii) is gravity sliding of sediments in foreland basins [71] and on passive margins [70,72]. We focus here on describing the geometry development of multi-layer folding in these three settings and analysing the variation of effective viscosity owing to the development of structures and changes in the strain-rate-sensitive rheological parameters.

## 2. Mathematical model

The force balance equations for slow, steady-state fluid dynamics (no inertial terms) are [31]
2.1
where *x* and *y* are the horizontal and vertical Cartesian coordinates, respectively; *σ*_{xx}, *σ*_{yy} and *σ*_{xy} are the three total stresses, *ρ* is the density and *g* is the gravitational acceleration. The constitutive equations for the incompressible power-law viscous rheology are given by [31]
2.2
where *τ*_{xx}, *τ*_{yy} and *τ*_{xy} are the three deviatoric stresses, *p* is the pressure, *v*_{x} and *v*_{y} are the horizontal and vertical velocities, respectively, and *η*_{e} is the effective viscosity given by [32]
2.3
where *η*_{0} is a reference viscosity, *D* is a reference deformation rate and *E* is the square root of the second invariant of the rate of deformation tensor and is given by
2.4
For *n*=1, the rheology is linear viscous, i.e. Newtonian; for *n*≠1, the rheology is nonlinear, i.e. power-law viscous. While shear thickening fluids (*n*<1) exist, rock usually exhibits shear thinning behaviour (*n*>1) [73]. The reference strain rate, *D*, could be included in the definition of *η*_{0}. However, *D* is introduced because (i) *η*_{0} is then a normal reference viscosity in terms of units (i.e. Pa s) and can be used to define a reference viscosity ratio, *R*, between strong and weak layers, and (ii) one can better distinguish if variations in the ratio of *η*_{e} between strong and weak layers are due to variations in material properties (i.e. *η*_{0}) or due to variations in the flow field (i.e. *E*/*D*≠1). Consequently, *D* is set to the corresponding value of *E* for bulk pure or simple shear, depending on the setting. If flow perturbations lead to *E*>*D*, then *η*_{e} decreases with increasing strain rate with respect to the reference viscosity (i.e. viscosity strain rate softening), but the stress is still increasing with increasing strain rate (i.e. strain rate hardening). If *E*<*D*, then *η*_{e} increases with decreasing strain rate with respect to the reference viscosity. In our simulations, the applied power-law viscous flow law does not describe a ‘true’ strain rate softening for which the stress decreases with increasing strain rate [74]. Such strain rate softening could also be investigated with power-law viscous flow laws by setting *n*<0 [75]. The nonlinear power-law rheology is notable only if the deformation inside the model domain deviates from the applied bulk deformation and power-law effects do not exist if the strain rate in the model domain is constant. This is not the case for the simulations of gravity sliding because, for these simulations, the bulk deformation exhibits a heterogeneous strain rate field that also varies with progressive deformation. For gravity sliding, we use a value of *D* that corresponds roughly to the average value of *E* in the area where the folds later develop. However, values of *E* vary by up to two orders of magnitude even for the first time step and, therefore, the choice of *D* has a considerable impact on the final fold shapes for given parameters. Alternatively, one could use experimentally determined rheological parameters, but such parameters are either non-existent or highly uncertain for weakly consolidated sedimentary rocks such as sandstone, limestone and shale that are usually involved in slump folding.

We solve the above equations with the finite-element method [76–78] using a mixed formulation with Crouzeix–Raviart triangles with quadratic velocity shape functions enhanced by a cubic bubble function and discontinuous linear interpolation for the pressure field (so-called P2^{+}–P1 element). We consider incompressible fluids and apply the penalty method in combination with Uzawa iterations to enforce a divergence-free velocity. The nonlinearity owing to the power-law rheology is treated with Picard iterations. A Lagrangian approach with re-meshing is used to handle the large strains. The model geometry is described through external and internal interfaces with up to 2000 nodes on a single-layer interface. The mesh is generated with Triangle [79]. Both the finite-element mesh and the interfaces are moved with the resulting velocity fields. Once the quality of the finite-element mesh becomes unsatisfactory, it is regenerated using the evolved model geometry. Typical simulations require between 200 and 600 time steps. The applied finite-element algorithm has been described previously [48,80] and relies on the fast finite-element method solver Milamin [81].

## 3. Numerical simulations

### (a) Pure shear

In the first set of numerical simulations, all model boundaries are free slip and kept straight, and moved with pure shear velocities corresponding to a constant rate of deformation (figure 2). The horizontal velocities at the lateral boundaries are therefore decreased after every time step to maintain a constant ratio of shortening velocity to model width. The reference viscosity ratio, *R* (i.e. the ratio of *η*_{0} between strong and weak layers), is varied between 25 and 100 and the power-law stress exponent of the strong layers, *n*, and weak layers, *n*_{1}, varies between 1 and 5. The strong and weak layers have initially the same thickness for most simulations. Additionally, two simulations have been performed in which the initial thickness of the weak layers is one-half and one-tenth of the initial thickness of the strong layers, respectively, and one simulation has been performed in which the thickness of each strong and weak layer is different. Another two simulations have been performed in which the top and bottom weak layer are initially 100 times thicker than the strong layers to represent a half-space above and below the multi-layer.

Two initial geometrical perturbations have been applied. The first is a bell-shaped perturbation only in the bottom layer given by [82]
3.5
where *a*=5*H* and *b*=0.1*H* with *H* being the initial thickness of individual strong layers. Folding of single layers with such an initial perturbation has been studied analytically [82], experimentally [83] and numerically [84], and the finite-amplitude evolution for the applied velocity boundary conditions can also be investigated with analytical finite-amplitude solutions [15,16]. For an easier later discussion of the multi-layer results, we shortly repeat and visualize the fundamental impact of this bell-shaped initial perturbation on single-layer folding with two numerical simulations (figure 2); one simulation for *R*=25, *n*=1 and *n*_{1}=1 and another one for *R*=25, *n*=3 and *n*_{1}=3. For both simulations, the ratio of initial layer length to thickness is 200 and the matrix above and below the layer is five times thicker than the layer. Laboratory experiments [83,85], numerical simulations [84] and results of fluid dynamics analytical solutions [36] for folding with a bell-shaped (or similar) initial perturbation agree well and predict lateral fold propagation (or serial folding or cellular buckling) in viscous layers. The fold at the position of the initial perturbation grows first while folds to the left and right start amplifying later sequentially. The fold in the middle stops growing at large amplitude stages (i.e. the fold shape locks) and the other folds grow to more or less the same amplitude. Final fold amplitudes are more or less similar; however, the thickness of folds closer to the lateral model boundaries gets larger because these folds remained in the phase of uniform layer thickening longer than the fold in the middle, which amplified from the onset of bulk shortening. For 60 per cent bulk shortening, the simulation with a linear viscous flow law shows localized fold geometry while the simulation with a power-law viscous flow law shows a more or less regular fold geometry (exhibiting localized fold geometry only for less than 40% bulk shortening). The reason is that fold amplification rates in the simulation with power-law viscous rheology are faster than for a linear viscous rheology and, therefore, folds also propagate faster laterally with respect to the bulk shortening rate. The simulations demonstrate a well-known feature of viscous single-layer folding [84,85], namely that, for both linear and power-law viscous flow laws and for velocity boundary conditions, localized fold geometries can develop from a localized initial geometrical perturbation.

The second type of initial perturbation is a random perturbation of all internal interfaces with a red noise distribution approximating natural rough interfaces [86]. The perturbations on the top and bottom layer interfaces are different, which causes lateral variations in the layer thickness. Single-layer folding in viscous and viscoelastic layers with random initial perturbations has been studied intensely numerically for two-dimensional [34,53,84,86,87] and three-dimensional [17,88] folding. The largest initial amplitude of these perturbations is always set to be equal to one-twentieth of the initial layer thickness. For all pure shear simulations, the initial ratio of length to thickness of individual strong layers is 200 and all simulations were performed with 15 strong and 16 weak layers (except for one simulation; see below).

Figure 3 shows folding of a multi-layer with an initial bell-shaped perturbation only in the bottom layer. Folds first propagate vertically from the perturbation in the bottom layer across the entire multi-layer and afterwards propagate laterally to the lateral boundaries of the model domain. This fold propagation (or serial folding) explains the larger fold amplitudes in the horizontal middle of the model domain compared with amplitudes at the left and right boundary. Also, fold amplitudes in the vertical middle are largest and smallest in the top and bottom strong layers because the confined model boundaries inhibit fold amplification rates [32,44,89].

Figure 4 displays folding of a multi-layer with initial random perturbations. In this case, all folds grow more or less simultaneously. The amplitudes and the ratios of arc length to thickness of individual folds vary widely and fold shapes are considerably irregular within individual layers. The strong irregularity makes the identification of a dominant wavelength, as predicted by analytical solutions [55], challenging (see §4). The amplitude of individual folds is generally largest in the middle (in vertical direction) of the model, whereas amplitudes are smallest in the top and bottom strong layers because of the confined boundaries. Such irregular, non-periodic fold geometry is termed ‘localized’ by Hobbs *et al.* ([74], fig. 33).

Figure 5 illustrates the results of four different simulations after 52 per cent bulk shortening. Figure 5*a*,*b* are for *R*=25 and *n*=3 and differ only in the value of *n*_{1}. The fold amplitudes and shapes are considerably different, indicating the control of *n*_{1} on the multi-layer fold evolution. Figure 5*c*,*d* are for *R*=100; however, figure 5*c* shows results for linear viscous layers, whereas figure 5*d* shows results for power-law viscous layers. Between these, fold amplitudes and shapes vary considerably. Therefore, even small values of *n* and *n*_{1} around 3 are sufficient to modify the fold evolution in multi-layer folding significantly compared with linear viscous multi-layer folding. The fold amplification rate is considerably larger for *R*=100 than for *R*=25. Therefore, the lateral fold propagation rate is also faster and folds reach the lateral model boundaries within a smaller amount of bulk shortening. After 52 per cent bulk shortening, folds have already been propagated to the lateral model boundaries for *R*=100, while, for *R*=25, the layers at the lateral boundaries are still more or less horizontal. This is in agreement with results of single-layer fold propagation (figure 2).

Figure 6 shows two results for a multi-layer with initial random perturbations where the weak layers are thinner than the strong layers. In both cases, high-amplitude folds develop in the vertical middle of the model. Fold hinges are generally less curved than for the simulations with thicker weak layers and fold shapes are more similar to chevron folds and kink bands. Despite the confinement and thin weak layers, high-amplitude folds developed in the middle of the model for moderate values of *R*, *n* and bulk shortening.

Figure 7 shows folding of 11 strong layers with different initial thickness, spacing and with initial random perturbations (figure 7*a*). The largest initial amplitude of these perturbations for all layers is always set to be equal to one-twentieth of the dimensionless layer thickness of 1. Values are *R*=25, *n*=5 and *n*_{1}=3. The non-dimensional thickness of the strong layers is, from bottom to top, 1, 0.4, 0.3, 0.7, 1.2, 0.8, 0.6, 0.4, 0.6, 0.4 and 0.3, and the non-dimensional thickness of the weak layers is, from bottom to top, 1, 0.5, 1, 1.5, 0.5, 0.75, 0.6, 2, 0.5, 0.4, 0.7 and 1. After 26 per cent shortening, some of the amplitudes in the multi-layer are considerably larger than others (figure 7*b*). Also, the amplitudes and limb dips of the thinner layers are larger than those of the thicker layers because the initial ratio of amplitude to thickness is larger in the thinner than in the thicker layers. Therefore, the initial phase of uniform layer thickening takes place within a smaller amount of bulk shortening in the thinner layers. In particular, the thin layers developed many individual folds with a more or less regular arc length; however, the thin layers are also deformed by growing folds in the thicker layers. After 51 per cent shortening, the fold shapes are considerably irregular (figure 7*c*). The thin layers developed larger scale folds owing to the influence of folding of the thicker layers. The orientation of fold axial planes of individual folds varies significantly and shapes of some individual folds are strongly asymmetric. Some of the individual folds can also be described as box-folds, i.e. between two limbs there is a more or less flat segment bounded by two hinges of the same type (e.g. both convex-upward).

The bulk normal viscosity, *η*_{N}, of the entire model can be calculated by dividing the average horizontal deviatoric stress, *τ*_{xx}, by the applied bulk pure shear shortening strain rate. Figure 8*a* illustrates the evolution of *η*_{N} with increasing bulk shortening. For all performed pure shear simulations, *η*_{N} decreases with increasing shortening after a certain amount of shortening. Such a decrease in *η*_{N} has already been demonstrated analytically for linear viscous single layers [90] and is referred to as structural softening during folding [91]. Structural softening results from the formation of structures, here folds, inside the numerical domain that remains rectangular. Depending on the configuration, the final value of *η*_{N} is between 2 and 20 times smaller than the initial value. We observe that (i) multi-layers with power-law viscous rheology start softening after a smaller amount of bulk shortening than the corresponding multi-layers with linear viscous rheology, (ii) multi-layers with larger values of *R* soften more than the corresponding multi-layers with smaller values of *R*, and (iii) the initial perturbation has no impact on the maximum softening but influences when softening starts with respect to bulk shortening. Figure 8*b* shows the evolution of *η*_{N} and of the area-averaged value of *η*_{e}, i.e. *η*_{ea}, with increasing shortening. While *η*_{N} is decreasing, *η*_{ea} is increasing. The reason is that strain rates in the fold limbs of the strong layers decrease and, therefore, the values of *η*_{ea} increase. However, as fold limbs mainly rotate without significant internal deformation, the overall resistance against shortening decreases, although values of *η*_{e} increase in the fold limbs (figure 9). If only the weak layers are power-law viscous, then *η*_{N} decreases, whereas *η*_{ea} stays more or less constant with shortening. This indicates that both *η*_{ea} and *η*_{N} are controlled by the strong layers.

Figure 9 shows the distribution of *τ*_{xx} (figure 9*a*) and *η*_{e} (figure 9*b*), for a bulk shortening of 50 per cent for the multi-layer simulation with an initial random perturbation and *R*=25, *n*=3, *n*_{1}=1. Values of *τ*_{xx} indicate extensive stresses in many of the more or less straight limbs and in the outer parts of the fold hinges. However, individual folds with smaller ratios of arc length to thickness still exhibit compressive values of *τ*_{xx} in their fold limbs. Values of *η*_{e} are significantly increased in the long, relatively straight fold limbs. The evolution of the corresponding value of *η*_{N} for this simulation is illustrated in figure 8*a*.

### (b) Simple shear

The left-hand, right-hand and inclined bottom model boundaries are kept straight in these simulations and the top boundary is a free surface (zero traction, figure 10). The right-hand vertical boundary is moved upwards with a constant velocity and the left-hand vertical boundary is kept fixed (no slip), resulting in a constant bulk simple shear strain rate. This simple shear rate is used to set the value for *D*, which is twice the vertical velocity at the right-hand model boundary divided by the model width (the factor of 2 is necessary so that the value of *E* is equal to *D* for homogeneous simple shear). The bottom boundary is initially oblique and deformed with the corresponding vertical simple shear velocities. The weak and strong layers have initially the same thickness and the uppermost weak layer is 150 times thicker than the other layers, representing a half-space above the multi-layer. The lowermost weak layer is five times thicker than the strong layers. In the lower half, the lowermost strong layer is shifted upwards by 50 per cent of its initial thickness and the corresponding thickness of the above weak layer is decreased accordingly. This step-like initial perturbation in the lowermost strong layer (inset in figure 10*a*) is asymmetric in order to favour the development of asymmetric fold shapes in this setting of multi-layer detachment folding under bulk simple shear.

Figure 10 shows the evolution of power-law viscous multi-layers for *R*=50, *n*=3 and *n*_{1}=1. Most developed fold shapes do not exhibit a significant asymmetry despite the simple shear and the asymmetric initial perturbation. For a shear strain, *γ* (i.e. vertical displacement at the right-hand boundary divided by the model width), of 1.2 individual folds on the two limbs of the larger scale fold in the middle of the multi-layer (inset in figure 10*c*) exhibit typical asymmetric shapes of so-called parasitic folds, i.e. s- and z-shaped [58]. Also, fold axial planes across the multi-layer folds are curved and not parallel. Individual folds in the middle of the multi-layer have higher amplitudes, straight limbs and sharp hinges while folds on the top of the multi-layer have smaller amplitudes and round shapes. Folds also develop at the two lateral boundaries of the model. This is a boundary effect because the layers are cut obliquely to their layering at the boundary to match the model boundaries. The application of a constant simple shear velocity leads to a (attempted) constant simple shear rate at the boundary. However, the layer geometry at the boundary is incompatible with a continuous tangential and normal stress across the layer interface, which is the reason for the deflection of the layer interface and folding at the boundary. Fold shapes for a simulation in a linear viscous multi-layer for *R*=50 (not shown here) are similar but with smaller amplitudes at the same amount of shear strain, *γ*, because the amplification rates in the linear viscous multi-layer are smaller.

Figure 11 documents results for three different simulations with only five layers for *γ*=1.8. Figure 11*a* shows results for *R*=50, *n*=5 and *n*_{1}=1. The multi-layer folds do not exhibit a significant asymmetry. However, the fold axial planes across the multi-layer are curved and not parallel, because the entire multi-layer in the middle of the model developed a larger scale fold. Figure 11*b* shows results for *R*=50, *n*=8 and *n*_{1}=3. The multi-layer folds are more localized around the initial perturbation and show asymmetric shapes. The strong layers are locally very close. Figure 11*c* shows results for *R*=50, *n*=5 and *n*_{1}=1 with a thick, strong (*R*=50) layer above the multi-layer. The multi-layer folds are more or less symmetric and fold axial planes are sub-parallel, because the overlying layer inhibits the formation of a larger scale fold such as the one developed in figure 11*a*.

### (c) Gravity sliding

The initial model geometry for gravity sliding has a slope in the left two-thirds of the model domain that becomes flat in the right one-third (figure 12). The lateral, vertical boundaries are kept straight with free slip, the bottom boundary is fixed with no slip and the top boundary is a free surface. Gravity acts vertically downward. The strong layers and the weak material in between have initially the same thickness. The strong layers do not touch the left-hand model boundary to allow downward sliding. The thickness of the weak material below the multi-layer stack is five times the initial thickness of the strong layers. The following parameters are used: the model width is 150 km; the thickness of the strong layers is 200 m and the lateral height difference is 21 km, corresponding to a slope of about 12°; the density for all layers is 2600 kg m^{−3}; the reference viscosity of the weak layers is 10^{20} Pa s and of the strong layers is 5×10^{21} Pa s; *g*=9.81 m s^{−2}; and *D*=1 Myr^{−1}.

Figure 12 displays the results for a linear viscous multi-layer with *R*=50. After reaching the gravitational equilibrium (top model boundary is more or less horizontal), the fold amplitudes are still relatively small. The folds develop to the right (i.e. down-slope direction) of the slope break. The multi-layer sequence exhibits a vergence towards the down-slope direction, i.e. the up-dip direction of the fold axial plane is towards the right. The folds are localized with respect to the horizontal direction.

Figure 13 shows results for the same model as in figure 12 in which the only difference is that the strong layers are power-law viscous with *n*=5. Fold amplitudes are significantly larger than for the linear viscous case and are localized without significant lateral propagation. Fold shapes are asymmetric, fold axial planes are strongly curved and exhibit a vergence towards the down-slope direction and individual fold limbs are even overturned. The results indicate that the activation of power-law viscous rheology in the strong layers has a significant impact on the fold amplification and fold shapes for folding driven by gravitational forces.

Figure 14 documents the results of three additional simulations for gravity sliding illustrating the variety of fold shapes for slump folding. For the simulation in figure 14*a*, the viscosity of the strong layers is 2.5×10^{21} Pa s and *n*=*n*_{1}=3. All other parameters are the same as for the simulations in figure 12. The initial ratio of *η*_{e} between the strong layers and the weak material at the position where the folds later develop is around 50. For the simulation in figure 14*b*, the viscosity of the strong layers is 1×10^{22} Pa s and both *n*, *n*_{1}=3 (all other parameters are also the same as for the simulations in figure 12). Here, the initial ratio of *η*_{e} between the strong layers and the weak material at the position where the folds later develop is around 200. The slope is halved in the simulation in figure 14*c*. The viscosity of the weak material is 1×10^{19} Pa s; the viscosity of the strong layers is 50×10^{19} Pa s, *n*=5 and *n*_{1}=1.5; the density of the strong layers is 2600 kg m^{−3} and of the weak material 2300 kg m^{−3}. All other parameters are the same as for the simulations in figure 12. The initial ratio of *η*_{e} between the strong layers and the weak material at the position where the folds later develop is around 40.

## 4. Discussion

All presented simulations for the three deformation settings exhibit significant growth of folds, even in cases of low viscosity ratios, confinement and thin weak interlayers. The reason for this is the increased growth rates in a multi-layer, compared with a single layer, as analytically determined [32,44,45,55]. The main difference between the simulations for pure and simple shear and the simulations for gravity sliding is the control of the layer shortening. Pure and simple shear are constrained by velocity boundary conditions that provide a constant rate of bulk shortening or shearing while gravity sliding is constrained by the gravity force that provides a constant body force but highly variable deformation rates in both space and time.

The fold shapes that developed during pure shearing in individual simulations exhibit ranges of arc length to thickness ratios, which makes it challenging to establish a link between the numerical multi-layer results and the single dominant wavelength yielded by analytical solutions [55]. The final fold shapes, in part, show significant variations in arc length and also amplitude of individual folds (figure 5). In our simulations, this variation can arise owing to fold propagation where folds do not amplify simultaneously but sequentially. Consequently, the phase of initial uniform layer thickening lasts longer in these areas of the layers that are furthest away from the initial, localized perturbation (figures 2 and 3). The process of fold propagation (or serial folding) has been suggested based on field observations [2] and has been studied with laboratory experiments [85,92] and mathematical models based on the elastic stability approach where fold propagation is also termed cellular buckling [22,93]. Furthermore, the irregularity in final fold shapes can be owing to random initial perturbation on all layer interfaces (figure 4). If the strong layers have different initial thickness and spacing, then the final fold shapes can be considerably irregular (figure 7). The multi-layer folds are frequently not periodic and may therefore be classified as localized according to the terminology suggested by Hobbs and co-workers ([74], fig. 33). The final multi-layer folds in figures 6*a* and 7 do not exhibit continuous vertical fold axial planes and, along some vertical lines, antiforms transform into synforms so that these fold shapes may even be classified as chaotic according to the study of Hobbs *et al.* [74]. These localized and chaotic fold shapes developed for power-law flow laws with viscosity ratios between 25 and 100 and standard power-law stress exponents between 1 and 5. The applied flow laws were always strain rate hardening and no additional strain or strain rate softening was applied. The irregular fold shapes are, therefore, a consequence of initial heterogeneities; here either geometrical perturbations in the layer interface or differences in initial layer thickness. Numerical simulations have shown that, for an initial localized geometrical perturbation, a multi-layer shortened above a detachment develops significantly more localized fold shapes for a power-law viscous rheology than for a linear viscous rheology [94]. Other heterogeneities, not considered in this study, are variations in material properties (e.g. *η*_{0}) within the layer, which could be caused by variations in chemical composition or grain size. One possibility to obtain information from irregular fold shapes on the initial amplification spectrum and the dominant wavelength would be to look at the bandwidth and distribution of arc length to thickness ratios, such as what Fletcher & Sherwin [38] did for single-layer folds.

The effective material behaviour of the multi-layer is rather complex despite the application of a simple power-law flow law with constant material parameters, i.e. *η*_{0} and *n*. We can make the following observations based on our pure shear simulations. Locally *η*_{e} varies owing to the differences in the rate of deformation. The area average of these viscosities, *η*_{ea}, shows, in the case of strong power-law layers, first a stiffening. The reason for this is that the strong layers are initially almost straight and parallel to the shortening direction, which means that the imposed rate of deformation prevails throughout the model. With the development of folds, the limbs become regions of essentially passive rotation and high *η*_{e}, while the hinges locally exhibit high rates of deformation and consequently low *η*_{e}. The weak power-law layers in between allow for slip and exhibit the highest deformation rates and lowest *η*_{e}. This has no substantial influence on *η*_{ea} as it is dominated by the strong layers. In the final stages of deformation, the limbs rotate into the field of extension and undergo stretching and thinning. However, depending on the multi-layer fold geometry, the values of *η*_{ea} can continue to increase or actually start to decrease (figure 8*b*). Parameter *η*_{ea} is misleading because the effective bulk rheology of the entire model is anisotropic. Given the initial layered configuration, the bulk normal viscosity *η*_{N} is highest here (Voigt bound); the corresponding effective shear viscosity would be lowest (Reuss bound) but is not considered here. The growth of folds strongly affects *η*_{N}. Owing to progressive deviation from the ideally layered configuration, *η*_{N} drops rapidly once the fold limbs start to rotate away from the shortening direction, which is opposite to *η*_{ea}.

The reduction of the bulk normal viscosity *η*_{N}, termed ‘structural softening’, is important during the deformation of layered rocks on many scales (figure 8). Structural softening results from the development of internal structures, here folds, inside the model domain, whereby the geometry of the model domain stays rectangular. Structural softening occurs when the internal structures are caused by instability, whereas no structural softening will occur if passive folds (no instability) develop, e.g. when different layers have the same material properties and the layer interfaces exhibit initial geometrical perturbations. The resulting passive folds will change the internal geometry of the numerical domain but will not cause a softening (this is why we do not use the term geometrical softening). All material properties remain constant during structural softening. The material properties applied here, i.e. *η*_{0} and *n*, are always constant during our simulations and the power-law viscous flow law is strain rate hardening, i.e. the stress always increases with increasing strain rate. Natural examples of structural softening include the formation of crenulation cleavage on the microscale [95], the formation of multi-layer folds in turbiditic sequences on the mesoscale [96] or the multi-layer folding of the entire lithosphere on the macroscale [97]. Numerical simulations of lithospheric shortening usually model lithospheric units, e.g. the upper crust, with isotropic flow laws for single minerals, e.g. quartz. As we demonstrate here this is not sufficient to capture the first-order behaviour. Crustal rocks are usually layered and develop folds on all scales during lithospheric shortening. Using an isotropic viscosity will take care of the area-averaged viscosity *η*_{ea}, which, however, can exhibit essentially opposite behaviour when compared with *η*_{N}. The presented numerical simulations are useful to calculate the strain-dependent average material properties of deforming anisotropic rock exhibiting nonlinear flow laws. A better understanding of these strain-dependent average properties of deforming anisotropic rock is essential for a better understanding of large-scale tectonic processes, such as crustal and lithospheric shortening, and these properties have been mainly studied with analytical methods [42,98]. Consequently, the evolution of *η*_{N} and the resulting structural softening and its feedback to the crustal and lithospheric shortening should be taken into account in further studies.

Fold asymmetries develop in our simulations, but, with the exception of gravity sliding in power-law materials, are not strongly accentuated and do not allow for inference of the dominant simple shear component. However, natural folds often show clear asymmetry development owing to simple shear such as that illustrated by the minor folds in the antiform of figure 1*c* that are due to flexural slip. One reason for the lack of significant asymmetry in the presented simple shear simulations is the constant thickness and homogeneous material properties in the strong layers. If significant lateral variations in thickness and material properties are present, asymmetric fold shapes are more likely to develop. In nature, such variations can also occur in the third dimension. The minor folded layers in figure 1*c* exhibit significant lateral variations in thickness. Asymmetric shapes of individual folds can be observed in figure 7 for bulk pure shearing with different thickness of the individual layers and in the inset in figure 10*c* on the limbs of the larger scale fold (parasitic s- and z-shaped folds).

The fold axial planes obtained are usually curved. This is important because this feature is frequently interpreted as evidence for multiple phases of tectonic deformation where a straight fold axial plane, developed during an earlier tectonic event, has been re-folded in a later tectonic event. However, especially in simulations for simple shear and gravity sliding, curved fold axial planes can easily develop during a single deformation event. Therefore, curved fold axial planes do not reliably indicate more than one tectonic deformation event and one must be cautious when inferring multiple tectonic events from observed fold shapes.

The simulations for gravity sliding are strongly sensitive to the power-law viscous rheology because the strain rates and therefore the values of *E* vary significantly in the model domain throughout the entire simulation (even for the first time step). This is the major difference when compared with the pure and simple shear setting, where values of *E* vary little in the model at the first time step and are close to the applied bulk pure or simple shear rate. Also, the density structure has a major impact on the fold development. In this study, only one density structure is considered and a detailed analysis of the impact of different densities on folding during gravity sliding is beyond the scope of this study.

In our simulations, the rocks are considered incompressible. However, in Nature, a significant amount of compressible deformation and also mass loss can occur, especially during folding in limestone and sandstone at shallow depth. The reason is that pressure solution could be the dominant deformation mechanism in the rock and a considerable amount of volume change may occur especially in the fold hinges. On the other hand, the solution and precipitation of material may occur locally on the small scale and the deformation of the entire fold on the large scale may be considered as incompressible. Also, brittle faulting can be active around the fold hinges. The simultaneous brittle fracturing in the extensional regions of fold hinges and fold amplification has been studied with a novel finite-element algorithm, as explained by Jager *et al*. [99]. However, a significant amount of research is still required to better understand the impact of compressibility, local volume change and brittle fracturing on rock folding.

From a technical point of view, several challenges must be handled to obtain the results presented here: (i) the sharp interfaces between the different phases must be adequately represented, (ii) the large strain evolution must be captured, and (iii) the strong variations in material properties must be resolved. The applied finite-element method with a deforming Lagrangian mesh in combination with re-meshing is well suited for this. Figure 15 shows a zoom into the finite-element mesh together with a colour plot of *η*_{e} for the simulation shown in figure 5*d*. The resolution of the mesh varies significantly and allows the thinning weak layers between fold limbs to be resolved. Values of *η*_{e} vary by more than five orders of magnitude but the sharp jump in *η*_{e} across layer interfaces is accurately resolved with the applied method. The large variation in values of *η*_{e} is mainly owing to the fact that the fold limbs rotate nearly passively and the deformation rate within the limbs decreases significantly causing a significant increase in *η*_{e}.

Mathematical models of rock folding provide a better understanding of the mechanical process of folding. They also allow for developing methods to estimate the amount of shortening and the material properties of folded rocks. Some methods have been developed and tested for single-layer folds [38,53,54,100]. For multi-layer folds, so-called dynamic retro-deformation (or reverse modelling or dynamic unfolding [101,102]) of observed fold profiles with numerical algorithms is a promising method to estimate and better constrain the deformation history of multi-layer folds.

## 5. Conclusions

The applied finite-element method using a deformable Lagrangian mesh with re-meshing is a suitable numerical method to simulate the large strain folding of a power-law viscous multi-layer exhibiting large effective viscosity ratios.

For bulk pure shearing, a multi-layer with power-law viscous flow law always develops finite-amplitude folds despite thin weak interlayers and rigid confinements for low-to-moderate viscosity ratios (*R*=25) and power-law stress exponents (*n*=3). The bulk normal viscosity of a multi-layer decreases (i.e. softens) significantly during folding under bulk pure shear for the applied strain rate hardening flow law. This structural softening is caused by active folding and the rotation of fold limbs without significant limb deformation. For the simulations presented here, the bulk normal viscosity decreased by a factor of 2–20. Some of the pure shear simulations developed multi-layer folds with curved fold axial planes. Fold shapes can be considerably irregular and the fold shapes can be classified as localized and chaotic. The irregularity of the fold shapes in our simulations results from initial geometrical heterogeneities such as localized and random perturbations of the layer interfaces and different initial thickness of individual layers.

For bulk simple shearing, multi-layer folding above a detachment with an asymmetric initial geometrical perturbation does not generate significantly asymmetric fold shapes. This indicates that natural fold shapes exhibiting a more or less symmetric geometry cannot be used as an indicator for bulk pure or simple shear. However, fold axial planes in multi-layer folds are often curved and non-parallel.

For gravity sliding, multi-layer slump folding can generate significantly asymmetric fold shapes, strongly curved fold axial planes and overturned fold limbs. Folds are localized along the multi-layer resulting from stress and strain rate concentrations owing to the gravity-controlled heterogeneous stress and strain rate field. The final fold shapes are strongly sensitive to the power-law rheology in the strong layers for folding driven by gravity.

## Acknowledgements

We thank Giles Hunt for inviting us to contribute to this special volume. We thank Bruce Hobbs for a critical and constructive review and acknowledge helpful comments by an anonymous reviewer. The work of S.M.S. has been supported by the University of Lausanne. D.W.S. acknowledges the Center of Excellence grant from the Norwegian Research Council to PGP (Physics of Geological Processes) at the University of Oslo and the access to the computational facilities provided by the Norwegian Metacenter for Computational Science (NOTUR).

## Footnotes

One contribution of 15 to a Theme Issue ‘Geometry and mechanics of layered structures and materials’.

- This journal is © 2012 The Royal Society