## Abstract

Geometrical heterogeneities along layer interfaces play a key role in determining the geometries of folds developed during shortening of competent layers. We present a series of numerical simulations to investigate the influence of initial sinusoidal perturbations on the folding of single layers. Models consist of a competent viscous single layer embedded in a softer matrix, with the layer orientated parallel to the shortening direction. We first generalize the wide spectrum of sinusoidal perturbations accounting for asymmetries along and across a competent single layer, using two parameters: transversal asymmetry (*A*′) and longitudinal asymmetry (*φ*). These two parameters allow the transition between classical fold shapes and pinch and swell geometries to be studied. The parameter *A*′ describes the development of fold hinges with different geometries between the upper and lower layer interfaces, and abnormal curvatures between the outer and inner arcs of fold hinges. The parameter *φ* induces a strong polarity on the folds, with a systematic preferred orientation of the pinch and swell regions of the layer, even if there is no shear component parallel to the layer. Our results demonstrate the importance of structural inheritance on the resulting fold geometries, and suggest that caution must be taken when using certain types of asymmetrical folds as strain markers and kinematic indicators.

Folds form when stiff or competent viscous layers get shortened. Folds are observable across all length scales in orogenic belts (from micrometre to kilometre scales), and are widely used as kinematic and mechanic indicators (e.g. Hudleston & Treagus 2010; Llorens *et al.* 2013*a*, *b*; Schmalholz & Mancktelow 2016). To commemorate the 50th anniversary of the publication of *Folding and Fracturing of Rocks* by John Ramsay (Ramsay 1967), we here focus on how pre-existing interface geometries of folding layers affect fold formation. Ramsay & Huber (1987) discussed in their book the so-called ‘fish-hook’ folds (Fig. 1) (see also fig. 17.11 of Ramsay & Huber 1987). Initially described by Sorby (1879), and discussed in Ramsay & Huber (1987), the train of fish-hook folds is observed in a bioclastic crinoidal limestone embedded in slates. The structure is characterized by a marked fold polarity (or asymmetry), with large differences in thickness between consecutive limbs and also with marked hinge geometry variations along the upper and lower layer interfaces (Fig. 1). Cleavage is well developed in this classical example and gets refracted in the limestone single layer, indicating a more competent behaviour than that of the host slates.

An issue when interpreting this type of structure is how to explain fold polarity. Ramsay & Huber (1987) wrote that they were not able to provide a confident explanation of this structure, but suggested a way of interpreting it based on a tri-layer system, in which the incompetent host rocks above and below the competent folding layer have different viscosities. They concluded that the only satisfactory explanation for the polarity of the structure was that the initial perturbations already presented this polarity, and suggested that such geometries would arise from competence differences of the two sides of the folding layer (Fig. 1b, stage 1). During non-coaxial deformation the viscosity contrast between the host rock and the folding layer induces a distinct wavelength selection along the upper and lower interfaces of the layer (Fig. 1b, stage 2), and the progressive rotation of the principal shortening axis with respect to the layer causes the fold asymmetry and the preferential thinning of one of the fold limbs (Fig. 1g, stage 3).

Ramsay & Huber (1987) assumed that the original thickness of the layer was uniform, and that thickness variations were developed during deformation. In terms of fold perturbation, the explanation they proposed implies an asymmetry on the amplitude of perturbations developed on the upper and lower interfaces of the competent layer, resulting in an asymmetry transversal to the layer. In addition, although not represented in the figures of Ramsay & Huber (1987), the selected wavelengths could be different along both interfaces, and the fold hinges at one of the interfaces did not match the location of the corresponding one on the other interface. This condition results in the development of a longitudinal asymmetry parallel to the layer.

Recent in-depth reviews have summarized the information that can be obtained from fold analysis (e.g. Hudleston & Treagus 2010), as well as the fundamental parameters controlling folding (e.g. Schmalholz & Mancktelow 2016). However, most of the existing studies on folding have only considered folds developed from perturbations defined by regular arrays of low-amplitude and symmetrical geometries, and the influence of irregular and/or inherited perturbations has not been addressed, although geometrical heterogeneous folds are abundant in the field (Fig. 2a–c). Assuming that certain fold types, such as ‘fish-hook’ folds, originate from inherited layer perturbations, the geometries of perturbations along and across the folding layer have to be clearly defined. This contribution presents an investigation of the influence of initial asymmetries along and across a competent single layer embedded in a ductile matrix during layer-parallel shortening. This is systematically studied by means of numerical simulations. We first define the parameters used to generalize the sinusoidal perturbations accounting for asymmetries along and across a competent single layer. These variations are described using two parameters that measure the degree of asymmetry along and across to the layer envelope. After a preliminary analysis of the problem using the analytical theory of folding, the finite evolution of the structures nucleated from these perturbations is systematically explored using a series of numerical simulations that consider viscous rheology. Finally, the main implications of our results on the use of folds structures as strain and kinematic indicators are discussed.

## Layer-interface geometry

The occurrence of small and regular/random irregularities on layer interfaces is assumed in all theories of single-layer folding (e.g. Biot 1957, 1961; Ramberg 1960; Fletcher 1974, 1991; Smith 1975, 1977, 1979; Schmalholz & Podladchikov 2000). Finite wavelength, amplitude and shape of evolved folds are strongly influenced by the initial geometry of these perturbations (Abbassi & Mancktelow 1990; Mancktelow 1999), and have been investigated by means of analogue experiments (Cobbold 1975; Abbassi & Mancktelow 1990, 1992) as well as numerical simulations (Williams *et al.* 1978; Zhang *et al.* 1996, 2000; Mancktelow 1999). However, these studies have only considered single layers of constant thickness defined by regular arrays of low-amplitude perturbations and symmetrical shapes. Results from experiments of Williams & Jiang (2001, see their figs 2 & 7) on single-layer folding using rock analogue materials show the strong influence of initial irregularities during folding. Large variations on fold structures were observed in their experiments, although the samples were deformed using the same configuration, materials and deformation conditions. Following a kinematic study, they interpreted that slight differences in the relative geometry between both layer interfaces can produce high deviations of the nucleated and amplified structures during layer-parallel shortening.

Smith (1975, 1977) defined four main types of structures that can grow from the deformation of a single layer embedded in a matrix, depending on: (1) the viscous contrast between the layer and matrix; and (2) whether the layer is subjected to shortening or lengthening. When a single layer experiences shortening, folds or mullions are expected to develop if the layer is stiffer (folds) or weaker (mullions) than the matrix. Additionally, when the layer is subjected to parallel extension, pinch and swell structures (or necking) are expected to form if the layer is stiffer than the matrix, while reverse folding structures develop if the layer is weaker than the matrix.

From a geometrical point of view, these four structures can arise from two basic initial layer-interface perturbation geometries (Fig. 3c): (a) symmetrical, when the perturbation is mirrored both along the layer and normal to it (i.e. with a ‘pinch and swell’ or ‘neck’ shape) and (b) antisymmetrical, when the perturbation is only mirrored normal to the plane but not parallel to it (i.e. resulting in a deflection or ‘fold’ shape). Although these perturbations can reflect different deformation conditions not considered here, from a geometrical point of view these geometries can be regarded as the end members of a progressive transition between ‘fold’ and ‘pinch and swell’ shapes (Fig. 3a). This has been explicitly taken into account in several numerical and experimental studies, in which mullion structures were modelled using both symmetrical and antisymmetrical perturbations (Sokoutis 1987, 1990; Lan & Hudleston 1991).

In general, layers in natural rocks are likely to deviate from perfectly planar geometries and tend to contain random irregularities along layer interfaces, as well as lateral variations of layer thickness (Fig. 2d–g). These irregularities can sometimes be systematic as a result of primary structures (e.g. sedimentary cross-bedding, lenticular bed geometries, irregular igneous dykes, etc.). At other times, a new deformation phase is superposed to previously deformed materials, resulting in the formation of arrays of periodic and/or random layer-interface irregularities, as, for example, pinch and swell structures or asymmetrical boudinage (Fig. 2d–f). However, most of the published experimental and numerical studies only focused on layers with planar geometries, initial symmetrical layer-interface perturbations or those characterized by very small random irregularities. However, non-systematic perturbations are expected across and along layer interfaces in the vast majority of natural rocks.

## Generalization of the perturbation geometry

Prior to presenting and discussing the series of numerical simulations, here we present a generalization of the perturbation geometry. Similarly to that in classical fold theories (e.g. Biot 1957, 1961; Ramberg 1960; Fletcher 1974, 1991; Smith 1975, 1977, 1979; Schmalholz & Podladchikov 2000), a regular and sinusoidal shape of perturbations is considered. Assuming an initial flat layer parallel to the reference *x*-axis, the geometry *ζ*(*x*) of perturbations along the upper interface of a layer can be described as:
*H* is the layer thickness, *A* is the amplitude and *L* the wavelength of the perturbation (Fig. 3b). A similar expression can be used to define the geometry of the lower interface:

If the amplitude of both interfaces has the same sign, the shape of the perturbation is antisymmetrical (Fig. 3c). Contrarily, if the amplitudes of the two interfaces have different signs, then the shape of the perturbation is symmetrical. The variation between both end members can be described using a normalized parameter *A*′, namely transversal asymmetry. *A*′ is defined as the ratio between the amplitude of the upper and lower interfaces, *A*′ = *A*_{upper}/*A*_{lower}. The geometry of the perturbation is antisymmetrical for *A*′ = 1, while *A*′ = −1 denotes the case of a symmetrical configuration (Fig. 3b). Values of *A*′ ranging from −1 to +1 define geometries in-between the two end members. *A*′ = 0 is a particular case in which one of the interfaces is planar.

In addition to the asymmetry of perturbations across the layer, asymmetries along the layer can also be defined. If we consider trains of sinusoidal perturbations with constant wavelength, an additional term can be added to describe the relative shift between hinges of the upper and lower interfaces. This parameter, namely longitudinal asymmetry (*φ*), can be expressed as:
*d*′ is the distance between consecutive maximum and minimum amplitudes of the lower and the upper interfaces, respectively (Fig. 3b). Note that this definition is similar to the one used for expressing the phase difference between two sinusoidal waves with the same wavelength. Using equation (3), *φ* = *π*/2 describes a train of sinusoidal perturbations with an antisymmetrical geometry (Fig. 3a, *A*′ = −1 and *φ* = *π*/2), while *φ* = 0 denotes a symmetrical geometry (Fig. 3a, *A*′ = −1 and *φ* = 0). Values between the two end members (*φ* ∈ [0, *π*/2]) are equivalent to perturbations that display geometries similar those of asymmetrical boudins or shear bands (e.g. Swanson 1992; Goscombe & Passchier 2003). Although the treatment is different, a value of *φ* ≈ *π*/2 yields similar geometries those used in the analogue experiments of Abbassi & Mancktelow (1990).

Using the transversal and longitudinal asymmetry parameters, equation (1) can be rewritten as:

Sinusoidal perturbations with asymmetries along and across the layer (i.e. in the *x*- and *y*-directions, respectively) can be defined using both parameters. The range of explored values is limited to *A*′ ∈ [−1,1] and *φ* ∈ [0, *π*/2] in this study.

In this contribution, the term ‘boudin’ is used to describe the swell region of the layer (i.e. where the layer thickness reaches a maximum), whereas the term ‘neck’ describes the pinch part of the layer. Note that these terms do not imply a mechanical interpretation of the origin of these perturbations and are only used to describe layer geometries. From a physical point of view, these ideal perturbations let us explore the influence of the layer-thickness variation and the effect of non-coincident hinge/inflection points on fold geometries.

## Analytical solution

The classical theory of the folding of linear viscous layers is based on the studies of Biot (1959, 1961) and Ramberg (1960). Later on, other authors modified the theory to include non-linear rheologies (Fletcher 1974; Smith 1975, 1977) or extended it for the three-dimensional case (Fletcher 1991, 1995; James & Watkinson 1994). However, these analytical studies are restricted to infinitesimal-amplitude solutions and perturbations with simple sinusoidal geometries. Remarkably, Schmalholz & Podladchikov (1999, 2000) proposed a new finite-amplitude theory of folding that overcomes these restrictions and provides results that match those obtained from numerical simulations (Schmalholz & Podladchikov 2001). A useful review of theoretical studies and analytical solutions for single-layer folding is provided by Schmalholz & Mancktelow (2016).

For a two-dimensional system, the rate of amplification of an infinitesimal perturbation in a linear viscous medium is given by (e.g. Schmalholz & Mancktelow 2016):
*A* is the amplitude, *t* is the time and *q*, corresponds to the amplification due to a passive growth rate. According to Johnson & Fletcher (1994), the infinitesimal amplitude solution for *q* can be expressed as:
*m* is the viscosity contrast between the layer and matrix (*η*_{l}/*η*_{m}), and *k* is the wavenumber (i.e. *πH/L*). For layer-parallel shortening, Smith (1975, 1977) demonstrated that the expression of *q* is similar for both symmetrical and antisymmetrical perturbation cases. The difference between the two of them is the change in sign of the terms in the quotient of equation (6). Upper signs are used for symmetrical cases, while lower signs refer to asymmetrical ones. The sum of the kinematic (i.e. passive amplification) and dynamic growth rates (1 + *q*) is plotted against *L*/*H* for both types of perturbations in Figure 3. For the antisymmetrical case, the dynamic growth rate is higher than the kinematic amplification and thus the perturbation is amplified. In such a scenario, the amplification rate follows an exponential growth rate (Ramberg 1960; Biot 1961) but quickly changes to a non-linear growth rate related to the finite amplitude of folding (see Schmalholz & Podladchikov 2000). For symmetrical cases, dynamic growth rates, *q*, are negative and, therefore, the predicted amplification is lower than the kinematic amplification due to homogenous deformation. For a constant viscosity contrast, the curves obtained for the antisymmetrical and symmetrical perturbation cases constrain the maximum and minimum values of the amplification fold growth-rate spectra. Note that the maximum and minimum amplification rates occur at the same *L*/*H*, and are coincident with the dominant wavelength to thickness ratio (*L*_{d}*/H*). For these cases, the dominant wavelength is *L*_{d} = 2*πH*(*m/*6)^{1/3}. General perturbations are expected to have growth rates ranging between these limits, as shown in Figure 4.

Burg *et al.* (2004) defined an analytical solution for the growth rate of a folding layer with a planar interface (i.e. *A*′ = 0). The dynamic growth rate, *q*, can be expressed as follows, based on equation (8) of Burg *et al.* (2004) and assuming folding in the absence of gravity:

The dominant wavelength is *L*_{d} = 2*πH*(*m/*3)^{1/3} in this situation. Note that for the case with a flat interface, the growth rate decreases, while the dominant wavelength increases with respect to the previous scenario with antisymmetrical perturbation (Fig. 4). For a hypothetical case with *m* = 100, the dominant wavelength is *L*_{d}/*H* = 16 for the antisymmetrical perturbation case, while for the example with a planar interface it is *L*_{d}/*H* = 20, slightly larger. For the cases with 1 < *A*′< 0, the growth rates are expected to be essentially very similar (Fig. 4). However, the expected range of growth rates is larger in models with *A*′< 0. An analytical solution of the growth rate using stability analysis is not available for the case of longitudinal asymmetries, *φ*.

The analytical solution of Burg *et al.* (2004) is based on a linear first-order approximation, and is only valid for gentle and infinitesimal-amplitude perturbations. Therefore, this solution cannot be used to predict the finite evolution of folds. As mentioned above, Schmalholz & Podladchikov (1999, 2000) proposed a new finite-amplitude theory of folding that overcomes these restrictions and provides results coherent with numerical simulations (Schmalholz & Podladchikov 2001). However, the present problem cannot be simplified to the study of only one of the interfaces, because there are differences in amplitudes and shifts of hinge/inflection points across the layer. This can potentially imply hinge migration and non-homogenous distribution of shear stress across and along the layer, at least at low-deformation stages. Schmalholz *et al.* (2008) obtained an analytical solution for finite-amplitude necking based on the assumption that plane sections remain plane (PSRP). This approach can potentially be used to predict the finite evolution of asymmetrical perturbations because it can handle large lateral variations in layer thickness. However, this is beyond the scope of this paper and the problem is treated here by means of numerical simulations.

## Methods and model set-up

### Numerical method

The numerical simulations were performed using the 2D finite-difference program FLAC (Cundall & Board 1988; Itasca 1998). The equations of motion in a continuum medium were discretized and solved using a dynamic relaxation scheme (Cundall & Board 1988). The medium was simulated using a structured mesh composed of four-node polyhedral elements, where a mixed discretization scheme was used for the handling of the volumetric constraints. Both schemes (i.e. dynamic relaxation and mixed discretization) provide an efficient strategy and robustness in modelling strain localization (Poliakov & Hermann 1994). The code has been used widely to simulate layer folding (e.g. Zhang *et al.* 1996, 2000; Toimil & Griera 2007) and other tectonic structures (Takeda & Griera 2006 and references therein).

### Model geometry and boundary conditions

The geometry of the model consisted of a competent layer embedded in a weaker matrix. The thickness of the layer was 2 units, the model length 80 units and the width 44 units. Following other studies based on single-layer folding simulation (e.g. Zhang *et al.* 1996, 2000; Mancktelow 1999), the contacts between layer and matrix were considered coherent or ‘welded’: that is, without the possibility of slipping and/or opening. Between 2880 and 3520 quadratic elements were utilized, and the mesh was progressively refined around the layer–matrix contact. A total of six elements were used to represent the layer width. This resolution was high enough for a correct distribution of stress and strain inside and around the layer. All the models were deformed under coaxial plane-strain conditions with a constant strain rate of 2.5 × 10^{−14} s^{−1}. Progressive shortening parallel to the layer was imposed by velocity boundary conditions applied to the sides of the model. Velocities normal to the sides were fixed, while parallel velocities were not constrained (i.e. free slip boundary conditions). The maximum shortening applied in all models was 55%, equivalent to a natural strain of *ε* = 0.8.

Two basic perturbation configurations were used: (a) models with periodic small perturbations (i.e. constant wavelength and amplitude); and (b) models with a single isolated perturbation located at the centre of the layer. These configurations can be considered equivalent, but not similar, to previous models used by other authors in their numerical simulations (e.g. Williams *et al.* 1978; Zhang *et al.* 1996, 2000; Mancktelow 1999). The geometries of these perturbations were defined using equations (1) and (3). For models with a periodic perturbation, the maximum amplitude (i.e. hinge) of the perturbations coincides with the model boundaries (i.e. hinge point at boundary). In addition, periodic boundary conditions in the *x*-direction were imposed to minimize boundary effects, using an approach similar than that of the ELLE platform (e.g. Jessell *et al.* 2009; Llorens *et al.* 2013*a*, *b*; Gomez-Rivas *et al.* 2017). For all models, the geometry of the upper interface was kept constant, and was defined by initial ratios of *L/H* = 20 and *A*/*H* = 0.2, while the geometry of the lower interface varied according the values of *A*′ and *φ*.

### Material properties

A linear elastoviscous model was used to describe the material rheology. This constitutive model is equivalent to a combination of a linear elastic element and a linear viscous element in series following a Maxwell model (e.g. see Ranalli 1987). The total strain rate is assumed to be the sum of the elastic and viscous strain rates. The elastoviscous properties of the material were defined using the shear modulus (*G*), the bulk (or incompressibility) modulus (*K*) and the shear viscosity (*η*) of the material. A Poisson's ratio (*ν*) value of 0.30 was assumed for both the folding layer and the matrix. In all of the simulations, volumetric strain was considered purely elastic. Material properties are summarized in Table 1. These values are within the range inferred from experiments of natural rocks in ductile conditions (Turcotte & Schubert 1982) and similar to those utilized by other authors (e.g. Mancktelow 1999; Zhang *et al.* 2000).

For elastoviscous materials, important material parameters are the contrasts in elasticity (*R*) and viscosity (*m*) between the folding layer and the matrix (Zhang *et al.* 2000). Here, we have opted to follow an approach where both parameters are independent, rather than using similar values for both parameters (i.e. *R* = *m*). This assumption is coherent because viscosity in natural rocks ranges by up to four orders of magnitude (e.g. Carter & Tsenn 1987), while the ratio of elastic properties is only up to one order (e.g. Turcotte & Schubert 1982; Mancktelow 1999). For all numerical simulations an elastic contrast of *R* = 2 and a viscosity contrast of *m* = 20 were used. The ratio between the applied strain rate and the relaxation time of the material was used to scale the elastoviscous rheology. The Deborah number (*De*) is the ratio of the stress magnitude to the elastic shear modulus (Reiner 1964). *De* was lower than 10^{−4} for the layer and matrix in these models, representing a mechanical response close to the expected flow for an ideal viscous material (Poliakov & Hermann 1994).

## Results

### Folding involving transversal asymmetry *(A′)*

The effect of transversal asymmetry on fold shapes for periodic perturbations is examined in Figure 5. All fold geometries correspond to a natural strain of *ε* = 0.8 and can be compared with their initial shape (Fig. 3a; configurations with *φ* = 0 and −1 *≤ A*′ *≤* 1). As expected, the fold shapes that develop are strongly controlled by the initial perturbation geometry. The increase in the parameter *A*′ tends to generate an asymmetry across the layer. In such cases, differences in the geometry and curvature of external and internal hinges at both interfaces are observed. Hinges remain fixed at their initial material points during deformation, and migrations or shifts of hinges during fold amplification are not observed. For the condition *A*′ = −1, layer-parallel shortening is mainly accommodated by layer thickening and slight amplification of the perturbation. This is in agreement with the analytical solution, where low rates of dynamic amplification are observed (Fig. 4). After 70% shortening, the symmetry of the structure across and along the layer is still conserved and fold hinges remain aligned (Fig. 5, *A*′ = −1). However, layer thickening is not homogenous at this point and is mainly accommodated at pinch regions. An increase in *A*′ progressively produces structures similar to typical fold shapes, presenting reinforced amplification on both interfaces (i.e. amplifications in the same direction) with an attenuation of layer-thickness variations (Fig. 4). For all cases (except *A*′ = −1) the final fold shapes are dominated by the lower interface. This clearly applies to the case of *A*′ = −0.5, a model that starts with perturbations on both the upper and lower interfaces with amplitudes in opposite directions, but ends with reinforced amplitudes (Fig. 5). In models with *A*′ ≤ 0, a marked difference in the fold shapes between hinges located in pinch and swell regions is observed. Outer and inner arcs show smooth curvatures in swell areas (Srivastava & Lisle 2004). In addition, pinch regions are characterized by high arc curvatures, and shapes ranging between sinusoidal and hyperbolic. Differences in the amplification between both interfaces are observed, with the lowest *A*/*H* values measured at the upper interface of the layer. Note that these geometries resemble arc-and-cusp structures (e.g. Ramsay & Huber 1987, p. 403), with strong hinge curvatures pointing to the lower amplified interface. The progressive amplification of the model *A*′ = 0 (with only one horizontal interface) is shown in Figure 6. In this model, hinge zones are fixed at the same material points during deformation, and the initial zones with maximum curvatures remain fixed at the material points. In this case fold amplification is associated with the migration of inflection points towards the inner hinges. The final geometry of the folding layer presents thickness variations, with pinch and swell regions being located at fold hinges and resulting in an asymmetry across the layer. In the early deformation stages, the differential stresses are higher at necking (or pinch) regions as a consequence of the thinner layer in these areas. However, progressive folding in the layer-parallel stress field displays a typical distribution with high tensional stresses located along external fold arcs and compressional stresses at inner hinges.

The evolution of the averaged horizontal stress transversally to the layer in the neck and pinch regions is shown in Figure 7. With the exception of the *A*′ = −1 case, all numerical simulations record a decrease in the averaged stress with increasing natural strain, in accordance with strain softening related to fold amplification (e.g. Schmalholz *et al.* 2005; Llorens *et al.* 2014). A marked stress relaxation occurs in the pinch regions rather than in swell regions (e.g. model *A*′ = −0.5). This larger structural softening of pinch regions in the models *A*′ = 0 and *A*′ = −0.5 results in a stronger stress reduction compared to that found in models with a constant layer thickness (*A*′ = 1). Contrarily, the *A*′ = −1 case maintains relatively large values of the average stresses and non-homogeneous evolution. Regardless of bulk softening or hardening, swell parts of the layer stiffen with progressive deformation while pinch portions soften. Therefore, a symmetrical configuration of the perturbation geometry is associated with a nearly constant bulk horizontal stress at all deformation stages.

Models with an isolated perturbation at the centre of a single layer show a similar evolution to those with periodic perturbations (Fig. 8). However, the influence of transverse asymmetry on fold geometry is more evident. A progressive transition from ‘pure’ fold shapes to folded swell regions is observed with decreasing *A*′. The model *A*′ = −1 is the end-member case in which all shortening is accommodated by layer thickening. Except for *A*′ = 1, fold hinges are located in the thickening part of the layer. There is a weak tendency to develop and propagate lateral folds at the edges of swell regions. This is more relevant in models with −1 < *A*′< 0, in which the maximum strain localization and fold amplification is observed on sideways folds (e.g. model *A*′ = −0.5 in Fig. 8). Additionally, in this range of models the maximum hinge curvatures are preferentially observed in the outer fold arcs compared to the inner fold arcs. This observation contradicts the general tendency of folds to show stronger curvatures along inner arcs (Hudleston & Lan 1994; Hudleston & Treagus 2010).

### Folding involving longitudinal asymmetry *(φ)*

The influence of the longitudinal asymmetry (*φ*) on fold shape was investigated using periodic perturbations (Fig. 9) and isolated perturbations (Fig. 10). The shapes of the resulting folds show a marked variation with the parameter *φ*. The addition of a relative shift between the hinges of both interfaces produces the development of asymmetrical folds (Figs 9 & 10). This also applies to cases in which *φ* is small. For example, the perturbation in the model *φ* = *π*/32 is strong enough to induce asymmetrical folding, even if the hinge shifts are smaller than *d′* = 0.03 units (Fig. 9). An increase in *φ* causes more intense fold amplification and a reduction in the deformation accommodated by layer-parallel thickening. Although the initial wavelength is predefined as *L* = 20, the initial hinge lateral shift produces a variation in layer thickness, thus causing differences in the initial wavelength to thickness ratio (*L*/*H*) between pinch and swell regions.

Figure 11 shows the evolution with shortening of the model with *φ* = *π*/16. During the first folding stages, hinges migrate and new wavelengths emerge (Fig. 11). The models tend to develop shorter wavelengths at pinch regions or thinned portions of the layer (*L*_{pinch} in Fig. 11 and *ε* = 0.8), while longer wavelengths are observed at thicker zones of the layer (*L*_{swell} in Fig. 11 and *ε* = 0.8). When limb dips attained values of 10°–20°, hinge migration stops and hinges remain at the same position during further fold amplification. Differences in dip, length and layer thickness can be observed between both fold limbs (Figs 10 & 11). The final geometry is that of folded pinch and swell structures with a preferred location of swell zones in the long limbs, while pinch regions are located in the short fold limbs. A decrease in the parameter *φ* (Fig. 9) defines the transition between folds with symmetrical (*φ* = *π*/2) and asymmetrical limbs (*π*/2 < *φ* < 0). The wavelength of short limbs progressively decreases up to the end-member case of *φ* = 0, where the short limb is only equivalent to the pinch region (Fig. 10).

The influence of parameter *φ* on the limb dip (*α*) of pinch and swell regions is shown in Figure 12. For a constant *φ* value, pinch regions attain larger dips than swell regions. The highest dip values at swell regions increase with *φ*. However, this is not observed for pinch regions in which higher dip values are observed in intermediate models (i.e. in-between *φ* = *π*/4 and *φ* = *π*/16). The first derivative of limb dip change with respect to the natural strain (∂*α*/∂*ε*: ‘rotation rate’) is shown in Figure 12b with respect to the limb dip. This diagram provides an alternative way to analyse the influence of *φ* on the fold growth rate. Limb rotation rates first increase and then decrease after attaining a maximum limb dip ranging between 20° and 40° (note that this applies for both positive and negative limb dips in Fig. 12). At low limb dips (−15 < *α* < 15°), the rotation rate is linear and independent of *φ*. This observation is coherent with equation (4), where the dynamic growth rate of a linear viscous material only depends on the viscosity contrast and wave number. However, when limb dips increase, rotation rates become dependent on the longitudinal asymmetry *φ*. In general, the rotation rates at pinch regions are up to two or three times higher than those observed at swell areas for the same *φ* configuration. For the case of *φ* = 0, there is a thickening of the pinch and swell structure without growth of fold structures (i.e. *α* = 0). Systematically, the values of rotation rates at swell zones are lower than in the *φ* = *π*/2 case (i.e. antisymmetrical case, where fold limbs are symmetrical). However, rotation rates at pinch regions are higher than in the *φ* = *π*/2 case. The influence of increasing *φ* on fold limb rotation rates is not straightforward. In swell zones, an increase in the longitudinal asymmetry *φ* raises the rotation rate and displaces the peak of maximum rotation towards lower limb dip values. Conversely, a direct relationship between *φ* and rotation rates is not observed in pinch regions (Fig. 12). The highest rotation rate values occur at intermediate values of *φ* in pinch regions. This is coherent with the peak displacement at higher limb dips.

These observations can be interpreted as being a consequence of the relative ratio of wavelength to thickness (*L*/*H*) of pinch and swell regions with respect to the theoretical dominant wavelength to thickness ratio (*L*_{d}/*H*). For the case of a viscosity contrast of *m* = 20, the initial *L*/*H* = 20 used in the models is higher than the dominant wavelength to thickness ratio (*L*_{d}/*H =* 10). The introduction of the *φ* parameter produces a modification of the initial *L*. If we assume a layer with a constant thickness, the wavelength increase at swell zones produces higher *L*/*H* ratios than the initial ones, resulting in lower amplification rates (sees Fig. 3). Contrarily, a reduction in *L* at pinch areas results in a shorter *L*/*H*, which is close to the *L*_{d}/*H* value, and thus produces higher amplification rates.

Geometrical differences between folds located at the centre and edges of the model are detected, even if we have used periodic boundaries along the *x*-axis to minimize boundary effects (Figs 9 & 11). Typically, folds in the centre of the model show less fold limb rotation but more thickening than folds at the edges. In models with longitudinal asymmetry (*π*/2 < *φ* < 0), a slight rotation of the layer envelope is observed, although the initial layer orientation was parallel to the shortening *x*-axis (Fig. 9). The sense of rotation of the layer envelope is the same as that for swell areas (or long limbs), and opposite to that of pinch zones (or short limbs), where larger dip limbs develop.

Models with isolated perturbations (Fig. 8) present similar tendencies to those observed from simulations with periodic perturbations, although marked differences on the final geometry can be observed. Deformation is strongly partitioned at the margin of the swell regions (i.e. overthickened parts of the layer) in models with isolated perturbations. In such cases, folds propagate on the sides of the swell. Deformation is preferentially accommodated by layer thickening in simulations with low longitudinal asymmetry *φ*, while increasing *φ* leads to a higher amplification and rotation of the swell zone. As for the case with periodic perturbations, the swell is located at the long and less rotated fold limb. However, the condition of free amplification and propagation of lateral folds reduces the final rotation of swells compared to the models with periodic perturbations. From a mechanical point of view, it seems more efficient to propagate folds along the thin layer than deforming and rotating the overthickened swell area. The latter case would require a higher stress build-up. Layer envelopes do not rotate in isolated perturbation models.

### Buckling involving transversal and longitudinal asymmetries

Figure 13 shows an example of fold amplification of a single layer with periodic perturbations defined by a combination of transversal (*A*′ = −0.5) and longitudinal (*φ* = *π*/16) asymmetries. As expected, according to previous results, differences in fold shape across and along the layer envelope clearly influence the resulting fold trend. Remarkably, the developed folds display polarity, in the way that thick and thin fold hinges are preferentially located. This geometry is not a consequence of heterogeneous layer deformation during folding, but reflects the initial perturbation geometry. At the final stage, the hinges of the upper and lower layer interfaces are not aligned and, hence, fold-axial planes show vergence with respect to the layer-envelope normal.

The influence of both types of asymmetries (*A*′ and *φ*) on the final fold geometry is also explored for the case of an isolated layer perturbation (Fig. 14). All the models presented in this figure were carried out using a constant value of *A*′ = −1/4, with the longitudinal asymmetry *φ* ranging between *π*/120 and *π*/10. The end-member case of *φ* = 0 is shown in Figure 8, illustrating the complexity of the finite fold geometries. A progressive geometrical transition is observed from a case in which the thick part of the layer corresponds to the hinge of an antiformal symmetrical fold (*φ* = *π*/120) to a situation in which the thick part of the layer is folded and defines the hinge and limb of the developed fold (*φ* = *π*/10). In general, strong disharmonies of fold geometries with significant differences of curvature of fold hinges at both interfaces are developed. Although not displayed, axial-plane orientations are highly variable in the matrix surrounding folds, with orientations varying from vertical attitudes in the case of pinched synformal hinges to subhorizontal ones in the case of external arcs of antiformal swell portions.

## Discussion

Modelling of natural structures normally requires the simplification of observations in order to capture the key elements for understanding rock structure evolution (e.g. Schmalholz *et al.* 2008; Schmalholz & Mancktelow 2016). For example, folds are generally described as periodic structures, and geologists tend to study examples where regular fold trains are observed (e.g. Ramsay & Huber 1987; Hudleston & Treagus 2010, figs 1 & 2). However, many natural folds are non-periodic structures with irregular distributions of wavelength and amplitude, and thus present variable fold geometries (e.g. Fletcher & Sherwin 1978; Hudleston & Holst 1984; Abbassi & Mancktelow 1990, fig. 3). Fold shapes are generally more variable in terranes affected by multiple deformation phases (Fig. 2a–c), in which there are layers with non-homogenous thickness and where pre-existing structures induce structural inheritance. These observations are coherent with experiments in which initial layer irregularities exert a strong influence on the final fold geometry (e.g. Abbassi & Mancktelow 1990, 1992). These authors observed that strongly asymmetrical folds could develop from a small asymmetrical perturbation despite the imposed layer-parallel shortening in pure-shear conditions. Our numerical results reaffirm this interpretation. Fold shapes in all our simulations reflect the geometry of the initial perturbations. This is valid at least for elastoviscous materials in which the elastic response is low, as in the present case. For high strain-rate conditions (d*ε*/d*t* = 1 × 10^{−10} s^{−1}), where build-up of stresses cannot be relaxed by viscous flow and where stresses can reach values of the order of GPa, it is possible that fold development is independent of the initial irregularities (Zhang *et al.* 2000, fig. 9). However, such stresses are in general uncommon considering the properties of crustal rocks, and thus the material would fail by brittle or ductile strain localization rather than being folded.

The issue of structural inheritance is not restricted to small-scale structures such as the folds simulated in this study, but represents a multiscale problem in structural geology and tectonics. This issue has attracted attention from the large-scale tectonics community; such as, for instance, the classical example of tectonic inversion of pre-rifted margins during a collision stage (e.g. Butler *et al.* 2006 and references therein). Frequently, the geometry of thrust sheets and folded nappes developed in orogens are related to geometrical features of pre-existing extensional basins, including their length, width or orientation with respect to the main reactivated faults (e.g. von Tscharner *et al.* 2016). A direct application of our results to larger-scale tectonic problems is not straightforward due to the more complex material rheologies and the requirement to include gravity effects. Therefore, more advanced mechanical models are required to address similar problems at the crustal scale. However, our results are useful for illustrating the influence of inherited structures on the geometries of subsequent ones.

One of the aims of this study is to propose a simple approach to generalize the geometry of sinusoidal perturbations and use it to explore their influence on the geometry of the resulting folds. The definition of the two proposed geometrical parameters (i.e. transversal asymmetry *A*′ and longitudinal asymmetry *φ*) allows a wide spectrum of geometrical perturbations, ranging between classical antisymmetrical (or ‘fold shape’) and symmetrical (or ‘boudin shape’) geometries, to be described. Any geometry in-between the two end members can be described with these two parameters, either in isolation or by combining them (Fig. 3). Therefore from a geometrical point of view, pinch and swell geometries result in a transition case between the two end members. The presence of transversal or longitudinal asymmetries along a layer results in the development of a directional polarity of the structure. For example, *A*′ produces finite folds displaying differences in fold shapes between the upper and layer lower interfaces. For conditions of *A*′ > 0, upper and lower interfaces reflect the initial perturbations. However, in cases with *A*′ < 0, the initial geometry of one of the interfaces is not directly observable from the finite geometry of folds and yields anomalous geometries in which the curvature of inner arcs is lower than that outer arcs (Figs 5 & 8). This contradicts the expected geometries of folds formed in viscous conditions from an initially flat layer (e.g. Hudleston & Lan 1994).

The folding of a layer with longitudinal asymmetry (summarized using the parameter *φ*) results in a systematic asymmetry of fold limbs and the formation of folds with vergence. Asymmetrical folds are generally abundant in natural examples and are typically utilized to infer the sense of shearing. However, numerical simulations demonstrate that the development of asymmetries during folding is a critical issue. Results of single-layer folding simulations in simple-shear conditions show that, in general, folds formed in such conditions present more or less symmetrical geometries and are relatively similar to those formed under layer-parallel compression (Viola & Mancktelow 2005; Llorens *et al.* 2013*b*). Numerical models also indicate that asymmetrical folds only develop, and can be preserved with progressive deformation, in systems with a multilayered matrix (Llorens *et al.* 2013*a*). Frehner & Schmalholz (2006) showed that the initial folds in multi-layered systems are symmetrical, but the longer fold limbs rotate and small folds become asymmetrical with progressive strain due to the relative shearing between fold limbs. Recently, Frehner & Schmid (2016) have shown that a pre-existing geometrical asymmetry can survive during multilayer folding if the pre-existing perturbation is relatively tight and exhibits relatively high amplitudes. Otherwise flexural flow and flattening related to the longer fold limbs will result in de-amplification and unfolding. Asymmetrical folds can probably also develop in systems with an anisotropic matrix and non-homogenous simple-shear conditions (Ran *et al.* 2018). In any case, one of the ways to explain asymmetrical folds is certainly based on the presence of initial systematic irregularities associated with primary structures (bedding and other stratigraphic surfaces, intrusive emplacement geometries, etc.), resulting in the formation of longitudinal asymmetries. In general, however, a systematic polarity of asymmetries is not expected in natural cases. Another option is that asymmetries arise from pre-existing deformation stages, such as in terranes deformed in multiple tectonic phases. For example, an early phase of layer-parallel extension can produce a systematic development of initial irregularities that cause the development of fold asymmetries of folds during layer shortening in a subsequent deformation phase.

Although transversal and longitudinal asymmetries have been treated here as pre-existing geometries, there are mechanical processes that can induce their development with progressive deformation. For example, Gardner *et al.* (2015, 2016) showed that during extension of a viscous layer embedded in a matrix with different viscosities on the two sides of it (i.e. in a tri-layer model), the pinch and swell or neck structures that develop present differences in amplitude between the upper and lower interfaces. Therefore, an effective way to develop a transversal anisotropy (*A*′) is by the presence of a viscosity contrast in the materials surrounding the layer. Additionally, Gardner *et al.* (2015, 2016) demonstrated that asymmetrical boudins can form by strain localization and shearing. As shown by Duretz & Schmalholz (2015), this asymmetrical shearing during layer-parallel extension is only effective for multilayers if the material is power-law viscous. Accordingly, this phenomenon can dynamically produce longitudinal asymmetries, such as those arising from numerical models, without requiring large amplitudes to induce asymmetries during subsequent folding. Finally, the emergence of systematic asymmetries under symmetrical deformation conditions (as in the case of pure-shear boundary conditions) has also recently been observed by Schöpfer *et al.* (2017) in numerical simulations of the faulting of layers subjected to layer-parallel extension. Domains with parallel faults that have a systematically similar displacement sense form when the matrix surrounding the brittle layer has low strength. The formation and slip along parallel faults induces a fault-bonded block rotation of the layer to balance the vorticity components. A similar effect was observed in our simulations, as indicated in Figure 9, with periodic sinusoidal perturbations with an initial longitudinal asymmetry. In our case, the layer envelope rotates with the same sense as that of swell areas (or long limbs), and opposite to that of pinch zones (or short limbs) where larger rotations developed. Even if this phenomenon is, in part, influenced by model boundary effects, the layer envelope rotates to counteract the unbalanced vorticity originated by the final asymmetry of the folded array, attending to the initial orientation of the layer (i.e. parallel to the shortening direction) and the coaxial boundary conditions, which present a theoretical vorticity of zero (e.g. Lister & Williams 1983).

Figure 2c, d show examples of folded quartz veins hosted in schists from the Cap de Creus Variscan Belt (Spain). Folds present complex shapes and display irregular layer-thickness distributions. These folds were interpreted to have formed by two different deformation phases (Carreras & Druguet 1994). An early extension parallel to quartz veins and schists produced the development of boudinage, pinch and swell structures and a heterogeneous layer thickness of quartz veins. During a subsequent deformation phase, quartz veins were folded. Some of the fold geometries resemble those obtained from the numerical models presented here, although natural patterns are more complex than those in our simulations. For example, the asymmetrical arrangement of boudins and preferred rotation of pinch parts with respect to boudins is in agreement with the presence of an initial longitudinal asymmetry (see Fig. 2a–c). Moreover, folded swell portions, the localization of folds at boudin margins and folds with smoothed inner arcs are indicative of the presence of initial layers that already had strong thickness variations. However, a systematic variation compatible with a transversal asymmetry (i.e. a transverse polarity in folds) was not observed in this field example. The early deformation phase may have produced an asymmetrical boudinage of the rocks. This boudinage would have predetermined the geometry of the amplified folds developed during the layer-parallel shortening phase. Alternatively, asymmetrical folds can be explained as a result of deformation with a sinistral sense of shear parallel to the layers during shortening (Druguet *et al.* 1997). Assessing the suitability of either a kinematic interpretation or providing an explanation based on the arrangement of initial layer-interface perturbations requires a more extensive study of the geometries of these structures in the field area. Although in this case, both deformation kinematics and initial perturbations probably influenced the final structure. However, numerical simulations demonstrate that during pure-shear deformation of a stiff layer parallel to the shortening direction, asymmetrical structures can develop as a consequence of the initial shape of layer irregularities, without implying non-coaxial deformation. The use of multiple kinematic indicators can help to constraint the sense of shear, but it is critical to take into account that discrepancies between kinematic indicators and fold asymmetry may be due to the systematic influence of the geometry of the initial perturbations.

The generalization of perturbation geometries illustrates well the problem of the folding of layers with a non-homogenous thickness. A simple and idealized case represents the folding of a pre-existing pinch and swell structure. It is worth noting that our analysis is restricted to a limited spectrum of the possible range of initial perturbations. Our numerical simulations were carried out using a constant wavelength between pinch and swell regions, and a linear elastoviscous behaviour. However, there is a much wider field of initial settings to be explored, including variable layer length ratios or different thickness ratios between pinch and swell regions. In addition, if a power-law viscous rheology is considered, the strain-localization pattern and fold-amplification behaviour would be more intense (e.g. see Llorens *et al.* 2013*b*; Ran *et al.* 2018). A systematic analysis is beyond the scope of this paper, but the results presented here are a starting point to promote further investigations on the influence of a non-homogenous layer thickness on fold development.

Finally, and returning to the ‘fish-hook’ folds of Ramsay & Huber (1987), our results can be used for reinterpreting non-periodic folds with irregular distributions of wavelength and amplitude. As described in the introduction, a unique feature of ‘fish-hook’ folds is that they present a marked geometrical polarity, with differences in fold-hinge geometry between the upper and lower layer interfaces. The latter authors proposed an interpretation of ‘fish-hook’ structures based on a tri-layer system, where the matrices above and below the competent layer have different viscosities. During non-coaxial deformation, the viscosity contrast between the layer and the matrix induces distinct wavelength selection along the upper and lower interfaces of the layer (see fig. 17.18 of Ramsay & Huber 1987). A more simple explanation is that the structure is derived from the folding of a pinch and swell structure with longitudinal and transversal asymmetry. Note the similarity between the models *ϕ* = *π*/16 shown in Figure 13 and *ϕ* = *π*/10 in Figure 14 with the geometry of the ‘fish-hook’ structures in Figure 1.

## Conclusions

We have investigated the effect of the initial layer-interface geometry on single-layer folding by means of numerical simulations. The geometry of a periodic or isolated sinusoidal perturbation has been generalized using two shape parameters, namely transversal asymmetry *A*′ and longitudinal asymmetry *φ*. In this way, we can study the gradual transition between the two classical end-member fold geometries (i.e. symmetrical and antisymmetrical). These parameters describe the initial variation in layer thickness, with thickened layer regions (i.e. swell) and thinned layer parts (i.e. pinch). During layer-parallel shortening, the initial perturbation geometry strongly controls the finite amplitude and geometry of resulting folds. These parameters can describe a directional polarity of the resulting folds, either along or across the layer. Increasing these parameters enhances fold amplification. In the case of transversal asymmetry, the results indicate that the upper and lower fold interfaces can exhibit different shapes, and curvature can thus differ along consecutive fold hinges. However, variations in layer thickness tend to be compensated. Models ranging between −1 < *A*′ < 0 are exceptional cases, because the initial geometry of one of the interfaces is not directly observable from finite folds and yields anomalous geometries where curvatures of inner arcs are lower than those of outer arcs. In the case of longitudinal asymmetries, a preferential arrangement of pinch and swell regions develops, with variations in wavelength and dip along consecutive fold limbs. In general, thin and short limbs (i.e. pinch portions) rotate faster than thick and long ones (i.e. swell portions), and produce a variation in wavelength and dip along consecutive limbs. The migration of hinge points with respect to the initial perturbation geometry takes place during the first folding stages and is required to accommodate the wavelength selection. Despite the limitations of numerical simulations, in which an idealized and sinusoidal perturbation is utilized, the results can be qualitatively used to better understand irregular and non-periodic natural fold patterns, and improve structural analysis. Our study illustrates the influence of inherited structures on the geometries of folds developed in successive tectonic phases, and demonstrates that caution must be taken when using asymmetrical folds for strain analysis and as kinematic indicators.

## Acknowledgements

Paul D. Bons is acknowledged for reviewing an early version of this manuscript. We thank M. Cooper and E. Moulas for thoughtful comments that helped to improve the manuscript.

## Funding

This work was supported by Spanish Ministerio de Economía y Competitividad (MINECO) (Funding agency ID: http://dx.doi.org/10.13039/501100003329 with grant CGL2014-54180-P to A. Griera and M.G. Llorens. E. Gomez-Rivas acknowledges the support of the Beatriu de Pinós programme of the Government of Catalonia's Secretariat for Universities and Research of the Department of Economy and Knowledge (2016 BP 00208).

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