A director theory for visco-elastic folding instabilities in multilayered rock

A model for finely layered visco-elastic rock proposed by us in previous papers is revisited and generalized to include couple stresses. We begin with an outline of the governing equations for the standard continuum case and apply a computational simulation scheme suitable for problems involving very large deformations. We then consider buckling instabilities in a finite, rectangular domain. Embedded within this domain, parallel to the longer dimension we consider a stiff, layered beam under compression. We analyse folding up to 40% shortening. The standard continuum solution becomes unstable for extreme values of the shear/normal viscosity ratio. The instability is a consequence of the neglect of the bending stiffness/viscosity in the standard continuum model. We suggest considering these effects within the framework of a couple stress theory. Couple stress theories involve second order spatial derivatives of the velocities/displacements in the virtual work principle. To avoid C-1 continuity in the finite element formulation we introduce the spin of the cross sections of the individual layers as an independent variable and enforce equality to the spin of the unit normal vector to the layers (-the director of the layer system-) by means of a penalty method. We illustrate the convergence of the penalty method by means of numerical solutions of simple shears of an infinite layer for increasing values of the penalty parameter. For the shear problem we present solutions assuming that the internal layering is oriented orthogonal to the surfaces of the shear layer initially. For high values of the ratio of the normal-to the shear viscosity the deformation concentrates in thin bands around to the layer surfaces. The effect of couple stresses on the evolution of folds in layered structures is also investigated. (C) 2002 Elsevier Science Ltd. All rights reserved.


Introduction
It has long been recognized that most sedimentary rocks were originally deposited in sequences of soft horizontal layers at the bottom of shallow seas and hardened over time, resulting in a characteristically layered structure. Layering may be also be caused by purely mechanical, hydro-mechanical or chemomechanical means (e.g. see Williams (1972) and Robin (1979) for geological evidence and conceptual models for the role dissolution and volume replacement on the present day appearance of fold structures; Ortoleva (1994) presents models for chemical drivers and metamorphic differentiation of various kinds).
A significant challenge in many instances is to explain the evolution of the enormous variety of shapes and forms observed in deformed sedimentary-and crustal rock structures. Observed structures include folds, shear-and kink bands and fractures. The strains under which these structures formed are typically high and the strain rates are low (10 À14 s À1 ), viscous, elastic, and brittle effects influence the observed structures. In this paper we concentrate upon folding as an example of a geological process which is at once simple to describe qualitatively, yet rich in complexity when a quantitative mechanical explanation is sought. There have been numerous steps taken on the road to such an explanation, including the seminal work of Biot (1937Biot ( , 1957Biot ( , 1964Biot ( , 1965a and Ramberg (1963), and the relentless consideration of greater levels of rheological complexity (Chapple, 1969;Fletcher, 1974Fletcher, , 1982Johnson and Fletcher, 1994;Hunt et al., 1997;Hobbs et al., 2001;M€ u uhlhaus et al., 2002) and sophisticated numerical modeling (Schmalholz and Podladchikov, 1999;Vasilyev et al., 1998). Different styles of instabilities as a function of structural and material parameters in a prestressed cohesive layer on a viscous substratum are rigorously analysed within the framework of a linear instability analysis by Leroy and Triantafyllidis (1996). The theory finds application in geological settings far removed from buckling of single rock strata including, for example, the Himalayan orogeny (Caporali, 2000) and possible buckling modes of the oceanic and continental lithosphere (McAdoo and Sandwell, 1985;Martinod and Davy, 1992).
We attempt to unify much of this work with a description of the buckling of layered visco-elastic materials, accounting for the influence of the microstructure, in both small deformation theoretical analysis and in large deformation numerical simulation. The benefit in developing a clear understanding of folding and other geological deformation processes is in being able to infer the physical conditions of the deformation and the constitutive behaviour of the rocks at the time of deformation. The difficulty in extrapolating small-scale, short-duration laboratory deformation experiments to geological time and length scales often means that it is difficult to say whether particular fold formations were produced by viscous, viscoelastic, or visco-elastic-plastic deformation styles, or whether the presence of layering influenced the constitutive behaviour. A computational model capable of simulating large-amplitude folding of a general class of layered, visco-elastic materials is therefore an important first step in being able to use fold characteristics to determine ambient conditions during deformation.
We develop a mechanical model including a large deformation formulation for multi-layered rock. The formulation intended for materials with fine internal layering, which can be described by a single director orientation. This constitutive model is specifically designed for geological deformation problems involving very large deformations. Although there are more general descriptions possible, this formulation is, in fact, very broadly applicable to crustal rocks, where the preponderance of layering arises from deposition of one rock type onto another under gravity.
We revisit the basic finite element formulation for viscous materials and demonstrate how the standard element vectors and matrices can be extended to include anisotropy. We next describe a computational method capable of following the evolution of macroscopic interfaces and the internal layering direction introduced in the constitutive relationship. The particle-in-cell (PIC) finite element method, as this technique is known, is a hybrid scheme which falls between the finite element method and a purely Lagrangian particle method such as the discrete element method. PIC is derived from standard finite elements but includes moving integration points to carry director orientation information and other history variables.
We explore scenarios from global to internal buckling in non-linear finite element studies. These show that buckling can be induced at much lower viscosity contrasts between the matrix and the embedded beam or plate than would be the case for isotropic materials. Numerical solutions based on the standard continuum formulation assumed initially become may become unstable if the contrast between the normal--and the shear viscosity becomes very severe. If the layered material is embedded in an incompressible medium this effect is not as pronounced as in the case of free boundaries or a strongly compressible embedding medium. Although the applications in this paper focus mainly on the former case we present, as an outlook and extension of the theories presented in sections one and two a couple stress formulation, which is numerically robust, independent of the nature of the embedding medium. The physical significance of couple stresses in the layered material is in the consideration of a finite thickness of the individual layers and the allowance for stress variation across the thickness of the layers. These stress variations are--in essence in the same way as in the classical beam bending theories--statically equivalent to couple stresses. We illustrate the couple stress model by means of numerical and an analytical solution (the latter is derived in the Appendix D) of simple, finite shear of an infinite layer and revisit the folding problem in the light of the couple stress theory.

Mathematical formulation
From a mechanical point of view, the salient feature of layered materials is that there exists a distinguished orientation given by the normal vector field n i ðx k ; tÞ of the layer planes, where (x 1 ; x 2 ; x 3 ) are Cartesian coordinates, and t is the time. Initially we assume linear viscous behaviour and designate with g the normal viscosity and g S the shear viscosity in the layer planes normal to n i . The orientation of the normal vector, or director as it is sometimes called in the literature on oriented materials, changes with deformation. A comprehensive account of director theory in the context of liquid crystals can be found in de Gennes and Prost (1995, p. 100ff and 150ff). Using a standard result of continuum mechanics, the evolution of the director of the layers is described by where L ¼ D þ W is the velocity gradient, D is the stretching and W is the spin. The superscripted n distinguishes the spin W n of the director n (the unit normal vector of the deformed layer surfaces) from the spin W of an infinitesimal volume element dV of the continuum. The 2D matrix representation of (1) as needed for our computational applications is represented in Appendix B for easy reference. We define a corotational stress rate as: Again the superscripted n distinguishes the stress rate _ r r n as observed by a material observer corotating with the director n from the material stress rate _ r r observed by a spatially fixed observer.

Specific viscous and visco-elastic constitutive relations
We consider layered, viscous and visco-elastic materials. The layering may be in the form of an alternating sequence of hard and soft materials or in the form of a superposition of layers of equal width of one and the same material, which are weakly bonded along the interfaces. We designate the normal shear modulus and the normal shear viscosity as l and g respectively; the shear modulus and the shear viscosity measured in simple, layer parallel shear we designate as l S and g S .
In the following simple model for a layered viscous material we correct the isotropic part 2gD 0 ij of the model by means of the K tensor (see Appendix A for derivation) to consider the mechanical effect of the layering; thus where a prime designates the deviator of the respective quantity, and þ n j n k d il þ n i n l d kj þ n j n l d ik Þ À 2n i n j n k n l Á : Similarly, a visco-elastic constitutive relationship for a layered medium may be written as: where _ r r n ij is the corotational stress rate introduced at the beginning of this section. We could have equally well used the Jaumann derivative of r which is obtained by replacing W n in (2) by the spin of an infinitesimal element of the continuum W. A remark on the notation: We use index notation which is less ambiguous than symbolic notation when vectors, second and fourth order tensors (such as K ijkl ) appear simultaneously in the equations.

Finite element formulation
The constitutive relationships derived in the previous section translate naturally into standard finite element matrix formulation for almost incompressible materials. The generic form of the finite element model for slow flows (zero Reynolds number) of an incompressible fluid reads: K is the so-called global effective stiffness matrix which contains all the material property and history parameters, G is the divergence expressed in matrix form, u and p are the unknown velocities and pressure respectively, and F is a vector of driving terms comprising body forces and in the case of non-linear problems so called out of balance force terms as well. A specific form of F will be specified below.
The matrixes K and G are global matrixes composed in the usual way of elemental matrixes; in the following we designate element matrixes and vectors by a superscripted E. The components of an element stiffness matrix may be written as The matrix B consists of the appropriate gradients of interpolation functions which transform nodal point velocity components to strain-rate pseudo-vectors at any point in the element domain. The constitutive operator corresponding to (3) is composed of two parts C ¼ C iso þ C layer representing the isotropic part of the constitutive model and a correction term considering the influence of layering. In two dimensions, in which D 0 ¼ 2n 2 1 n 2 2 and D 1 ¼ ðn 1 n 3 2 À n 3 1 n 2 Þ. A visco-elastic equivalent of the viscous equations can be obtained by inserting the corotational rate (2) into (5) and subsequently replace the time derivative by the corresponding first order difference quotient. After rearranging, the constitutive equations can be written in the form (3) with the viscosities replaced by effective viscosities and an additional term which we define below representing the stress history and the influence of the stress rotation (see Eqs. (10) and (11)). The result reads: Details of the derivation of (9) are represented in the Appendix B (Eqs. (C.1) and (C.2)). The effective viscosities read: where the superscripted t refers to the previous time step and Dt is the time increment. The parameters a ¼ g=l and a S ¼ g S =l S are the relaxation times for pure-and simple shear respectively. The time increment is limited by a local Courant condition ensuring that during one time step a particle does not travel further than a typical element dimension (this limitation ensures that changes to the stiffness matrix, K, are small during a timestep). The position of the particle is updated using a standard second or fourth order Runge-Kutta integration scheme to integrate the local particle velocity (interpolated from the nodal point velocities using element shape functions in the usual way) forward in time to the new position. A corrector step to improve particle positions, and revise the computed nodal point velocity solution can help improve accuracy.
In Moresi et al., the time increment Dt used for the advection of material points usually differs from the time increment used for the integration of the stress history. This is necessary to ensure that the timescale of interest for elastic processes, and hence the effective viscosity, is determined by the physics of the problem, and not by the mesh-dependent Courant condition.
The superscripted parameter 0 6 b 6 1 in the corotational term in (10) was introduced to provide flexibility in the numerical treatment of the problem: during a predictor step we put b ¼ 0 and during subsequent iterations we put b ¼ 1.
Furthermore, we need to modify the force vector F tþDt for elasticity and also correct for the observer rotation (see Section 2 Eqs. (1) and (2)). To achieve this we introduce the auxiliary vector s as defined by (10) and write: where F E ext is the external load vector. After each time step the incremental solution may be improved iteratively. During iteration the director spin, particle positions etc are replaced by their values at t þ Dt and are continuously updated until the increment of the velocities between two successive iterative steps are sufficiently small in the sense of a suitable norm. The above strategy allows one to modify existing codes for viscous materials without major interference with the rest of the code. There are many possibilities for refinement but for the purpose of this paper, for the examples presented in Section (5), this simple formulation is perfectly applicable.

The particle-in-cell finite element method
Some difficulties arise in devising a practical implementation of the finite element formulation described in Section 3 for the large deformation modeling of layer folding. In particular, since the C matrix is a continuously evolving function of position through its dependence on director orientation, it is necessary that we are able to track an evolving vector function of the material during deformation.
We have therefore developed a hybrid approach--a PIC finite element method that uses a standard Eulerian finite element mesh (for fast, implicit solution) and a Lagrangian particle framework for carrying details of interfaces, the stress history etc.
Our PIC finite element method is based closely on the standard finite element method, and is a direct development of the material point method of Sulsky et al. (1995). The standard mesh is used to discretize the domain into elements, and the shape functions interpolate node points in the mesh in the usual fashion. The problem is formulated in a weak form to give an integral equation, and the shape function expansion produces a discrete (matrix) equation. For the discretized problem, these integrals occur over sub-domains (elements) and are calculated by summation over a finite number of sample points within each element. For example, in order to integrate Eq. (7), over the element domain X E we replace the continuous integral by a summation In standard finite elements, the positions of the sample points, x p , and the weighting, w p are optimized in advance. In our scheme, the x p 's correspond precisely to the Lagrangian points embedded in the fluid, and w p must be recalculated at the end of a time step for the new configuration of particles. Constraints on the values of w p come from the need to integrate polynomials of a minimum degree related to the degree of the shape function interpolation, and the order of the underlying differential equation (e.g. Hughes, 1987). These Lagrangian points carry the history variables including the director orientation which are therefore directly available for the element integrals without the need to interpolate from nodal points. Moresi et al. (2001) give a full discussion of the implementation of the PIC finite element scheme used here including full details of the integration scheme and its assumptions.

Numerical simulations
We present an example of a simulation of folding of a layer of anisotropic visco-elastic material sandwiched between two isotropic layers of equal viscosity on each side (Fig. 1). To accommodate the shortening of the system, the outermost isotropic layers are compressible. In benchmarking, this sandwich of incompressible and compressible embedding material was found to give good agreement with analytic results assuming an infinite domain . The Deborah number, De (relaxation time/process time, e.g. length of block divided by axial velocity prescribed on the sides of the block, Fig. 1), is defined as De ¼ ðgV Þ=ðlLÞ. We consider four cases in which De ¼ 0, 0, 0.0, 0.5 and 5 respectively. The first two examples are purely viscous. In the last two, the visco-elastic cases, we assume that g=g S ¼ l=l S , and an unstressed initial condition. We trigger the folding process by assuming initial perturbations of the form dx 2 ¼ 0:05h cos ðqx 1 Þ, 0 6 x 1 6 2.
For case 1 (Fig. 2), we assume that g=g M ¼ 10 and g=g S ¼ 10. We observe that the growth rate is higher for smaller wavenumber in the infinitesimal deformation limit, and this persists with finite amplitude deformation. Growth at wavenumber q ¼ 2p is significantly more developed after 40% shortening than growth at q ¼ 10p. The shading indicates the extend of the local director rotation as measured by n 1 j j; blue ¼ small; red ¼ large. At q ¼ 2p the director rotations are concentrated around to the surfaces of layer 3 indicating a surface wrinkling or surface instability-like phenomenon. For isotropic beams buckling is insignificant at g=g M ¼ 10 (Biot, 1965a).
For case 2 (Fig. 3) we assume that g=g M ¼ 100 and g=g S ¼ 10 we observe the same overall trend as case 1: high wavenumber perturbations do not grow as fast as low wavenumber. However, we also observe that low wavenumber modes are excited in the finite deformation limit irrespective of the perturbation wavenumber. The perturbation causes a noticeable secondary variation in the interface deflection.
Turning to visco-elastic materials (cases 4 and 5) with De ¼ ðgV Þ=ðlLÞ ¼ 5 and 0.5 respectively we observe for the large wavelength perturbation no or at least no significant mode coupling: the initial perturbation is amplified in an unstable fashion. The growth coefficient of the homogeneous rectilinear ground solution is L=V ¼ 0:2 would produce an amplification of the initial perturbation of about 4.5 at 40% shortening which is much smaller than what we observe for the long wavelength perturbations in case of the viscous layers and for both short and long wavelength perturbations in the visco-elastic cases. A second, dominant mode emerges at q ¼ 10p in case 4 in a rather spectacular way between 25-35% shortening (Fig.  4). A more viscous case with De ¼ 0:5 is represented in Fig. 5. Fig. 1. Initial geometry for the folding experiment. Layer 1 is compressible, viscous (g M ) background material, layer 2 is identical to layer 1 but incompressible (see text for an explanation), layer 3 is the test sample: visco-elastic (l; l S ; g; g S ) with a director orientation (n). The anisotropic layer contains small perturbations to the otherwise horizontal internal layering. V ¼ 10 is constant during any given experiment and unchanged between different experiments. The length of the block is L ¼ 2 and the width of the central layer 3 is h ¼ 0:12.

Couple stresses
We conclude the main body of this paper with an excursion into the potential influence of couple stresses on the mechanical behaviour of layered materials. We restrict ourselves to viscous materials, mainly for algebraic convenience; the extension to more general constitutive relationships is reasonably straightforward however.
Couple stresses are significant in situations where the gradient of n i changes strongly over a short distance (limiting case: disclination).
In such cases we have to take the variations of the normal stresses across the layer cross sections into consideration (e.g. M€ u uhlhaus, 1993). The couple stress theories (see e.g. Mindlin and Thiersten, 1962;M€ u uhlhaus, 1993;M€ u uhlhaus and Aifantis, 1991a,b) provide a convenient framework for the consideration of stress fluctuations on the layer-scale without having to abandon the homogeneity properties of the anisotropic, standard continuum approach. Generalised continuum models--so called gradient plasticity models--motivated by dislocation structures in single polycrystalline metals were proposed by Aifantis (1984Aifantis ( , 1987. In these models the stress tensor is symmetric and non-local hardening effects e.g. due nonlocal interactions between dislocation sources and obstacles or immobile dislocations are considered to the lowest order in terms of a Laplacian of the equivalent accumulated plastic strain. In the present case the couple stress enhancement leads naturally to the superposition of visco-elastic bending stiffness on our standard continuum model (6). In connection with layered materials the internal length scale introduced by the couple stresses is proportional to the layer thickness (ranging from microns to kilometers in geological applications) and to the differences between the viscosities and shear moduli governing pure and simple shear respectively (see e.g. M€ u uhlhaus, 1993). In layered materials the explanation why the stress tensor is non-symmetric in couple stress materials is straight forward: in a continuum description the stresses represent average values over multiples of the layer thickness. In bending the shear stress obtained by averaging normal to the layering is different in general from the shear stress parallel to the layering. The latter may even be zero--for instance, in the case of a stack of perfectly smooth playing cards (a standard continuum model would break down in this case). Within the framework of a couple stress theory one considers the variation of the normal stress across the layer thickness (in much the same way as in the standard engineering beam and plate theories), introduces statically equivalent couple stresses balancing the difference between the shear stresses (Fig. 6).
The couple stress tensor l (moment per unit area) is conjugate in the rate of energy to the rate of curvature j. In the context of layered materials, a natural choice for the rate of curvature reads: Fig. 6. Evolution of folding in viscous layer. Isotropic embedding material has viscosity 1, layer has, normal viscosity 1000. Results are shown for perturbation to the director orientation with wavenumber q ¼ 10p at 40% shortening. (a) layered beam, g=g S ¼ 1000 (b) couple stresses included, same as (a) otherwise. (Layer thickness)/(beam thickness) ¼ 0.2. In this example the internal length was chosen relatively large hence the Cosserat terms act on the larger wavelength, increase the larger length scale. The effect of couple stresses in connection with higher values of the shear viscosity and in particular in the context of visco-elasticity requires further investigation.
In 2D deformations in the (x 1 ; x 2 ) plane if the director is parallel to the x 2 -axis, (14) reduces to j 31 ¼ v 2;11 and j 32 ¼ v 2;12 : The power balance for a couple stress medium reads: where the superscripted dot means differentiation with respect to time, _ d d i is the angular momentum, b i and m V i are volume forces and couples respectively, are the surface stress and couple stress tractions and N i is the unit outward normal vector to the surface S of the body. There are a number of complications associated with the application of couple stress theories in finite element analyses. For instance rotations cannot independently varied on surfaces where the normal component of the displacements or velocities are prescribed. In this case the velocity gradients on the surface have to be decomposed into normal and surface parallel components. The surface parallel part will--after application of the surface divergence theorem--produce a contribution to the stress traction (see e.g. M€ u uhlhaus and Aifantis (1991a, b) and Fleck and Hutchinson (1997) for details of the analysis within the context of a strain gradient--and a couple stress theory respectively). Another difficulty arises from the fact that the volume integrals in (16) contain second order velocity gradients so that the shape functions in a finite element model must be continuously differentiable across element boundaries (C 1 continuity). Standard finite element programs mostly support C 0 continuity only. Both problems, the non-independence of variations of surface gradients and the C 1 continuity, can be circumvented by relaxing the constraint in the spirit of a penalty approach (see (1) for definition of W n ij ) by adding the term to the lhs of (16). In (19) P is the so-called penalty parameter. If the energy supply is bounded then we expect that the original constraint (18) will be satisfied in the limit as P ! 1. In the relaxed form of the governing equations we have recovered the equations of the full (unconstrained) Cosserat continuum where the rotations u i are independent degrees of freedom. In our case this is true as long as P is finite. The independence of the rotations means that the non-independence of surface gradients and the C 1 continuity problem have disappeared. In finite element calculations the velocities and rotations are approximated independently and since the highest order derivative of both velocities and rotations are of the first order (in the power balance) both may be approximated by using the same type of shape function, which needs to be C 0 continuous only. From the expanded form of the power balance (obtained by adding (19) to the lhs of (16)) and invariance of the power balance with respect to superimposed rigid body rotations we conclude in the usual way that where the anti-symmetric Cosserat stress is defined as As usual in Cosserat theory we may express the rate of deformation measures D ij and W ij À W c ij in terms of a single measure, the relative velocity gradient We note that The relative rotation on the lhs of (21), expressed in terms of c ij reads: where In the derivation of (25) we have used the relations (1). The power balance including the penalty term can be written as where r ij ¼ r ji and r c ij ¼ Àr c ji . The following relationships are used in the finite element implementation: Also in connection with finite element calculations it is customary to arrange the components of the relative velocity gradients in a pseudo vector c T ¼ ðc 11 ; c 22 ; c 33 ; c 12 ; c 13; c 21 ; c 23 ; c 31 ; c 32 Þ. In 2D the matrix representations of A ijkl and D ijkl read: A visco-elastic constitutive relationship inspired by the usual plate theory relating j to the couple stress l and the corotational couple stress rate c (defined in analogy to the corotational stress rate (2)) is derived in the Appendix C. In 2D, for zero relaxation time (purely viscous case) the result reads: where C ½ ¼ n 6 2 þ n 2 1 n 2 2 ð1 þ n 2 2 Þ n 1 n 2 n 1 n 2 n 6 1 þ n 2 1 n 2 2 ð1 þ n 2 1 Þ where t is the thickness of the individual layer. In the following example we assume that t ¼ const:; this can easily be relaxed by assuming _ t t=t ¼ n T ðDnÞ for example. For illustration of the model we consider the simple shearing of an infinite, viscous layer (x 1 ; 0 6 x 2 6 h). On the surface x 2 ¼ h the shear stress r 12 ¼ s ¼ const: is prescribed and we assume that r 22 ¼ 0. Also at x 2 ¼ h we assume that u 3 ¼ 0 and we designate the local, unknown velocity as V. At x 2 ¼ 0 we assume that v 1 ¼ v 2 ¼ u 3 ¼ 0. The initial director orientation is n T ¼ ð1; 0Þ i.e. initially the material layer surfaces are oriented orthogonal to the surfaces of the shear layer. First we consider the convergence of the penalty scheme for increasing values of the penalty parameter P. Fig. 7 shows the dimensionless velocity V ¼ V g S =ðshÞ as a function of the dimensionless parameter P ¼ P =g S for three cases: analytic solution (see Appendix D), finite element solution assuming full numerical integration and a finite element solution using one point integration for the integration of the penalty stiffness terms. The analytic-and the one-point numerical solution coincide to the first three digits, however the finite element solution based on full integration diverges for increasing values of the penalty parameter: the velocity on top tends to zero. The reason for the divergence is that for the present choice of finite elements (four noded quadrilaterals, periodic boundary conditions on the sides) W n is constant within an element whereas W c is linearily variable. The latter produces to a positive definite contribution to the argument of the penalty stiffness depending quadratically on x 2 .
In Figs. 8 and 9 we show the evolution of V as a function of the director orientation as described by the angle U 0 ¼ U 3 (x 2 ¼ 0) between the x 2 -axis and n and the evolution of the dimensionless Cosserat rotation (see Figure caption for definition) at the center of the shear layer respectively. Fig. 7. Dimensionless velocity V ¼ V g S =ðshÞ versus dimensionless penalty parameter P ¼ P =g S analytical solution (crossed); numerical, full integration (broken line); numerical, one point integration (solid line). Finite element model: eight by twelve four noded quadrilaterals; sixteen particles per element. Periodic boundaries on the sides, i.e. velocities and rotations are the same on both sides; if one particle leaves the domain on one side it re-appears on the other side. g=g S ¼ 2, t=h ¼ 0:2. During the calculation the director orientation is fixed at n ¼ ð1; 0Þ, i.e. the internal layering is always orthogonal to the x 2 ¼ const: Fig. 8. V as a function of sinðU 0 Þ, U 0 ¼ U 3 ðx 2 ¼ 0Þ Á t=h ¼ 0:2.

Conclusions
We have presented a simple formulation for the simulation of large, visco-elastic deformations in layered systems. The influence of the bending stiffness of the individual layers is considered within the framework of a couple stress theory. The combination of the basic model with a large deformation, PIC finite element method allows the simulation of a diverse range of crustal deformation problems. By way of examples we have given a realistic treatment of folding and simple shear processes, which includes the mechanical influence of fine-scale. The model is relatively simple in its present form but still gives a useful insight into the physical processes involved in certain types of folding processes involving simple shear.
One of the most interesting results occurs for purely viscous, layered simulations where low-wavenumber folding is induced even for very low viscosity contrasts between embedded and embedding media. In the past, the very large viscosity contrasts required to produce Biot-type folding in purely viscous media have led people to discount the possibility that viscous buckling occurs at all in geology.
In folding couple stresses have an appreciable effect only if the ratios of the parameters (layer material/ embedding material; normal viscosity/shear viscosity etc.) are very extreme, e.g. equal to 1000 or more. In fact in the case of the shear layer the influence of couple stresses appreciable already at relatively low parameter ratios as long as the ratio of the thickness of an individual layer to the width of the shear layer is of the order of 1/10 or more. If couple stresses are considered the deformation of the layer is no longer homogeneous; we observe, depending on the viscosity ratio the formation of a slowly deforming core, rapid shearing concentrates around the layer surfaces. Corresponding results are shown in Figs. 8 and 9.