Large Amplitude Folding in Finely Layered Viscoelastic Rock Structures

Abstract—We analyze folding phenomena in finely layered viscoelastic rock. Fine is meant in the sense that the thickness of each layer is considerably smaller than characteristic structural dimensions. For this purpose we derive constitutive relations and apply a computational simulation scheme (a finite-element based particle advection scheme; see Moresiet al., 2001) suitable for problems involving very large deformations of layered viscous and viscoelastic rocks. An algorithm for the time integration of the governing equations as well as details of the finite-element implementation is also given. We then consider buckling instabilities in a finite, rectangular domain. Embedded within this domain, parallel to the longer dimension we consider a stiff, layered plate. The domain is compressed along the layer axis by prescribing velocities along the sides. First, for the viscous limit we consider the response to a series of harmonic perturbations of the director orientation. The Fourier spectra of the initial folding velocity are compared for different viscosity ratios.¶Turning to the nonlinear regime we analyze viscoelastic folding histories up to 40% shortening. The effect of layering manifests itself in that appreciable buckling instabilities are obtained at much lower viscosity ratios (1:10) as is required for the buckling of isotropic plates (1:500). The wavelength induced by the initial harmonic perturbation of the director orientation seems to be persistent. In the section of the parameter space considered here elasticity seems to delay or inhibit the occurrence of a second, larger wavelength.¶Finally, in a linear instability analysis we undertake a brief excursion into the potential role of couple stresses on the folding process. The linear instability analysis also provides insight into the expected modes of deformation at the onset of instability, and the different regimes of behavior one might expect to observe.


Introduction
A common feature of mechanical systems producing infrequent, catastrophic releases of free energy such as earthquakes or rockbursts is a shared capacity to store energy over prolonged periods. This capacity is absent from models with purely dissipative, viscous rheologies which have conventionally been used to simulate large-deformation geodynamic processes such as mantle convection, large amplitude folding, necking and so on. On the largest scale, catastrophic transitions in plate motion may occur when accumulated gravitational potential energy in the cool thermal boundary layer is released when the yields strength of the lithosphere is exceeded. On a smaller scale, the vigor and the geometric nature of the buckling event depend crucially on the capacity of the folding structure to store and release elastic energy.

Overview
In the following section we develop a mechanical model including a large deformation formulation for viscoelastic, multi-layered rock. The formulation is devised with the goal of describing materials with fine internal layering, which can be described by a single director orientation. This constitutive model is new; 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 the inescapable control of gravity. Indeed, one of the most enduring tenets of geology is the existence of an organized stratigraphy in the rock record. Another generalization we might make about the mechanical origins of geological formations is that deformation is almost certain to involve very high strain. For example, Figure 1 shows the largeamplitude folding of multi-layered rock on the scale of a few tens of centimeters. The macroscopic deformation is assumed volume preserving for convenience.
In the following two sections we turn to a computational method capable of following the evolution of elastic stresses, macroscopic interfaces and the internal layering direction introduced in the constitutive relationship of section 2. The method needs to be able to deal with very large strains associated with largeamplitude folding, while faithfully tracking material history and interfaces. Versatility and robustness are usually associated with the various formulations of the Finite Element Method (FEM). The need to track material history strongly suggests a Lagrangian formulation, which provides a reference frame, locked with the material itself. Unfortunately, large deformations are quite difficult to handle within the FEM because mesh distortion and remeshing are required to maintain optimal element configurations.
In section 3 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 and elasticity. In section 4 we describe an extension to standard finite element methods, which incorporates moving integration points to carry director orientation and other history variables.
The Particle-In-Cell (PIC) finite element method, as this technique is known, is a hybrid scheme which falls somewhere between the Finite Element Method (FEM) and a purely Lagrangian particle method such as the Discrete Element Method (DEM). The PIC scheme attempts to combine the versatility of the continuum FEM with the geometrical flexibility of DEM.
In Section 5, we explore scenarios from global to internal buckling and from elastic to purely viscous folding in nonlinear finite element studies. It turns out that the characteristic length scale of the emerging folding pattern tends to zero with increasing relaxation time. An explanation for this is attempted in section 6 within the framework of a linear stability analysis.

Mathematical Formulation
Layered materials are ubiquitous in geological formations. Accordingly, virtually every mathematically-minded structural geologist has, at some stage, contributed to this field (CHAPPLE, 1969;FLETCHER, 1974FLETCHER, , 1982JOHNSON and FLETCHER, 1994;SCHMALHOLZ and PODLADCHIKOV, 1999;HUNT et al., 1997;VASILYEV et al., 1998). Layering may be caused by purely mechanical, hydro-mechanical or chemomechanical means (e.g., WILLIAMS, 1972;ROBIN, 1979;ORTOLEVA, 1994). From a mechanical point of view, the salient feature of such 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 behavior 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. 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 2-D matrix representation of (1) as needed for our computational applications is represented in Appendix A.2 for easy reference. We define a corotational stress rate as: Again the superscripted n distinguishes the stress rate _ r r n as observed by an material observer corotating with the director n from the material stress rate _ r r observed by a spatially fixed observer.

Specific Viscous and Viscoelastic Constitutive Relations
We consider layered, viscous and viscoelastic 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.1 for derivation) to consider the mechanical effect of the layering; thus where a prime designates the deviator of the respective quantity, and Similarly, a viscoelastic 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. The present choice seems more natural in the context of layered continua and also leads to algebraically simpler expressions in the linear stability analysis in section 6. 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 as follows: K is the so-called global stiffness matrix which contains all the material property 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 surface tractions. The matrices K and G are global matrices composed in the usual way of elemental matrices; in the following we designate element matrices 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 an. 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 viscoelastic 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 equations (10) and (11)). The result reads: Details of the derivation of (9) are represented in the Appendix 3. The effective viscosities read: where the superscripted t refers to the previous time step and Dt is the time increment. The time increment is limited by a 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 time step). In MORESI et al. 2001 the time increment for the advection of the material points usually differs from the time increment used for the integration of the stress history. This is necessary to ensure time scale of interest for elastic processes is determined by the physics of the problem and not by the mesh dependent Courant conditon. The parameters a ¼ g=l and a S ¼ g S =l S are the relaxation times for pure and simple shear, respectively. The superscripted parameter 0 b 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 equations (1, 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. In this case, for t fixed, we replace s t by s tþDt , which during iteration is defined by s tþDt r tþDt ij as defined by (9). 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.

Possible Numerical Schemes
In fluid dynamics, where strains are generally very large, but do not appear in the constitutive relationship of the material, it is common to transform the equations to an Eulerian mesh and deal with convective terms explicitly. Problems arise whenever advection becomes strongly dominant over diffusion since an erroneous numerical diffusion dominates. In our case, the advection of material boundaries and the stress tensor are particularly susceptible to this numerical diffusion problem. Mesh-based Lagrangian formulations alleviate this difficulty, but at the expense of remeshing and the eventual development of a less-than optimal mesh configuration. This increases complexity and can hinder highly efficient solution methods such as multigrid iteration. The Natural Element Method eliminates remeshing difficulties but is associated with considerable complexity of implementation, particularly in 3-D. A number of alternatives are available which dispense with a mesh entirely: smooth particle hydro dynamics and discrete element methods are common examples from the fluid and solid mechanics fields, respectively. These methods are extremely good at simulating the detailed behavior of highly deforming materials with complicated geometries (e.g., free surfaces, fracture development), and highly dynamic systems. They are, in general, formulated to calculate explicitly the interactions between individual particles which ultimately means that a great many time steps would be required to study creeping flow where the time scales associated with inertial effects are very many orders of magnitude smaller than typical flow times. We have therefore developed a hybrid approach -a particle-in-cell 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.

The Particle-in-Cell Approach
Our particle-in-cell 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 subdomains (elements) and are calculated by summation over a finite number of sample points within each element. For example, in order to integrate equation (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 and corotational stress rate, which are therefore directly available for the element integrals without the need to interpolate from nodal points to fixed integration points. We therefore store an initial set of w p 's based on a measure of local volume and adjust the weights slightly to improve the integration scheme. MORESI et al. (2000) give a full discussion of the implementation of the particle-in-cell finite element scheme used here including full details of the integration scheme and its assumptions. They also discuss the specific modifications to the material point method required to handle a convecting fluid.

Figure 2
Schematic of Particle-in-Cell Method for representing large deformation in materials with interfaces and material history (including storing/transport of tensorial information such as stress). Mesh points remain fixed; particles move relative to the mesh and carry interface information via their relative positions and directional information directly on the particles.

Numerical Simulations
We present an example of a simulation of folding of a layer of anisotropic viscoelastic material sandwiched between two isotropic viscous layers of equal viscosity. To accommodate the shortening of the system, one of the isotropic layers is compressible. In benchmarking of viscous folding, this sandwich of incompressible and compressible embedding material was found to give good agreement with analytic results assuming an infinite domain (MORESI et al., 2000). Throughout, we deal with a special case where the relaxation time for the normal and shear components of the rheology are always identical (i.e., a ¼ g=l ¼ a s ¼ g s =l s ).
We first examine the viscous limit (a ! 0, a S ! 0) for infinitesimal deformation of an embedded layer with a normal viscosity contrast of g=g M ¼ 1, 10 or 100 to the embedding medium. The boundary conditions in this case are slippery, undeforming boundaries on the vertical sides and the base. No density variations are assumed. For reference, in an isotropic sample with a viscosity contrast of 1 or 10 between the layer and the embedding material no appreciable folding occurs during shortening by 50%, as predicted by linear instability analysis (e.g., BIOT (1965a); see also section 6 in which we explore the potential influence of couple stresses within the framework of a linear instability analysis).
The initial orientation of the internal layering is approximately parallel to the macroscopic layering of the system with a small harmonic perturbation, introduced particle-by-particle: Figure 4 plots the dominant components of the Fourier transform of the vertical velocity along the mid-line of the embedded layer in each of three cases at a time t ¼ 0. Case 1 has no normal viscosity contrast to the embedding medium (g ¼ g M ), and a shear-to-normal viscosity contrast of 100 (g ¼ 100g s ). The growth rate is strongly peaked for a perturbation wavenumber q ¼ p. Case 2 has a normal viscosity contrast of 10 between the embedding and embedded materials (g ¼ 10g M ) and g ¼ 10g s . There are two strong signals in the Fourier transform of the vertical velocity: one at the wavenumber of the director perturbation and another, stronger signal at q ¼ p=2. The low wavenumber signal has a growth rate which is almost independent of the perturbation wavenumber. Growth at the wavenumber of the director perturbation falls off with increasing q, but reaches a plateau by around q ¼ 8p. At a higher contrast between embedded and embedding materials, Case 3, g ¼ 100g M , g ¼ 10g s shows a similar pattern to case 2, except that the growth rates are amplified. The initial growth rates described in Figure 4 correspond to infinitesimal deformation of the layer. At finite shortening, the initial deformation pattern may be modified. For each of the three cases above, we plot the finite amplitude response to shortening in Figure 5.
For case 1 (Fig. 5a), the small growth rates observed at the outset produce only very small deflection of the layer interfaces after 40% shortening. Of interest, however, is the fact that a perturbation with wavenumber q ¼ 10p produced a corresponding layer deflection of larger magnitude than a perturbation at q ¼ 2p, suggesting that the growth rate remains flat as a function of q even at finite amplitude deformation.

Figure 3
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: viscoelastic ( g; g s ; l; l 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.

Figure 4
Rate of growth at wavenumbers introduced through perturbation to the director orientation expressed as Fourier coefficients of vertical velocity at the mid-line at time t ¼ 0 for purely viscous cases (g=l ! 0). Isotropic embedding material has viscosity 1. Case 1, no contrast in normal viscosity (g=g s ¼ 100; g=g M ¼ 1). Case 2, g=g s ¼ 10; g=g M ¼ 10 Case 3, g=g s ¼ 10; g=g M ¼ 100.
For case 2 (Fig. 5b), 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.
For case 3 (Fig. 5c), we observe the same overall trend as case 2: 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 secondary variation in the interface deflection.
Shortening by 40% is examined in Figure 6 for a viscoelastic layer. This case corresponds to case 3 of the purely viscous simulations, but with finite relaxation time: g ¼ 100g M , g ¼ 10g s , l ¼ 100-100000. The introduction of elasticity strongly reduces the tendency to generate long-wavelength buckling modes in the layer. For example, when l ¼ 1000, the perturbation wavenumber is dominant even for q ¼ 10p, and the amplitude of the deformation is considerably larger for small wavenumber than the corresponding deformation for the viscous case.
The linear instability analysis predicts a strong tendency to generate very high wavenumber (short wavelength) folds. The numerical simulations uphold this prediction, amplifying the deformation at the finest available wavelength: the one provided by the finite element mesh. The anisotropic layer itself shortens almost uniformly, but the internal layering direction develops an extremely strong periodicity in the shortening direction. Couple stresses (neglected in our numerical analysis) would stabilize the solution at a finite wavenumber and, thus, at least ameliorate the mesh-sensitivity of the result (MÜ HLHAUS, 1993). In the viscous simulations, the natural development of low-wavenumber buckling of the anisotropic layer from initially fine-scale perturbations suggests that the role of couple stresses is considerably less important. We conclude the main body of this paper with a brief excursion into the potential role of couple stresses within the framework of a linear instability analysis.

Couple Stresses and Linear Instability Analysis
Here we restrict ourselves to an incrementally elastic constitutive relation for the couple stress, partially for algebraic convenience but also because the numerical difficulties mentioned above seem to occur only in connection with significant elastic deformations. Couple stresses may be significant in situations where the gradient of n i changes strongly over a short distance (limiting case: disinclination).
In such cases we have to take the variations of the normal stresses across the layer cross sections into consideration (e.g., MÜ HLHAUS, 1993). The couple stress theories (see e.g., MINDLIN and THIERSTEN, 1962;MÜ HLHAUS, 1993;MÜ HLHAUS and Figure 5a Evolution of folding in anisotropic viscous layer. Isotropic embedding material has viscosity 1, layer has shear viscosity 0.01, normal viscosity 1. Results are shown for perturbation to the director orientation with wavenumber q ¼ 2p and q ¼ 10p.

Figure 5b
Evolution of folding in anisotropic viscous layer. Isotropic embedding material has viscosity 1, layer has shear viscosity 1, normal viscosity 10. Results are shown for perturbation to the director orientation with wavenumber q ¼ 2p and q ¼ 10p.
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. In the present case the couple stress enhancement leads naturally to the superposition of viscoelastic bending stiffness on our standard continuum model. 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., MU¨HLHAUS, 1993). In layered materials the explanation why the stress tensor is nonsymmetric in couple stress materials is straightforward: 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 is even zero for instance in the case of a stack of perfectly smooth 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. As usual in plate theories the couple stresses (moment per unit area) are conjugate in energy to a suitable measure for the rate of curvature. We assume that in the reference configuration the layer normal is parallel to the Cartesian x 3 axis. Without going into details we mention the following definition of the curvature rate, which is suitable in the context of layered materials which upon linearization, and nothing else is required for our linear instability analysis, reduces to the familiar expression All other components of j vanish in the linear case. In (14) symm(.) means the symmetric part of the argument and k ¼ nn T (see equation (1)). We assume the simplest possible constitutive relationship for the couple stress by putting _ m m n 11 ¼ _ m m 11 ¼ Hlh 2 j 11 . The corotational rate is equal to the material rate because m 11 ¼ 0 in the ground state; H is a dimensionless scalar typically smaller than unity and h is the width of a periodic cell of the layered material. In the linear instability analysis (see Appendix 4 for details) we assume one-dimensional modes of Figure 5c Evolution of folding in anisotropic viscous layer. Isotropic embedding material has viscosity 1, layer has shear viscosity 10, normal viscosity 100. Results are shown for perturbation to the director orientation with wavenumber q ¼ 2p and q ¼ 10p.
b Figure 6 Layered viscoelastic case. Isotropic embedding material has viscosity 1, layer has shear viscosity 10, normal viscosity 100. Results are shown for a range, l and for perturbation to the director orientation with wavenumber q ¼ 4p and q ¼ 10p. The ratio of elastic moduli is constant throughout: l=l s ¼ 10.
the form reflecting the fact that typical folding scenarios are usually one-dimensional (Fig. 2), at the onset at least, depending only on the coordinate parallel to the layer surfaces, x 1 (e.g., BIOT, 1964BIOT, , 1965aBIOT, ,b, 1967HILL and HUTCHINSON, 1975). The onedimensionality itself does not represent an assumption per se; it does however in the context of our specific problem: Here and in the computational examples we consider a multilayered structure embedded in an isotropic, incompressible viscous medium of infinite extent. Again, within the context of this specific problem, the assumption of one-dimensionality expressed in (16) precludes gross changes in layer thickness during folding as is observed in many natural folds. Next in our stability analysis we insert (16) into the constitutive relations and the result into the incremental equilibrium conditions (see Appendix 4). Because of the one-dimensionality of our problem equilibrium needs to be considered in the x 2 direction only. The resistance of the embedding medium against the folding (see Fig. 7) of the structure is considered BIOT-style (1965a) by the traction 2Q ¼ À4g M xq expðxtÞU cosðqx 1 Þ. Again, details of this derivation are included in the Appendix.
Insertion of (A.18) into (A.23) and nondimensionalising leads to the following characteristic equation for the dimensionless growth coefficient of the folding instability: Figure 7 Moment equilibrium diagram for a slice of the layered structure. The macroscopic bending moment M (a); couple stress m 11 (b); the traction 2Q represent the resistance of the embedding medium against the deformation of the layered structure. x where d is the thickness of the folding plate (Fig. 7) and For h=d ! 0 we recover the standard continuum case (equation (A.19)). However, unlike in the case of the standard continuum ( Fig. 9) the positive branch of (17) tends to zero as q ! 1 i.e., the ill-posed nature of the folding problem is removed.
A maximum of xðqÞ exists as long as h > 0. For illustration we consider the extreme case g S ¼ 0. From (17) we obtain: Unstable modes are obtained for q ðr=ðHlh 2 ÞÞ 1=2 and the maximum growth rate x max ðq max Þ is obtained as: The wavelength predicted by (20) for the extreme case g S ¼ 0 ranges between 15 to 20 times the thickness h (which is well within the ranges of wavelength observed in real rocks) if the stress ratio r=l ranges between 1/10 and 1/20. The maximum shifts to the longer wavelength side with increasing magnitude of the coefficient ofq q 3 in (20). If h ¼ 0 the maximum degenerates to a boundary maximum at q ! 1. In this case the perturbations with the shortest wavelength grow fastest, are dominant, and we have to expect a strong dependency of finite element solutions on the mesh size. A similar situation occurs in connection with strain localisation and strain softening (DE BORST et al., 1993). However, in many cases the maximum occurs at wavenumbers, which are beyond the applicability of the concept of equivalent (or homogenized) continua and prior to the maximum, the situation is usually satisfactorily described as a standard continuum.

Conclusions and Future Research
We have presented a simple formulation for the consideration of viscoelasticity in deforming layered systems. The combination of the basic model with a large deformation, particle-in-cell finite element method allows the simulation of a diverse range of crustal deformation problems.
Our demonstration examples include a realistic treatment of folding which includes the mechanical influence of fine-scale laminations and viscoelasticity. 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.
While we have noted the importance of considering couple stresses for the elasticdominated simulations, the story will not be complete unless we also consider that Dispersion function for viscoelastic, layered material with bending resistance, g M = viscosity of embedding material, g S and l S are the shear viscosity and the weak shear modulus, t and h are the width of the structure and the individual layer, respectively, L = wavelength, x = growth coefficient. The growth coefficient x is bounded if r=l S < 1 (see Appendix 3).
the large couple stresses produced by layer bending could also lead to failure of the rock. The consideration of layered materials with couple stresses and a yield criterion in large deformation is the subject of ongoing research, however, we would predict that yielding within the layering will produce a localized band structure with a strong discontinuity in the director orientation. 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.
Viscoelastic layered buckling tends to emphasize the finest scale imposed by explicitly assumed perturbations as well as discretization artifacts. The latter implies mesh-dependency not at all uncommon in numerical simulation but unwelcome all the same. The possible cure would be to extend the numerical representation to include terms which, in nature, become dominant at large wavenumber (in this case consideration of the bending stiffness of the individual layers).
The particle-in-cell formulation for including couple stresses is conceptually no different from the one presented here though more complicated in detail.
Linear instability analysis gives a good insight into the expected modes of deformation at the onset of instability, and the different regimes of behavior one might expect to observe.