1. Introduction
Folds are one of the most common geological features in ductile compressional regimes. It is a wavy and continuous structure that provides valuable insights into the deformation history and rheological properties of deformed rocks (Schmalholz et al., Reference Schmalholz, Podladchikov and Schmid2001; Bastida et al., Reference Bastida, Aller, Bobillo-Ares and Toimil2005; Schmalholz, Reference Schmalholz2006; Druguet et al., Reference Druguet, Alsop and Carreras2009; Hudleston and Treagus, Reference Hudleston and Treagus2010; Yakovlev, Reference Yakovlev2012a, Reference Yakovlev2012b; Llorens et al., Reference Llorens, Bons, Griera, Gomez-Rivas and Evans2013). Over the past few decades, scientists have worked on the formation of various fold geometries using analogue modelling (e.g., Currie et al., Reference Currie, Patnode and Trump1962; Ramberg, Reference Ramberg1963; Hudleston, Reference Hudleston1973; Abbassi and Mancktelow, Reference Abbassi and Mancktelow1992), analytical techniques (e.g., Biot et al., Reference Biot, Ode and Roever1961; Johnson and Fletcher, Reference Johnson and Fletcher1994; Hunt et al., Reference Hunt, Muhlhaus, Hobbs and Ord1996) and numerical simulations (e.g., Dieterich, Reference Dieterich1970; Cobbold, Reference Cobbold1977; Zhang et al., Reference Zhang, Hobbs and Ord1996; Mancktelow, Reference Mancktelow1999; Schmalholz et al., Reference Schmalholz, Podladchikov and Schmid2001).
Depending on the number of layers involved in the folding event, folds can be categorized as either single-layer or multilayer structures. Ramberg (Reference Ramberg1960, 1961) first clarified a layered system as a true multilayer or as independent single layers based on the spacing in relation to the dominant wavelength. Later, Schmid and Podladchikov (Reference Schmid and Podladchikov2006) redefined true multilayers (1/λ b < h/s < λ b ) and effective single layer (h/s < 1/λ b ) on the basis of the ratio of thicknesses of soft to hard layers (s/h) and viscous dominant wavelength (λ b ). Hudleston and Treagus (Reference Hudleston and Treagus2010) grouped the classical layered system into three subsystems: effective single layer, true multilayer and independent layer.
In nature, independent single-layer folding is rare, as the ideal assumption of an infinite surrounding medium without influence from adjacent layers is seldom met. Folding of multilayer is relatively more complex, as it involves more number of controlling parameters than independent layer folding. These complexities give rise to various fold patterns such as chevron folds (e.g., Williams, Reference Williams1980; Ramsay, Reference Ramsay1974; Bastida et al., Reference Bastida, Aller, Toimil, Lisle and Bobillo-Ares2007), parallel and similar folds (e.g., Johnson and Pfaff, Reference Johnson and Pfaff1989), kink folds (e.g., Dewey, Reference Dewey1965) and parasitic folds (e.g., Frehner and Schmalholz, Reference Frehner and Schmalholz2006; Frehner and Schmid, Reference Frehner and Schmid2016; Liu et al., Reference Liu, Eckert, Connolly and Thornton2020).
Such layered systems produce harmonic, polyharmonic and disharmonic folds (Fig. 1). Typically, harmonic folds refer to different folded layers with similar wavelength and axial planes, which are in vicinity with each other (Fig. 1a). Typical harmonic folds, completely free of higher-order folds, are formed in true multilayer system. In independent layered system, as the spacing is large enough for the layers to behave independently within an effectively infinite medium, disharmonic folds of different orders develop depending on the thickness, spacing, and viscosity contrast between the layers and the embedded medium (Fig. 1b). When layer thicknesses differ and the spacing lies between the above two scenarios, polyharmonic folds develop (Fig. 1c). These are characterized by small-scale folds with short wavelengths and low amplitudes within larger folds having greater amplitudes and wavelengths (Ramsay and Huber, Reference Ramsay and Huber1987; Price and Cosgrove, Reference Price and Cosgrove1990).
Schematic diagram of (a) Harmonic folds, (b) Disharmonic folds and (c) Polyharmonic folds.
Unlike classical concept, harmonic folds can also be developed in layered mediums, which belong to neither ‘true multilayer’ nor ‘independent layer’. The present study is based on a model in which competent layers of different thicknesses are embedded within an infinitely extended incompetent matrix, as opposed to being sandwiched between thin incompetent units, as typically observed in multilayer systems. Consequently, the notion of a systematic alternation of competent and incompetent layers does not align with the conceptual framework of this research. Furthermore, the system under investigation involves individual competent layers that influence each other’s folding behaviour within their respective contact strain zones (Ramberg, Reference Ramberg1963). While the setup may resemble a multilayer system, the underlying conceptual distinctions are significant, and the current study specifically addresses a different class of layered system.
Biot (Reference Biot1957) introduced the concept of the ‘dominant wavelength’ of a folded structure, governed by the viscosity contrast and layer thickness. Previous studies (Biot, Reference Biot1957; Ramberg, Reference Ramberg1962, Reference Ramberg1964) also reveals that, in a system with similar rheological contrasts, thinner layers are folded first during deformation and subsequent folding of thicker layers play a crucial role in determining the geometry of large-scale folds (Fig. 1). During progressive deformation, higher-order folds in thinner layers are rarely preserved, rather mimic the geometry of thicker layers with longer wavelengths (lower-order folds), resulting in secondary folding. Using Biot’s theory of dominant wavelength, we often estimate viscosity contrast from folded layers. However, natural systems often deviate from this ideal scenario, leading to potential misinterpretations of viscosity contrasts based on observed fold geometries, as adjacent layers can influence local strain distribution.
In deformed terrains composed of metasedimentary rocks, alternating layers of varying competence and thickness are common. If both thinner and thicker layers have the same rheology, the thinner layer is typically perturbed first during early deformation stages (Ramberg, Reference Ramberg1962). According to Ramberg (Reference Ramberg1960), the mutual, layer-perpendicular influence of a folded layer extends significantly only up to a distance equal to its initial wavelength on either side. Beyond this distance, the influence decreases exponentially and becomes negligible if the spacing exceeds twice the initial wavelength – a region referred to as the ‘contact strain zone’ (Ramberg, Reference Ramberg1963). During progressive deformation, folding of the thinner layer may influence the instability patterns and development of folds in adjacent thicker layers within the contact strain zone.
Since the primary objective of this study was to investigate the effect of mutual interaction on fold development in a layered system, all thicker layers were placed within the contact strain zone of the central thinner layer and aligned parallel to the shortening direction. In this study, we use finite element modelling (FEM) to investigate fold development in a layered system, focusing on mutual interactions between layers. We also examine deviations in fold wavelength from the dominant wavelength as defined by Biot (Reference Biot1957).
2. Numerical modelling
2.a. Model considerations
2.a.1. Material
To understand the effect of mutual interaction on fold geometry, we considered a three-layer system of similar rheology embedded within a relatively softer medium. Numerical simulations were performed using the commercial FEM software ABAQUS (Liu et al., Reference Liu, Eckert and Connolly2016, Reference Liu, Eckert, Connolly and Thornton2020). All simulations were conducted in two dimensions under plane strain conditions (Ramsay and Lisle, Reference Ramsay and Lisle2000; Arslan et al., Reference Arslan, Passchier and Koehn2008; Samanta and Deb, Reference Samanta and Deb2014; Samanta et al., Reference Samanta, Basu Majumder and Sarkar2017; Basu Majumder and Samanta, Reference Basu Majumder and Samanta2023). Both the layers and the embedding medium were modelled with Maxwell viscoelastic rheology (cf. Zhang et al., Reference Zhang, Hobbs and Ord1996; Mancktelow, Reference Mancktelow1999; Passchier and Druguet, Reference Passchier and Druguet2002; Samanta and Deb, Reference Samanta and Deb2014; Basu Majumder and Samanta, Reference Basu Majumder and Samanta2023) to simulate high-pressure and high-temperature lower crustal conditions (Turcotte and Schubert, Reference Turcotte and Schubert1982).
The Poisson’s ratio and relaxation time were kept constant at 0.25 and 8 × 108 s, respectively (Larsen et al., Reference Larsen, Motyka, Freymueller, Echelmeyer and Ivins2005; Samanta et al., Reference Samanta, Basu Majumder and Sarkar2017; Basu Majumder and Samanta, Reference Basu Majumder and Samanta2023), for all materials used. The rheological behaviour is defined by the following constitutive equation:
$${{{d\varepsilon }}\over{{dt}}} = \;{{1}\over{G}}\;{{{d\sigma }}\over{{dt}}} + \;{{\sigma }\over{\eta }}$$
where ε, σ and t are the instantaneous strain (finite under steady-state deformation), stress, and time, respectively; G and η are the Maxwell shear modulus and viscosity (Passchier and Druguet, Reference Passchier and Druguet2002, Arslan et al., Reference Arslan, Passchier and Koehn2008, Arslan et al., Reference Arslan, Koehn, Passchier and Sachau2012, Samanta and Deb, Reference Samanta and Deb2014, Samanta et al., Reference Samanta, Basu Majumder and Sarkar2017; Basu Majumder and Samanta, Reference Basu Majumder and Samanta2023).
2.a.2. Model specification
To obtain optimal results within reasonable computational limits, a rectangular block of 5000 × 3000 units was used (e.g., Passchier et al., Reference Passchier, Mancktelow and Grasemann2005; Eckert et al., Reference Eckert, Connolly and Liu2014; Liu et al., Reference Liu, Eckert and Connolly2016, Reference Liu, Eckert, Connolly and Thornton2020; Damasceno et al., Reference Damasceno, Eckert and Liu2017), as shown in Fig. 2. A central thinner layer of 5 units thickness (t h ) was placed at the centre. Two thicker layers of equal thickness (2nd, 4t h or 6t h ) were symmetrically placed on either side of the thinner layer, with variable initial spacing (s) between adjacent layers.
Schematic diagram of model. The thickness of the central thin (dark) layer is 5 units (t h ) and the thicknesses of adjacent thick (dark) layers are 10 units (2t h ) which varied in different model. The initial spacing (s) between two adjacent layers is defined by a spacing factor, n (s = nλ dv , where spacing factor, n = 0.5, 1.0, 1.5 and 2.0) and wavelength (λ dv ) of initial sinusoidal perturbation of the central thin layer.
A bulk layer-parallel shortening of 30% was applied throughout the simulation (e.g., Biot, Reference Biot1957; Ramberg, Reference Ramberg1964). All simulations were carried out using the designated ‘VISCO’ step in ABAQUS (Liu et al., Reference Liu, Eckert and Connolly2016, Reference Liu, Eckert, Connolly and Thornton2020), employing a full Newtonian solution technique and a direct equation solver. The model was discretized using structured quadrilateral meshing with a global mesh size of 2.5, and element type CPE4 (4-node bilinear plane strain quadrilateral). The total number of elements was approximately 2,400,000. A constant natural strain rate of 10⁻¹3 s⁻¹ was maintained (e.g., Jeng et al., Reference Jeng, Lai and Teng2002; Barraud et al., Reference Barraud, Gardien, Allemand and Grandjean2004; Hobbs et al., Reference Hobbs, Regenauer-Lieb and Ord2008; Ord and Hobbs, Reference Ord and Hobbs2013) under steady-state conditions.
According to Ramberg (Reference Ramberg1964), layer buckling requires the presence of flaws or heterogeneities either within the layers or in the surrounding medium, especially in the absence of geometric irregularities. In the present study, sinusoidal as were applied only to the central thinner layer (e.g., Biot, Reference Biot1957, Reference Biot1959a, Reference Biot1959b; Schmid and Podladchikov, Reference Schmid and Podladchikov2006; Frehner and Schmid, Reference Frehner and Schmid2016). Since the primary objective of this study was to investigate the effect of mutual interaction on fold development in a layered system, all thicker layers were placed within the contact strain zone of the central thinner layer and aligned parallel to the shortening direction (Fig. 2).
The nature of the initial perturbation depends on the R factor, which determines whether a layer behaves viscously (R < 1) or elastically (R > 1) during folding (Schmalholz and Podladchikov, Reference Schmalholz and Podladchikov1999; Schmalholz et al., Reference Schmalholz, Podladchikov and Schmid2001; Liu et al., Reference Liu, Eckert and Connolly2016, Reference Liu, Eckert, Connolly and Thornton2020). The R factor is defined as the ratio between the viscous dominant wavelength (λ dv ) and the elastic dominant wavelength (λ de ):
$$R = \;{{{{\lambda _{dv}}}}\over{{{\lambda _{de}}}}} = \;\scriptstyle{}\root{3}\of{{{{1}\over{6}}{{{{\eta _M}}}\over{{{\eta _L}}}}}}\sqrt {{{{{P_0}}}\over{G}}} = \scriptstyle{}\root{3}\of{{{{1}\over{6}}{\eta _R}}}\sqrt {{{1}\over{G}}4{\eta _M}\dot \varepsilon}$$
Where η R = η L /η M is the viscosity ratio of the competent layer (η L ) to the embedding medium (η M ), G is the shear modulus, P 0 is the initial layer-parallel stress, and εis the strain rate.
In our simulations, the R-value ranged between 0.01 and 0.10, indicating deformation occurred within the viscous regime (Table 1). Therefore, a sinusoidal perturbation with an initial amplitude of 0.005 λ dv (i.e., 0.5% of the dominant wavelength) was imposed on the central thinner layer. The dominant wavelength for a viscous layer (Biot, Reference Biot1957) is given by:
$${\lambda _{dv}} = 2{t_h}\pi \;\;\scriptstyle{}\root{3}\of{{{{1}\over{6}}\;{\eta _R}}}$$
Material properties used in numerical modelling
Where t h is the initial thickness of the competent thinner layer. No perturbation was applied to the thicker layers, which remained straight and aligned parallel to the shortening direction. The central thinner layer’s thickness was fixed at 5 units, while the thicker layers had thicknesses of 10 (2t h ), 20 (4t h ) and 30 (6t h ) units in three different models (Fig. 2). In our simulation, the spacing among layers was varied according to:
$$S = n{\lambda _{dv}} = 2\pi n{t_{h\;\;}}\scriptstyle {}\root{3}\of{{{{1}\over{6}}\;{\eta _R}}}$$
Where the spacing factor n = 0.5, 1.0, 1.5 and 2.0.
Throughout all simulations, the viscosities of the thinner and thicker layers were kept the same, while the viscosity ratio η R between the layers and the surrounding medium was varied across different models (Table 1).
2.a.3. Boundary conditions
The entire model was compressed using layer-parallel shortening under pure shear conditions. Bulk shortening was fixed at 30% to maintain the overall thickness of the individual layers throughout their length. The equations for the applied shortening, expressed in terms of displacement at all four boundaries, are as follows:
$$U_x = - k.x $$
$$U_y = {{k}\over{{1 - k}}}.y$$
In our case, since bulk shortening was fixed at 30% for all simulations, k = 0.3.
2.b. Model results
Total 36 number of numerical simulations were run by varying three fundamental parameter of buckle fold. They are: a) viscosity ratio (η R = η L /η M ) of layer to embedded medium, b) thickness ratio of thick to thin layers and c) spacing between two adjacent layers (s). The aim of this study is to find out how the layers of same competency, with different thicknesses and their mutual spacing, affect the final fold geometry. To achieve the objective, we introduce W-factor (λ ob /λ f ), the ratio of wavelength obtained from numerical simulation (λ ob ) to the wavelength derived after 30% shortening of the dominant wavelength (λ f = 0.7λ dv ) of the respective layer. If W = 1, it means that the dominant wavelength is preserved after deformation. W > 1 indicates the deviation from the dominant wavelength.
2.b.1. Progressive development of fold geometry
The focus of this study is to examine how minute instabilities in a thinner layer can perturb adjacent, initially straight, thicker layers embedded in a homogeneous medium. It is well established that a homogeneous layer embedded in a uniformly distributed medium will not buckle if it is oriented parallel to the shortening direction, even if the viscosity contrast exceeds the critical threshold for buckling (Ramberg, Reference Ramberg1962; Reference Ramberg1964).
To investigate this, we placed thicker layers on either side of a minutely perturbed thinner layer (see Section 2.a.2) with identical rheology and observed the resulting fold geometry during progressive layer-parallel shortening. In the first setup, the layer-to-medium viscosity ratio (η R ) and the thickness of the thicker layers were 100 and 20 units, respectively, with a spacing of 1.5λ dv (as defined in Equation 4).
At 5% shortening, there was no visible instability in the thicker layers (Fig. 3a). At 10% shortening, the amplitude of the thinner layer increased, and mild instabilities appeared in the thicker layers (Fig. 3b). Our numerical results reveal that higher-order folds in thinner layer exhibit asymmetry along the limbs of lower-order folds, while maintaining symmetry in their hinge zones, which are similar to the findings of Ramberg (Reference Ramberg1963). At 15% shortening, well-developed folds became visible in the thicker layers (Fig. 3c), and the thinner folded layer began to mimic the large-scale folds, forming a polyharmonic pattern. With further shortening (30%), the fold amplitude of the thicker layers increased (Figs. 3d–3f), while the folds in the thinner layer were modified, resembling chevron-like fold in a few cases (Figs. 4 and 6).
Progressive development of folds. Viscosity contrast (η R) is 100. Thicknesses of the thin and thick layers are 5 and 20 units, respectively. The spacing between adjacent thick and thin layers is 1.5 λ dv . The percentage of shortening of the models (a to f) ranges from 5% to 30% with an interval of 5%. Note that initial sinusoidal perturbation with wavelength λ dv and amplitude 0.005λ dv (Eqn. 3) was imposed only on the central thin layer. Adjacent thick layers was initially straight and folded later by the heterogeneous strain produced by the central thin layer and the fold geometry gradually transforms from disharmonic to polyharmonic in progressive deformation.
Patterns of harmonic folds generated by varying viscosity contrast. (a) η R = 100, (b) η R = 500 and (c) η R = 1000. Thicknesses of the thin and thick layers are 5 and 30 units, respectively. Spacing (1.5λ dv ) between layers increases with increasing viscosity contrast (η R) following Eqn. 4. Finite shortening is 30%.
Graphical comparison of fold geometry of thick (solid line with open marker) and thin (dashed line with solid marker) layers for different layer thickness (t h ). (a) η R = 100, (b) η R = 500 and (c) η R = 1000. Hollow diamond, triangle and square markers represent 10, 20 and 30 units thickness of the adjacent thick layer, respectively. Similarly, solid markers represent thin layers under the influence of corresponding thick layers. Finite shortening is 30%. Note that deviation from deformed dominant wavelength (λ f ) is more for thicker layer and it increases with the decreasing spacing.
Overall, homogeneous thickening at the hinge of the thicker layer was only 0.95% after 30% shortening, which further decreased with increasing viscosity contrast. For example, with η R = 1000, the thickening was only 0.3%, which was deemed negligible in this study. These results suggest that minimal perturbation in the central thinner layer can both initiate and control instability in adjacent thicker layers during early deformation stages.
2.b.2. Effect of viscosity ratio
In our simulations, the viscosity ratio (η R ) of the layers to embedded medium was taken 100, 500 and 1000 (Fig. 4). It was observed that increasing viscosity contrast reduces the preservation of higher-order folds in the central thinner layer; these are mostly retained at hinge zones. While visual inspection shows minimal impact of viscosity ratio on fold geometry, the graphs provide more precise insight (Fig. 5). W-factor decreases with increasing spacing factor (n), regardless of viscosity ratio. However, W increases significantly with higher viscosity ratio at lower n-values. Comparatively, the W-factor for the thinner layer shows limited variation relative to the thicker layer, indicating that higher-order folds in thinner layers effectively represent rheological contrasts.
Fold patterns generated by varying thickness of the adjacent thick layers. (a) 10 units, (b) 20 units and (c) 30 units. t h = 5 units, η R = 1000, n = 1.5. Finite shortening is 30%.
2.b.3. Effect of layer thickness
To analyse the influence of layer thickness, we fixed the thinner layer’s thickness at 5 units and varied the thicker layers’ thickness to 10, 20, and 30 units. The graph of W-factor versus spacing factor (n) shows a pattern similar to viscosity ratio changes (Fig. 7).
Graphical comparison of fold geometry of thick (solid line with open marker) and thin (dashed line with solid marker) layers for different viscosity contrast (η R ). (a) 2t h , (b) 4t h and (c) 6t h . Hollow diamond, triangle and square markers represents thick layers with viscosity contrast (η R) of 100, 500 and 1000, respectively. Similarly, solid markers represent thin layers under the influence of corresponding thick layers. Finite shortening is 30%. Note that deviation from deformed dominant wavelength (λ f ) increases with the decreasing spacing.
As thicker layer’s thickness increases, the fold geometry shifts progressively from disharmonic to polyharmonic, and the number of small-scale folds diminishes (Fig. 6). Unlike viscosity ratio effects (Section 2.b.2), the impact of layer thickness on fold harmonicity is more prominent.
2.b.4. Effect of spacing
As the spacing factor (n) increases, the discretization of layers is preserved, leading to a smaller deviation from the single-layer results, where ‘deviation’ refers to the W-factor (λ ob /λ f ), indicating how the observed wavelength differs from the dominant wavelength in a deformed state. This deviation is demonstrated by the shift in slope from steep to gradual at the transition point of n = 1. Beyond this point, the deviation slowly decreases as the spacing factor increases. This pattern indicates that the discrete nature of the layers diminishes, resulting in interactions within the contact strain zone, as depicted in Figures 5 and 7.
With increasing spacing, the fold geometry transitions from harmonic to polyharmonic, and eventually to disharmonic (Fig. 8). As spacing increases, fold harmonicity decreases, and the number of small-scale folds in the thinner layer increases.
Fold patterns generated by varying spacing (nλ dv ) in between adjacent thin and thick layers. (a) n = 0.5, (b) n = 1, (c) n = 1.5 and (d) n = 2. Viscosity contrast (η R) is set as 500, and thickness of thin and thick layers are 5 and 20 units. Finite shortening is 30%.
From the W-factor versus layer thickness graph, it is evident that increasing spacing reduces the W-factor, suggesting a diminishing effect from the contact strain zone (Fig. 9).
Graphical comparison of fold patterns of thick (solid lines open marker) and thin (dashed lines with solid marker) layers for different spacing (nλ dv ) in between two adjacent thick and thin layers. (a) n = 0.5, (b) n = 1, (c) n = 1.5 and (d) n = 2. Hollow diamond, triangle and square markers represents thick layers with viscosity contrast (η R) of 100, 500 and 1000 respectively. Similarly, solid markers represent thin layers under the influence of corresponding thick layers. Finite shortening is 30%. Note that deviation from deformed dominant wavelength (λ f ) increases with the thickness of the adjacent layer.
3. Discussions
3.a. Influence of contact strain zone on the harmonicity of thinner-layered folds
From the results presented in Section 2, it is evident that the ‘contact strain zone’ (Ramberg, Reference Ramberg1960) plays a pivotal role in determining the harmonicity of fold geometry. One of the most striking findings of our study is that even a slight instability in a thinner competent layer can induce perturbation in a thicker layer placed at a distance up to twice the initial wavelength of the thinner layer. This observation reinforces Ramberg’s theory and analogue models (1960, 1962, 1963 and 1964), confirming the significance of the contact strain zone. In our case, higher-order folds in thinner layers are modified and imitate the geometry of lower-order folds of thicker layers, resulting in secondary folding in later stage of deformation. This gives rise to a wide variety of fold geometries, ranging from harmonic to disharmonic.
Interestingly, the overall fold pattern becomes disharmonic when the small-scale primary folds of the thinner layer are preserved. These higher-order folds typically vanish when the thicker layer is positioned too close to the thinner layer. This obliteration is most pronounced when the spacing factor (n) is 0.5, as also documented by Ramberg (Reference Ramberg1964). Within this close spacing, the individual identities of discrete layers are lost, and the system behaves as a single package. With increasing spacing, this loss of identity gradually reduces. In other words, the mutual interaction of adjacent layers within their respective contact strain zones dictates the final fold geometry. The degree of influence is largely governed by the initial wavelength: larger wavelengths dominate the final geometry, but their impact diminishes as spacing between layers increases. A thinner layer with higher-order folds can mirror a thicker layer’s large-wavelength folds only if the thicker layer is adequately close; otherwise, the thinner layer’s original folds tend to be preserved as spacing increases (Fig. 8). In overall, our numerical results (Fig. 4 to 9) align with Ramberg’s (Reference Ramberg1960) theory of folding within the contact strain zone.
3.b. Estimation of viscosity ratio
In Section 3.a, we attempted to identify the possible causes behind folds with varying harmonicity. Thicker layers, due to their longer wavelengths, often appear to govern the overall fold pattern in a layered medium. It is commonly believed that the dominant wavelength is preserved most clearly in the thicker layer. However, this notion may lead us up the garden path – it could result in inaccurate estimations of finite shortening and mislead us in calculating the viscosity ratio between layer and medium. This effect was recognized by Sherwin and Chapple (Reference Sherwin and Chapple1968) for folds in single layer, showing that thicker layers undergo significant uniform layer–parallel shortening during folding – causing their observed wavelength–to–thickness ratios to fall below the prediction of Biot et al. (Reference Biot, Ode and Roever1961) – whereas thinner layers, which better satisfy Biot’s small–shortening assumptions, preserve the dominant wavelength more accurately.
From Figures 5, 7 and 9, it is evident that the W-factor is sensitive to both the thickness of the layers and the spacing between adjacent layers, even under a constant rheological contrast. As our results suggest, the dominant wavelength of a thinner layer is preserved more clearly in the higher-order folds located at the hinges of its larger, secondary folds. This is because the thinner layer is typically the first to undergo folding without significant interference from neighbouring layers.
In the gneissic rocks of Purulia, India, quartzo-feldspathic layers of varying thickness, embedded within a biotite-rich matrix, exhibit folded structures with different spacing (Figs. 10 and 11). Natural data collected from the broad folds of thicker layers and the residual small folds of thinner layers in the hinges of secondary folds (Table 2) show striking similarity with our numerical findings. We observed that the wavelength-to-thickness ratio decreases as the layer thickness reduces and as the spacing between two adjacent layers increases. The study reveals that this ratio for thinner unit are more close to the result of single-layer model (horizontal line in Fig. 12). Therefore, if the bulk shortening is known, plotting the graph of the wavelength-to-thickness ratio against interlayer spacing allows us to estimate the viscosity ratio between the layer and the embedding medium (Fig. 12).
Example of natural folding in gneissic rock, generated by the mutual interaction among the layers. Photographed from Purulia, India. (23.16oE, 86.28oN). Scale bar is 4 cm. Details of measurement are mentioned in Table 2.
Example of natural folding in gneissic rock, generated by the mutual interaction among the layers. Photographed from Purulia, India. (23.16oE, 86.28oN) Scale bar is 4 cm. Details of measurement are mentioned in Table 2.
Measurements from natural folds
Estimation of viscosity contrast from folded layers of different thicknesses. Horizontal line represents an ideal single layer. Square and triangular markers denote thick and thin layers, respectively. (a) η R = 100, (b) η R = 500 and (c) η R = 1000. Finite shortening is 30%. Note that the competency contrast can be estimated more precisely from the folds of thinner unit.
The data presented in the graph indicates that in layered rock system, higher-order folds in thinner units also provide a more precise estimation of rheological contrast compared to thicker units. In summary, the thinner layers are critical for determining the more accurate viscosity ratio in a layered system, which is in conjunction with the work of Sherwin and Chapple (Reference Sherwin and Chapple1968), done in single-layer fold system.
4. Limitations
In this research, we examined a simplified model that employs homogeneous materials and a three-layered configuration consisting of layers with two distinct thicknesses, yet exhibiting the same rheological contrast in relation to the surrounding medium. It is important to note that natural geological environments frequently feature layers with variable thickness, spacing, and viscosity, which can lead to significantly more complex and intricate fold geometries. Nevertheless, our study establishes a foundational understanding by concentrating on the evolution of fold geometries as a result of the interactions among layers within this specific system.
5. Conclusions
Based on the above observations and results, the following conclusions can be drawn:
-
i) In a layered rock system, the overall fold geometry is primarily influenced by the folding behaviour of the thinner layer, which perturbs the surrounding strain field and subsequently affects the folding of adjacent thicker layers.
-
ii) During progressive deformation, the thicker layers influence the folds in thinner layers, causing transitions from harmonic to polyharmonic, and eventually to disharmonic folds as finite strain increases.
-
iii) The contact strain zone has a significant impact on the development of fold geometries within a layered system.
-
iv) Apart from the contact strain zone, parameters such as viscosity ratio, relative layer thickness, and interlayer spacing are major determinants of the final fold geometry.
-
v) The geometry of higher-order folds in thinner layers provides the most reliable estimates for determining the viscosity contrast between the layer and the embedded medium.
Acknowledgements
We are grateful to Prof. Bruce Hobbs and an anonymous reviewer for their scientific comments that greatly contributed to the improvement of the paper. We sincerely acknowledge the constructive suggestions by Professor Olivier Lacombe. Masud Rana helped in the field study. The present work is funded by Ministry of Earth Science, Govt. of India through the project (MOES/P.O.(Geo.)/185/2018). Infrastructural facilities provided by the Department of Geological Sciences, Jadavpur University, is gratefully acknowledged.
Author statement
DBM: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Validation; Visualization; Roles/Writing – original draft; Writing – review & editing. SD: Data curation; Formal analysis; Investigation; Methodology; Validation; Visualization; Roles/Writing – original draft; Writing – review & editing. SKS: Conceptualization; Data curation; Formal analysis; Funding acquisition; Investigation; Methodology; Project administration; Supervision; Validation; Visualization; Roles/Writing – original draft; Writing – review & editing.
Competing interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open access
