# The interplay of material and geometric instabilities in large deformation s of viscous rock

H. B. Muhlhaus, L. Moresi, Frédéric Dufour, B. Hobbs

► To cite this version:

H. B. Muhlhaus, L. Moresi, Frédéric Dufour, B. Hobbs. The interplay of material and geometric instabilities in large deformation s of viscous rock. IUTAM Symposium on Analytical and Computational Fracture Mechanics of Non-Homogeneous Materials, 2001, Cardiff, United Kingdom. pp.105-115, 10.1007/978-94-017-0081-8_14 . hal-01009021

HAL Id: hal-01009021

https://hal.science/hal-01009021

Submitted on 21 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
THE INTERPLAY OF MATERIAL AND GEOMETRIC
INSTABILITIES IN LARGE DEFORMATIONS OF VISCOUS
ROCK

H-B MÜHLHAUS, L.MORESI, F. DUFOUR, B. HOBBS

CSIRO Division of Exploration and Mining
39 Fairway, Nedlands WA 6009, Australia

hans@ned.dem.csiro.au

**Abstract**

We analyse folding phenomena in layered viscous rock in which we assume that the thickness of each layer is much smaller than characteristic structural dimensions. We derive constitutive relations and apply a computational simulation scheme suitable for problems involving very large deformations of layered viscous materials. 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 effect of layering manifests itself in that appreciable buckling instabilities are obtained at much lower viscosity ratios (1:10) than required for the buckling of isotropic beams (roughly 1:500). The wavelength induced by the initial harmonic perturbation of the layer orientation seems to be persistent. A linear instability analysis is conducted to clarify the numerical results and also to examine the potential role of couple stresses on the folding process

**1 Introduction**

The founding fathers of modern geology established that most sedimentary rocks were originally deposited in sequences of soft horizontal layers at the bottom of shallow seas and hardened over time. The challenge remains in many instances to explain the evolution of the enormous variety of shapes and forms observed in 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⁻¹⁴ s⁻¹), viscous, elastic, and brittle effects influence the observed structures. The example in figure 1 illustrates another reason for the popularity of folding as a subject of geological interest: the frequent association of fluid flow and mineralization with tightly folded rocks.

We develop a mechanical model including a large deformation formulation for multi-
layered rock. The formulation describing materials which also have a 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.

<figure><img src="image_2.png" /><figcaption>Figure 1. A plaque in Castlemaine commemorating gold discoveries in the Bendigo Gold Fields: "ANTICLINAL FOLD: This fine exhibit was disclosed when Lyttleton Street East was constructed in 1874. Saddle reefs occur in similar folds of the sandstones and slates in lower geological horizons".</figcaption></figure>

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 nonlinear finite element studies. These show that buckling can be induced at much lower strength contrasts between the matrix and the embedded beam than would be the case for isotropic materials. An explanation for this is attempted within the framework of linear stability analysis which also provides insight into the expected modes of deformation at the onset of instability, and the different regimes of behaviour one might expect to observe.

## 2. 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, 1974, 1986; Johnson & Fletcher, 1994; Schmalholz & Podladchikov, 1999; Hunt et al., 1997; Vasilyev et al., 1998). Layering may be caused by purely mechanical, hydro-mechanical or chemo-mechanical means (eg Williams, 1972; Robin, 1979; Ortoleva, 1994). From a mechanical point of view, the salient feature of such materials is a distinguishing 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. We assume linear viscous behaviour and designate with $\eta$ the normal viscosity and $\eta_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

$$ \dot{\mathbf{H}}_i = W_{ij}^n n_j \text{ where } W_{ij}^n = W_{ij} - (D_{ki}\lambda_{kj} - D_{kj}\lambda_{ki}) \text{ and } \lambda_{ij} = n_i n_j \quad (1) $$

where **L**=**D**+**W** is the velocity gradient, **D** is the stretching and **W** is the spin. The superscribed *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 material 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 μ and η respectively; the shear modulus and the shear viscosity measured in simple, layer parallel shear we designate as μ_S and η_S.

In the following simple model for a layered viscous material we correct the isotropic part $2ηD'_{ij}$ of the model by means of the Λ tensor (see appendix Mühlhaus et al, 2001, for the derivation) to consider the mechanical effect of the layering; thus

$$ \sigma_{ij} = 2\eta D'_{ij} - 2(\eta - \eta_S)\Lambda_{ijlm} D'_{lm} - p\delta_{ij} \quad (3) $$

where a prime designates the deviator of the respective quantity, and

$$ \Lambda_{ijkl} = \left( \frac{1}{2} (n_i n_k \delta_{lj} + n_j n_k \delta_{il} + n_i n_l \delta_{kj} + n_j n_l \delta_{ik}) - 2n_i n_j n_k n_l \right). \quad (4) $$

### 3. 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:

$$ \begin{pmatrix} \mathbf{K} & \mathbf{G} \\ \mathbf{G}^{\mathrm{T}} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \mathbf{u} \\ \mathbf{p} \end{pmatrix} = \begin{pmatrix} \mathbf{F} \\ \mathbf{0} \end{pmatrix} \quad (6) $$

$\mathbf{K}$ is the so-called global stiffness matrix which contains all the material property parameters, $\mathbf{G}$ is the divergence expressed in matrix form, $\mathbf{u}$ and $\mathbf{p}$ are the unknown velocities and pressure respectively, and $\mathbf{F}$ is a vector of driving terms comprising body forces and surface tractions (e.g. see Hughes, 1987).

The matrixes $\mathbf{K}$ and $\mathbf{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
$$ \mathbf{K}^E = \int_{\Omega^E} \mathbf{B}^T(\mathbf{x})\mathbf{C}(\mathbf{x})\mathbf{B}(\mathbf{x})d\Omega \qquad (7) $$

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,

$$ C_{iso} = 2\eta \begin{bmatrix} 1 \\ \frac{1}{2} \end{bmatrix}; \quad C_{layer} = -2\eta \left(\frac{\eta_s}{\eta} - 1\right) \begin{bmatrix} -\Delta_0 & \Delta_0 & -\Delta_1 \\ \Delta_0 & -\Delta_0 & \Delta_1 \\ -\Delta_1 & \Delta_1 & -\frac{1}{2} + \Delta_0 \end{bmatrix} \qquad (8) $$

in which $\Delta_0 = 2n_1^2 n_2^2$ and $\Delta_1 = (n_1 n_2^3 - n_1^3 n_2)$.

## 4. 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 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.

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 sub-domains (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 $\Omega^E$ we replace the continuous integral by a summation

$$ \mathbf{K}^E = \sum_p w_p \mathbf{B}^T(\mathbf{x}_p) \mathbf{C}_p(\mathbf{x}_p) \mathbf{B}(\mathbf{x}_p) \qquad (12) $$

In standard finite elements, the positions of the sample points, $\mathbf{x}_p$, and the weighting, $w_p$, are optimized in advance. In our scheme, the $\mathbf{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 particle-in-cell finite element scheme used here including full details of the integration scheme and its assumptions.

## 5. Numerical simulations

<figure><img src="image_3.png" /><figcaption>Figure 2. Initial geometry for the folding experiment. Layer 1 is compressible, viscous ($\eta_M$) background material, layer 2 is identical to layer 1 but incompressible (see text for an explanation), layer 3 is the test sample: viscoelastic ($\eta, \eta_s$) with a director orientation ($\mathbf{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.</figcaption></figure>

We present an example of a simulation of folding of a layer of anisotropic viscous material sandwiched between two isotropic layers of equal viscosity (Figure 2). To accommodate the shortening of the system, one of the isotropic layers is compressible. In benchmarking this sandwich of incompressible and compressible embedding material was found to give good agreement with analytic results assuming an infinite domain (Moresi et al, 2001). We present three cases which have normal viscosity contrasts between the embedding material and the beam ($\eta_M/\eta$) of 1 (case 1), 1/10 (case 2) and 1/100 (case 3). For case 1, the ratio $\eta/\eta_s = 100$, whereas for cases 2 and 3 $\eta/\eta_s = 10$.
<figure><img src="image_4.png" /><figcaption>Figure 3. Case 1: 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 = 2\pi$ and $q = 10\pi$.</figcaption></figure>

For case 1 (Figure 3), 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 = 10\pi$ produced a corresponding layer deflection of larger magnitude than a perturbation at $q = 2\pi$, suggesting that the growth rate remains flat as a function of $q$ even at finite amplitude deformation.

For case 2 (Figure 4), the growth rate is higher for smaller wavenumber in the infinitesimal deformation limit, and this persists with finite amplitude deformation. Growth at wavenumber $q = 2\pi$ is significantly more developed after 40% shortening than growth at $q = 10\pi$.

For case 3 (Figure 5), 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.

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Case 3: 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 = 2\pi$ and $q = 10\pi$.</figcaption></figure>

# 6. Linear instability analysis

The numerical simulations show a propensity to, amplify 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 (Muhlhaus, 1993). We conclude the main body of this paper with a linear instability analysis in which we consider the potential effect of couple stresses as well. We present only the key assumptions and results of the analysis; for details of the derivation we refer to (Muhlhaus et al 2001).

We consider a plate of thickness $h$ consisting of a layered viscous material. The layering is assumed to be parallel to the plate surfaces. The plate is embedded in an infinite, isotropic material. The viscosity of the embedding material is $\eta_M$. The thickness of the individual layers is $t$. The resistance of the embedding medium against the folding of the structure is considered Biot-style (1965a) by the traction $2\sigma_n = -4\eta_M \omega q \exp(\omega t) U \cos(qx_1)$. 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. Muhlhaus, 1993). The couple stress theories (see e.g. Mindlin and Thiersten, 1962; Muhlhaus, 1993; Mühlhaus & 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 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. Mühlhaus, 1993). In layered materials the explanation for non-symmetry of the stress tensor 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 may be zero 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 (cf. standard engineering beam and plate theories), and introduces statically equivalent couple stresses to balance 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, in our 2D linear instability analysis, the only non-zero component of the curvature rate is: $\kappa_{11} = v_{2,11}$ (Muhlhaus et al 2001). The simplest possible constitutive relationship for the couple stress is obtained by setting $\dot{m}_{11} = 1/3(\mu - \mu_S)^2 \kappa_{11}$. Note that the corotational rate is equal to the material rate because $m_{11}=0$ in the ground state. Assuming modes of the type

$$ \delta u_2 = U \exp(\omega t) \cos \frac{2\pi}{L} x_1 = U \exp(\omega t) \cos qx_1, \quad \delta u_1 = 0, \qquad (13) $$

the critical condition is obtained as:

$$ \omega = \frac{\sigma k + \sigma \frac{\eta - \eta_S}{3\eta_S} k^3}{\frac{\eta k^2}{3\eta_S} \left( \eta_S k + \frac{(\eta - \eta_S)}{3} \left(\frac{t}{h}\right)^2 k^3 + 4\eta_M \right) + \frac{(\eta - \eta_S)}{3} \left(\frac{t}{h}\right)^2 k^3 + 4\eta_M} \qquad (14) $$

where $k = 2\pi h / L = qh$, $L$ is the wavelength of the perturbation (13) and $\omega$ is its growth coefficient (dimension 1/time). In the limit $\eta_S \to 0$ we obtain

$$ \omega = \frac{\sigma k}{4\eta_M + \frac{1}{3}\eta\left(\frac{t}{h}\right)^2 k^3} \qquad (15) $$

Next we consider the limit $\eta_S = \eta$. In this case no micro lamination exists. The structure behaves like a monolithic beam of thickness h.:

$$ \omega = \frac{\sigma k}{4\eta_M \left(1 + \frac{1}{3}k^2\right) + \frac{\eta}{3}k^3} \qquad (16) $$
The key results can be summarized as follows (Figures 6 and 7): For $t/h=0$ we obtain the standard continuum case (Figures 6a and 7a). Structural buckling does not occur if the ratio of $\eta/\eta_M$ is low (Figure 6a). In all standard continuum cases the short wavelength limit (obtained as $k \to \infty$) of the growth coefficient $\omega$ is obtained as $\omega = \sigma/\eta_s$. If structural buckling is suppressed due to a low ratio of $\eta/\eta_M$ (see Figure 6a) the shortest wavelength grow fastest ($\omega = \sigma/\eta_s$), 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). The same holds true if structural buckling occurs (Figure 7a) if the ratio of $\eta/\eta_s$ sufficiently high. For moderate ratios of $\eta/\eta_s$, a dominant wavelength associated with a maximum of the graph $\omega(k)$ still exists. In contrast to a monolithic plate we have $\omega(k \to \infty) = \sigma/\eta_s \neq 0$ in general. If couple stresses are considered ($t/h>0$, Figures 6b and 7b) the graph $\omega(k)$ always has a maximum, and $\omega(k \to \infty) \to 0$ i.e. the mesh sensitivity of finite element solutions for decreasing element size should, at least in principle, not occur. However the maximum may occur at such small wavelength that they are beyond the resolution of practical finite element models i.e. the effect of micro-couples will not be "felt" by the model.

<figure><img src="image_7.png" /><figcaption>Figure 6. Disperion relationship: (a) $\eta/\eta_M = 1, t/h=0$ (no couple stresses); (b) $\eta/\eta_M = 1, t/h=1/100$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Dispersion relationship: (a) $\eta/\eta_M = 1000, t/h=0$ (no couple stresses); (b) $\eta/\eta_M = 1000, t/h=1/100$</figcaption></figure>
## 7. Conclusions

We have presented a simple formulation for the consideration of viscous deformation in 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. By way of examples we have given a realistic treatment of folding 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.

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.

Linear instability analysis gives a good insight into the expected modes of deformation at the onset of instability, and the different regimes of behaviour one might expect to observe. Consideration of the finite thickness of the individual layers within the framework of a couple stress theory sometimes does not suffice to remove the mesh sensitivity of the finite element solution. A possible remedy may be to introduce a suitable scalar measure for the director fluctuations as a dynamic quantity in order to separate the effect of small wavelength contributions from the mean field response. Such a strategy would not be unlike the treatment of the kinetic energy of the velocity fluctuations in the turbulence theories of fluid dynamics.

## 8. References

* Biot, M.A. 1964. Theory of Internal Buckling of a confined multi layered Structure. Geol. Soc. America Bulletin, Vol 75. 563-568

* Biot, M.A. 1965a. The Mechanics of Incremental Deformations. John Wiley, New York

* Biot, M.A. 1965b. Theory of similar folding of the first and second kind. Geol. Soc. America Bull. Vol.76. 251-258

* Biot, M.A. 1967. Rheological stability with couple stresses and its application to geological folding. Proc. Royal Soc. London. v. A 2298. 402-2423

* de Borst, R., Sluys, L.J., Mühlhaus, H.-B. and Pamin, J. 1993. Fundamental issues finite element analyses of localization of deformation. Engineering Computations 10, 99-121.

* Chapple, W.N. 1969. Fold Shape and Rheology: The Folding of an Isolated Viscous-Plastic Layer. Tectonophysics. Vol 7. 97-116

* Fletcher, R.C. 1974. Wavelength Selection in the Folding of a Single Layer with Power Law Rheology. Am. J. Sci. Vol 274. 1029-1043

* Fletcher, R.C. 1982. Coupling of diffusional mass transport and deformation in tight rock. Tectonophysics, Vol. 83. 275-291

* Green, A.E. and Zerna, W. 1968. Theoretical elasticity. 2nd Edition. Oxford at the Clarendon Press.

* Hill, R. and Hutchinson, J.W. 1975. Bifurcation Phenomena in Plane Tension Test. J. Mech. Phys. Solids. Vol. 23. 239-264
Hobbs, B.E., Muhlhaus, H-B, Ord, A and Moresi, L. (2001) The Influence of Chemical migration upon Fold Evolution in Multi-layered Materials. Vol. 11, Yearbook of Self Organisation. Eds H.J. Krug and J.H. Kruhl.

Hobbs, B.E., Means, W.D. and Williams, P.F. 1976. An Outline of Structural Geology. Wiley

Hunt, G.W., Muhlhaus, H-B and Whiting, A.I.M. (1997). Folding processes and solitary waves in structural geology. Phil. Trans. of R. Soc. London, 355, 2197-2213.

Johnson, A.M. and Fletcher, R.C. (1994) Folding of Viscous Layers. Columbia University Press, New York

Mindlin, R.D. and Thiersten, H.F. 1962. Effects of couple stresses in linear elasticity. Arch. Rat. Mech. Analysis. Vol. 11. 415-448

Moresi, L., Muhlhaus, H.-B., Dufour, F. 2001. Particle- in-Cell Solutions for Creeping Viscous Flows with Internal Interfaces . In: Proceedings of the 5th International Workshop on Bifurcation and Localization in Geomechanics (IWBL’99), Perth, W. A., Australia. Ed’s: H-B Muhlhaus, A Dyskin and E Pasternak Australia, Balkema, Rotterdam.To appear: Dec2000

Mühlhaus, H.-B., Aifantis, E.C. 1991a. A variational principle for gradient plasticity. Int. J. Solids and Structures. Vol. 28. 845-857

Mühlhaus, H.-B., Aifantis, E.C. 1991b. The influence of microstructure-induced gradients on the localization of deformation in viscoplastic materials. Acta Mechanica Vol. 89. 217-231.

Mühlhaus, H.-B. 1993. Continuum models for layered and blocky rock. In: Comprehensive Rock Engineering. Invited Chapter for Vol. II: Analysis and Design Methods. Pergamon Press. 209-230.

Muhlhaus, H-B, Moresi, L, Hobbs, B. E., Dufour, 2001, Large amplitude folding in finely layered viscoelastic rock structures, Pure Applied Geophys, In press. http://www.ned.dem.csiro.au/research/solidmech/Publications

Ortoleva, P. J. 1994. Geochemical Self-Organization. Oxford University Press.

Robin, P.-Y.F. 1979. The Theory of Metamorphic Segregation and Related Processes. Geochim. Cosmochim. Vol 43. 1587-1600

Schmalholz, S.M. and Podladchikov, Y. 1999. Buckling versus folding: Importance of Viscoelasticity. Geophysical Research Letters, Vol. 26, 17, 2641-2644

Sulsky, D., Zhou, S.-J., Schreyer, H. L. 1995. Application of a Particle-in-Cell Method to Solid Mechanics, Comput.Phys. Commun. 87, 236-252.

Vasilyev, O.V., Podladchikov, Y.Y. and Yuen, D. A. 1998. Modelling of compaction driven flow in poro-viscoelastic medium using adaptive wavelet collocation method, Geophysical Research Letter, 17, 32-39

Williams, P., F. 1972. Development of metamorphic layering and cleavage in low-grade metamorphic rocks at Barmagui, Australia. Am. J. Sci., 272. 1-47