## Abstract

Previous multi-layered folding models have struggled to describe the geometry of the interfaces between layers. In the level set method (LSM), a single function *ϕ*(*x*, *y*, *t*) can be used to encode all the information of a multi-layer material in which all layers are in contact and thus is a very natural way of dealing with the geometry. This paper shows the potential for the LSM (in multi-layer problems) by demonstrating that it can describe the geometry of multi-layer folding patterns including those with singularities. The method is then applied to describe the mechanics of the modelling of parallel folding in multi-layered structures.

## 1. Introduction

Deformation patterns of multi-layered materials under compression are strongly influenced by the way that the layers interact and slide over each other, and hence are quite different from the folding patterns of homogeneous materials. Such folding arises naturally when layered materials are compressed in their plane, in a medium that allows them to slide over each other but not to separate. An important example occurs in the buckling of layered rocks under tectonic compression. The same mechanisms also appear in the compression of layers of paper (for which we provide some experimental results) and in certain types of composite materials. Structural geologists have classified folds according to their profile and specific names have been given to those that frequently appear (Price & Cosgrove 1990). Their geometries range from almost sinusoidal patterns to the straight limbs and sharp corners of kink banding and chevron folding. However, in all cases, the resulting deformation is strongly influenced by the constraints of the multi-layered geometry; a subtle mix of *geometrical* constraints imposed by the need for the layers to fit together, and *mechanical* constraints of bending stiffness and interlayer friction.

Here, we are mostly interested in *parallel folding* (folds where the orthogonal thickness remains constant). This occurs when a finite number of layers, loaded in their plane, deform into a softer surrounding medium (foundation or matrix) while slipping at their interfaces (Edmunds *et al*. 2006). Parallel folds formed under large overburden pressures would be expected to limit voids between the layers (carrying a large energy penalization), and we therefore take the process to be one of buckling in the complete absence of voids (Budd *et al*. 2003; Edmunds *et al*. 2006).

Figure 1*a* shows an example of parallel folding in rocks with a layer thickness of around 10–25 cm. Since the layers fit together without voids, each has a slightly different geometry from its neighbour and as a consequence, a singularity can appear on one interface, as seen in the figure. Here, the layers appear smooth (almost sinusoidal) at the bottom of the picture, but as we move towards the top they become increasingly constrained by the geometry, until at a particular interface a singularity occurs. Past this singularity, the layers take a non-differentiable ‘V’ shape which propagates without change as we move further up; this V shape is also seen in kink bands and in chevron folds (figure 1*b*).

A series of experiments showing regular parallel folding in layers of paper has been undertaken by two of the current authors and figure 2 shows a typical outcome. A method for modelling the fold geometry, developed by Budd *et al*. (2003) for two flexible layers and extended to *n* layers by Edmunds *et al*. (2006), considered a sinusoidal Galerkin approximation for the shape of the central layer and simply extended it to the outer layers. This technique is unable to deal with the geometrical constraints of multi-layer folding leading to singularities. This means that the associated problems of kink banding and chevron folding (layers bend with infinite curvature and straight limbs such that no voids form; figure 1*b*) require a separate treatment for what is essentially a similar problem (e.g. Hunt *et al*. 2000; Wadee *et al*. 2004).

Here, we present a method for describing multi-layer parallel folding based on the level set approach. This naturally copes with the *geometry* of smooth parallel folding and also allows a consistent description of the geometry associated with singularity formation and kink banding (§2). The level set approach propagates the shape of the central layer at a prescribed rate in the normal direction. By doing this, the shapes of all outer layers can be determined (and controlled) rather than being just approximated. The normal propagation is achieved using the level set method (LSM; Dervieux & Thomasset (1979, 1981); Osher & Sethian 1988) in which a function *ϕ*(*x*, *y*, *t*) is introduced such that the shape of the layer indexed by *t* is given by the (*x*, *y*) curve satisfying *ϕ*(*x*, *y*, *t*)=0 (the zero level set). An application of this method to the rock layers pictured in figure 1*a* is shown in figure 3. A spline function has been fitted to a central layer (dashed line) and the resulting curve propagated by calculating *ϕ*. We see that all of the geometrical features of the layers are resolved, including the cusp and singularity formation as we propagate in the upwards direction and the simple geometry in the downwards direction. To apply this method to find the actual form of a deformed material requires information only about the shape of a single reference layer. The LSM then gives the position of all the other layers and encodes all the information of the geometry in terms of *ϕ*. Additional information about the *mechanical* behaviour of the system, in particular the bending energy, the work done by friction and the work done in compressing the embedding medium is then determined for the individual layers. This provides a total potential energy function for the whole system of the deformed layers, which when coupled to the geometrical description of the layers given by the LSM contains a significant (stiffening) nonlinearity. We then look for states which are stationary points of this energy functional. As the energy can be expressed in terms of the single function *ϕ* this calculation becomes relatively straightforward.

The paper first describes the geometry of multi-layer folding, showing that it naturally leads to the formation of singularities. We then describe the LSM for calculating this geometry, and discuss its performance and errors when singularities are present. Following this, the LSM is used to construct a potential energy function for the deformed material, up to the point of a singularity which we take as the limit of the folding pattern. (Beyond this, elasticity assumptions would be expected to break down.) From this energy functional, we deduce the profile of the deformed multi-layer material.

## 2. The geometry of multi-layer folding

We now consider the basic processes of multi-layer folding, first using a Lagrangian view to describe the geometry of the expected patterns. A natural consequence is the observation that singular behaviour arises naturally in such systems.

### (a) The basic folding process

Sedimentary rocks can be formed under the sea bed as loose particles (sediments), laid down layer upon layer, and forced together by the overburden pressure caused by the material above. This process forms horizontal layers, which fold when tectonic plate movement produces an axial load (figure 4). In this figure, we see identical layers of length *L* and width Δ*t*, under an overburden pressure *q*, undergoing an end-shortening due to the compressive load.

Consider a layered material to be characterized by the position of each layer. We can then consider curves *Γ*_{t}, indexed by *t*, which describe the interface between layers. To describe such a layered material, let the upper surface of top layer be described by the curve *Γ*_{nΔt} and the lower surface of the bottom layer by the curve *Γ*_{−nΔt} with the interfaces between the other layers described by the curves *Γ*_{iΔt}, −*n*<*i*<*n*. In *parallel folding,* we assume the layers are always in contact and the curves *Γ*_{iΔt} separated by a *constant* normal distance Δ*t*. The geometry of the set of curves can then be described by taking a *reference curve* (without loss of generality *Γ*_{0}) and propagating this in a normal direction at a distance *i*Δ*t* to give the curve *Γ*_{iΔt}. Generally, the shape of *Γ*_{iΔt} will differ from *Γ*_{0} (although the chevron folding phenomenon of figure 1*b* gives an example where an identical family of curves can coexist). In particular, an initially smooth curve *Γ*_{0} may give rise to curves *Γ*_{iΔt} with singularities and sharp corners as seen in figure 1. To propagate past such a singularity, great care needs to be taken.

### (b) The Lagrangian representation of parallel folding

We now give a Lagrangian calculation of the process of parallel folding and show that the length of the interface is conserved, provided: (i) the surface remains smooth (i.e. no singularities develop) and (ii) the angle of the interface with the horizontal is the same at either end. We consider the set of curves given parametrically by *Γ*_{t} : → so that *Γ*_{t}={(*x*(*s*, *t*), *y*(*s*, *t*)) : *s*∈[0, 1]}. The parallel folding assumption is that the *normal separation* between two curves parameterized by *t*=*t*_{1} and *t*=*t*_{2} is given by |*t*_{2}−*t*_{1}| and does not depend on *s*. This leads to a simple construction of the *entire* set of curves for all *t*, given a reference curve for *t*=0. The exact position of *Γ*_{i} at time *t* may be constructed by advancing each point of *Γ*_{0} in its unit normal direction, ** n**, at a distance

*t*=

*i*Δ

*t*(Sethian 1999) where(2.1)Here, the superscript 0 refers to the reference curve. It follows from standard analysis that for each fixed

*t*=

*i*Δ

*t*, the Lagrangian description of

*Γ*

_{t}takes the form(2.2)Equation (2.2) is valid for all

*t*=

*i*Δ

*t*, although the resulting curve may have points of singularity characterized by a lack of a well-defined normal vector. To find

*x*

_{s}(

*s*,

*t*) and

*y*

_{s}(

*s*,

*t*), we differentiate equation (2.2) with respect to

*s*and simplify by considering the curvature

*κ*(2.3)We make the significant observation that this vector vanishes (and hence we cannot define a normal vector) when

*t*=1/

*κ*(

*s*, 0) which is the radius of curvature of the reference curve at this point. By moving

*t*forward and then backwards by the same amount, it is immediately clear that (1−

*κ*(

*s*, 0)

*t*)(1+

*κ*(

*s*,

*t*)

*t*)=1, so that(2.4)

Note that the curvature *κ*(*s*, *t*) becomes infinite at the first value of *t*, which occurs on the normal through the point on the reference curve with the *smallest radius of curvature*. If *κ*(*s*, 0) is *negative* then the curve *Γ*_{t} has a well-defined normal for all positive *t*. However, if *κ*(*s*, 0) is *positive*, then there will be a first value of *t*==min(1/*κ*(*s*, 0), at which the normal vector first vanishes. The value of gives the maximum width of the layers before a singularity occurs. For *t*>, the curve *Γ*_{t} is multi-valued and has a swallowtail singularity with points of infinite curvature (Sethian 1999). This situation is illustrated later in figure 5 where we take a cosine to be the reference curve.

As the layers are assumed to be incompressible in their local tangent direction, it follows that the *arc-length* of each layer must remain constant. This imposes a constraint on the allowable deformations. To consider the effects of this constraint, we calculate the arc-length. Suppose that d*σ*_{t} is the infinitesimal arc-length on the curve *Γ*_{t}. If the total arc-length of *Γ*_{t} is *σ*_{t}, thenNowTwo cases then arise: either *t*< in which case 1−*κ*_{0,s}*t* has a constant sign over the length of *Γ*_{t} or it changes sign at certain points along the curve. If *θ*(*s*, *t*) is the angle of the curve *Γ*_{t} at the point (*x*(*s*, *t*), *y*(*s*, *t*)), then tan(*θ*)=d*y*/d*x* and *κ*=∂*θ*/∂*σ*. Thus, if *t*< we haveWe arrive at the result:

*While equation* (2.2) *holds and t*< *total arc-length of Γ*_{t} *is given by*:(2.5)

*If t*< *and* , *the total arc-length of each curve Γ*_{t} *remains constant*.

If *t*>, the normal ceases to exist and the direction of d*σ*/d*s*(*t*) changes sign at certain points so that |1−*tκ*(*s*, 0)| changes from 1−*tκ*(*s*, 0) to −1+*tκ*(*s*, 0). Our reasoning breaks down, and there is no guarantee that arc-length is preserved.

### (c) Singularities

If *t*>, the curve *Γ*_{t} has a singularity resulting in a self-intersecting curve with infinite curvature. We now explore what happens as *t*→^{−}. To study this situation, we consider what happens when a parabola is propagated forward at a constant speed. We take a reference curve *Γ*_{0} defined for *s*∈[0, 1] byIf we set *S*=(*s*−1/2) then the unit normal ** n** is given byand hence,The radius of curvature of the reference curve

*Γ*

_{0}takes its smallest value of 1/

*K*when

*S*=0 and we deduce that the singularity occurs when

*t*==1/

*K*. If

*t*>1/

*K*then the curve self-intersects when (

*x*(

*S*,

*t*),

*y*(

*S*,

*t*))=(

*x*(−

*S*,

*t*),

*y*(−

*S*,

*t*))=(0,

*Y*), which occurs whenAt the point of self-intersection, the curve

*Γ*

_{t}has a locally V-shaped form. The gradient of the curve is given by d

*y*/d

*x*=

*y*

_{s}/

*x*

_{s}. Substituting

*K*

^{2}

*S*

^{2}=

*K*

^{2}

*t*

^{2}−1 in the expressions for

*y*

_{s}and

*x*

_{s}we have that at the point of self-intersection(2.6)

To examine the form of the curves as we approach the singular point, we set *t*=1/*K*−*δ* and let *S* and *δ*>0 be small. To leading order we then haveIt follows thatIt follows immediately that when *x*=0 we have d^{2}*y*/d*x*^{2}=∞ if *δ*=0 and d^{2}*y*/d*x*^{2}=1/*δ* if 0<*δ*≪1. This is significant in the context of the mechanics of buckling. The bending energy associated with a buckled layer is proportional to the integral of the square of this derivative which approaches infinity as *δ*→0. The buckling energy is thus a strongly nonlinear function of the displacement of the reference layer, which has a significant effect on the resulting buckling profiles. It is difficult to know exactly what happens to the rock layers physically as a singularity is approached and possibly transgressed. However, the assumed geometry forces this to happen and a comparison of figures 1 and 5 shows clear similarities. It is clear, however, that the rocks do not pass through themselves. Indeed, it seems that the rock layers of figure 3 (shown as a dashed line) take the geometry of the *non-self-intersecting* part of the curve. Hence, to fully realize all possible geometries, we are motivated to look at a construction method which is different from the Lagrangian approach and allows us naturally to resolve and pass through the singularity, enabling a natural description of the non self-intersecting part of the curves when *t*>1. Such a procedure is given by the LSM, which relies on an Eulerian representation of the set of layered surfaces and a weak formulation of the equation of propagation of the layers. Sethian (1985) uses an *entropy condition* to overcome the problem of self-intersection. Informally, this strategy deletes the swallowtail, a process that causes a loss in length. The process is equivalent to using a numerical method to solve the level set equation described in §3 to give a regularized (viscosity) solution.

## 3. The Eulerian representation, the theory of level sets and the viscosity solution

### (a) The Eulerian representation and the level set method

In contrast to the previous description, the Eulerian representation aims to find the (*x*, *y*) equation of the curves *Γ*_{t} thinking of *t* as a continuous variable and looking at the differential equations governing the propagation of *Γ*(*t*). Following Sethian (1999), let ** x**(

*t*)=(

*x*(

*t*),

*y*(

*t*)) be a point on the curve

*Γ*

_{t}≡{(

*x*,

*y*) :

*y*=

*w*(

*x*,

*t*)} and

**(**

*v**x*,

*y*) be the speed in the normal direction. Starting from the identity

*ϕ*(

*x*,

*y*,

*t*)=0 it follows that

*ϕ*satisfies the following hyperbolic partial differential equation (the

*level set equation*)(3.1)Solving this equation for

*ϕ*(

*x*,

*y*,

*t*) together with a prescribed initial function

*ϕ*(

*x*,

*y*, 0) determines the curves

*Γ*

_{t}={(

*x*,

*y*)|

*ϕ*(

*x*,

*y*,

*t*)=0}.

In the level set formulation, the condition for *parallel folding* is that the speed function, ** v**(

*x*,

*y*) should be

*constant*. However, it would be straightforward to include other forms for the function

**to allow, for example, for situations in which the rocks can compress laterally (Wadee**

*v**et al*. 2004). Other examples of speed functions include the cases where

**depends on the curvature**

*v**κ*(which arises in flame front propagation (Rhee

*et al*. 1995)) and where it depends on |∇

*ϕ*| (which arises in electro-machining). With this formulation, the important properties (figure 4) of the geometry of the layers of rock can be easily calculated in terms of

*ϕ*.

### (b) Numerical implementation

To propagate the interface using the LSM, we calculate the whole function *ϕ* using a time-stepping method. Assuming that we have a good approximation to *ϕ*(*x*, *y*, *t*) for some time *t*=*i*Δ*t*, we first use a numerical method to find an approximation to the solution *ϕ* of (3.1) for the time level *t*=(*i*+1)Δ*t*. A second numerical method is then used to locate the zero contour of *ϕ*(*x*, *y*, (*i*+1)Δ*t*) to find an approximation to *Γ*_{(i+1)Δt}. The algorithm we use for implementing this strategy is as described in Osher & Sethian (1988). To implement this, we consider an approximation to the function *ϕ*(*x*, *y*, *t*) at the point (*x*, *y*, *t*)=(*jh*, *kh*, *n*Δ*t*). Here *h* and Δ*t* are small and constant. An explicit discretization of (3.1) for a constant ** v** takes the form(3.2)where(3.3)Here the terms , etc. are shorthand for the difference operators defined by(3.4)with similar expressions for . Local errors in this scheme are of (

*h*, Δ

*t*). The LSM as implemented is stable provided that the Courant–Friedrichs–Lewy condition (Courant

*et al*. 1928) is met. There are three sources of error when numerically solving the level set equation (Sethian 1999). These are the

*initialization error*which is error obtained when finding the signed distance function on a discrete mesh from an initial curve; the

*update error*which is error in finding an approximation to

*ϕ*(

*x*,

*y*, (

*n*+1)Δ

*t*) given

*ϕ*(

*x*,

*y*,

*n*Δ

*t*); and the

*measurement error*which is the error associated with finding the position of the zero level set using, for example, the Matlab contour routine. Note that in order to calculate

*Φ*

^{n+1}we must calculate

*Φ*

^{n}at

*all*spatial points, not just those close to the interface. This is an inefficiency in the method. We note that there are more efficient schemes such as the ‘narrow band method’ (Chopp 1993), in which the only entries to be updated are those in a narrow band around the interface. Once the interface reaches the edge of the band the signed distance function is recalculated, the narrow band is re-formed around the current position of the interface and the problem is solved until the interface hits the edge of the band again and so on. Hence, a balance is struck between saving time for updating only a small number of elements and the cost of recalculating the signed distance function.

### (c) Convergence of the LSM and the viscosity solutions

The level set equation (3.1) considered in the current paper is an example of a wider class of equations called *Hamilton–Jacobi* equations, in this case given by(3.5)In general, there can be no smooth solution to equation (3.5) lasting for all times *t*>0. This corresponds precisely to the swallowtail singularity observed in §2*b*. However, we may regularize the Hamilton–Jacobi equation by adding artificial viscosity, *ϵ*Δ*ϕ*^{ϵ}, with *ϵ*>0 to give(3.6)In general, provided *ϵ*>0, equation (3.6) admits smooth solutions. If the solutions *ϕ*^{ϵ} of (3.6) converge weakly to a function *ϕ* as *ϵ*→0, then this is a weak solution of (3.5). It is called the *viscosity solution* and can have gradient discontinuities (Crandall & Lions 1983; Crandall *et al*. 1984). The viscosity solution picks out the correct weak solution when no classical solution exists. Informally, it ‘deletes’ the multi-valued part of the swallowtail curve. In the present numerical implementation, the regularization is automatically provided by the numerical truncation error which introduces an artificial viscosity term. The numerical method is effectively solving a problem close to (3.6) with *ϵ*≈(*h*) (Enright & Fedkiw 2003). Consequently, the numerical method will give solutions close to the viscosity solution, automatically smoothing out the singularities and regularizing the sharp corners.

### (d) Singularities found in rock folding

We conclude this section by looking at some examples in which we explore the applications and limitations of the LSM for calculating the geometry of certain problems given a known reference layer *Γ*_{0}.

#### (i) Examples with re-entrant corners

As a first quantitative calculation, which has a close link to a rock folding problem which we will consider further in §4, we consider a reference curve given by *Γ*_{0}≡{(*x*, *y*) : *x*=*s*, *y*=cos(2*πs*)}, 0≤*s*≤1. To apply the LSM, we take an initial function *ϕ*(*x*, *y*, 0)=*y*−cos(2*πx*). This is not a *signed distance function*, meaning that gradients are steeper and therefore harder to approximate accurately numerically, but it does give an accurate and easy to implement initial zero level set. We now compute the resulting layers *Γ*_{t} and compare these with the solutions predicted by the Lagrangian formulation. Using this, the exact parametric equation of the layer at the time *t* is given byThe curvature of the reference curve takes its maximum value of 4*π*^{2} at the point *s*=1/2 and hence a singularity occurs when *t*=1/4*π*^{2}, *x*=1/2, *y*=−1+1/4*π*^{2}. A close-up of the singularity of the exact (multi-valued) solution arising from the Lagrangian description is plotted in figure 5*a*. We now compare the Lagrangian solution with that derived by using the LSM. A calculation using the method for *h*=0.01 and Δ*t*=0.005 is presented in figure 5*b* (here the corresponding close-up of the singularity is shown). Observe that, in contrast to the Lagrangian description, the LSM has deleted the self-intersecting part of the curves and the resulting curves have an apparent gradient discontinuity at the centre. The local V-shaped nature of these curves is very similar to that of the layers in the chevron folding pattern illustrated in figure 1*b*. Indeed, if we take *Γ*_{0} to be the V-shaped curve *Γ*_{0}≡{(*x*, *y*): *x*=*s*, *y*=|*s*−1/2|} then the resulting calculation of the layers *Γ*_{t} using the LSM is given in figure 6*a* and a close-up in figure 6*b*. We see that the LSM has successfully coped with the gradient singularity, reproducing the self-replicating feature of parallel folding in this case where all layers *Γ*_{t} have exactly the same shape and the same arc-length. We note that in figure 6*b*, the corner is slightly smoothed due to the error. This effect can be reduced by refining the mesh.

As a further measure of the accuracy of the calculation of the propagating cosine reference curve *Γ*_{0}≡{(*x*, *y*): *x*=*s*, *y*=cos(2*πs*)}, 0≤*s*≤1, we plot the total arc-length of the resulting curves. The choice of reference curve ensures that [*θ*_{t}]=0. Hence, from theorem 2.1, the total arc-length of the curve *Γ*_{t} stays constant up to the point of singularity formation when *t*=1/4*π*^{2}. For larger values of *t*, the curve self-intersects and the total arc-length of the curve, omitting the self-intersecting part, decreases until a steady state is reached. This arc-length can be calculated exactly. If *s*=*s*^{*} is the parameter value on the reference curve for which the normal intersects with the line *x*=1/2, thenAs the normal to the reference curve is also normal to *Γ*_{t}, the angle *θ* is given by *θ*=tan^{−1}(−2*π* sin 2*πs*). The arc-length *σ*(*t*) is then given by theorem 2.1 as(3.7)This is a standard elliptic integral and can be evaluated analytically. Similarly, the arc-length of *Γ*_{t} can be found numerically by applying quadrature to the sets found by the LSM. The two values are compared in figure 7.

As a second example of an application of the LSM, we consider again the calculation illustrated in figure 3. A reference layer *Γ*_{0} was obtained directly from the photograph by sampling the photographed curve at regular points at a roughly middle layer and then fitting a spline through the resulting data. The resulting layers propagated (through the singularity) by solving the level set equation (3.1). Here, we took *h*≈5 and Δ*t*≈2.5, and chose to plot the zero level set at every 15 time-steps (there are 640 pixels in the horizontal direction). The qualitative agreement between the calculations and the observed geometry appears very good including, apparently, the representation of the singularity and of the rock layers formed beyond the singularity.

#### (ii) Examples with salient corners

The viscosity solution does not always agree with the observed patterns of folded rock. In the photograph of figure 1*b*, we see an example of *chevron folding* in which rock layers fold in a zigzag with both re-entrant and salient corners, each of which has a gradient discontinuity. Instead of propagating salient corners without change, the artificial viscosity in the LSM picks out a rarefaction fan of solutions (Sethian 1999), which is not consistent with that observed in practice. To overcome this difficulty, we could follow the procedure outlined in Russo & Smereka (2000) in which the LSM is used to describe crystal growth with sharp corners. In such a method, additional slip is introduced between the layers. This has the advantage of preserving the total arc-length of the solution, but we shall not consider this approach in the present paper.

## 4. The mechanics of folding

Having determined a procedure for finding the overall geometry of the layered system from a single reference layer *Γ*_{0}, we now examine how the form of *Γ*_{0} itself can be calculated by using properties of the function *ϕ*. This section shows how we can incorporate mechanical features into the geometrical description given by the LSM to find the shape of the reference curve. The profile of the compressed layered material is determined by the interplay of several mechanisms namely the effect of the applied force, the bending properties of the layers, the work done into the external medium, the effects of the geometry constraints and the frictional forces acting to oppose sliding between the layers. These are best described by calculating the total energy *V* which is a combination of the bending energy, frictional energy and the work done by the external forces. Crucially, this total energy *V* can be calculated in terms of the level set function *ϕ* depending on *Γ*_{0}. We can then find the profile of the material by finding stationary points of *V* with respect to (appropriate) variations in *Γ*_{0}.

### (a) Total energy in terms of *ϕ*

We consider the geometry illustrated in figure 4 in which an initially undeformed multi-layer material enclosed within a foundation material is compressed and buckled by an external load. This load is resisted by the stiffness of the layers, the frictional force as the layers slide in contact and the resistance of the matrix in which the layers are embedded. In Budd *et al*. (2003), the potential energy of two axially and transversely incompressible layers of small thickness *t* in contact along an interface line was considered. In this, the two layers each comprised a material of bending stiffness *EI*, embedded in a soft foundation of transverse stiffness *k* per unit length and compressed tangentially by a load *q* per unit length, and moved longitudinally a horizontal distance (the *end-shortening*) by a load *P*. The total energy was decomposed into four components: *V*=Bending Energy+Foundation Energy−Work Done by Load+Work Done against Friction, or(4.1)The function *Χ*=±1 is included in (4.1) as a *friction indicator* which ensures that the friction makes a positive contribution to the total potential energy when friction opposes the external force, and conversely, makes a negative contribution when friction acts in the same sense as the external force. For details on the interaction between *Χ* and *μ* see Budd *et al*. (2003). In this two-layer model, *w* is the vertical deflection of the interface between the layers, *x* the horizontal distance and *σ* the arc-length. The total length of each of the layers is a constant *L*. Expressions for each of the terms making up *V* are derived in Budd *et al*. (2003). To extend this two-layer model to a multi-layer formulation, we assume as before that there are 2*n*+1 layers each of width Δ*t* and length *L* and that *w*_{iΔt}(*x*), *i*=−*n*, …, *n* is the vertical deflection of the interface between two successive layers. In the experiments conducted, each of the layers is compressed by the load by *the same horizontal amount* . The simplest way to ensure this occurs with the level set formulation is to set . The total energy of the multi-layer material is then(4.2)(4.3)(4.4)(4.5)where *b* is the breadth of the layers and dots denote differentiation with respect to *σ*. In this expression, the external foundation is assumed to have a Winkler force law with the deflection of the top and bottom layers making a contribution to the energy. We also assume that friction acts only between the layers and not between the foundation and the boundary layers. It is important to observe that the bending energy has a component which varies as the integral of . As we have already seen, this term becomes unbounded as we approach a singularity, and the expression (4.2) is not defined at this point. For the present, we will assume that no singularity has formed so that this issue does not arise and will return to consider the singular case. We can now combine the geometric description of the layers with the mechanical description in the previous subsection by expressing *V* in terms of the single function *ϕ*. Before we reformulate the energy contribution in terms of the level set function *ϕ*, it is useful to state a lemma that expresses various terms which arise in the energy expression in terms of *ϕ*.

*The derivatives of x and y with respect to σ are as follows*:(4.6)*The curvature κ can be expressed as*(4.7)*The angle of the tangent θ is given by*

(4.8)

*Proof*. See Sethian (1999).

We now use these expressions to reformulate the various components of the energy contributions (4.2)–(4.5) in terms of *ϕ*(*x*, *y*, *t*) where, to simplify the expression we set *ϕ*^{i}=*ϕ*(*x*, *y*, *i*Δ*t*).

*In the case of a* non-singular *deformation the energy terms are given by*(4.9)(4.10)(4.11)(4.12)The delta function is included in equations (4.9)–(4.12) so that it is the position of the *zero level set*, *Γ*_{iΔt}, that makes the contributions to the energy and not the part of *ϕ* defined over the rest of *Ω*. To determine *V* from this expression, the delta function would have to be numerically approximated and this could potentially cause the evaluation of the integrals to be inaccurate. However, having located the zero level set of *ϕ* using the contour plotting algorithm and given the form of the function *ϕ*^{i}(*x*, *y*), a more direct approach can be used. To do this, we formulate the equations over each zero level set *Γ*_{iΔt} to give the following result.

*Integrating over each level set we have*(4.13)

(4.14)

(4.15)

(4.16)

We prove both the theorem and the lemma by considering each of the contributions to the energy in turn.

*The bending energy*. The total bending energy of the multi-layer system is given by equation (4.2). Now, the curvature, *κ*, is defined as ∂*θ/*∂*σ* and the tangent angle by . Hence,(4.17)Therefore, if *κ*_{i} is the curvature on the curve *Γ*_{iΔt}, equation (4.2) becomes(4.18)Using the expression (4.7), it follows that the bending energy is given by

*The foundation energy*. By definition, the vertical displacement *w* is the same as *y* such that *ϕ*(*x*,*y*,*t*)=0. So that(4.19)

*The work done by the load*. We have assumed that each layer undergoes the same end-shortening under the axial load *P* and therefore the work done by the load is the same for one layer as for multiple layers. Thus, typically for *i*=0,

*The work done against friction*. Substituting (4.6) into (4.5) gives

### (b) The energy close to a singularity

To refine this calculation, we briefly consider the form taken by the bending energy close to a singularity where the curve *Γ*_{t} develops an infinite curvature when *t*=. In this calculation, we will assume that the reference layer *Γ*_{0} is a smooth curve, with maximum curvature of *K*>0 occurring at a minimum point where it has a locally parabolic form. As *t*→=1/*K*, the maximum curvature of the curve *Γ*_{t} increases, and hence the bending energy also increases. To determine the asymptotic behaviour of the bending energy, we will approximate *Γ*_{0} locally by the parabola *y*=*Kx*^{2}/2, with −1<*x*<1. It then follows from the analysis presented in §2 that if *t*=1/*K*−*δ* with *δ*>0 be small, then the maximum curvature of the folded layer is given by . Similarly, if *t*=1/*K* then close to the point *x*=0 the second derivative of the profile of the singular layer is given by d^{2}*y*/d*x*^{2}=*K*^{1/3}*Ax*^{−2/3} where *A*=2^{1/3}/3. Note that 1/*δ*=*K*^{1/3}*Ax*^{−2/3}if *x*=*l*≡*A*^{3/2}*K*^{1/2}*δ*^{3/2}. The function d^{2}*y*/d*x*^{2} obtained by a numerical differentiation of the Lagrangian formulation is illustrated in figure 8*a*, where we take *K*=1 and *δ*=0.01. In this figure, we compare this second derivative both with the function *K*^{1/3}*Ax*^{−2/3} and with the maximum estimate of 1/*δ*. When *δ* is small, it is clear from this figure that we can estimate the function d^{2}*y*/d*x*^{2} close to the point *x*=0 byThis estimate allows us to make an asymptotic estimate for the bending energy in terms of *δ*. As the dominant contribution to the bending energy occurs close to the singular region at *x*=0 where d*y*/d*x* is small, we may closely approximate the curvature by the second derivative so that in the limit of *δ*→0Hence, we have the asymptotic estimate that as *t* approaches 1/*K* there is a constant *B* so that *U*_{B}∼*B*(1/*K*−*t*)^{−1/2}. In figure 8*b*, we compare a numerically computed integral of the bending energy with the asymptotic estimate and see good agreement.

### (c) Deformation profile for a particular reference layer profile

To determine the configuration of the multi-layer material, we must determine the stationary values of the energy functional over all possible configurations of the reference layer *Γ*_{0}. This is difficult, in general, as *Γ*_{0} can take many possible forms, leading to non-unique solutions. A study of some of these for a two-layer material, obtained by approximating *Γ*_{0} by cubic splines, is described in Hunt *et al*. (2006). To make progress in the multi-layer problem, we will, for the present, restrict the class of solutions to those for which the reference deflection *Γ*_{0} has a sinusoidal profile. Of course, this is a restrictive class of solutions, but an inspection of figure 2 shows this to be not unreasonable for certain configurations. A fundamental difference between this calculation and earlier approaches is that we only make the assumption that the *reference* layer is a sinusoid but allow other layers to take different profiles. In other calculations, it has been assumed that *all* of the layers have a sinusoid form. This has led to significant underestimates of the total bending energy. Hence, for this calculation, we make the approximation (Edmunds *et al*. 2006) that *Γ*_{0} is given in terms of the arc-length *σ* by(4.20)and proceed to calculate the other layer profiles using the LSM. For convenience, we restrict the calculation to the half-wavelength of total arc-length *L* given by taking 0≤*σ*<*L*. This profile has a local minimum with maximum curvature *K* given by *K*=*Qπ*^{2}/*L*^{2}. We assume, further, that the value of *L* is known *a priori* from the linear buckling length-scale (Peletier 2001; Budd *et al*. 2003) and that the value of *Q* is to be determined. For a given value of *Q*, the horizontal shortening of the layers is given by(4.21)To determine *Q* numerically, we calculate the level set function *ϕ*, evaluate the energy using the formulae in lemma 4.3, and then find for which values of *Q* this is stationary. To calculate *ϕ*, a computational grid is set up over a rectangle by dividing the horizontal region [0, *L*− ] into *M*_{1}−1 equal intervals of size *h* such that *x*_{j}=*jh*, *j*=1, …, *M*_{1}. A similar division is made of the vertical region into *M*_{2} regions of width *h* so that *y*_{k}=*kh*. The value of *h* is chosen to be consistent with the CFL condition, given that the time-step is determined by the layer thickness Δ*t*. In practice, it is often convenient to calculate the level sets at more refined intervals than the layers themselves, so we may take a timestep Δ for calculating *ϕ* given by Δ=Δ*t*/*m* where *m* is a suitable integer. The level set problem is then initialized by taking(4.22)The approximate values of *ϕ* are then determined at time-intervals of Δ by using the methods outlined in §3 and the level sets *Γ* found by using the Matlab contour function. The integrals (4.13)–(4.16) are evaluated using a 16 point Gaussian quadrature rule at every *m* time-steps (corresponding to each layer) and finally the variation of the total energy *V* with respect to small changes in *Q* is found using central differences. The linear dependence of the total energy on the (fixed) load *P* allows us to find *P* as a function of the single variable *Q* without solving nonlinear equations. In particular, setting ∂*V*/∂*Q*=0, we have simply(4.23)A plot of the resulting (, *P*) curve in the case of a material with *n*=220 layers of width 0.1 mm and with *L*=35.9 mm is given in figure 9. Here, we have taken *M*=180 which corresponds to *h*≈0.22.

Using the previous theory, we can understand the qualitative form of this figure. If *Q* is small then each of the terms *U*_{B}, *U*_{F} and are locally quadratic functions of *Q* and hence their derivatives vary linearly with *Q*. However, *U*_{μ} is a locally linear function of *Q* and its derivative is constant. Thus, as *Q*→0 the load *P*→∞ (Budd *et al*. 2003). In contrast, as *Q* increases further, the deflection of the layers furthest from the reference layer becomes more nearly singular. As the maximum curvature of the reference layer is given by *K*=*Q*(*π*/*L*)^{2} and the half thickness *t* of the total is given by *t*=*n*Δ*t*, it follows from §4*b* that the bending energy varies asIn particular, if *Q*_{S}=*L*^{2}/(*π*^{2}*n*Δ*t*) the bending energy *U*_{B} becomes unbounded as *Q*→*Q*_{S} and there is a constant *C* so that . Hence, as the end-shortening only varies slowly with *Q* as *Q*→*Q*_{S}, we conclude from (4.23) that as *Q*→*Q*_{S} we have *P*∼(*Q*_{S}−*Q*)^{−3/2}. Mechanically, the formation of the singularity in the material is equivalent to a dramatic re-stiffening of the whole system, seen also in other folding configurations such as kink banding (Wadee *et al*. 2004). Of course we do not see an infinite load in practice, and for the rock formations considered in this paper we would expect to see some sort of fracturing or plastic behaviour in this limit or alternatively the formation of multiple folds.

## 5. Conclusions

This paper has shown the potential of the LSM when modelling multi-layer problems. The ability to model parallel folding has been presented and it is anticipated that it can be adapted to model other folding scenarios. The LSM naturally encodes the nonlinear terms. The model is designed only to model the formation of the first hump. If extra humps are to be formed it will be important to understand what happens as the layers lockup, and the role that the singularities have in this (Hunt *et al*. 2006). It is interesting that the formation of singularities is purely due to the geometry and yet these singularities then impact on the system in a way which is presently unclear.

## Footnotes

- Received November 1, 2006.
- Accepted January 24, 2007.

- © 2007 The Royal Society