## Abstract

Most natural fold systems are not sinusoidal in profile. A widely held view is that such irregularity derives solely from inherited initial geometrical perturbations. Although, undoubtedly, initial perturbations can contribute to irregularity, we explore a different (but complementary) view in which the irregular geometry results from some material or system softening process. This arises because the buckling response of a layer (or layers) embedded in a weaker matrix is controlled in a sensitive manner by the nature of the reaction forces exerted by the deforming matrix on the layer. In many theoretical treatments of the folding problem, this reaction force is assumed to be a linear function of some measure of the deformation or deformation rate. This paper is concerned with the influence of nonlinear reaction forces such as arise from nonlinear elasticity or viscosity. Localized folds arising from nonlinearity form in a fundamentally different way than the Biot wavelength selection process. As a particular example of nonlinear behaviour, we examine the influence of axial plane structures made up of layers of different mineralogy formed by chemical differentiation processes accompanying the deformation; they are referred to as *metamorphic layering*. The alternating mineralogical composition in the metamorphic layers means that the embedding matrix exerts a reaction force on the folded layers that varies not only with the deflection or the velocity of deflection of the layer, but also in a periodic manner along the length of the folded layers. The influence of this spatially periodic reaction force on the development of localized and chaotic folding is explored numerically.

## 1. Introduction

The theory of buckling of a layer, or layers, embedded in another material has been extensively studied in geology for almost 200 years since the classical work of Hall [1], although the concepts involved clearly go back to Coulomb and Euler. Interestingly, the folds produced by Hall in layers of cloth (figure 1*a*) are *localized* in the sense that they are not sinusoidal, and such irregularity has been widely recognized in the geological literature (figure 1) [2–12]. This paper is concerned with such departures from sinusoidal behaviour. Most studies of fold irregularity in the geological literature derive from the classic work of Biot in the period 1937–1984, some of which is brought together in Biot [13]. That body of work has had a profound impact on the way in which geologists view the folding process, but one should appreciate that it is based on small deflections of thin, inextensible layers and linear stability analyses employing linear constitutive relations and geometries. The results are strictly valid for infinitesimal deflections, but there is a tendency in the geological literature to extend the results of the theory to describe finite deformations of thick, extensible layers with nonlinear matrix materials [12] and to consider the folding process independently of processes operating in the softer matrix between the buckled layers, in particular (except for Kocher *et al.* [14]), the development of structures such as axial plane cleavages and foliations.

One would expect from the linearity of the Biot problem that fold systems at finite strains would be sinusoidal. However, individual natural fold systems show a range of wavelength-to-thickness ratios, ranging from about two to somewhere in the range 15–20 (table 1 in Hudleston & Treagus [12]). This observation is commonly interpreted as resulting from the Biot theory, but representing a relatively flat dispersion relation, so that the wavelength amplification process is not very efficient at selecting a dominant wavelength and is influenced by the statistics of initial geometrical perturbations. There is no doubting that initial geometrical imperfections of sufficiently large wavelength are capable of inducing an irregular buckling response in linear materials with no need to appeal to nonlinear material behaviour of some kind [5–10]. An analytical solution to this problem has been given for finite deformations by Mühlhaus *et al.* [15], pp. 228–231 for Newtonian viscous materials with constant-velocity boundary conditions. The analysis shows that initial perturbations with a wavelength much smaller than the Biot wavelength will grow slowly compared with initial imperfections with wavelengths close to or larger than the Biot wavelength. The larger wavelength imperfections result in localized or irregular fold systems. The issue is this: Are all irregularities in natural fold systems derived from initial imperfections or are there other ways of inducing such irregularity? Other mechanisms of inducing an irregular or localized buckling response involve the development of some form of nonlinearity in the geometry or constitutive behaviour of the materials involved, and this process is distinct from the Biot model [16,17].

An important question is this: In the absence of suitable initial perturbations, does nonlinear behaviour emerge at large deflections even in systems that are linear at small deflections? Biot specified boundary conditions of constant force parallel to the layer so that the layered system is loaded by the equivalent of a dead weight. Such a loading system leads to exponential amplification, but because the Biot result is applicable only for the instant that folds are nucleated, it is not clear that such an explosive growth rate holds for sinusoidal folding or whether some form of geometrical nonlinearity arises at large deflections, as indicated by Wadee [18]. Schmalholz & Podladchikov [19], using constant-strain-rate boundary conditions, place a constraint on the classical Biot theory and propose that the exponential growth is limited by softening within the layer due to stretching of the layer as the limb dip increases. Even so, although a bifurcation occurs in this system [20], it appears that no nonlinearity arises that leads to localization; we explore this further in §2. There is a possibility that, among the many variables and assumptions inherent in the small-deflection, thin-layer theories, there remains an aspect that induces non-sinusoidal behaviour in the absence of initial perturbations once thick layers and/or large deflections are taken into account, and we expand upon this possibility below.

Another fundamental aspect of the Biot theory is that it considers the force exerted on the folding layer by the embedding matrix to be a linear function of the deflection of the layer (for elastic materials) or, for viscous materials, a linear function of the rate of deflection [13,17,21], and it is this relationship that is the focus here. The purpose of this paper is to consider nonlinear reaction forces exerted by the embedding material on the folding layer, including reaction forces that arise from a periodic or quasi-periodic distribution of constitutive parameters within the embedding matrix. In particular, the aim is to consider the influence of new layers that grow in the matrix as a result of metamorphic mineral reactions coupled with deformation and that develop parallel or approximately parallel to the axial plane of the folding layer. Some examples of such metamorphic layering are given in figures 1*d* and 2 and towards the end of §5.

Thus, departures from sinusoidal folding can conceivably have two origins other than the growth of initial imperfections:

— The Biot theory, when extended to include large deflections and thick plates, may introduce nonlinearities that cause the system to localize or the dispersion curve may broaden with increased strain as proposed by Mühlhaus [22] and Schmalholz & Podladchikov [19].

— Nonlinear behaviour of the embedding matrix may be responsible for localization as has been well documented in elastic, plastic and viscous materials by Tvergaard & Needleman [16] and in a large number of papers by Hunt, Budd and co-workers [21,23–31].

We explore both possibilities here. We first discuss the concept of localized folding (§2) and draw out the distinctions between Biot folding and localized folding. We then revisit (§3) the classical linear formulation of Biot and summarize the results of subsequent analytical analyses, involving no approximations, of finite deformations of thick layers as a reminder that the present formulations of the problem, in the absence of long-wavelength initial irregularities, give rise to strictly periodic solutions. We then extend the treatment to numerical examples (§4) involving nonlinear reaction forces exerted by the embedding material, in particular, some of those softening relations considered extensively by Hunt [23]. Then follows a brief discussion of the development of metamorphic layering from the point of view of reaction–diffusion theory (§5). This sets the scene for a numerical examination of reaction forces that vary periodically along the length of the folded layer (§6). As expected from the work of Hunt and co-workers [23,24], the result is a rich set of localized and chaotic responses. The notation used in this paper is summarized in table 1.

## 2. Localization of folding

Consider a single layer of linear elastic, viscous or visco-elastic material embedded in a matrix of weaker material. Biot [13] showed that, if this system is shortened, then the initial flat (but with small perturbations) configuration of the layer is unstable and folds develop provided the contrast in mechanical properties between the layer and matrix is large enough. Our aim in this discussion is to understand the form of deflections that ultimately arise in this layer when shortened parallel to its initial orientation by a constant force or constant velocity at its ends. We make no excuse for considering only a single layer. Our approach is to understand the simplest of geometrical situations before proceeding to more complicated multi-layered geometries. Our main goal is to investigate the influence of nonlinear behaviour in the matrix on the development of folds. We also want to understand the influence of any initial geometrical irregularities on the form of the folds that develop with a nonlinear matrix. Again, in order to advance to this level of complication, we need first to understand the situation with well-prescribed, simple geometrical imperfections.

### (a) Types of softening

Before proceeding, we need to be a little more specific about what is meant by a nonlinear softening response to deformation. Many systems show nonlinear softening effects as they deform; such effects may be material or geometrical in origin or arise from the boundary conditions. The material nonlinearity arises because of the effective or real degradation of some constitutive parameter and may or may not be associated with a strain-softening response. Thus, in the elastic range, many materials show an initial linear stress–strain relation followed by a nonlinear but still hardening response [33]. This arises through a degradation of the elastic moduli and is sometimes expressed in one dimension as a relation of the form [17]
2.1where *σ* is the axial stress, is a material constitutive parameter such as the shear or bulk modulus, *n* is a number commonly in the range 1–4 and *ε* is the axial strain. Such a relation is equivalent to a degradation of the material parameter from an initial value according to
2.2Discussions of the origin of such nonlinearity are widespread; a thermodynamics-based discussion based on damage accumulation with comparisons with experimental work is given by Lyakhovsky *et al.* [33] and Lyakhovsky & Ben-Zion [34].

Nonlinear viscous materials can show an effective degradation of the viscosity. A common example is a power-law viscous material represented in one dimension by a constitutive law of the form
2.3where is the axial strain rate, *N* is a number greater than 1 and is a material parameter. These materials are effectively anisotropic in that they can be described by two different viscosities defined by
2.4where *η*_{T} is the tangent viscosity and *η*_{S} is the secant viscosity. These are referred to as the normal viscosity and shear viscosity by Smith [35], who shows that *η*_{T} is relevant to pure shearing deformation histories, while *η*_{S} is more relevant to simple shearing deformation histories. Thus, the materials behave as anisotropic materials, while always remaining isotropic at any instant during the deformation [5]. For *N*>1, *η*_{S}>*η*_{T} and both viscosities soften with increasing strain rate. These materials are strain-rate hardening, but viscosity strain-rate softening.

A second form of material softening arises through strain softening, where the stress supported by the material decreases as the material continues to deform. This behaviour has been extensively studied in the field of nonlinear elasticity, with special reference to martensitic transformations [36,37], although, in the past decade, the applications have extended to plastic behaviour [38–40]. This is the behaviour studied in the classic paper by Tvergaard & Needleman [16], where it is shown that sinusoidal buckling is the response at peak load. As the system softens, a bifurcation occurs *after* peak load that corresponds to a transition from sinusoidal to localized folding.

There is also a class of materials in which material strain-rate softening occurs [41–44]. This behaviour occurs when some diffusive process (thermal or chemical) competes or interacts with other deformation mechanisms such as strain hardening arising from dislocation interaction to produce a decrease in the effective viscosity with an increase in strain rate. Examples include the influence of heat dissipated by deformation [45–47] or by chemical reactions on the viscosity [41]. Such strain-rate weakening processes lead to an effective degradation of material properties coupled with the deformation and consequent localization [43,44].

A common form of geometrical softening relevant to structural geology arises from the rotation of a plane of slip or of anisotropy during deformation so that its normal tends to align with the principal axis of compression. This form of weakening has been extensively studied with respect to kink and chevron fold formation by Hunt *et al.* [27] and to sequential fold development by Budd *et al.* [31]. Geometrical softening has also been studied with respect to crystal plasticity by Ortiz & Repetto [38] and Ord & Hobbs [48]. An expression for the stress supported by an anisotropic viscous material with respect to the orientation of the anisotropy relative to the deformation-rate tensor is given by Johnson & Fletcher [5], p. 422 and shows that the stress supported by the system softens or hardens depending on orientation as it rotates during deformation.

Thus, we recognize three different types of nonlinearity (table 2):

— Material nonlinear behaviour. We emphasize that nonlinear in this context does not necessarily imply strain or strain-rate weakening. For instance, the nonlinear elastic behaviour employed by some workers [17] involves experimentally well-documented softening of the elastic moduli, but the stress–strain curve is hardening [17], fig. 9.

— System nonlinearity arising from geometrical and structural effects or internally imposed constraints.

— External nonlinearity arising from the externally imposed boundary conditions.

### (b) Localized folding

In the literature on folding of layered materials, two different types of deformation responses are identified. We call these two kinds of behaviour *Biot folding* and *localized folding*. The first originates from the work of Biot [13,55], where layers of a material with linear constitutive behaviour (elastic or Newtonian viscous or combinations of the two) are inter-layered with material (which we call the matrix or embedding medium) with different linear constitutive behaviour. In these models, the force parallel to the initial orientation of the layers remains constant so that, at least for small deflections, *w*, the whole system remains linear. No softening behaviour arising from constitutive or geometrical origins or from boundary conditions is considered. The result is that the governing equation, for small deflections (by small is meant limb dips of less than *ca* 10^{°}, see §3*a*), is
2.5where *α* is a function of the constitutive behaviour of the layers and *F*(*w*) is a function that represents the reaction force on the layer from the matrix as the layer deflects; *x* is the distance measured along the layer from some reference point. The layer deflection, *w*, is a function of both *x* and time, *t*. For the linear problem described above, *F* is always a linear function of *w* and depends also on the wavelength of the deflection [55]. One can readily confirm that, if *F* is linear, one solution to (2.5) is
2.6where *k* is the wavenumber for that deflection, *k*=2*π*/*λ*, where *λ* is the wavelength of the deflection, and *A*(*t*) is a function that describes the way in which the amplitude evolves with time. Biot [13] shows that , where *w*^{0} is some initial deflection of the layer and *ω* is the amplification rate of the deflection that arises. Because the problem is linear, the actual solution to (2.5) can be expressed as a Fourier series that is a linear combination of an infinite number of solutions of the form of (2.6) with different values of *w*^{0} and *ω* corresponding to each wavenumber:
2.7Biot [13] showed that the dispersion curve (the plot of *ω*_{n} against *k*_{n}) is bell-shaped, skewed towards smaller wavenumbers, and sharply peaked so that one wavelength grows faster than all others; this wavelength is called the dominant wavelength. This wavelength selection process, together with (2.6), means that the folds that grow according to this theory are always sinusoidal as long as (2.5) holds. We refer to this first type of folding as *Biot folding*. Fletcher & Sherwin [2] proposed that departures from a sinusoidal form that are commonly observed in natural folds arise from a broad shape for the dispersion curve, meaning that the spread of growing wavelengths is large and the wavelength selection process is not very efficient. We explore this proposition below.

The second type of folding behaviour identified in the literature follows naturally as a generalization of Biot folding and arises if the reaction of the matrix to deflection of the layer is no longer linear, as occurs if the matrix softens in some manner. Then *F* is no longer a linear function of *w* and (2.6) is not a solution to (2.5). The result is that a linear combination of terms such as (2.7) can no longer be written as a solution to (2.5). This means that, in general, a wavelength selection process based on preferred growth of Fourier components does not operate, and the process is quite different from Biot folding. The solutions to (2.5) are now non-sinusoidal or localized along the layer, as shown in figure 3, and very sensitive to the magnitude of *α* in (2.5) and the form of *F*; many examples are given in the book devoted to the solutions of (2.5) by Peletier & Troy [56]. This is known as *localized folding*. As an example, if *F*(*w*)=*w*−*w*^{2}, then Champneys & Toland [57] show that there are an infinite number of localized solutions to (2.5). If *F*(*w*)=*w*−*w*^{3} [21], then some of these localized solutions are given by
2.8and an example of the time evolution of this solution together with the time evolution of the axial force for constant-velocity boundary conditions is given in figure 3. The values of the parameters *ς*_{1},*ς*_{2},*ς*_{3},*ς*_{4},*ϑ*_{1},*ϑ*_{2} (all of which vary with time) are taken from fig. 5 of Hunt *et al.* [21]. Notice that *F*(*w*)=*w*−*w*^{3}/6 is an approximation to .

If packets of folds interfere with one another, then *chaotic folding* can arise [25]. An important distinction between the processes of Biot and of localized fold formation is that, in the Biot theory, all folds form at the same time and all folds with the dominant wavelength grow at the same time at the same rate. In contrast, localized folds form at different times so that there is a sequential development of folds; some folds may form only to collapse later in the deformation history or folds can grow in what were formerly planar regions [29–31]. The Biot wavelength selection process does not occur.

Although there is a fundamental difference in the processes that operate during the development of Biot and localized fold systems, Biot folding can be (but is not always) the precursor for localized folding. Thus, Tvergaard & Needleman [16] show that, for softening plastic materials under constant loading rate, the stress in the layer first rises to a maximum where sinusoidal folds develop; this we identify with Biot folding. Further, Tvergaard & Needleman showed that, as the stress within the plastic layer decreases due to softening, the system behaviour bifurcates and a localized fold system develops, only to be replaced later in the deformation history by a sinusoidal fold system with a wavelength different from that at peak load. An example where a localized system evolves towards a sinusoidal system is given by Knobloch [58]. Many examples of this sequential development of localized folds are given by Budd & Peletier [29] and Budd *et al.* [30,31]. This kind of behaviour is explained, in the spirit of Thompson & Hunt [59], by the system seeking, at any particular time, the lowest energy configuration. In linear systems there is but one energy minimum and this corresponds to a sinusoidal solution. For nonlinear systems the energy landscape has many or perhaps an infinite number of local minima and the system jumps from one localized solution to another as deformation proceeds [60]. The energy landscape may even be fractal [25]. Occasionally, the local minimum also corresponds to a new sinusoidal solution. At present, the mathematics that enables analytical descriptions of these processes has not been developed [58]. Although analytical solutions may sometimes be obtained, they may not be unique or stable and an infinite number of other solutions may exist. This is particularly true for solutions obtained by some form of perturbation analysis. A common approach [58,60] is to explore the phase space appropriate to the system and search for bifurcations using software such as AUTO [61] or XPPAUT [62]. Examples of localized folding are presented in §4.

## 3. Biot's theory: a linear matrix with sinusoidal folds

This section explores some extensions to the Biot solution to finite deflections of linear materials and systems to see if some form of nonlinearity develops that could initiate fold localization. It is clear that, for linear systems, initial geometrical perturbations with wavelengths greater than the Biot dominant wavelength can initiate localized fold growth [6–9,15]. The question we address here is this: Is there an intrinsic nonlinearity in the folding of linear materials, in the absence of these kinds of perturbations, that results in fold localization? Although there is a considerable literature that proposes ‘exact’ analytical solutions to the problem of ‘large’ deflections of ‘thick’ plates (see Hudleston & Treagus [12] for a review of some of this literature), such literature needs to be examined with considerable care. For instance, Johnson & Fletcher [5], pp. 196–207 and 224–236 report *exact* solutions for large deflections but with the assumptions that the functions cos(*x*) and exp(*x*) are truncated at the second or third terms. Thus, the *exact* theory is for *approximate* formulations. Unfortunately, however, some of these approximations only hold for small deflections (that is, (2.5) is assumed either explicitly or implicitly), so their use for large deflections is not valid. Johnson & Fletcher [5] express the boundary conditions in terms of deformation rate, instead of the constant-load boundary conditions of Biot, with the (unstated) assumption that the load is constant to large deflections. However, they fail to include the kinematic constraint (3.11) for constant velocity or deformation rate in the governing equations. We will see that the inclusion of this constraint introduces a geometrical nonlinearity that makes a considerable difference to the outcome. Any attempt to treat large deflections using the ‘biharmonic equation’ (2.5) to express the layer displacement in terms of distance along the layer automatically restricts the problem to low amplitudes (see §3*a*). Any attempt to treat the ‘thick plate’ theory must specify whether Timoshenko, Kirchhoff, Reissner or some other constraint for thick plates is assumed so that the strain tensor and hence the stress state within the layer are constrained. Hudleston & Treagus [12], p. 3 claim that ‘Thick plate theory does not lead to an explicit expression for the dominant wavelength’, whereas such an exact expression, with no approximations, for thick Kirchhoff and Reissner plates at large deflection for any (three-dimensional) stress state and for power-law viscous layers can be obtained from the available literature [63], eqns (4.22) and (6.2). Any attempt to use (2.5) with a linear form of *F*(*w*) with a nonlinear matrix material needs to be justified since, in general, *F*(*w*) is nonlinear. The situation is quite confused in that no distinction is made between linear and nonlinear systems and nonlinear systems are treated as though they are linear. Although a general solution to the problem of large deflections of thick plates of linear viscous materials is still not available, we hope that the following discussion may add some clarity.

The classical formulation of Biot is extended below to include the results of analytical analyses involving *finite deformations with no approximations* in order to address the following two specific sets of questions:

— What is the nature of the buckling dispersion relation at large deflections? Is it significantly different from that derived for small or approximate deflection theory? Does the dispersion relation for large deflections suggest that a broad range of wavelengths and amplitudes should grow or is the result strictly sinusoidal waveforms? Do geometrical nonlinearities emerge that lead to non-periodic behaviour?

— For thick plates and constant-force boundary conditions, does the growth rate remain exponential to large deflections? And what implications does this have for the applicability of these boundary conditions for natural deformation conditions? How is the resulting geometry of buckling influenced by adopting constant-strain-rate or constant-velocity boundary conditions?

Although aspects of these questions have been explored by previous work [15,19,20,22,63], some of that work is difficult to access; the intent of this section is to bring such work together into one place and to present analytical results for fold evolution involving thick plates to large deflections with different boundary conditions as a foundation for discussions of localization in natural fold systems. We restrict the discussion to isotropic materials, although, as Mühlhaus *et al.* [63] point out, the treatment of power-law viscous plates automatically solves the folding problem for orthotropic plates. Anisotropic matrix materials immediately introduce a geometrical softening and so are beyond the scope of the linear Biot theory. A discussion of the folding of anisotropic linear viscous materials to finite deflections together with an analytical solution for the amplification rate is given by Hobbs *et al.* [64]. Numerical solutions show that localized folding is a common response to deformation for both anisotropic and power-law viscous matrix and layer materials [14,49,64].

### (a) Geometry and kinematics of the folding problem at large deflections

In discussions of the finite deflection of thick plates, the geometry shown in figure 4 is relevant. Since we are only dealing with geometry, the discussion is independent of the constitutive behaviour and is relevant to elastic, viscous and plastic materials or any combinations of these. In the undeformed state with coordinates (*X*_{1},*X*_{2}), an element of the layer, of length d*S*, with unit normal to the mid-surface, , and thickness, *H*, is shown. In the most general situation relevant to geological deformations, this element in the deformed state with coordinates (*x*_{1},*x*_{2}) has a new thickness, *h*, and length, d*s*, and is rotated through an angle, *φ*, while the normal is deformed to become the vector inclined at *ϕ* to *x*_{2} and with a deformed length .

The general expression for the evolution of a buckling elastic plate is given [17,59] as
3.1where , *s* is the distance measured along the layer, , *P* is the force driving the buckling, *R* is a quantity related to the rheology and *F*(*w*_{2}) is related to the matrix reaction force. This equation is also applicable to the buckling of Maxwell visco-elastic materials as discussed by Hunt *et al.* [26].

The folding theory developed by Biot in the period 1937–1984 [13,55,65] has been incredibly influential within the structural geology community. Biot's theory, as spelled out by himself and a host of subsequent workers, employs a number of simplifying assumptions. The theory is formulated for small deflections of a *thin*, *inextensible* layer so that , d*S*=d*s*, *H*=*h* and *φ* is small, so that if *w*_{2} is the deflection normal to the initial orientation of the layer, then [59] the quantities and nonlinear terms of similar magnitude can be neglected. This means that the general equation for a deforming plate (3.1) reduces to the well-known equation [5,13,66,67]:
3.2where *α* is a function of both the rheological properties of the layer and the compressive force in the layer and *F*(*w*_{2}) is related to the reaction force exerted on the layer by the matrix as the layer deflects. Since , the assumption that can be neglected with respect to 1 restricts the use of (3.2) to quite small deflections (say *φ*≈10^{°}). The small deflection theories also neglect higher-order terms for the strain relative to the deformed state. For finite deflections, the axial strain is given by Fung [68] as

This is the expression used by Mühlhaus *et al.* [15] to explore the influence of constant-velocity boundary conditions on fold amplification. Wadee [18] points out that inclusion of high-order terms past ∂*w*_{1}/∂*x*_{1} for *ε*_{11} can be responsible for localized behaviour, but the analysis of Mühlhaus *et al.* [15] does not lead to this conclusion for linear viscous plates. Other assumptions are that the distortion of the plate is plane and that the thickness of the layer is small compared with the wavelength of buckling. Note that (3.2) is a stationary state for the Swift–Hohenberg equation [56,58,60,69], which is a special form of a reaction–diffusion equation and well documented as having sinusoidal, modulated, localized and chaotic solutions [56,70] depending on the value of *α* and the nature of *F*(*w*_{2}).

The Biot theory assumes that the loading corresponds to a dead weight, meaning that the force initially parallel to the layer and that drives the buckling instability remains constant throughout the deformation. Biot [55] also derived a relation for the reaction force, ** q**(

**), exerted by a linear embedding matrix on the deforming layer. For a linear elastic embedding medium with shear modulus,**

*w**μ*

^{e}, this force is a linear function of

**and of the wavenumber,**

*w**k*

_{n}, for the deflection, where

*k*

_{n}=2

*π*/

*λ*

_{n}and

*λ*

_{n}is the corresponding wavelength: 3.3For a linear viscous embedding medium with viscosity

*η*

^{e},

**is given by 3.4where the overdot denotes differentiation with respect to time,**

*q**t*, and

*w*

_{1},

*w*

_{2}are the components of the layer displacement,

**(figure 4). Note that (3.3) and (3.4) are not applicable if**

*w**μ*

^{e}and

*η*

^{e}are nonlinear functions of

**or , respectively.**

*w*There are three common approaches [68,71] to the buckling of thick layers (figure 5). Specification of the relevant plate model is essential to constrain the stress state in the layer. Although the special forms of plates discussed below had their origins in describing the deformation of elastic plates, the discussion is purely geometrical and can be applied to plates composed of any materials. The important point is that each model defines a deformation from which the stresses may be derived once the constitutive behaviour is defined. We use these models below to describe the deformation of elastic, plastic and viscous plates. The *Timoshenko* model of buckling assumes that the beam is extensible and thick (the thickness remains constant). The normal to the mid-surface does not remain normal to this surface and is extensible. This means that *φ*≠*ϕ*, , *H*=*h* and d*S*≠d*s*. A special case of this theory where *φ*=*ϕ* and corresponds to the *Kirchhoff* model. The case *φ*≠*ϕ*, , *H*≠*h* and d*S*≠d*s* corresponds to the *Reissner* model. In all three theories, initially plane sections through the layer remain plane after deformation. The Timoshenko and Kirchhoff models produce concentric flexural shear folds and concentric buckle folds, respectively, whereas the Reissner model is capable of producing class 1A, 1C and 3 folds of Ramsay [66] and also similar folds as a special case. In the most general case, and that most relevant to geological folding, the deformation also involves the non-affine distortion of plane sections initially normal to the central layer together with inhomogeneous extension of and changes in thickness, so that we have *φ*≠*ϕ* and *ϕ* changes across a profile of the fold, , *H*≠*h* and d*S*≠d*s*; this clearly demands constitutive relations that define the evolution of and *h*. To date, we know of no treatment of this general case. Mühlhaus [22] develops the Kirchhoff model for large values of *φ* for an elastic layer embedded in a Newtonian viscous medium and assumes Biot's linear relations (3.3) and (3.4) for the reaction force. Mühlhaus *et al.* [15,63] develop the Kirchhoff and Reissner models for finite *φ* for elastic and viscous plates embedded in a linear viscous medium and for viscous materials in three dimensions. Although Mühlhaus [22] focuses on an elastic plate in a viscous medium, much of the development for finite deflections of thick plates in that paper is general and independent of constitutive relations. Schmalholz & Podladchikov [19] use a Kirchhoff model, but with an approximation to a thin plate.

Various authors have made steps along the path to generalizing the strict Biot theory. These studies constitute a progression that involves more and more constraints on the strain tensor and its variation within the layer and hence the stress state once the constitutive behaviour of the layer is defined. The strict Biot theory proposes that the layer is inextensible. Sherwin & Chapple [72] continued the thin-plate assumption, but dropped the inextensible assumption. Fletcher [32], and many other papers that are summarized in Hudleston & Treagus [12], together with Johnson & Fletcher [5], approached the thick-plate, large-deflection theory with a series of approximations. Those studies resulted in definition of the dispersion curve for the approximations involved and some understanding of the ways in which the dominant wavelength and amplification are influenced by increasing the layer thickness and by considering large values of *φ*. One outcome, for instance, is that the amplification remains exponential but is slower than the strict Biot result [12]. Schmalholz & Podladchikov [19] extended the process of generalization a step further but reverted to a thin plate, where the normal stress in the plane of the layer and parallel to the mid-surface decreases as the layer deflection increases. This has the effect of decreasing the growth rate below exponential for constant-strain-rate boundary conditions. Although the dispersion curve broadens at high deflections, there is no indication that a nonlinearity develops that results in fold localization. This type of constraint differs from the above models, which constrain only the deformation in the layer, in that it places a constraint on how the stress evolves within the layer.

Many of these studies assume constant-force boundary conditions. Mühlhaus [22] and Mühlhaus *et al.* [15,63] extended the studies to thick Kirchhoff and Reissner plates with both force and velocity boundary conditions. This series of studies starting with Sherwin & Chapple [72] represents a steady progression to define a general strain state within a layer (that is, the components of the Green tensor defined in (3.5) below), but there is still some way to go before a general theory is available. The outcome has been an understanding of why naturally observed folds have *λ*/*h* values different from that predicted by the strict Biot theory, an understanding of controls on amplification by the stress state in the layer and by boundary conditions, and the controls on fold shape by initial perturbations at least for constant-velocity boundary conditions. However, no irregularity develops and no nonlinear behaviour emerges, so that the folds that develop for linear matrix responses are sinusoidal with no localization. The only influence on the development of non-sinusoidal folding for these models is that arising from the influence of geometrical perturbations.

Although there is a hint [22] that coupling between the rates of stress diffusion and fold growth rates might lead to nonlinear behaviour in visco-elastic thick-layer systems, it turns out that the different time scales involved in these two processes result in decoupling of the processes. Strictly, periodic behaviour remains throughout the folding history, although there is some broadening of the dispersion curve. For some combinations of elasticity and viscosity, two wavelengths of folds develop [63,73] and the dominant wavelength may differ from the Biot wavelength.

In general, the reaction force exerted by the matrix on the deforming layer is expressed as a relationship (not necessarily linear as in (3.3) and (3.4)) between the traction, ** q**, on the surface of the layer and the corresponding surface displacement rates, , for a viscous matrix or the displacements, (

*w*

_{1},

*w*

_{2}), for an elastic matrix. Mühlhaus [22] assumes quite generally that the deformation of the layer is given as a function of time by

**=**

*x***(**

*x***,**

*X**t*) and that the plate deformation can be completely described in terms of two geometric quantities, namely, the plate mid-surface with coordinates (

*ξ*

_{1},

*ξ*

_{2}) and the field of the deformed initial normals, , to this surface, where . As we have seen,

*ϕ*=

*φ*corresponds to the Kirchhoff constraint, while

*ϕ*≠

*φ*corresponds to the Timoshenko constraint. The deformation of the plate is given by and so the deformation gradient is given in terms of generalized coordinates by where the primes denote differentiation with respect to

*x*

_{1}. Mühlhaus [22] shows that the Green–Lagrange strain tensor,

*γ*

_{IJ}, is then given by 3.5where . Thus, for large deformations,

*w*

_{2}(or its derivative with respect to time) is a suitable measure of the deformation to define the softening of the reaction forces so long as the plate is inextensible. Thus, for more general deformations, a dependence of

**on measures of the deformation other than**

*q**w*

_{2}is necessary. In this paper, we restrict attention solely to a dependence of

*q*

_{2}on

*w*

_{2}or , but future work should explore other approaches to defining

**(**

*q**D*

_{ij}), where

*D*

_{ij}is the stretching tensor.

### (b) Finite deformations of a thick power-law viscous layer

Mühlhaus *et al.* [63] derive analytical solutions with no approximations for the finite deformation, in three dimensions, of a viscous power-law plate of thickness, *h*, embedded in a linear material; this means that the reaction force from the embedding material is linear as derived by Biot [55]. Note that, for this treatment, the matrix is a linear viscous material; if the matrix is nonlinear viscous, then *F*(*w*) is nonlinear and localized folding is expected. This is outside the scope of the linear Biot theory.

The system is described by a Cartesian coordinate system (*x*_{1},*x*_{2},*x*_{3}) with *x*_{3} initially normal to the deforming plate. The shear strain rate, , is given in terms of the stretching tensor, *D*_{ij}, as and the viscosity, *η*, of the layer is , where *η*_{0} is a reference viscosity at a reference strain rate, , and *N* is a power-law exponent; *N*=1 corresponds to Newtonian viscosity; *η*^{e} is the corresponding viscosity in the embedding material. The normal tractions in the plane of the plate are given in terms of the Cauchy stress, *σ*_{ij}, as
Both Kirchhoff and Reissner models are examined, with broadly similar results for both. Mühlhaus *et al.* [63], p. 3141 enumerate the departures from the strict Biot results. Using an incremental approach, they arrive at (3.6), which governs the evolution of the finite deflection *w*(*x*_{1},*x*_{2},*t*) of the plate. Note that this differs considerably from (3.2) used in the Biot theory (see, however, [21,25,74]) in that the time derivative of the deflection appears:
3.6where *q* is the reaction force of the embedding material on the layer, and
The *c*_{ij} are material parameters and are also functions of the deformation-rate tensor:
A Fourier transform of (3.6) yields
3.7where
and *k*_{i} are the wavenumbers in the *x*_{1} and *x*_{2} directions, respectively. A linear stability analysis of (3.7) gives the general solution of (3.7) as
in which the dispersion relation is
3.8Introducing the dimensionless quantities
where
the dimensionless form of the dispersion relation becomes
3.9with
where
Note that, for a Newtonian viscous material (*N*=1), any dependence on the nature of the deformation-rate field vanishes although there remains a dependence on the tractions *T*_{1} and *T*_{2} through . For power-law materials, the nature of the deformation-rate field corresponding to pure stretching, plane stretching or some other form of stretching becomes important.

These equations enable the dispersion relations and hence the dominant wavelength at large deflections to be calculated for any power-law exponent, *N*, and stretching field, *D*_{ij}. In order to compare directly with the other results of the Biot theory, we assume Newtonian behaviour (*N*=1) and plane straining, *D*_{22}=0, *T*_{2}=*T*_{1}/2. Then the dispersion relation (3.9) becomes, for ,
3.10Note that, for this case, . If one puts , then the dominant wavelength for these conditions is given as or, in other words, the dominant wavelength is the Biot wavelength. Similar relations can be obtained for any stress or stretching state and value of *N* from (3.9). The dimensionless dispersion relation (3.10) is plotted in figure 6*a*. Also tabulated on that diagram are the wavelength-to-thickness ratios, *λ*/*h*, for a given and value of . The range of values of *λ*/*h* for natural folds (2 to 20) occurs on the left-hand side of the dimensionless dispersion curve and some of these values for are shaded. Thus, the large-deflection thick-plate theory indicates that, for the range of *λ*/*h* observed for natural folds, amplification factors are large only for values of , say, 100 to 1000. This is much more restrictive than the results obtained by Johnson & Fletcher [5], fig. 5.5, where values of *λ*/*h* in the range 10–30 for have large amplification factors. This is emphasized in figure 6*b*, which is a plot for Newtonian viscosity and of the relevant measure of the amplification against *λ*/*h* from Hudleston & Treagus [12], fig. 5a, based on Fletcher [32], compared with the results for finite buckling of a thick layer with no approximations [63]. For the approximate theories, the quantity *q*_{F}, which is the label we give the measure of amplification adopted by Hudleston & Treagus [12], is plotted for comparison with the quantity , where *ω*_{k} is the amplification of Mühlhaus *et al.* [63]. It appears from (7) of Hudleston & Treagus [12] that is the measure of amplification equivalent to (*q*_{F}+1). Figure 6*b* indicates that the amplification rates for the thick-plate finite-amplitude theory with no approximations [63] are much smaller than those for the approximate theories [12]. Importantly, the range of *λ*/*h* spanning the regions of maximum amplification are restricted in the finite theory relative to the approximate theory.

### (c) Influence of boundary conditions

Mühlhaus *et al.* [15] have calculated the dispersion relation for large deflections for a Newtonian viscous layer embedded in similar material in two dimensions for constant-velocity boundary conditions. Different from most treatments of the buckling problem, the kinematic constraint expressing the constant velocity is solved in conjunction with the governing equation describing the deflection of the layer. If *L* is the undeformed length of the layer, and *L** is its component parallel to *x*_{1} in the deformed state, then [59] and constant-velocity boundary conditions are expressed as
3.11Given this constraint, the dispersion relation at the initiation of folding is obtained as
3.12where *P* is the initial value of the compressive force parallel to the layer. Note that (3.12) is identical to the solution for constant-force boundary conditions at finite deflections and three-dimensional thick plates for *k*_{2}=0 (3.9). However, for constant-velocity boundary conditions, *ω*_{k} decreases with time since the kinematic constraint (3.11) means that *P* must decrease during deformation. The time history of the deflection *w*_{2n} can be obtained from
3.13A plot of against dimensionless time, , is given in figure 7*a*, where . The dimensionless quantities in (3.13) are given by and , where *λ* is the wavelength, and *λ*_{Biot} is the Biot wavelength given by *λ*_{Biot}=2*π*/*k*_{dom}. Note that, in figure 7*a*, the curve for , the Biot wavenumber, coincides with the Biot solution for small values of , as it must.

Equation (3.11) places a constraint on the way *φ* can evolve with time. This is another form of constraint analogous to the constraint derived by Schmalholz & Podladchikov [19] for the evolution of stress within the layer and also results in a growth rate that is non-exponential.

### (d) Influence of initial perturbations

Mühlhaus *et al.* [15] present an analytical solution for the finite growth of any initial geometrical perturbation of any wavelength and amplitude under constant-velocity boundary conditions in Newtonian materials. Results are given for an initial spike in the layer but any initial perturbation geometry could be investigated. The influence of an initial perturbation in producing a localized fold increases as its wavelength becomes greater than the Biot wavelength. The dependence on viscosity ratio is illustrated in figure 7*b*. The dispersion function for the growth of a perturbation in a viscous layer undergoing finite deflections with constant-velocity boundary conditions is (from eqn. (41) of Mühlhaus *et al.* [15]):
The results are plotted in figure 7*b* for together with the range (2–20) of *λ*/*h* values tabulated by Hudleston & Treagus [12] for natural folds. Each curve has a maximum at the Biot *λ*/*h* value. For the range of plotted, the maximum growth rate of perturbations coincides with the expected range of *λ*/*h* for , although the growth rate is always small. Large amplifications exist in this range even for large values of . This analytical solution indicates that for constant-velocity boundary conditions fold irregularity due to the presence of initial perturbations is expected if the *λ*/*h* of the perturbation is close to the Biot *λ*/*h* value for , but the effect is present to higher values of . One should emphasize that these conclusions hold for linear viscous materials. There is no analytical solution for power-law materials.

In answer to the questions raised near the beginning of §3, we can state the following:

— The dispersion relation at large deflections for power-law viscous thick plates in a general deformation-rate field,

*D*_{ij}, and for a Kirchhoff model is given in dimensionless form in (3.9) and plotted for plane straining in figure 6*a*. Figure 6*b*shows that the amplification is significantly smaller for the finite-amplitude solution and linear viscosity [63] than for the approximate theories [12], as is also the wavelength-to-thickness spread at high amplifications. Only for large values of viscosity ratio, say , would one expect a significant spread of*λ*/*h*and, for more reasonable values (say , the dispersion curve suggests the buckling response to be strictly sinusoidal. In the buckling theories developed to date involving both Kirchhoff and Reissner models, at large deflections the behaviour remains linear and geometrical nonlinearities or mode broadening due to nonlinear behaviour do not arise. Schmalholz & Podladchikov [19] have proposed a thin-plate model where the stress within the layer decreases as the limb dip increases. This results in a decrease in the amplification rate but does not result in mode broadening except at large deflections.— For Kirchhoff and Reissner models, the growth rate remains exponential to high deflections for constant-force boundary conditions. The introduction of constraints such as those proposed by Mühlhaus

*et al.*[15] or by Schmalholz & Podladchikov [19] will modify all of these behaviours in that the growth rate becomes non-exponential. Of course, the question arises as to how relevant constant-force boundary conditions are for natural deformations. The answer is open but constant-velocity boundary conditions are easier to envisage for natural deformations than dead-weight types of boundary conditions. The Schmalholz & Podladchikov [19] constraint leads to non-exponential growth because the force in the layer decreases with increased deflection. Under constant-velocity (or constant-strain-rate) conditions, the dominant wavelength gradually increases as the deformation continues, the force in the layer asymptotically approaching zero with continued deformation even for small limb dips [15,21]. The deformation approaches homogeneous shortening, so that early formed folds are passively amplified. This is the type of deformation behaviour necessary to produce class 1C folds. Such a scenario is in contrast to the Biot constant-force scenarios. The question remains as to the thick-plate model of greatest relevance to deformed rocks. The important point is that, of the models so far developed assuming linear reaction forces,, nonlinear behaviour leading to fold localization does not arise.*q*

The outcome of this brief review of thick-plate, large-deflection theories for layers made up of elastic and viscous (including power-law) materials, for linear matrix materials and for constant-stress and constant-strain-rate boundary conditions is that a sinusoidal response remains as the only response unless there are large-wavelength initial geometric perturbations in the layer. Although there are hints [22] that nonlinear behaviour will arise, no such behaviour emerges and we are left with exploring a nonlinear response of the embedding medium in order to see constitutive behaviour appear as an influence on fold geometry. The possibility remains that more general thick-plate models might lead to non-sinusoidal solutions for a linear matrix response but such work remains to be done.

## 4. Localized folds with a softening matrix

As discussed in §2, if the reaction force, ** q**, is a nonlinear function of the displacement or displacement rate and, in particular, if the behaviour is softening, then it is well documented [16,23,29–31,56] that the buckling response is generally non-sinusoidal, so that the folds localize at one or more places along the layer. If individual packages of folds interfere with each other, the fold system can become spatially chaotic. Such systems are characterized by a wide range of amplitudes and wavelength-to-thickness ratios; in some instances, the systems are fractal [25,28]. Fold amplification arising from some form of matrix softening is intrinsically different from the Biot type of folding behaviour in that the solutions to (3.2) are localized rather than sinusoidal and no wavelength selection process based on amplification of specific parts of a Fourier series operates [60]. Localized folding is the response for any matrix composed of nonlinear elastic, power-law viscous, plastic materials such as Mohr–Coulomb material or anisotropic materials.

Hunt *et al.* [21,25] and Whiting & Hunt [17] explore fold formation for an embedding matrix with Maxwell constitutive properties and where the elastic element exhibits softening of the form given in (2.2); the viscous element shows no softening. They show that localized folds develop due to the degradation of the elastic properties. In what follows, we present examples of the behaviours described above. We first show the development to large deflections of sinusoidal folds in a linear system to illustrate how such a system behaves with constant-velocity boundary conditions rather than with the constant-force boundary conditions of Biot. By large deflections, we mean situations where (2.5) no longer holds. This leads to system softening as the force driving the folding decreases. However, in this case involving linear constitutive behaviour, the system nonlinearity does not lead to localization of folding and the folds remain sinusoidal to large deflections. These and subsequent models have been produced using the commercially available finite-difference computer code FLAC [75].

Secondly, we show an example with constant-velocity boundary conditions where constitutive softening ultimately leads to fold localization. The system bifurcates below peak load to form a localized fold system. After peak load, some of these folds amplify rapidly; others collapse.

As indicated by previous work [59], imperfections reduce the peak load and enhance the possibility that folds will form at the site of the imperfection. For linear systems, this enhancement exists but tends to be smeared out. Nonlinear systems are highly sensitive to initial imperfections [25].

### (a) Response of linear system to softening arising from boundary conditions

Before we explore material softening behaviour, we show in figure 8 the folding response for Maxwell materials with the linear Biot reaction forces. The conditions and material parameters are given in the legend. The boundary conditions are constant velocity, so that ultimately the force at the end of the layer decreases with increased shortening. The response of the system remains sinusoidal despite this nonlinearity. All folds grow at the same rate at the same time, with visually perceptible buckles beginning to grow at peak load, although the folds nucleate while the load is increasing. At peak load, the growth is exponential but soon drops below exponential as the load drops in a manner described by Mühlhaus *et al.* [15]. Figure 9 shows a sinusoidal response for a three-layer system with the same material and physical set-up as figure 8.

### (b) Response of system to material softening

In figure 10, the same materials as in figures 8 and 9 are deformed with the same deformation conditions except that various forms of weakening are introduced. In figure 10*a*, both elastic moduli and the viscosity are softened according to
4.1where is the property involved and *n*=3. The resultant viscosity as a logarithmic plot is shown in figure 10*b*. Strong localization results and is to be contrasted with the strictly sinusoidal response in figure 8. Figure 10*a* is typical of the response with other forms of weakening investigated. These include: (i) weakening only of the elastic moduli according to (4.1), with no weakening of viscosity; (ii) weakening only of the viscosity according to (4.1); (iii) weakening of the elastic moduli, with weakening of the viscosity according to (4.1) but with substituted for *w*_{2}; and (iv) weakening of the viscosity only, with substituted for *w*_{2}. The conclusion is that weakening of both elastic moduli and viscosity as in (4.1) in single-layer Maxwell materials results in localized folding.

Figure 10*c* shows the development of localized folds in the same three-layer system as in figure 9 except that weakening of elastic moduli and viscosity is included. The outer two layers behave in a very similar fashion, while the central layer is close to sinusoidal but slightly modulated. The plot of viscosity close to the central layer (figure 10*d*) shows strong localization.

The detailed development of localized folds in a single layer is shown in figure 11. The materials are homogeneous and only one geometrical perturbation, as a spike in the form , is introduced. Here, we use *a*=1, *b*=−100, *c*=2600, so that the perturbation is centred at the first synform from the left (*x*_{1}=50 m) that grows with continued deformation. The initial amplitude of the deflection is 0.01 m. Folding deflections begin as the load increases and one well-defined localized buckle (not at the imposed perturbation) is developed at peak load. At the site of the perturbation, a less-defined localized buckle also develops. It ultimately begins to grow exponentially but growth slows under the influence of the constant-velocity boundary conditions. Just after peak load, deflections occur all along the layer in a periodic manner. By the time the load has dropped significantly, other localized buckles have developed and most of these continue to amplify. The force–deflection history of four of these buckles is shown in figure 12. Some folds continue to amplify as the load decreases, others begin to grow and then collapse to begin growing in the opposite sense as the load slowly decreases. This overall behaviour contrasts strongly with the sinusoidal fold development in figure 8.

## 5. The development of axial plane structures

We now turn to the influence of a special kind of anisotropy that is common in deformed metamorphic rocks, namely metamorphic layering. This layering grows in association with the development of the associated fold system, and in view of the preceding discussion such anisotropy is expected to influence fold growth. We present below a brief discussion of the development of metamorphic layering in order to understand the controls on the orientation of such layering relative to the fold geometry and the characteristics of the anisotropy that develops.

Deforming rocks at elevated temperatures and pressures undergo chemical reactions to produce new mineral phases stable under the ambient conditions. These reactions are commonly coupled so that they do not proceed independently of each other, but the products of one reaction act as the reactants for another reaction [76–78]. Such chemical reaction systems are called *networked* or *cyclic systems* [77,79,80]. Commonly, phases such as quartz or mica act as catalysts in that they appear both as reactants and products in the coupled mineral reaction system as a whole. In such systems, the opportunity exists for both spatial and temporal instabilities to develop [69,79,81,82]. The simplest of such systems is described by coupled reaction–diffusion equations of the form:
5.1These equations are written for a networked reacting system with one spatial dimension, *x*, with concentrations, written as *A* and *B*, of two chemical components A and B, and respective diffusion coefficients *D*_{A} and *D*_{B}; *F* and *G* are functions, normally nonlinear, of *A* and *B*, that describe the kinetics of the coupled reactions involved. Under suitable conditions, these systems of equations can become unstable [69,79] to produce temporal composition oscillations, compositional travelling waves [83] and spatial compositional patterns [84]. These spatial compositional patterns, particularly those that develop as layers, interest us in this paper.

As an example, the simplest networked mineral reactions of interest here are those involving redox reactions such as: *Fe*^{2+}→*Fe*^{3+}+*electron*. For a system where phase A is converted to phase E, *A*→*E*, with intermediate products B and D that involve redox reactions, we can write:
5.2where the *k*′_{i} are the rate constants for the various reactions. The set of equations that describes the evolution of [Fe^{2+}] and [Fe^{3+}] is:
5.3Square brackets have been used for clarity to denote concentrations.

The system (5.3) is the classical Brusselator reaction that has been widely studied [69,79,85]. The condition that the system passes through a Hopf bifurcation from a state that oscillates in time to one that forms spatial Turing patterns [84] is *b*>*a*^{2}+1, where *a*=[*Fe*^{2+}]_{steady state} and *b*=[*Fe*^{2+}]_{steady state}[Fe^{3+}]_{steady state}. Fisher & Lasaga [85] treat this system in depth and the reader is referred to their paper for details. In particular, they derive the conditions under which spatial compositional differentiation will arise. They also derive the dispersion relation for the rate of growth of such layering in one spatial dimension, *x*, from an initial compositional distribution that consists of small perturbations, *g*(*x*), above a background steady-state composition, *c*_{0}, defined by a Fourier series:
5.4In (5.4), *L* is some length scale associated with the system and *a*_{n} is the initial compositional fluctuation associated with the *n*th wavenumber. The spatial evolution of the composition, *c*, of the system is given by
5.5where *ω*_{n} is the amplification of the *n*th mode and is discussed by Fisher & Lasaga [85]. Most values of *ω*_{n} for the spatial evolution of the reactions (5.2) are negative except for *n*=2 and 3. In the example investigated by Fisher & Lasaga [85], the layering grows with an amplitude that varies approximately as .

The orientation of the compositional layers that form in the Brusselator reaction is controlled by the boundary conditions and the geometry of the reacting system. Of more interest in the folding context are reactions (of which the Brusselator can be one) that undergo instability to produce spatial patterning due to coupling with deformation.

Any process additional to chemical reaction that is coupled to changes in concentration of A and B in (5.1) can, in principle, lead to instability in the reacting system [83,86]. Diffusion is one such process but alternative or additional processes are advection of concentration through the system in a moving fluid or by deformation [86]. The reason for this is that these processes are capable of moving the stationary states defined by *F*=0 and *G*=0 in (5.1) relative to each other in compositional *A*–*B* phase space so that what is a stable state in the absence of these influences becomes an unstable state. Examples are presented by Ortoleva & Ross [83], Ortoleva [86], Ortoleva *et al.* [87], Rusinov & Zhukov [88], Hobbs & Ord [81] and Hobbs *et al.* [82]. The effect is expressed by writing (5.1) in the form:
5.6where the functions *Q* and *S* represent the influence of advection and/or deformation on the concentrations of A and B and *ε*_{A}, *ε*_{B} represent the relative importance of these two processes. If *ε*_{A}≠*ε*_{B}, then instability can arise even for *D*_{A}=*D*_{B} and for stationary states defined by *F*=0 and *G*=0 that are stable. In particular, Dewars & Ortoleva [89] consider the simple system:
5.7

This system in the absence of deformation is stable [81] but Dewars & Ortoleva [89] show that layering of quartz-rich and muscovite-rich regions can emerge from an initial homogeneous mixture of quartz and muscovite if the reaction (5.7) is coupled to deformation. Moreover, Ortoleva *et al.* [87] show that the maximum growth rate for layering exists for a wavevector oriented parallel to the compressive principal stress direction. This means that, if the layering forms early in the folding history, it grows approximately parallel to the fold axial plane. Examples of the compositional layering that develop by these coupled processes are given in figure 13 together with one possible mechanical model arising from the layering.

As noted above, the metamorphic layering that develops by these processes is commonly parallel to the axial planes of folds or is at a high angle to the folded layers (figures 1*d*, 2 and 13*b*), so that the folded layer is supported during folding by a periodic or a quasi-periodic distribution of mechanical properties along its length. Thus, if one were to consider the classical Winkler support [90] shown in figure 13*c*, then the metamorphic layering would modify the support, as shown in figure 13*d*.

## 6. Localized folds with metamorphic layering

Figure 14 shows a single-layer system in which metamorphic layering has been imposed normal to the initial orientation of the buckling layer according to
6.1where is the differentiated value of the elastic moduli or the viscosity, is the value of this quantity in figure 8, and *m*=1 and 3 in figure 14*a*,*b*. In addition, the constitutive parameters weaken as in (4.1) with *n*=3. The resulting localized fold pattern is not sensitive to whether the viscosity weakens according to *w*_{2} or . The resulting fold pattern shows a weak sensitivity to values of *m*, but the introduction of metamorphic layering does produce strong localization when compared with figure 8. The metamorphic layering is distorted on the limbs of folds so that it remains approximately normal to the folded layer. This effect is seen in natural fold systems associated with metamorphic layering (figure 2*b*).

In figure 15, a three-layer system is shown with the same weakening conditions as in figures 8 and 9; *m*=1/*π* and 1 in figure 15*a*,*b*, respectively, and *n*=1. It is notable that the central layer undergoes strong localization, in contrast to the situation without layering depicted in figure 10*c*, where the central layer, although localized, is close to sinusoidal in form. The individual folds developed in the central layer of figure 9 are more or less duplicated although distorted. The outer layers are strongly localized independently of each other, contrasting again with the situation with no metamorphic layering (figure 10*c*). The number of folds in these outer layers increases as the wavelength of the metamorphic layering decreases. We have not undertaken a systematic study of the influence of *m* and *n*, but clearly there is an effect that needs further investigation.

## 7. Discussion

In this paper, we have investigated the influence of softening behaviour of the embedding matrix on the development of folds. In keeping with a wide range of investigations, including Tvergaard & Needleman [16], Hunt *et al.* [24] and Hunt [23], the matrix softening leads to some form of localization, which is manifest as localized folding where packets of folds of different amplitudes develop rather than a sinusoidal wave train. These arise simply because sinusoidal solutions such as (2.6) are not solutions to (2.5) and are replaced by localized solutions such as (2.8). As indicated by Tvergaard & Needleman [16], these localized folds can be nucleated as sinusoidal fold systems; at this early stage, the fold development is described by the Biot theory of folding. As the folds grow and weak zones develop as the embedding matrix softens, localized fold growth emerges, illustrating a strong sensitivity to initial geometrical imperfections. There are as yet no general rules for the nucleation and evolution of an initial localized system [58]. Tvergaard & Needleman [16] point out that an initial sinusoidal solution can localize only to become sinusoidal with a larger wavelength at higher strains. This transition from localized to sinusoidal has also been described by Knobloch [58]. Another possibility is that packages of localized folding interfere to produce chaotic fold systems [21,25,28]. The picture that emerges from these studies is that the folding systems adopt configurations that minimize the local energy in the system (cf. [59]). For linear systems, there is one energy minimum and that corresponds to sinusoidal solutions to (3.2). For nonlinear systems, there are a large (even infinite) number of local minima and the system jumps from one minimum to another during deformation [60]. In some instances [28], these minima may have fractal geometries.

We have been particularly interested in examining the influence of periodic distributions of softening response along the length of the folded layer imposed by the growth of metamorphic layering. The analyses of Ortoleva *et al.* [87] and Hobbs & Ord [81] show that those forms of metamorphic layering that form through the coupling of mineral reactions with deformation develop normal to a principal axis of compression or, equivalently, for coaxial constitutive relations, parallel to a principal plane of deformation rate. Here, we have assumed that the metamorphic layering forms relatively early in the folding history when the principal axis of compression is everywhere approximately parallel to the layer undergoing compression. This means that the metamorphic layering grows approximately parallel to the axial plane for incipient folds. One can define a dimensionless number, , which is the ratio of the growth rate for metamorphic layering, *ω*_{metlayer}, to the growth rate for folding, *ω*_{folding}:
*ω*_{metlayer} is considered for the Brusselator reaction by Fisher & Lasaga [85], eqn. 46 and an example of *ω*_{folding} is given in (3.8). If , the layering forms early and is likely to be distorted locally by the folding, as shown in figures 14 and 15. If , then the layering grows late in the folding history and is controlled in orientation by the local stress field; this means that it must develop a fan-like distribution in the hinges of folds. Depending on the deformation history, the axis of this fan may or may not be parallel to the fold axis.

Although the discussions of Fisher & Lasaga [85] and of Ortoleva *et al.* [87] arrive at a sinusoidal distribution of metamorphic layering, it should be borne in mind that many reaction–diffusion equations can be expressed as Swift–Hohenberg equations [56,60,69], so that non-sinusoidal solutions to the spatial distribution of metamorphic layering can arise in the same way that non-sinusoidal solutions can arise from the identical class of equations that describe initial fold evolution.

In fact, the relationship between buckling instabilities and the formation of metamorphic layering may be closer than one first expects. Elsewhere [82], we have shown that a combination of the first and second laws of thermodynamics says that the evolution of a quantity *w* within a deforming reacting body is given in a one-dimensional spatial coordinate, *x*, by an equation of the form
7.1where *F* is commonly a nonlinear function of *w* and *α* is a parameter with the dimensions of diffusivity. In many instances, the evolution of a system requires ℵ coupled equations of the form of (7.1) involving different quantities *w*_{1},*w*_{2},…,*w*_{ℵ}. Equation (7.1) is a standard reaction–diffusion equation that has been widely studied over the past 50 years [69]. A special form of Equation (7.1) is the Swift–Hohenberg equation [69], p. 873:
7.2Both (7.1) and (7.2) appear in various forms as the basis for describing many of the structures that we observe in deforming metamorphic rocks. In particular, a stationary form of (7.2), given by ∂*w*/∂*t*=0, is
7.3Solutions to (7.3) take on many forms depending on the value of *α* and the nature of *F*(*w*). These solutions can be sinusoidal, modulated, localized, structured kinks and pulses, and chaotic [56]. Some special forms of (7.1) and (7.3) are as follows:

— If

*F*(*w*) is linear then (7.3) reduces to the standard biharmonic equation discussed widely in structural geology [13,66,67,91] for the initial buckling of a strong layer embedded in a weaker material. In this case,*F*(*w*) expresses the reaction forces exerted by the embedding material on the layer as the layer buckles. The solutions to the biharmonic equation for a linear*F*(*w*) are strictly sinusoidal and play a prominent role in the classical Biot theory of folding, where*w*is the displacement of the layer normal to its initial orientation.— If

*F*(*w*) is a nonlinear softening function then the application of (7.3) to the low-amplitude buckling of layered materials results in fold systems more complicated than those arising from the Biot theory. A sinusoidal wave form is no longer a solution to (7.3) and the buckles localize sequentially along the deforming layer, forming packets of folds. A wavelength selection process no longer operates. Depending on the values of*α*and the nature of*F*[56], the buckles are ‘multi-bumped’ but periodic or localized and can be chaotic [17,21,25,27]. Examples of all of these forms of folds are illustrated by Hudleston & Treagus [12]. An important remaining problem is to establish the conditions under which the structures that develop in the embedding material are also localized to form micro-lithons and crenulation cleavages.— If (7.1) is thought of as applying to ℵ coupled chemical reactions then chemically unstable behaviour can lead to metamorphic differentiation and/or compositional zoning in growing mineral grains [81].

— Metamorphic differentiation resulting from (7.1) leads to a sinusoidal reaction response of the embedding material on a buckling layer, which in turn, through (7.3), produces localized folding with the metamorphic layering axial plane to the folds. Equations (7.1) and (7.3) have profound implications for structural-metamorphic geology. Understanding the rich assemblage of deformed metamorphic fabrics observed in Nature depends on developing physically realistic and geologically relevant descriptions of

*α*and of*F*and the ways in which folding and the simultaneous development of metamorphic layering may be expressed by the same set of coupled equations.

The behaviour reported in this paper suggests a simple mechanical model for the folded system with metamorphic layering (figure 16). Although individual layers have their own intrinsic buckling modes, they are connected to each other by a system of Maxwell units (figure 16*a*) that weaken in response to the folding deflections. The behaviour of one part of the system is coupled to all other parts of the system. The introduction of metamorphic layering (figure 16*b*) imposes another spatial periodic or quasi-periodic response. The model suggests that some insight into the behaviour of folded multi-layered systems with metamorphic layering might benefit from work on the buckling of nonlinear struts.

## 8. Concluding remarks

We have used the term *Biot folding* or *Biot folding theory* to mean the selective amplification of initial perturbations in a deforming layer embedded in a matrix of weaker material. For folding theory, in general, two different cases can be distinguished based on whether the response of the matrix to deflections of the folding layer is linear or nonlinear and the characteristics of the two are summarized below.

*Linear matrix*.

— Layer(s) may be linear or nonlinear.

— Biot folding theory holds together with many subsequent modifications involving linear matrix materials.

— Analytical solutions exist to finite deformations for a range of plate models and boundary conditions.

— An amplification spectrum exists based on selection of fast-growing Fourier components. All folds form and grow synchronously.

— All solutions are sinusoidal unless imperfections are present. The energy landscape, in the absence of initial geometrical perturbations, contains just one smooth minimum and solutions rapidly evolve to this attractor.

— Amplification rate, and its evolution, and

*λ*/*h*depend on assumptions concerning the constitutive relation for the plate, boundary conditions and strain distribution in the layer.— The only origin for non-sinusoidal solutions is imperfections (either initial or ones that grow during folding); imperfections that have their own growth spectrum wavelengths close to the Biot wavelength grow fastest (at least for Newtonian materials).

*Nonlinear matrix*.

— Layer(s) may be linear or nonlinear.

— Biot folding theory does not hold.

— Analytical solutions may sometimes be obtained for special cases, but these are not necessarily unique and commonly an infinite number of solutions exist. If one looks for sinusoidal solutions, one might discover some especially using a linear stability analysis. However, one needs to show that such a solution is unique and stable during fold growth. In general, solutions are found by exploring the phase space of the system and searching for bifurcations using programs such as AUTO [61] and XPPAUT [62].

— No growth selection process exists based on the fastest-growing Fourier modes. Folds nucleate along a layer at different times and grow (or collapse) at different times. Sequential development (or collapse) of folds is common.

— Localized solutions are the norm. The energy landscape commonly contains an infinite number of energy minima, some of which may be fractal in structure. The result is that folds evolve by the system jumping from one local minimum to another depending on the local conditions of deformation and geometry. Solutions can consist of sinusoidal waves, modulated waves, periodic and non-periodic localized structures, pulses, kinks and chaotic wave forms. The system can jump from any of these forms to another nearby form.

— There need not be any well-defined value for

*λ*/*h*but clearly for periodic solutions this is possible. The value of*λ*/*h*at any time at a particular part of the layer is controlled by the time evolution of parameters in the evolution equations. At present, there is no general understanding of how these are controlled by constitutive behaviour or geometry.— Non-sinusoidal solutions are common and the precise wave form is strongly dependent on initial imperfections or on imperfections, such as metamorphic layering, that grow during deformation.

Although initial geometrical perturbations of large wavelength undoubtedly have an influence in promoting the formation of some localized folds, an additional mechanism involves the localized response to weakening in the matrix between the buckling layers. Nonlinearities arising from large deflections and thick-beam theory using the linear non-softening matrix response of the Biot theory have not so far emerged as a mechanism for localization.

The structures developed in deformed metamorphic rocks are commonly localized but quasi-periodic over large distances compared with the layer thickness, and this suggests that a control arising from nonlinear matrix response is important in many instances. As such, these structures add yet another layer of richness to the complexity of the folding process observed in nonlinear layered materials. The Biot theory may appear as the analysis that describes the initial sinusoidal deflections in a layered system before softening of the embedding matrix appears. After the embryonic sinusoidal fold system forms, localization of the fold system can occur and chaotic systems develop if localized packets of folds interfere. The growth of metamorphic layering parallel or approximately parallel to the axial planes of folds introduces a periodicity in matrix response along the length of the folding layer, thus introducing the opportunity for even greater fold localization.

Some workers may regard the identification of localized folding behaviour based on a nonlinear matrix response as a retrograde step because it negates the proposal, derived from the Biot theory, that fold geometries offer information on viscosity ratios and fold development. In fact, the opposite is true. Identifying the mechanisms of localization in geological materials offers a vast new field of study that concentrates on the processes that operate during folding and that lead to nonlinear behaviour. Natural folds display a bewildering array of geometries in fold profile. Strictly sinusoidal wave forms are rare, and periodic localized, non-periodic localized and chaotic systems are the norm. The challenge for the future is to document field examples where the progressive development of localized systems can be observed and to identify the processes that have operated during folding. Attempts to fit localized fold systems into theories based on the Biot approach are not only misleading, but also neglect the rich array of opportunities these nonlinear systems offer to understand the processes that operate during deformation and metamorphism.

## Acknowledgements

We are indebted to Giles Hunt, Chris Budd, Tim Dodwell, Jill Hammond and Hans Mühlhaus for many discussions on fold localization. We thank Stefan Schmalholz and an anonymous reviewer who supplied critical and helpful reviews that greatly improved the paper.

## Footnotes

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

- This journal is © 2012 The Royal Society