Convoluted accommodation structures in folded rocks

Introduction

Accommodation structures in multilayered folds

Multilayered structures found in the Earth's crust may arise from a sequence of sedimentary deposits, subsequently held together by additional overlying layers. Tectonic plate movement causes in-plane compression of these layers of rocks that can fold them into convoluted patterns. The resulting deformation is not only dependent on the complex mechanical and material properties of each layer, but on the highly nonlinear geometric constraints of them fitting together.

It is widely accepted that a layered system buckles to a wavelength determined primarily by the bulk material properties of the stack [2,3], and as a result less competent layers can be forced to accommodate the overall geometry imposed. In some cases layers can be obliged to fit into very tight or even singular geometries [1,[3][4][5], and as a result various patterns, named accommodation structures [6], may form on a localized length scale. These structures include accommodation thrusts, keel-like hinges [3,Ch. 12], void or saddle reef formation [1,5] and in particular higher-order localised buckles. Two examples of such localised buckles are seen in Fig. 1. Developing a simplified structural model to understand this particular type of accommodation structure is the primary purpose of this paper.

Dodwell et al. [1] model the behaviour of an elastic layer forced by a constant overburden pressure into a simplified V-shaped geometry, where the layer is freely allowed to slip on its limbs. Constrained energy minimization shows that any such solution is convex, symmetric and unique [1]. Since the model gives only convex solutions we see that pressure, by itself, cannot initiate the formation of higherorder folds like those in Fig. 1. Any valid model must therefore embrace not only the complex topological effects of layers fitting together, but also those of slip relative to each other [7, pp. 167-168]. In the inevitable presence of friction, axial loads and/or corresponding constraints hence become part of the story, in the formation of such convoluted accommodation structures.

The model

To encapsulate the formulation of these structures at singular geometries, we construct the simple mechanical model of Fig. 2. A knee joint is made up of two rigid bars of length L/2, hinged by a plastic rotational spring with the static response characteristic shown at the bottom left; no rotation θ can take place at the hinge until the moment in the spring reaches the plastic moment M p , whereupon it rotates against this constant resistive moment. This knee joint represents the mechanics and geometry of a larger fold or competent layer. A weaker accommodating beam of bending stiffness EI is pinned to it at each end, where in-line springs of stiffness k 2 , which are unstressed in the flat state, undergo deflections Q representing relative slip between the components. The accommodating structure is taken as inextensional but is allowed to buckle, as seen at the bottom right, with displacements described by the angle ψ(s) where s represents measures along the buckled length as shown. The system is loaded horizontally by a load P , displacing by an amount ∆, and there is a further spring k 1 at the point of loading to give an initial, finite (pre-buckled) stiffness. The entire system is placed in a bath of pressure q, against which work needs to be done to create voids.

Dodwell et al [1] simultaneously define layer geometry in cartesian coordinates, characterised by the vertical displacement u(x) with respect to fixed horizontal coordinate x, both being measured from the center of the hinge. We note that the ability to switch freely between coordinate frames is an important part of these analytical developments, each description having its benefits. Intrinsic coordinates generally give concise governing equations, but a common fixed coordinate system using cartesian coordinates proves useful in describing voiding between layers [5].

The linear elastic springs of stiffness k 2 (see Fig. 2), control the amount of slip between the layers, and are taken to mimic axial compression in extended lengths of the accommodating layer which are able to unload into the buckle. The modelling of such lengths to put a value to k 2 would depend heavily on frictional properties at the interface and would certainly be non-trivial in practice, mirroring that of the selection of so-called active lengths in pipeline buckling problems [8].

We define u(x) over the set X = [−L/2 cos θ, L/2 cos θ], as well as the contact set Γ = {x ∈ X : u(x) = f (x)} , together with its complement the non-contact set Γ c . Throughout this work we restrict ourselves to consider both convex and nonconvex solutions of u, which are symmetric and have a single interval non-contact set. Finally the parameter = inf{x > 0 : u(x) = f (x)} characterises the length of the separation.

Variational toolkit and comments on modelling

Using variational techniques developed in [1,5], we construct a potential energy function for the system described above, and derive necessary equilibrium conditions for each of the free parameters u, ∆, θ and Q. This leads to representation of u as a four-order nonlinear free boundary problem, which once recast in intrinsic coordinates can be solved as an initial value problem with 3 unknown shooting parameters.

These solutions not only have clear similarities with the geological examples presented in Fig. 1, but demonstrate some characteristic features of such nonlinear systems, infinite critical loads and localization together with interactive buckling on different scales for instance. Similar interactive behaviour has also been identified with various forms of failure in manufactured materials, for example mode interaction in the buckling of sandwich structures [9], local delamination in composites [10] and local and global mode interaction in I-beams under uniform bending [11]. A recent interdisciplinary special issue of Philosophical Transactions of the Royal Society A has been devoted to highlighting similarities and differences between layered structural geometry and behaviour in geology, carbon-fibre composites, and paperboard used for packaging [12].

In the work presented, we see a departure from the trend of modelling advanced geological features with large complex finite element simulations (see [13] for a recent review). As well as having clear visual appeal, such studies have also been geared towards answering a host of important geological questions. However, with the vast range of parameters often involved in such formulations, it is sometimes unclear how their geological counterparts might accurately be represented. Given a sufficiently complex setup and the ability to tweek parameters, it would seem churlish to suggest that any given result might be achieved, but there could nevertheless be some degree of subjectivity in the choices of parameter values. While we fully appreciate the value of such approaches, our objective here is to complement them with a simplified toy model in which the parameters influencing the localised bucking process can be closely defined and carefully monitored. The model is a part of a sequence of such systems [1,5], moving towards a full description of multilayered buckling under the highly nonlinear constraints imposed by layers being required to fit together [14]. We believe that such simplifications are allowing new interpretations of the development of localized buckles in singular geometries, where voiding and related depressurization effects may be important.

Total potential energy and equilibrium solutions

We now construct each contributing energy in turn for the system described in section 1.2, before deriving equilibrium conditions with respect to perturbations of each of the free parameters. The potential energy stored in the layer-parallel spring U IS , the axial springs U A , and the work done in the plastic rotational spring U KJ , follow from standard elasticity and plasticity theories and simple geometric arguments, so that

The elastic layer does work in bending, U B , from a flat configuration and against overburden pressure in forming a void U V . Classical bending theory [15,Ch.1] gives the bending energy over a small segment of the beam ds as dU B = 1 2 Bκ(s) 2 ds, where κ is the curvature and B is the elastic stiffness of the beam. Integrating over all s we find that

If an overburden pressure q per unit length acts on the system, the work done in creating a void is given by q(u − f )dx, integrating over all x gives,

If q is large then U V imposes a severe energy penalty.

We assume the elastic beam u is inextensible. This imposes a length condition on the solution u, which in turn implies compression or extension, Q, of the axial springs. Therefore the inextensibility condition can be written as

The layer u is under hard unilateral constraint from the arms of the hinge; therefore the contact condition u ≥ f , is imposed. The total potential energy for the system is the sum of all contributing potential energies derived above, minus the work done by the load in end shortening, P ∆, therefore

Here W (u, θ, Q) defines the potential energy functional of the layer u, given by

which is constrained by the contact condition u ≥ f and the inextensibility constraint (2).

Equilibrium solutions

We now derive a series of equilibrium equations, by finding stationary conditions for the potential energy with respect to perturbations of ∆, θ, Q and u.

Stationarity with respect to ∆ and θ

Two equilibrium equations are derived by considering stationarity of V with respect to ∆ and θ, i.e.

and

Solving these two equations simultaneously provides two solution paths. Firstly the trivial solution for which θ = 0,

and secondly various non-trivial equilibrium solutions (θ = 0)

We note that ∂W/∂θ depends on the solutions profiles u, for which we now derive equilibrium conditions.

The Euler-Lagrange and Hamiltonian equations

We now apply Kuhn-Tucker theory [16, pp.249] and standard Lagrange multiplier theory [16], to derive necessary conditions for a stationary solutions of V with respect to perturbations of u, subject to the contact constraint u ≥ f and the inextensibility constraint J (u) = Q.

A standard variational approach considers stationarity of V with respect to vertical perturbations, i.e. u(x) → u(x)+εϕ(x), for some small constant ε and suitable test function ϕ. As a result the corresponding Euler-Lagrange equations gives the condition for a vertical force balance of the system [1]. Equally, one may insist on stationarity with respect to horizontal perturbations, i.e u(x) → u(x + εϕ). In this case the resulting differential equation is the static Hamiltonian for the system [1], and can be physically interpreted as a horizontal force balance.

Firstly, applying Kuhn-Tucker theory, stationarity with respect vertical perturbations leads to the Euler-Lagrange equation

where µ, a non negative measure, and λ ∈ R are Lagrange multipliers, arising from the contact constraint u ≥ f and the inextensibility constraint (2) respectively. Physically µ can be interpreted as the vertical component of the reactive contact load with the boundary f . Kuhn-Tucker theory represents this by the complementarity condition,

Therefore µ is zero wherever there is no contact, and equals a discontinuous shear force at points of delamination. The constant Lagrange multiplier λ can be interpreted as the in-plane load required to enforce the inextensibility constraint (2), and as with µ is to be found as part of the solution.

There is also a first integral of (9),

where

which follows by imposing continuity of first and second derivatives at x = .

Secondly we derive the Hamiltonian for the system by finding the stationary solutions with respect to horizontal perturbations, i.e. u(x) → u(x + εϕ). This yields a second fourth order nonlinear differential equation,

which describes the horizontal force balance on the beam. This also has a first integral, given by

(13) Here the right hand side follows from applying the fixed symmetric boundary conditions,

where u 0 and Λ ∈ R are free parameters to be found as part of the solution.

Intrinsic representation

Equations (11) and (13) can be written in terms of intrinsic coordinates, characterized by the arc length s, measured from the point of symmetry, and the tangent angle ψ with the horizontal. Recalling the relationship between both sets of coordinates [14],

equation (9) can be rewritten as

which is equivalent to

Similarly we rewrite (12), so that

Multiplying (17) by cos ψ and (18) by sin ψ and adding together, yields the equilibrium equation for u normal to the layer, given by

Similarly multiplying (17) by sin ψ and (18) by cos ψ and subtracting, one obtains an equilibrium equation for u parallel to the layer,

Substituting the free boundary condition at x = , u( ) = tan θ, u x ( ) = tan θ and u xx ( ) = 0, into the above we derive the following condition

rearranging for u 0

This condition now implies an additional fixed boundary condition at x = 0.

Stationarity with respect to Q

The condition for stationarity with respect to perturbations of Q, is obtained in a similar manner, so that

Physical interpretation of equilibrium equations

The equilibrium equations give conditions under which the energy V is stationary with respect to perturbations or displacements of each of the free parameters ∆, θ, u and Q. Each of these equations describes the corresponding force or moment balance for equilibrium with respect to displacements of each of these parameters. Mathematically they represent the dual elements associated with each of these free parameters, loads with displacements and moments with rotations. The equilibrium equation

describes the force balance associated with perturbations of ∆, and is interpreted as a balance of the reactive force exerted by the horizontal spring against the in-plane force P . The Euler-Lagrange equations and Hamiltonian equation describe force balances of u vertically and horizontally respectively. We can identify each term in both of these equations as follows.

Noting that Q is also written as a constraint equation on u, (2), it follows that the above equilibrium equation is equivalent to

This makes physical sense as λ(Q) is the axial load in the beam so k 2 Q is the force in each spring undergoing a displacement Q.

Finally the equilibrium equation (8),

gives the stationary condition with respect to rotational perturbations about the hinge θ; it can therefore be interpreted as the moment balance about the hinge. Additionally ∂W/∂θ is the contribution to the total moment made by the elastic beam. This observation is important, as it allows us to derive an analytical expression for ∂W/∂θ, which substantially simplifies the following analysis.

θ Figure 3. Free body diagram of whole system plus block of pressurized fluid. Figure 3 shows the moments acting about the hinge. The total moments in the anticlockwise direction are given by,

whereas moments in the clockwise direction are

Comparing each of these forces and moments (24) and (25), we note that P 2 = 1 2 BΛ 2 + λ and M 2 = Bκ(x = 0) = BΛ; therefore it follows that

Comparing this moment balance with equation (8), we observe that

3. Numerical solutions

Shooting for stationary solutions of u

We now obtain numerical solutions of u, for a range of values of θ and physical parameters k 2 , q and B. We note that for particularly convoluted solutions, u x becomes infinite, which means such solutions are unattainable in a cartesian framework. This problem is eliminated by solving the equilibrium equation normal to the layer (19). This emphasises the importance of calculating both the Euler-Lagrange and Hamiltonian equations.

Numerical results can be obtained by reducing this free boundary problem to an initial value shooting problem, and by rewriting (19) as a four-dimensional first order system,

A suitable numerical integrator (e.g. MATLAB's ode45) can then be used to shoot from the symmetric section

so that the length constraint (2) holds. Here u 0 is derived from the condition (22) and Λ, and λ are free shooting parameters to be found numerically. A three dimensional search routine was constructed [14] by nesting MATLAB's built-in function fminsearch. This proved efficient, given a reasonable initial guess.

Several different outcomes of this shooting method, from different initial guesses, are shown in Fig. 4. Unlike the earlier model [1] solutions are no longer unique, and up-buckled shapes of increasing complexity start to appear. Fig. 4(a), and intersects the red straight line just once, at point a. There is thus a unique convex deflection profile for every value of θ. However the same is not true for the up-buckled solution, seen as the dot-dash curve in Fig. 5. For the given value of θ the straight line −k 2 Q intercepts λ(Q) twice, at b and c, meaning that two different up-buckled shapes exist for the same value of θ.

The shapes of the curves also suggest that, whereas the convex (solid blue) curve is always going to be crossed once by a straight line passing through the origin at Q = λ = 0, the up-buckled curve typically is crossed either twice or not at all. This is seen in the plots of Fig. 6, showing the straight and curved loci for the up-buckled shape when θ is (a) small; (b) moderate, and (c) relatively large. It is clear that up-buckled shapes only exist for θ values greater than the threshold value θ min of Fig. 6(b). This significant point is highlighted by the load-deflection plots of the next section.

Once equilibrium solutions u to equation (26) are attained for each value of θ, an approximation for ∂W/∂θ can be constructed. Numerically this is done using the trapezoidal rule to calculate W for each θ, and then a finite difference approximation for ∂W/∂θ. Equation (8) can then be used to find the corresponding value of P and the following response curves constructed.

Overall system response

The differences between the convex and up-buckled equilibrium solutions can be seen most clearly in plots of load P against end-shortening ∆, as in Fig. 7, which embrace all possible values of θ as the displacement evolves. Fig. 7(a) shows the equilibrium paths for the convex shape of Fig. 4(a). Two possible equilibrium states exist, the fundamental path seen as a red solid line where the arms are completely aligned (θ = 0), and buckled path in the bent configuration (θ = 0) shown as

c.

d. blue/dashed. With increasing P , θ → 0 in the buckled state and these paths approach each other asymptotically, but never meet. As the load falls from this position, θ grows, until at some positive value of load the path reverses direction and follows a route with P increasing and θ continuing to grow. While falling the buckled path is unstable, but under conditions of controlled end-shortening it would restabilize at the point where ∆ is a minimum. Thus the possibility arises of a dynamic jump at constant ∆, from the fundamental state to the now-stable rising post-buckled path. To initiate this an energy hump needs to be overcome, but at high loads this would be expected to be relatively small. Because an elastic layer can only bend to finite curvature, to adopt this deflected state there must be significant slip between the layers taken up by the springs k 2 . The equilibrium path for the up-buckled shape, shown in Fig. 7(b), is similar but has important differences. In particular, although it is gets close to the fundamental path over some of its length, this is only over a finite load range. The solution reaches a limit load (seen here as a cusp) at the θ min of Fig. 6(b), where it separates into the simple up-buckled shape of Fig. 5(b) and its counterpart 5(c). The simple up-buckled shape eventually restabilizes with increasing θ much like the convex mode, whereas the more open configuration is on the way to contacting with itself and becoming the teardrop of Fig. 5(d). Here the path is stopped before this state is reached, because the delamination has extended to the ends of the rigid arms. It should be noted that, although the paths of Fig. 7(a) and (b) occupy the same region of P − ∆ space, they would in fact be far from each other in the full displacement space, separated in particular by different values of Q. In the absence of any further parametric variation such as a changeable imperfection there is, we think, no likelihood of any bifurcation linking the two sets of paths.

Maxwell stability criterion

For load versus corresponding deflection plots like those of Fig. 7, the total energy stored in the system can be inferred from the area under the appropriate graph [15]. The schematic of Fig. 8 shows relevant areas for three possible equilibrium states that can arise at a typical, constant, value of applied end-shortening ∆. Energy on the fundamental path, at A, simply involves the triangular area shown in 8(b). Energy at B adds to this the area E 1 between the fundamental path and the high load/small θ region of the buckled path (see 8(c)), and for the energy at C area E 2 is similarly subtracted, as 8(d). Note that, even when the paths approach asymptotically but never meet as in Fig. 7(a), these areas can often be found by direct integration (see for example [17]). E 1 represents the input required to get (a) Convex profile: the buckled (blue/dashed) solution and the fundamental (red/solid) paths approach each other asymptotically as P → ∞ but never meet. (b) First-order up-buckled profile: the buckled solution reaches a limit point (seen here as a cusp) at a finite load value where the response separates into shapes with growing and diminishing void sizes. over the energy hump at B, while E 2 is the difference between this maximum and C, which under constrained end-shortening would be an energy minimum.

There is of course no bifurcation point on the fundamental path, and so A is an energy minimum at all load values and therefore technically always stable. High loads however, would push the system into a dangerous state of metastability; as P increases, E 1 becomes increasingly thin and shrinks in area. In such situations it is common to adopt the Maxwell stability criterion [18, p. 53] where stability rests with the global energy minimum rather than local minima. This marks out the Maxwell displacement ∆ M -the point at which the energy at A equals the energy at C and E 1 = E 2 -as a lower bound for the point of buckling. For applied ∆ < ∆ M , the global energy minimum rests with the fundamental equilibrium path: for ∆ > ∆ M it rests with the buckled state. Depending perhaps on the amount of available energy from outside disturbances relative to the size of the energy hump E 1 , this may or may not be a good estimate of the actual buckling load; however in the absence of any other ballpark estimate it provides some kind of framework within which to assess the buckling process. The Maxwell displacement itself can rest with the convex or the up-buckled solution, depending on the choice of parameters. In practical terms the actual buckling mode is also likely to be heavily influenced by knock-on effects of initial imperfections [12].

Parametric variations

The plots of Figs. 7 and 8 are of the load P against its corresponding deflection (end-shortening ∆), and give a good overview of the system response in just two dimensions. Whether P or ∆ is chosen as the parameter to be controlled the same equilibrium curves are produced, although their stability is likely to depend on this choice [15]. Other parameters of the system could also be varied, and we now review the effect on the system response of two such variations. Consider for example the plots of Fig. 9, showing the response of the up-buckled shape of Fig. 4(b) at three different values at M p . While the convex solution always takes the same form whatever the chosen M p , this is clearly not the case here. Focussing on the limit point (seen here as a cusp) where the load is locally a maximum and θ = θ min , increasing M p has the effect of significantly increasing its corresponding load in comparison with the remainder of the response. We note that ∂W/∂θ of equation (8) is independent of M p , and therefore comes up with the same contribution whatever its chosen value; θ min is thus the same for all three cases. But the response for the higher values of M p begins to look much more like that of the convex shape than for low values; the relative energy hump to be overcome to initiate buckling clearly drops as M p is increased. We might therefore expect the Maxwell stability criterion to carry more relevance for high than for low values of increase in k 2 raises the load, but this time simultaneously reduces the value of θ min at the limit point. It also causes the post-buckling curve to stiffen more quickly, quickly becoming almost parallel to the fundamental path as shown. Movement in the k 2 springs is a necessary component of the up-buckled shape, equivalent of slip between the layers in the geological context, and as these springs stiffen up more work is being done in them to accommodate the same shape. For a given θ, the higher the value of k 2 the more quickly a buckling load will develop in the accommodating beam and hence the corresponding value of θ min is lowered.

Concluding remarks

This paper presents a simplified structural model to describe the formation of convoluted accommodation structures in tight geometries that arise naturally in layered systems. It has been formulated specifically to highlight the important nonlinear geometric constraints imposed by the need for layers to fit together, and demonstrates mechanisms by which the less competent layers accommodate such tight geometries imposed by their more competent counterparts. The work extends previous work by the authors [1], to include the more general formation of the V-shaped singularity, as well as important effects of axial constraint. The pin joints at the ends of the structure, along with the introduction of springs k 2 and the inextensibility, are useful simplifications but remain some way from reality. Relative slip is concentrated into deflection Q at the ends of the beam rather than being spread out over a finite, extensible, contact length as would be the case in more practical circumstances. The situation is not unlike that of a straight long pipeline on a flat bed carrying hot oil or gas, which can release heat-induced axial compression by up-buckling against gravity over a finite length [8]. Lengths of pipe at each end outside the buckled region stretch and slip along the bed to feed the buckle. In the absence of friction, these so called active lengths could be infinitely long and feed ad infinitum into the buckle. With realistic friction between the bed and the pipe this mechanism is restricted to finite lengths, over which the drop in compressive load from the unbuckled to buckled region is offset by friction from the bed. For pipelines these active lengths can be long in comparison with the buckled length, but in the geological setting friction would be expected to dominate and corresponding active lengths be relatively short. We note that the inclusion of inter layer friction is important for future models for multilayered structures; however it does offer a difficult challenge, because the normal contact forces change as the response develops.

One parametric change not pursued here is that of pressure q. For the present parameter values the teardrop shape of Fig. 5(d) never fully develops, as the separated region grows to the full length of the model beforehand. However, with relative increases in q this limitation may well disappear. Pressure forces acting normal to the surface of layer should tend to reduce the separation length involved, and could perhaps even generate tension in the accommodating beam. We leave this development for future work.

The paper draws on a philosophy found also in a companion paper in this volume [19], in selecting wherever possible linear elastic constitutive laws to express material behaviour. Thus a small element of the accommodating beam has linear bending stiffness B, and frictional effects are modelled by linear elastic springs k 2 . The exception is the rigid-plastic rotational spring M p . This could have been taken as a stiff linear elastic spring; bifurcation at a high rather than infinite value of P would then occur, but the post-buckling response remains much the same [14].

Linear constitutive laws may need replacing if good experimental comparisons are to be achieved, but they do have certain modelling advantages over more complex rheologies. First, they are usually more straightforward to analyse. But most importantly, they highlight the roles that nonlinear geometry and the associated potential for many alternative solutions can play in such systems. Modelling has sometimes failed in the past not because of the formulation, nor the resulting gov-erning equations, but because of poor choice of solution. Extra criteria for selection are sometimes necessary or useful. The buckling of the long axially compressed cylindrical shell gives a good example of just such a circumstance [20].

Because nonlinear geometry is central to the formulation, the system is capable of undergoing large deflections and rotations. Governing equations for the layer are derived in both cartesian and intrinsic coordinates, the latter meaning that for convoluted profiles as seen in Fig. 4, numerical solutions do not blow up. We again recognise the purely elastic assumptions are not realistic, but more complex rheological models might likewise be subject to error. Real geological systems would involve an intricate balance of nonlinear constitutive laws, strain history and temperature, pressure, chemical coupling etc. [21], such that precise analysis would appear an impossible dream. Real progress is perhaps best made by a combination of modelling on many levels, with proper interaction and cross-fertilization of ideas across all elements of the spectrum.