Abstract
Viscous folding in mechanically heterogeneous layers is modelled numerically in two dimensions for linear and power-law viscous fluids. Viscosity heterogeneities are expressed as circular-shaped variations of the effective viscosity inside and outside the layers. The layers are initially perfectly flat and are shortened in the layer-parallel direction. The viscosity heterogeneities cause a perturbation of the velocity field from the applied bulk pure shear, which perturb geometrically the initially flat-layer interfaces from the first numerical time step. This geometrical perturbation triggers interfacial instabilities, resulting in high-amplitude folding. We compare simulations with heterogeneities with corresponding simulations in which the heterogeneities are removed after the first time step, and, hence, only the initial small geometrical perturbations control wavelength selection and high-amplitude folding. Results for folding in heterogeneous and homogeneous layers are similar, showing that viscosity heterogeneities have a minor to moderate impact on fold wavelength selection and high-amplitude folding. Our results indicate that the interfacial instability is the controlling process for the generation of buckle folds in heterogeneous rock layers. Therefore, existing analytical and numerical solutions for folding in homogeneous layers, in which folding was triggered by geometrical perturbations, are useful and applicable to study folding in natural, heterogeneous rock layers.
Folds are one of the most common structures in deformed rocks (e.g. Fig. 1) and occur on all geological scales (e.g. Ramsay 1967; Ramsay & Huber 1987; Hudleston & Treagus 2010; Schmalholz & Mancktelow 2016). Folds often form during layer-parallel shortening of mechanically layered rocks (Fig. 1). Folding of mechanically layered rock under layer-parallel shortening is commonly explained as being a result of an interfacial instability (Biot 1961; Fletcher 1974; Smith 1977; Johnson & Fletcher 1994). Interfacial instability means that small geometrical perturbations on the layer interfaces, which always exist on natural rock surfaces, amplify with rates that are faster than the applied rate of horizontal deformation. Therefore, the layer interface locally displaces faster as it would do under homogenous thickening, thus leading to folding. Most analytical and numerical solutions quantify the interfacial instability for viscous folding in mechanically homogeneous materials: that is, the material properties inside the stiff layer and inside the weaker layers above and below the stiff layer, referred to here as matrix, are constant. The solutions considering mechanically homogeneous viscous materials provide analytical expressions for the dominant wavelength and the maximal fold amplification rate (Biot 1961; Ramberg 1962; Fletcher 1974; Smith 1975; Johnson & Fletcher 1994), which are useful for: (i) quantifying fold amplification, (ii) estimating the viscosity ratio between the layers; and (iii) assessing the bulk shortening associated with observed fold shapes (Hudleston & Treagus 2010; Adamuszek et al. 2016; Schmalholz & Mancktelow 2016).
Natural single-layer fold made from a folded quartz vein in sillimanite schists (Cap de Creus, Spain).
Most natural folded rock layers are, however, polymineralic or exhibit variations in grain size, so that the effective viscosity is never homogeneous in natural layers and in the surrounding rock matrix. Also, the viscous flow during natural rock folding is frequently accompanied by additional processes such as porous fluid flow (Eckert et al. 2016), thermomechanical coupling (Hobbs et al. 2008) or grain-size reduction (Gairola & Kern 1984). A common feature of such additional processes is that they can vary locally the effective viscosity within the folding layer or surrounding rock, because effective rock viscosity can depend on fluid pressure, temperature and grain size (e.g. Karato 2008). Therefore, such processes may cause effective viscosity variations within and around the rock layers during folding (e.g. Adamuszek et al. 2013). However, the impact of mechanical heterogeneities, such as viscosity variations, in viscous layers on the interfacial instability and subsequent fold amplification has not been investigated in detail. Therefore, we compare here numerical simulations of viscous folding with and without viscosity heterogeneities in layer and matrix to evaluate whether the interfacial instability controls high-amplitude folding or whether the heterogeneities control this folding.
We quantify the impact of viscosity variations inside the layer and weaker matrix on folding with two-dimensional (2D) numerical simulations of linear and power-law viscous flow. The viscosity variations are modelled with either weaker or stronger circular inclusions either inside the layer or in the matrix. We compare the results for folding with viscosity variations in the layers to the corresponding results for the folding of homogeneous layers. We focus on the first-order characteristics of the high-amplitude fold evolution, such as the amplification direction, fold amplitude and wavelength. The aims of this study are: (i) to test whether viscosity heterogeneities in the layers and matrix significantly affect the interfacial instability, fold wavelength selection and high-amplitude folding, and consequently; (ii) to evaluate whether existing analytical solutions and numerical simulations considering a homogeneous layer and matrix are applicable to natural, heterogeneous folded rocks.
Model and configuration
We applied the equations of continuum mechanics to describe 2D incompressible slow viscous flow without gravity (Johnson & Fletcher 1994, p. 403; Pollard & Fletcher 2005, p. 280). We used the open-access M2Di solver (Räss et al. 2017) to perform the numerical simulations for both linear and power-law viscous flows. M2Di utilizes the finite-difference method, and has been benchmarked with analytical solutions for linear and power-law viscous folding and is, hence, well suited to perform the folding simulations (Räss et al. 2017). To track the deformation of the layer and inclusion interfaces during high-amplitude folding, M2Di was elaborated by a level-surface approach based on the semi-Lagrangian backwards characteristics method (e.g. Beuchert & Podladchikov 2010). Similar to the level-set method (e.g. Sethian 1998), the level-surface approach relies on the usage of a smooth and continuous function to track the matrix-layer-inclusion boundaries. After each time step, the velocity field, calculated by M2Di, was used to advect the level-surface.
The model configuration included a viscous layer embedded in a less viscous matrix (Fig. 2a). All parameters were dimensionless, and were scaled by the initial layer thickness (H), the applied rate of horizontal shortening (Db) and the layer viscosity (μr). Far-field pure-shear deformation was applied with velocity boundary conditions, and horizontal shortening was parallel to the initial layer interfaces (x-direction). Since Db is negative in compression, we considered its absolute value for scaling. Free-slip was imposed on all model boundaries. The layer interfaces were initially perfectly straight. The effective viscosity ratio between layer and matrix (μr/μm) was 50. For power-law viscous flow, the effective viscosity was strain-rate dependent (e.g. Johnson & Fletcher 1994), and the effective viscosity ratio for each simulation was quantified using the applied rate of horizontal shortening (Db). The power-law exponent was always 3 (n = 3) and was chosen to be representative of many natural cases where dislocation creep was the dominant creep mechanism in rock (e.g. Bürgmann & Dresen 2008).
(a) Model configuration. A viscous layer rests in a less viscous matrix. The layer–matrix interfaces are initially perfectly flat and are indicated by black solid lines. Dashed lines indicate the region shown in (b) & (c) where the initial viscosity heterogeneities are displayed for model configurations 1 (M1) and 2 (M2). (b) A circular inclusion placed in the matrix below the layer (M1). (c) A circular inclusion placed within the layer (M2). H indicates initial layer thickness. See the text for details on the model configuration.
In order to investigate the effect of viscosity heterogeneities on the interfacial instabilities, we considered three model configurations, namely M1, M2 and M3. In the first two configurations we considered the effect of individual circular inclusions, outside (M1) or inside (M2) the layer. This allowed us to investigate the velocity perturbation from the background pure-shear velocity field due to the presence of the viscosity heterogeneity. As a more complex configuration, we considered many circular inclusions that were randomly distributed outside the layer (M3). For model configuration 1 (M1), an inclusion with a viscosity of μr/10 (i.e. five times more viscous than the matrix) was placed below the layer at the horizontal centre (Fig. 2b). The inclusion had a radius of 0.75H. The distance between the inclusion centre and the central horizontal axis of the layer was 3H (Fig. 2b). For model configuration 2 (M2), a strong inclusion with a viscosity of 10μr was placed within the layer and its radius was H/5. Vertically, the inclusion was displaced by a distance equal to three-quarters of its radius (0.75 × H/5) towards the base of the layer (Fig. 2c). For model configuration 3 (M3), the matrix around the initially perfectly flat layer contained 200 inclusions with a radius of H/4 which had a viscosity of μr/10. The distribution of the inclusions was random (Fig. 3). The initial height for models M1 and M2 was 40H, and their initial width was 20H. Model M3 had a larger initial height (80H) and width (40H). All of the applied viscosity variations were significant as they generated an effective viscosity variation that ranged from a factor of 5 to 10 between the inclusion and its surrounding material.
Model configuration 3 (M3). A viscous layer rests in a less viscous matrix. The two solid straight lines indicate the initially perfectly flat layer–matrix interfaces. The small circles indicate the viscous inclusions that are randomly distributed within the matrix. H indicates the initial layer thickness. See the text for details on the model configuration.
Perturbation of the background velocity field and the interface deflection
For each of the three configurations, two types of simulations were performed. For the first type, the simulation was run with the inclusions. For the second type, the inclusions were removed after the first numerical time step (<2% of bulk shortening) and the simulations were continued without inclusions. The removal of the inclusions did not change the numerical mesh since the applied Eulerian mesh is orthogonal with constant grid-spacing. After the first time step, the inclusions caused a deviation from the applied bulk pure shear locally by perturbing the pure-shear velocity field (Fig. 4). The perturbation of the velocity field depends on the inclusion size, and on the viscosity ratio between the inclusion and its surrounding material (Figs 4 & 5). In order to quantify the perturbation velocity induced by a circular inclusion, we used an analytical solution for incompressible viscous flow of a fluid with a circular inclusion (Schmid & Podladchikov 2003). When the inclusion experiences horizontal shortening, the vertical component of the perturbation velocity 
above the horizontal centre of the circular inclusion is given by:
(1)
where |Db| is the absolute value of the applied rate of horizontal shortening, R is the viscosity ratio between the inclusion and its surrounding material, y is the coordinate in the vertical direction and ri is the radius of the inclusion (Fig. 5). The solution provides the vertical perturbation velocity in the material surrounding the inclusion: that is, y ≥ ri. Equation (1) shows that the absolute magnitude of the perturbation velocity generally decreases with increasing distance from the inclusion, proportional to 1/y. Equation (1) also indicates that for no viscosity contrast, R = 1, the perturbation velocity is zero. Furthermore, the sign of the perturbation velocity is different for R > 1 and R < 1. For very large or very small viscosity ratios, R > 1000 and R < 1/1000, the impact of the viscosity ratio becomes negligible (Fig. 5c) because the ratio −(R − 1)/(R + 1) approaches −1 for R > 1000 and 1 for R < 1/1000.
Colour plot of the vertical component of the perturbation velocity 
at the beginning of shortening for configuration M1 (Fig. 2b). 
has been normalized by the initial layer thickness (H) and the absolute value of the horizontal rate of shortening (|Db|). The position of the layer is indicated with solid black lines. (a) The inclusion viscosity (μi) is 10 times larger than its surrounding matrix (μm). (b) The inclusion viscosity (μi) is 10 times smaller than its surrounding matrix (μm).
Vertical components of the perturbation velocity 
on top of a circular inclusion in an infinite matrix. (a) General geometrical configuration. The inclusion is centred at the origin, and the distance from the origin is normalized by the inclusion radius (ri). The dashed line indicates the location of the perturbation velocity that is plotted in (b) (ϕ = π/2). (b) Analytical solution for the vertical component of the perturbation velocity 
as a function of the distance from the inclusion. 
has been normalized by the inclusion radius (ri) and the absolute value of the horizontal rate of deformation (|Db|). Different lines correspond to different viscosity ratios between the inclusion and its surrounding matrix (R = μi/μm). (c) Influence of the viscosity ratio (R) on the magnitude of the perturbation velocity at the inclusion–matrix interface right on top of the inclusion (ϕ = π/2).
The perturbed velocity field causes a small heterogeneous deformation of the layer–matrix interfaces. Therefore, after the first time step, the layer interfaces are no longer perfectly flat but exhibit small geometrical perturbations. The maximal amplitudes of these geometrical perturbations after the first time step are always less than 0.01H (1% of the initial layer thickness). According to the analytical solutions for folding the maximal dimensionless amplification rate, q, for R = 50 and n = 1 is c. 16 and for R = 50 and n = 3 is c. 34 (Fletcher 1974). The amplitude to thickness ratio at which dynamic folding starts to dominate over kinematic layer thickening is c. 1/(2q) and termed nucleation amplitude (Schmalholz & Podladchikov 2001). The initial amplitude of 0.01H in the simulations with removed inclusions is, hence, always smaller than the nucleation amplitude. Therefore, the numerical simulations cover the full range of dynamic fold amplification and the associated wavelength selection process. After removal of the inclusions, it is, hence, only the small geometrical perturbations and the resulting interfacial instability which cause the subsequent high-amplitude folding and control the wavelength selection. The results of the two types of simulations were compared to quantify the impact of the inclusions, representing any kind of strength variation, on high-amplitude folding.
Results
Without any viscosity heterogeneity the initially flat layer is, as expected, homogeneously shortening and thickening without any folding, and its state of stress is homogeneous (results are not shown). The simulations show that any small variation in viscosity causes small geometrical perturbations of the layer interface, which trigger an interfacial folding instability (Figs 6, 7, 8). This can be deduced by observing the vertical component of the velocity perturbation near the geometrical centre of the model domain (Fig. 4a). The direction of the velocity perturbation above, within and below the layer is downwards, and, hence, locally deflects the layer interfaces downwards within the first small time step. Therefore, the layer interfaces are also slightly bent downwards after the first time step. This explains why fold amplification in the model centre, above the inclusion in the matrix, is directed downwards (Fig. 6a).
Colour plots of pressure (P) normalized over the layer viscosity (μr) and the absolute value of the horizontal rate of shortening (|Db|) for model configuration 1 (M1). The numbers shown in the left-hand part of the diagrams correspond to the per cent of bulk shortening 
. The dashed line in (a) indicates the line in the middle of the layer that is used to calculate the apparent fold amplitude (A). (a) & (b) Simulation with an inclusion. (c) & (d) The same model as in (a) & (b) but the inclusion was removed after the first time step.
Colour plots of pressure (P) normalized over the layer viscosity (μr) and the absolute value of the horizontal rate of shortening (|Db|) for model configuration 2 (M2). The numbers shown in the left-hand part of the diagrams correspond to the per cent of bulk shortening 
. (a) & (b) Simulation with an inclusion. (c) & (d) The same model as in (a) & (b) but the inclusion was removed after the first time step.
Colour plots of pressure (P) normalized over the layer viscosity (μr) and the absolute value of the horizontal rate of shortening (|Db|) for model configuration 3 (M3). The numbers shown in the left-hand part of the diagrams correspond to the per cent of bulk shortening 
. (a) & (b) Simulation with inclusions. (c) & (d) The same model as in (a) & (b) but the inclusions were removed after the first time step.
The first-order geometrical fold characteristics, such as deflection direction, wavelength, arc length and amplitude, are similar for folds generated with viscosity heterogeneities and without (Figs 6, 7, 8). Wavelength in this context refers to the horizontal distance between two fold hinges of the same style (i.e. synform or antiform). With respect to amplitude (A), we refer to the apparent amplitude: that is, the vertical distance of a material point that was originally placed in the middle of the layer at y = 0. The coordinates of the points that are used to calculate the apparent amplitude are highlighted with a dashed line in Figure 6a. In addition, the pressure field generally follows the classical patterns observed in earlier folding studies (Schmalholz & Podladchikov 1999). By pressure, we refer to the mean stress defined as positive in compression. The major difference is that the folds which developed in the presence of viscosity heterogeneities (Figs 6a, b, 7a, b & 8a, b) amplify with slightly larger rates, and, hence, these folds are amplified more for the same amount of bulk shortening than the folds in layers without viscosity heterogeneities (Figs 6c, d, 7c, d & 8c, d). In addition, the layers that folded without inclusions can be slightly thicker as they experience the bulk pure-shear shortening and thickening for a longer period of time before they dynamically amplify significantly (Fig. 7c, d). For example, for M1 a bulk shortening of c. 47% in the simulation without inclusions (Fig. 6d), it was necessary to achieve a similar fold amplitude to that which was reached already after 40% bulk shortening in the simulation with inclusion (Fig. 6a). However, for the stage with similar fold amplitude, the fold shapes for simulations with and without inclusions are similar, indicating that the inclusion has a minor impact on wavelength selection and the resulting final fold shape, but a moderate impact on fold amplification. The results are similar to M2 with a viscosity variation inside the stiff layer (cf. Fig. 7a and d).
The presence of a large number of inclusions (M3) in the matrix (Fig. 8a) created random perturbations in the velocity field and, consequently, random geometrical perturbations of the layer interfaces. The fold shapes and associated pressure fields in simulations with and without viscous heterogeneities are similar (Fig. 8b, d). For a more quantitative description of the model differences we plotted the normalized apparent amplitude (A/H) of the folds as a function of their x-coordinate and bulk shortening (
: Fig. 9). The results show that the fold wavelengths, the amplification direction and the location of fold hinges are similar for the runs with and without inclusions once the folds have developed (>40% of strain). The major difference between the models with and without the viscosity perturbations is, hence, the bulk shortening that is needed in order to develop the final fold shape.
Colour plots of the normalized apparent amplitude (A) of the central line in the folding layer (see Fig. 6a) v. bulk shortening 
and normalized horizontal distance (x/H). The apparent amplitude (A) is calculated from the horizontal reference at y = 0. (a) M1 configuration. (b) Same as in M1 but the inclusion has been removed after the first time step. (c) M2 configuration. (d) Same as in M2 but the inclusion has been removed after the first time step. (e) M3 configuration. (f) Same as in M3 but the inclusions have been removed after the first time step.
Simulations for M1, M2 and M3 for power-law viscous rheology show that the importance of the interfacial instability is not limited to linear viscous flow since the results with and without the viscosity heterogeneities are also similar for power-law viscous flow (Figs 10 & 11). Using a more quantitative comparison as in Figure 9, the power-law model results confirm the model results that employed a linear viscous rheology (Fig. 12). The major difference in the models that utilized power-law rheology is that they amplify earlier and therefore the effect of the inclusions on the high-amplitude amplification of the folds is minor. To demonstrate the major differences in the amplification of all the models we plotted the maximum apparent amplitude of each fold train v. the amount of bulk shortening for all model configurations (Fig. 13). These results show clearly that the impact of the inclusions on high-amplitude fold amplification is larger in the models with linear viscous rheology. In addition Figure 13 shows that all the runs significantly amplified after c. 10% of bulk shortening (much later than for the inclusion removal at c. 2%), which shows that the wavelength selection was not affected by a too large initial interface deflection.
Colour plots of pressure (P) normalized over the layer viscosity (μr) and the absolute value of the horizontal rate of shortening (|Db|) for power-law viscous flow in both the layer and matrix (n = 3). (a) Model configuration 1 (M1) after 25% bulk shortening with the inclusion indicated by the small black circle. (b) Same model configuration (M1) after 26% bulk shortening but the inclusion was removed after the first time step. (c) & (d) same as in (a) & (b) but for the model configuration 2 (M2).
Colour plots of pressure (P) normalized over the layer viscosity (μr) and the absolute value of the horizontal rate of shortening (|Db|) for power-law viscous flow in both the layer and matrix (n = 3). (a) Model configuration 3 (M3) after 30% bulk shortening with inclusions indicated by small black circles. (b) Same model configuration (M3) after 31% bulk shortening but the inclusions were removed after the first time step.
Same as in Figure 9 but for power-law rheology (n = 3).
Maximum normalized apparent amplitude (A) plotted v. bulk shortening 
for all models. The maximum value of A was calculated by monitoring the maximal vertical distance of the points that were initially in the middle of the layers (e.g. dashed line in Fig. 6a). The maximum value of A has been normalized over the initial layer thickness (H). Solid lines indicate the amplitude evolution for the runs with the inclusions and the dashed lines indicate the runs where the inclusions have been removed after the first time step. (a)–(c) Models M1, M2 and M3 for linear viscous rheology. (d)–(f) Models M1, M2 and M3 for power-law viscosity (n = 3).
In all the simulations, small viscous inclusions caused a perturbation from the applied bulk pure-shear deformation and, hence, generated small geometrical perturbations of the layer interfaces. The position of the viscous inclusions, and their relative viscosities, have a strong impact on the orientation of the perturbed velocity field. By inspection of Figure 5 and equation (1), we can see that the viscosity of the inclusion relative to its surrounding matrix can change the sign of the velocity perturbation, and therefore the effective viscosity of a heterogeneity can be responsible for the deflection of an initially parallel interface. For example, for M1, an inclusion with a higher viscosity than that of the matrix causes a downwards-directed perturbed velocity field (Fig. 4a), whereas the same inclusion with a viscosity smaller than that of the matrix causes an upwards-directed perturbed velocity field (Fig. 4b). In addition, the magnitude of the perturbation velocity diminishes when the non-dimensional distance (distance over the inclusion radius) from the inclusion grows (Fig. 5b). This means that large inclusions that are far from the layer interface will cause similar-in-magnitude velocity perturbations to those of smaller inclusions which are closer to the interface. Changing the initial position of the heterogeneities and/or their relative viscosity, therefore, affects the shape of the initial interface perturbation (i.e. upwards or downwards). However, the general fold characteristics, such as the fold wavelength or the amplitude, remain similar. In the case where the heterogeneities are placed in a perfectly symmetrical manner with respect to the layer interface, interfacial instability and folding of the layer do not occur (Fig. 14). The latter confirms that the dominant factor in triggering the interfacial instability and, consequently, high-amplitude folding is the asymmetry of the velocity field that causes an asymmetrical, with respect to the vertical direction, geometrical interface perturbation.
Colour plots of pressure (P) normalized over the layer viscosity (μr) and the absolute value of the horizontal rate of shortening (|Db|) for two models which assume that the distribution of the viscous inclusion(s) is perfectly symmetrical with respect to the initially flat central line of the layer at y/H = 0. Both models are shown for 40% bulk shortening. The viscosity of the inclusion(s) is μr/10. (a) Two weak inclusions located above and below the layer interfaces. (b) One weak inclusion initially located exactly at the origin (0,0).
Discussion
Different processes can be responsible for the development of heterogeneous effective viscosities in rock layers. Variations in grain size, temperature, fluid content or mineralogy all can lead to the effectively heterogeneous distribution of viscosity, stress and strain rate (Adamuszek et al. 2013). We have shown that heterogeneous viscosity distribution may cause perturbations from the applied pure-shear velocity field. In the theoretical limit in which the viscous heterogeneities are perfectly symmetrical and the associated perturbation velocities cancel each other out, the flow is stable and pure shear dominates. Consequently, any slight change from this perfectly symmetrical case will perturb the velocity field, causing asymmetrical geometrical perturbations of the layer interfaces and thus giving rise to interfacial instabilities. The effect of viscosity heterogeneities is more profound in the linear viscous models compared to power-law viscous models (Fig. 13).
Our simulations with a single viscosity variation (Figs 6, 7, 8) show that localized fold shapes can develop without any softening mechanisms, simply due to initially localized geometrical perturbations related to localized viscosity heterogeneities. Despite the different mechanisms that can be responsible for the viscosity variations during folding, our results suggest that the interfacial instability has a first-order control on folding in heterogeneous rock layers. Viscosity heterogeneities can trigger folding since they can cause asymmetrical geometrical perturbations of initially flat or symmetrical layer interfaces. Such perturbed layer interfaces trigger an interfacial folding instability. However, during the associated dynamic folding and wavelength selection these heterogeneities play only a minor role. We, therefore, argue that if natural folds exhibit typical features of buckle folds, such as continuous layers without prominent faults or shear zone (Fig. 1), then the controlling process of fold formation was most likely to have been interfacial instability.
Analytical solutions for folding that are based on interfacial instability and on homogeneous materials are, hence, also applicable to heterogeneous materials. However, if the number of the heterogeneities is very large, then the effective viscosity of the material including the heterogeneities can be different, which will affect the folding characteristics (Adamuszek et al. 2013). In that case, an appropriate effective viscosity should be used instead of the ‘matrix’ viscosity. Consequently, the application of these analytical solutions to strain and competence contrast estimation from fold shapes (Schmalholz & Podladchikov 2001), to finite amplitude folding (Muhlhaus et al. 1994; Schmalholz & Podladchikov 2000; Schmalholz 2006), and to folding of the crust and lithosphere (Ramberg & Stephansson 1964; Burov et al. 1993; Schmalholz et al. 2002) is also appropriate and useful for quantifying folding in heterogeneous rock units.
Conclusions
Numerical models of viscous folding in mechanically homogeneous layers are good approximations to folding in heterogeneous layers. Viscosity heterogeneities in initially perfectly flat layers cause geometrical perturbations of the layer interfaces during the initial stages of layer-parallel shortening. These geometrical perturbations trigger an interfacial folding instability. The impact of the studied viscosity heterogeneities on dynamic high-amplitude fold amplification and wavelength selection is minor to moderate, as fold shapes developing during high-amplitude folding are similar for corresponding simulations with and without heterogeneities. Heterogeneities for linear viscous folding have a moderate impact on high-amplitude folding, whereas for power-law viscous folding they have only a minor impact.
Our results indicate that for natural folded rock layers with typical characteristics of buckle folds, such as the absence of prominent faults or approximately constant layer thickness, the interfacial instability was the first-order mechanism that caused the natural folds. Hence, viscosity heterogeneities, or any other kind of strength heterogeneities, is likely to have had only a minor impact on high-amplitude folding. This may also hold true for material heterogeneities developing during folding due to any kind of material softening, which was not significant enough to generate considerable strain localization and associated shear-zone and/or fault development. Existing analytical and numerical solutions based on mechanically homogeneous layers without softening are, hence, useful and applicable to quantifying the generation and evolution of natural buckle folds.
Acknowledgements
We thank Albert Griera and Marta Adamuszek for their comments which helped us to improve our manuscript. The editorial handling by Clare Bond is greatly appreciated.
Funding
This work was supported by the University of Lausanne.
- © 2018 The Author(s). Published by The Geological Society of London. All rights reserved
References
Current volume
The importance of interfacial instability for viscous folding in mechanically heterogeneous layers
