Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation

Tectonophysics 484 (2010) 36–47

Contents lists available at ScienceDirect

Tectonophysics j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / t e c t o

Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation Boris J.P. Kaus ⁎ Geophysical Fluid Dynamics, Department of Earth Sciences, ETH Zurich, Switzerland University of Southern California, Los Angeles, USA

a r t i c l e

i n f o

Article history: Received 18 February 2009 Received in revised form 24 June 2009 Accepted 31 August 2009 Available online 12 September 2009 Keywords: Lithospheric deformation Rheology Plastic deformation Viscoelastoplasticity Viscoplasticity Numerical modelling Long-term tectonics

a b s t r a c t Numerical models of brittle deformation on geological timescales typically use a pressure-dependent (Mohr– Coulomb or Drucker–Prager) plastic flow law to simulate plastic failure. Despite its widespread usage in geodynamic models of lithospheric deformation, however, certain aspects of such plasticity models remain poorly understood. One of the most prominent questions in this respect is: what are the factors that control the angle of the resulting shear bands? Recent theoretical work suggest that both Roscoe (45°), Coulomb angles (45 +/− ϕ/2, where ϕ is the angle of internal friction) and Arthur angles (45° +/− (ϕ + Ψ/4) where Ψ is the dilation angle), as well as all intermediate angles are possible. Published numerical models, however, show a large range of shear band angles with some codes favoring Arthur angles, whereas others yield Coulomb angles. In order to understand what causes the differences between the various numerical models, here I perform systematic numerical simulations of shear localization around an inclusion of given length scale. Both numerical (element type), geometrical and rheological (viscoplastic versus viscoelastoplastic) effects are studied. Results indicate that the main factor, controlling shear band angle, is the non-dimensional ratio between the length scale of the heterogeneity d and the size of the numerical mesh Δx. Coulomb angles are observed only in cases where the inclusion is resolved well (d/Δx > 5–10), and in which it is located sufficiently far from the boundary of the box. In most other cases, either Arthur or Roscoe orientations are observed. If heterogeneities are one element in size, Coulomb angles are thus unlikely to develop irrespective of the employed numerical resolution. Whereas differences in element types and rheology do have consequences for the maximum obtainable strain rates inside the shear bands, they only have a minor effect on shear band angles. Shear bands, initiated from random noise or from interactions of shear bands with model boundaries or other shear bands, result in stress heterogeneities with dimensionless length scales d/Δx ~ 1–2. Such shear bands are thus expected to form Roscoe or Arthur orientations, consistent with the findings in previous numerical models. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Numerical modelling of the deformation of brittle crust on geological timescales is an important but challenging topic that has received an increasing amount of attention over the last two decades (Braun et al., 2008; Buck et al., 2005; Buiter et al., 2006; Fullsack, 1995; Gerbault et al., 1999; Hobbs and Ord, 1989; e.g., 2008; Lyakhovsky et al., 1993; Moresi et al., 2007; Poliakov et al., 1996; Poliakov et al., 1993; e.g., Poliakov et al., 1994; Popov and Sobolev, 2008). Experimentally, it has been demonstrated that brittle rock failure can be described by a pressure-dependent yield stress (‘Byerlees law’),

⁎ Geophysical Fluid Dynamics, Department of Earth Sciences, ETH Zentrum, Sonneggstrasse 5, 8092 Zurich, Switzerland. Tel.: +41 44 633 7539; fax: +41 44 633 1065. E-mail address: [email protected]. 0040-1951/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2009.08.042

irrespective of rock type (e.g., Byerlee, 1978; Scholz, 1990; Ranalli, 1995). In numerical models, this is frequently approximated by a (phenomological) Drucker–Prager or Mohr–Coulomb yield criteria in which the maximum differential stress of rocks is a function of its normal stress, the internal angle of friction ϕ and the cohesion C. In addition to this, modelling of plasticity requires the specification of a (phenomological) plastic flow rule. Typically, a flow rule is used which has the dilation angle Ψ as material parameter. Rocks have Ψ ≪ϕ or even Ψ ~ 0° (so called non-associated plasticity) which results in spontaneous localization of deformation into narrow bands (Rudnicki and Rice, 1975; Vermeer and de Borst, 1984). Besides a plastic yield criteria (and corresponding plastic flow potential), one also has to specify the rheology of non-yielded materials. Here, differences exist between various modelling groups, where some employ viscous rheologies, whereas others treat non-yielded materials as viscoelastic. Despite the widespread implementation of plasticity in geodynamic numerical models, some seemingly basic questions about

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

plasticity remain unclear. One of the most prominent is: which shear band angle (θ) should we expect in a brittle Mohr Coulomb material? From a theoretical point of view, this question was studied by Vardoulakis (1980) and Vermeer (1990) for elastoplastic rheologies. Vermeer (1990) demonstrated with a bifurcation analysis and with a post-failure analysis that a range of angles are mechanically stable: θ = 45-F

ψ 2

ð1aÞ

(Roscoe angle — Roscoe, 1970) θ = 45-F

ϕ 2

ð1bÞ

(Coulomb angle — Coulomb, 1773) ϕ+ψ θ = 45-F 4

37

Coulomb angles can be observed in fine grained sands, independent of the boundary conditions (Wolf et al., 2003). Note, however that the reason for the different shear band angles in experiments is likely to be caused by differences in material properties, whereas in numerical models they occur due to the non-uniqueness of the plasticity models. This brief (and necessarily incomplete) summary thus shows that there is currently no consensus on what the stable shear band orientation is, and how this orientation is affected by the numerical method and the employed rheologies. The goal of this work is therefore to systematically test the effects of numerical parameters, initial conditions and rheology on the angle of shear bands obtained in a simple setup of compression and extension of a brittle material. These insights might than help to interpret and understand the results in more complex numerical models. 2. Mathematical and numerical formulation

ð1cÞ

(Arthur or intermediate angle — Arthur et al., 1977) The Coulomb angle results in the fastest rate of pressure reduction inside the shear band. The Arthur (or intermediate) angle is related with (elastic) unloading of material outside the shear band and hence causes the fastest rate of macroscopic softening upon shear band formation. The Roscoe angle, on the other hand, is admissible but does not result in softening. Angles outside the Roscoe–Coulomb range are not expected to occur, since formation of such angles would result in a macroscopic strain hardening of the material. Yet, all angles between the Roscoe and Coulomb angle are feasible. A more recent analysis by Lemiale et al. (2008) focused on a bifurcation analysis of a viscoplastic solid (without post-failure deformation) and found that both the Roscoe angle and the Coulomb angle are possible solutions. Numerically, Moresi et al. (2007) found the Roscoe angle to be stable, using a viscoplastic finite element approach with linear elements and moving integration points. Lemiale et al. (2008), however, found the Coulomb angle to be stable using the same numerical code. The differences between the two studies were attributed to differences in numerical resolution (which was considerably higher in the Coulomb angle case), although the authors noted that it could not be excluded that parameters such as type of integration scheme might affect results. Popov and Sobolev (2008) employed an arbitrary Lagrangian– Eulerian FEM code with a viscoelastoplastic rheology and found that shear bands, initiated from random noise, initially develop at a range of orientations (including the Coulomb angle), but that the Arthur angle dominates after several percents of strain. Hansen (2003) employed a meshless finite element formulation with an elastoplastic rheology and obtained shear band angles close to the Coulomb orientation. Mancktelow (2006) performed viscoelastoplastic simulations using a finite element formulation and obtained Arthur shear band orientations. Gerya and Yuen (2007) employ a viscoelastoplastic finite difference methodology and obtain Arthur shear band orientations in extensional experiments both in incompressible cases and in cases with a nonzero dilation angle. Buiter et al. (2006) performed numerical experiments in order to reproduce laboratory experiments of brittle failure under both extension and compression. Although a qualitative agreement between numerical and analogue models was found, differences exist between numerical models with respect to location, spacing and dip angle of shear zones, which ranged from Roscoe to Coulomb angles. Interestingly, the wide spread in shear band orientations can also be observed in laboratory experiments in which sand is employed to simulate brittle materials. A literature summary given in Vermeer (1990) indicates that Roscoe angles form in coarse grained sands whereas

2.1. Mathematical formulation An incompressible Boussinesq approximation is employed which is given by ε˙ ii = 0;

ð2Þ

∂τij ∂P + = ρgi ∂xj ∂xi
i + where ε˙ ij = 12 ∂v ∂x

ð3Þ

−

j

∂vj ∂xi



indicates strain rate, vi velocity, P = − σ3ii

pressure, σij stress, τij = σij + P deviatoric stress, xi spatial coordinates

(in the 2D case employed here: x1 = x, x2 = z) ρ density and gi the gravitational acceleration (the sign convention taken here assumes that compressive stresses are negative). The rheology of rocks is assumed to be Maxwell viscoelastoplastic, which (for an incompressible rheology) is given by ε˙ ij =

1 1 Dτij ˙ ∂Q ; τ + +λ 2η ij 2G Dt ∂σij

ð4Þ

where η is the effective viscosity, G the elastic shear module, t time, λ̇ the plastic multiplier and Q the plastic flow potential (e.g., Moresi et al., 2007). The Jaumann objective derivative of the deviatoric stress tensor is given by Dτij ∂τij ∂τij = −Wik τkj + τik Wkj ; + vk Dt ∂t ∂xk where Wij =

1 2




∂vi ∂xj

ð5Þ

 ∂v − ∂xj is the vorticity. i

Rocks fail plastically if differential stresses exceed the yield stress. A commonly used, phenomenological, yield stress formulation under low temperatures is Mohr Coulomb plasticity, for which (in 2D) the yield function F and plastic flow potential can be written as (e.g., Vermeer and de Borst, 1984) F = τ⁎−σ ⁎ sinðϕÞ−c cosðϕÞ; Q = τ⁎−σ ⁎ sinðψÞ; where τ⁎ =

ð6Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σxx −σzz 2 2 + σ xz ; σ ⁎ = −0:5ðσxx + σzz Þ, c is cohesion, 2

ϕ friction angle, ψ the dilation angle (zero in incompressible cases). ϕ is typically 30° for dry and less for wet rocks, such as serpentine (e.g., Ranalli, 1995). Plasticity is activated only if stresses are above the yield stress: ˙ ≥ 0; F ≤ 0; λ ˙ F = 0: λ

ð7Þ

38

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

It should be emphasized at this point that it is crucial to employ the dynamic pressure for the evaluation of the yield function (Eq. (6)). If this is not done, resulting the shear bands will form at 45°. 2.2. Numerical formulation In order to solve the equations numerically, the time derivative in Eq. (4) is discretized in an implicit manner, which results in τij = 2μ ve ðε˙ ij − ε˙ ij Þ + χ τ˜ ij ; pl

old

ð8Þ

  γpl c = c0 + ðc∞ −c0 Þ min 1; ; γ0

with μ ve =

χ=

1 η

of the code, this is done by Picard iterations. More efficient strategies exist, in particular Newton–Rhapson iterations with a consistent tangent stiffness matrix (de Souza Neto et al., 2009). This is not yet implemented in the current version of the numerical code. A recent comparison of Newton–Rhapson versus Picard iterations suggests that the width of the shear bands is smaller and have larger strain rates in the case of Newton–Rhapson iterations (Popov and Sobolev, 2008). Shear band angles, however, seem to be unaffected by the type of iterations and the main conclusions made here are thus likely to remain valid. Strain weakening is applied to cohesion, according to

1 ; 1 + GΔt

pl where γpl = ∫ε˙ 2nd dt indicates plastic strain and γ0 = 0.1 is a critical value.

1 ; 1 + GΔt η

2.3. Numerical method

old old old τ˜ ij = τ ij + vk

∂τold ij ∂xk

old old

old

Spelled out in two dimensions, the strong form of the governing force balance equations is

old

−Wik τkj + τik Wkj ;

where τijold are deviatoric stresses from the last time step, and Δt is the employed time step (Kaus, 2005; Kaus & Podladchikov, 2006). If stresses are below the yield stress (F < 0), plastic strain rates ε̇iplj are zero and rheology is viscoelastic. Note that the viscoelastic viscosities μ νe are dependent on the employed numerical time step. A purely viscous rheology (τij = 2ηε̇ij) is recovered if G → ∞ (since then χ = 0, μνe =η). In a typical time step one first assumes that plastic strain rates are zero and the rheology is fully viscoelastic. If initial stresses are above the yield stress (F(τij) > 0, plastic strain rates should be computed in such a manner that stresses are brought back to the yield surface (F = 0). The rheological expression (Eq. (8)) is than given by τy = 2μ ve ðε˙ 2nd − ε˙ 2nd Þ + χ τ˜ 2nd ; pl

old

ð9Þ

∂vx ∂vz p + =− ; κ ∂x ∂z

−

−

which can also be written as τy = 2μ vep ε˙ 2nd + χ τ˜ 2nd ; old

where ε˙ 2nd =

ð10Þ

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r
ffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
  ˜ old 2 ˜ old 2 old ε˙ xx − ε˙ zz 2 τ old 2 xx − τzz ˙ ˜ ˜ + εxz ; τ2nd = + ðτxz Þ ; τy = 2 2

σ ⁎ sinðϕÞ + c cosðϕÞ is the yield stress and μ vep the viscoelastoplastic effective viscosity, given by μ vep =

old τy −χτ˜ 2nd : 2ε˙ 2nd

ð11Þ

The full viscoelastoplastic rheology employed here can thus be written as old

τij = 2μ vep ε˙ ij + χ τ˜ ij ;

ð12Þ

with

i

μ vep

ð14Þ

8 > ˜ old > < τy −χτ2nd ; ˙ 2 ε2nd = > > : μve;

9 > > if F ðτ initial Þ > 0 = ij

if F ðτijinitial Þ≤0

> > ;

ð13Þ

The present formulation can deal both with viscoelastoplastic as well as with viscoplastic rheologies (if G → ∞). Typically, iterations (indicated by the superscript i in Eq. (13)) are required to compute plastic viscosity since changing the effective viscosity locally affects the stress distribution in surrounding points. In the present version

      ∂P ∂ ∂vx ∂ ∂vz ∂vx i i + + + 2μ vep μ vep ∂x ∂x ∂z ∂x ∂x ∂z ! old ∂τ˜ xx ∂τ˜ old xz ; + =− ∂x ∂z       ∂P ∂ ∂νz ∂ ∂νz ∂νx i i + + + 2μ vep μ vep ∂z ∂x ∂z ∂z ∂z ∂x ! old old ∂ τ˜ xz ∂ τ˜ zz + ; = ρg− ∂x ∂z

ð15aÞ

ð15bÞ

ð15cÞ

where κ is a penalty parameter used to enforce incompressibility. Eq. (15a–c) are solved numerically with a recently developed finite element code MILAMIN_VEP, which solves the unknowns using a velocity pressure formulation. MILAMIN_VEP is based on the efficient MATLAB-based Stokes solver MILAMIN, which employs an iterative penalty method to enforce incompressibility (Cuvelier et al., 1986) and which is described in detail in Dabrowski et al. (2008). The main modifications to the MILAMIN solver are τ̃ old ij -related terms in Eq. (15b,c) which have to be added to the right hand side, as well as an iterative scheme that is required to determine μ iνep in case of plastic failure. The error of plasticity iterations ep are computed at every iteration step i with i

ep =

∥vi −vi−1 ∥∞ ∥vi ∥∞

where vi indicates nodal velocities at iteration step i. Iterations are terminated once ep < 10− 4 or once more than 25 iterations have been performed. Typically, the maximum number of iterations is reached during the first few time steps, but after that, iterations converge more rapidly (presumably due to the effects of geometric and strain softening). Tests with smaller cutoff values and more iterations revealed that the shear zones obtain a slightly larger strain rates. This, however, does not change the conclusions with respect to shear band angles. In all simulations, a constant time step of 30,000 years was employed. Compared to the original MILAMIN code, I have also added additional element types and shape functions, thermo-mechanical

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

coupling, particle-based advection of material properties, remeshing, phase transitions, sedimentation. Both quadrilateral and triangular elements are implemented. The quadrilateral elements use both bilinear (4-node) shape functions for velocity and a constant shape function for pressure (Q1P0), or (9-node) biquadratic shape functions for velocity and linear discontinuous shape functions for pressure (Q 2P− 1). Triangular elements are identical to the MILAMIN implementation, which uses quadratic shape functions for velocity and linear, discontinuous shape functions for pressure (P2P− 1). Lower and upper cutoff viscosities of 1019 and 1025 Pa s are employed. The penalty factor κ (Eq. (15a)) was taken to be 1000 times the maximum viscosity in the domain, which results in an incompressible solution in a few iteration steps. Viscoelastoplastic simulations have zero initial differential stresses. In all simulations, the code is employed in a purely Lagrangian manner and simulations are performed until the elements become too deformed (as indicated by negative jacobians for the elements), or until an overall strain of 2% is reached. In the case of viscoplastic simulations, it is important to use an initial stress field that is in agreement with the physics of the problem, since otherwise the plasticity algorithm fails to return to the (correct part of the) yield surface. One way to do this, is to compute the initial homogeneous state of the problem (without heterogeneity), and use this state of stress as an initial guess (Lemiale et al., 2008). A disadvantage of this method, however, is that it can be cumbersome to derive the correct initial state of stress for more complicated boundary conditions or initial geometries. For this reason, I use a different approach, in which the boundary velocity vb is incrementally increased during the first time step, according to vb =

 3 k 0 vb nk

where v0b is the desired boundary velocity, k the iteration step and nk the total number of iterations. It was found that shear band orientations are sensitive to the value of nk during the first time step (nk is taken to be one in all subsequent steps). Values employed in most simulations are nk = 25 for extension and nk = 1 for compression. These values are selected, as they maximize the potential for Coulomb angles, while simultaneously ensuring physically realistic stress states at the end of the first time step. Employing significantly smaller nk values during extension results in unphysical initial stress states, whereas the use of larger nk values during compression increases shear band orientation towards the Arthur angle (see Fig. 9 and related discussion). No boundary velocity iterations are required for viscoelastoplastic simulations, as the initial stress evolution in these simulations increases with time as a consequence of viscoelasticity. The time step in these simulations, however, should be chosen sufficiently small to prevent failure of the whole model domain during the first time step. With the values employed here, it takes 5–10 time steps until the system fully yields. Larger values of elastic shear modulus G thus require a proportionally smaller time step. 3. Results Shear bands do not initiate in numerical models with a perfectly homogeneous setup. To localize deformation, one therefore has to introduce some form or heterogeneity. This can consist of random variations in material properties, a heterogeneity of specified size, or random noise that develops due to round-off errors in the numerical approach (as is typically for explicit geodynamic numerical models). Here, it is demonstrated that the length scale of such heterogeneities has a large impact on the resulting shear band angles. For this reason, the model setup considered here consists of a single viscous heterogeneity of length scale d, embedded in a viscoplastic or a viscoelastoplastic crust (Fig. 1). The box is subjected to either compression or extension with

39

Fig. 1. Model setup. See Table 1 for material parameters.

constant background strain rate. Material parameters are given in Table 1 and are chosen in such a manner that shear localization occurs, starting at the heterogeneity. In all simulations, the angle of the resulting shear bands was determined in an automated manner by computing the location of maximum strain rate invariant at a distance of d/4 and d/4 + 2000 m from the sides of the inclusion. The average value from the left and right shear band is reported. Experiments in which the sampling location and width was varied revealed that the angle can be reproduced within +/− 3°. Unless stated otherwise, Q1P0 elements are employed. 3.1. Effect of friction angle The effect of the employed numerical resolution on the shear band angle is shown on Fig. 2 under both extensional and compressional boundary conditions. At small numerical resolutions, resulting shear band angles are close to the Arthur angle (52.5° respectively 37.5°), whereas at larger numerical resolutions they are closer to the Coulomb angle (60° respectively 30°). Maximum obtained strain rates inside the shear zone are larger for larger numerical resolutions. Convergence of the results with respect to shear band angle is obtained at a resolution of ~ 400 by 100 elements. The aspect ratio of the elements does not appear to play a major role, suggesting that mesh alignment effects (frequently occurring in strain-softening problems; see e.g. De Borst, 1991) are of minor importance here. A more systematic analysis on the effect of friction angle on the resulting shear band angle (Fig. 3), confirms that shear band angles are closer to 45° at smaller numerical resolutions. Small numerical resolutions and friction angles smaller than ~ 20° result in Roscoe

Table 1 Parameters employed in the models. Parameter

Symbol

Dimensional value

Dimensionless value

Width Height Weak inclusion size Background strain rate Density times grav. acceleration Viscosity inclusion Viscosity of viscoplastic layer Initial cohesion Cohesion after weakening Friction angle inclusion Friction angle visco(elasto)plastic layer Minimum yield stress (cutoff value) Maximum yield stress (cutoff value) Elastic shear module

W H d ε̇bg ρg

40 km 10 km 400 – 3200 [800] (+/−)10− 15 s− 1 27,000 kg m− 2 s- 2

4 1 0.04–0.32 [0.08] (+/−)1 2700

μi μ

1020 Pa s 1025 Pa s

1 105

c c∞ ϕi ϕ

40 × 106 Pa 5 × 106 Pa 0° 0–45° [30°]

400 20 0° 0–45° [30°]

σy,min

20 × 106 Pa

200

σy,max

1000 × 106 Pa

104

G

5 × 1010 Pa or ∞

5 × 105 or ∞

Characteristic values for length, viscosity and time are L⁎ = 10 km, μ = 1020 Pa s, and t = 1015 s− 1, respectively. Cohesion is weakened in a linear manner between plastic strain = 0 and 0.1.

40

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

Fig. 2. Effect of numerical resolution on shear band angle θ for both compressional (left) and extensional (right) cases for a viscoplastic rheology with ϕ = 30°. The inclusion length size d = 800 meter, and Q1P0 elements are employed. Thin colored lines indicate the Roscoe, Arthur and Coulomb orientation. Thick white lines indicate the measured shear band orientation. At small numerical resolutions (for which the heterogeneity is numerically less well resolved) the orientations are closer to the Arthur angle, whereas at larger resolutions they are closer to the Coulomb angle. Changing the aspect ratio of the elements has a minor effect on shear band angles (bottom row).

angles, whereas Arthur shear band angles are observed for larger values of the internal friction angle. At larger numerical resolutions, the shear band angles are between Arthur and Coulomb angles.

These results thus confirm the findings of Lemiale et al. (2008) that large numerical resolutions are required for Coulomb angles to develop.

Fig. 3. Effect of friction angle on resulting shear band angle for two different numerical resolutions and viscoplastic rheology. Lines indicate the Roscoe, Arthur and Coulomb angles, and symbols indicate measured shear band angles from numerical simulations. The inclusion length scale d = 800 m. All other parameters as in Table 1.

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

41

Fig. 4. a) Effect of normalized inclusion size on shear band angle for a resolution of 400 × 100 nodes, Q1P0 elements and viscoplastic rheology with ϕ = 30°. Simulations are performed with a resolution of 400 × 100 nodes and with increasing inclusion length scale. b) Effect of increasing the numerical resolution while keeping the inclusion length scale constant. Numbers in braces indicate the employed numerical resolutions. Coulomb angles are selected if the initial heterogeneity is numerically well resolved (d/Δx > ~ 5–10).

3.2. Effects of resolution and length scales Although numerical resolution clearly matters, it remains unclear whether this is because the heterogeneity is better resolved numerically or whether another effect plays a role. Therefore three additional sets of experiments have been performed. In the first experiment, the numerical resolution was kept constant but the inclusion size d was increased such that the non-dimensional parameter d/Δx increases (where Δx is the spacing between two horizontal grid points). Results show a weak dependence of shear band angle on d/Δx with small values (i.e. small numerical resolutions) favoring Arthur angles and larger values favoring Coulomb angles (Fig. 4a). A potential problem with this experiment is that with increasing inclusion size, the results might be influenced by the proximity to the boundaries of the model box. For this reason, a second experiment

was performed in which the inclusion size was kept constant (d = 800 m), but the numerical resolution was increased (such that d/Δx increases). Results show a clear dependence of shear band angle on d/Δx with Coulomb angles for values of d/Δx > 10–15, and a trend towards Arthur or Roscoe angles at smaller values of d/Δx (Fig. 4b). These results thus suggest that shear band angles are very sensitive to the manner in which the heterogeneity, from which the shear band initiates, is resolved numerically. In order to study effects of the proximity of the heterogeneity to the upper boundary, a third series of experiments was performed, in which d/Δx and d were kept constant but in which the size of the model box was increased. Results indicate that the proximity to the upper boundary indeed influences the shear band angle (Fig. 5). The clearest results are obtained for compressional boundary conditions, which show that the model box height should be at least 10 times larger than the inclusion size for Coulomb angles to

Fig. 5. Effect of model box size on shear band orientation. In all simulations, the inclusion length scale is maintained constant (d = 800 m). Numerical resolution is increased with increasing box size, such that the inclusion is resolved with the same number of elements in all simulations (d/Δx = 8). Material is viscoplastic and ϕ = 30°. a) strain rate plots that illustrate the steepening of the shear band angle with increasing box size. b) results of systematic simulations as a function of box size under both compression and extension. Results indicate a rotation of the shear band angle towards the Arthur angle for a model height < 10 km.

42

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

develop. Simulations with the same setup but without twice the gravitational acceleration (not shown here) gave a similar trend, which suggests that the effect observed here is not due to an increase of confining pressure, but can be solely attributed to the proximity to the boundaries. The results obtained here thus suggest that Coulomb shear band angles can be obtained if (1) the heterogeneity from which the shear band initiates is resolved numerically with at least 5–20 elements, (2) the heterogeneity is sufficiently far from the boundaries of the model box (at least 10 times the heterogeneity size). 3.3. Effect of element type So far results have been obtained with a single element type (Q1P0). Since this is a non-admissible element (see e.g., Cuvelier et al., 1986; Hughes, 1987), it is interesting to see how results change if higher-order elements are employed. A comparison of linear and quadratic shape functions for quadrilateral elements (Fig. 6), reveals that higher order shape functions result in shear band angles that are slightly closer to the Coulomb angle (34° instead of 36°). Triangular elements can be either fitted to the viscous inclusion (‘fitted P2P− 1’), or can be employed with a more or less regular distributed mesh (‘non-fitted P2P− 1’). Both cases result in shear band angles that are comparable to results obtained with Q2P− 1 elements. It can thus be concluded that the type of element (triangular versus quadrilateral)

has little effect on the shear band orientation, but that the type of shape function (linear versus quadratic) has a minor effect. 3.4. Viscoplastic versus viscoelastoplastic rheologies As outlined in the introduction, there are significant differences between various geodynamic modelling groups when it comes to the treatment of non-plastic rheologies. Whereas some use a viscous rheology, others use a (Maxwell) viscoelastic rheology. Rock deformation experiments clearly indicate that rocks behave elastically at small strains (e.g., Ranalli, 1995) which argues in favor of viscoelastic models. Yet, if one is interested in modelling brittle (plastic) behavior only, it can be argued that non-plastic behavior is irrelevant. Furthermore, implementing elasticity in viscoplastic numerical codes requires additional work, since stress tensors have to be advected along with the numerical grid. Moreover, a previous study found little differences between viscoplastic and viscoelastoplastic rheologies (Buiter et al., 2006), although the setup employed in this comparison was relatively complicated. Yet, it remains unclear to which extend elasticity influences results in a simpler setup. For this reason, a number of additional experiments have been performed in which the two rheologies are compared. Results for an extensional and compressional setup with a friction angle of 25°, reveals that there are little differences in the angle of the resulting shear bands (Fig. 7). There are however

Fig. 6. Effect of element type on shear band angle for a viscoplastic rheology with ϕ = 25° and a compressional setup. A resolution of ~ 40,000 nodes is used in all cases. For this setup, the Coulomb angle is 32.5° and the Arthur angle 39°. Higher order elements thus result in slightly steeper shear band angles as well as larger strain rates inside the shear band. The element shape (triangular of quadrilateral) has little effect on the shear band angle.

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

43

Fig. 7. Comparison of viscoplastic versus viscoelastoplastic simulations under extension and compression (after 2% strain, with ϕ = 30° and a numerical resolution of 800 × 200 nodes). In both cases, viscoelastoplastic shear band orientations are slightly closer to the Coulomb angle. Under compression, elasticity results in larger strain rates inside the shear band, and in smaller matrix strain rates. Under compression, the situation is reversed and shear bands are more pronounced in the case with viscoplastic rheology.

differences in the sharpness of the shear zones (which more pronounced in viscoplastic cases), as well as in the distribution of strain rates. Under compression, maximum obtained strain rates are an order of magnitude larger in viscoelastoplastic cases, whereas the opposite effect is observed under extension. A more systematic study on the dependence of shear band angle on fiction angle for viscoelastoplastic rheologies (Fig. 8a) gives very similar results as for viscoplastic cases (Fig. 3b). Similar trends are also observed for simulations in which the inclusion length scale is kept constant (d = 800 m) and the numerical resolution is increased (Fig 8b). Simulations, in which the inclusion is poorly resolved numerically, favor Arthur (or Roscoe) angles, whereas simulations with larger numerical resolutions result in Coulomb angles. The main differences between viscoelastoplastic and viscoplastic cases are

observed in compressional setups, where the tendency for Arthur angles is stronger in the presence of elasticity. This is consistent with the findings of Vermeer (1990), who showed that the Arthur angle is related to unloading of stresses outside the shear band. If elasticity is present, this unloading depends on the available elastic stresses (which are larger for compressional setups). Viscous deformation is irrecoverable and hence unloading should be less pronounced. 3.5. Effect of initial iterations in viscoplastic simulations A (computational) advantage of viscoelastoplastic over viscoplastic simulations is that stresses build up gradually in during the course of a numerical simulation. Plastic failure in the vicinity of the viscous heterogeneity thus only occurs after a finite amount of deformation,

Fig. 8. Simulations which address the effect of elasticity on shear band orientation. a) Friction angle versus shear band orientation for a viscoelastoplastic rheology for a numerical resolution of 400 × 100 nodes. b) Shear band orientation versus normalized numerical resolution. Simulations have been performed with a constant inclusion length scale of d = 800 m and with increasing numerical resolution. As in viscoplastic cases, low resolutions result in Roscoe or Arthur shear band angles, whereas Coulomb angles develop at large resolutions.

44

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

rheologies under extension and compression (for d = 800 m and ϕ = 30°). In these simulations, the maximum strain rate was computed at a height of 1200 m. Results show that under compression, viscoelastoplastic cases always result in larger strain rates inside the shear band (Fig 10). Interestingly, both viscous and viscoelastoplastic simulations result in a strain rate which saturates with increasing resolution. Extensional cases, however, show such saturation behavior only for viscoplastic rheologies (which have larger strain rates than viscoelastoplastic simulations at low resolutions). These results thus suggest that the viscoplastic simulations performed here do not have a strong mesh dependency. More work is required, however, to fully understand these effects. 3.7. Multiple heterogeneities

Fig. 9. Effect of number of initial iterations nk on shear band orientation for viscoplastic rheologies and Q1P0 elements. Extensional cases with nk < 10 resulted in numerical artifacts and are not shown. The manner of initializing the stress-state thus has an impact on the results, with compressional cases favoring small numbers of iterations, whereas extensional cases require 10 or more.

which is approximately what one would expect in physical experiments. In viscoplastic simulations, however, the complete model fails in a plastic manner during the first time step, which might cause difficulties for the plasticity algorithm. For this reason, I employ an algorithm in which the applied boundary velocity is slowly increased during the first time step. It is interesting to study the dependence of the resulting shear band angle on the number of iterations that are performed. Results (Fig. 9) indicate that at least 7 iterations are required for extensional cases and that results converge after ~15 iterations. Compressional cases, on the other hand, can be performed with a single iteration, and results in slightly steeper angles if more than ~ 10 iterations are employed. 3.6. Consequences for strain rate inside shear bands The effect of numerical resolution on the maximum strain rate inside shear bands was studied for both viscoelastoplastic and viscoplastic

All simulations performed in here focused on a relatively simple setup in which a single heterogeneity was responsible for strain localization. It is instructive to see how such simulations can be compared to more realistic cases in which several heterogeneities are present. Therefore an experiment was performed in which 4 heterogeneities with length scale d = 800 m were randomly distributed. Results for a viscoelastoplastic rheology as a function of strain (Fig. 11) show that shear bands with Arthur to Coulomb angles are initiated at the corners of the heterogeneities (somewhat depending on where the inclusion is located with respect to the upper boundary). Interactions between two shear bands that initiated at two heterogeneities, however, might even induce shear bands at an angle that is significantly shallower than the Coulomb angle. Moreover, shear bands are in many cases curved. An interesting observation is also that shear bands, which are reflected at the lower boundary, increase their shear band angle. In the light of the simulations presented here, this can be understood in the following manner: Once a shear band reaches a boundary, it serves as a location to initiate a new shear band. The heterogeneity in the stress field, however, is only a few elements in size and thus significantly smaller than d, with as a result shear band angles that are closer to the Arthur or Roscoe angle (e.g. Fig. 4). 4. Discussion The simulations performed here thus highlight the importance of resolving the initial viscous heterogeneity sufficiently well numerically. Further insight in this effect can be obtained by plotting the second invariant of the deviatoric stress field as a function of

Fig. 10. Effect of numerical resolution on maximum strain rate at a distance d above the inclusion. Vertical resolution is proportional to the horizontal resolution.

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

45

Fig. 11. Shear localization as a function of bulk strain for a model with 4 randomly distributed viscous inclusions with d = 800 m. Rheology is viscoelastoplastic, ϕ = 30° and a resolution of 1200 × 300 nodes is employed with Q1P0 elements. Magenta, black and yellow lines indicate Roscoe, Arthur and Coulomb angles respectively. Shear bands initiated from inclusions have an Arthur to Coulomb orientation. Once they reflect from side boundaries, however, their orientation rotates towards the Roscoe angle (see for example white arrows at 2.0% strain).

Fig. 12. Second invariant of deviatoric stress tensor around viscous heterogeneities computed with different numerical resolutions, Q1P0 elements and a viscous rheology. Small d/Δx values result in smaller absolute stress values.

46

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47

numerical resolution for a given heterogeneity size for a purely viscous rheology (Fig. 12). For small values of d/Δx, the two high stress ‘lobes’ at the corners of the inclusion are not well resolved. As a result, maximum stresses are smaller than in cases with larger d/Δx. In addition, stress differences (maximum minus minimum stresses) are smaller. Thus, the initial (pre-plastic) stress field is more homogeneous at smaller numerical resolutions, which appears to favor Roscoe or Arthur angles. The findings of this work can be used to reconcile previous numerical experiments on shear localization, which appear to have given contradictory results. Most numerical simulations reported shear band angles from simulations in which localization was induced by random noise (in either the viscosity field or in the plastic parameters). Doing such, induces stress heterogeneities that are on the order of one or several elements in size. The simulations performed here, suggest that in this case, resulting shear bands should have the Arthur orientation unless friction angles are 20° or smaller, where Roscoe angles dominate. Both Lemiale et al. (2008) and Hansen (2003) obtained Coulomb orientations in their simulations, but in both cases a heterogeneity was present that was numerically resolved sufficiently (d/Δx ~ 5 in Lemiale's model, whereas d/Δx ~ 5 in Hansen's models). Several parameters, however, remain to be studied. Firstly, I have employed a very specific algorithm to enforce incompressibility (an iterated penalty method). The code used by Lemiale et al. (2008) uses a different algorithm to compute incompressibility and it is unclear to which extend this affects model results. Secondly, some of the results obtained here seem to indicate that, at least in viscoplastic cases, the maximum obtainable strain rate inside shear zones saturates to a given value with increasing numerical resolution. This indicates that viscoplastic shear localization is not mesh sensitive, as opposed to non-regularized elastoplastic localization (e.g., Vermeer and de Borst, 1984; de Souza Neto et al., 2009). In the computational mechanics literature there is a rich tradition of numerically solving of strain-softening problems. Here, it was realized several decades ago that non-uniqueness of the solutions results in mesh-sensitivity of the results particularly for rateindependent elastoplastic materials (for which no internal length scales exists). Several solutions have been proposed to overcome this sensitivity. One way, called the Cosserat continua, is to add additional rotational degrees of freedom in the system (e.g., Muhlhaus and Vardoulakis, 1987; De Borst, 1991; Deborst and Sluys, 1991). Another way is to add viscosity (which effectively makes the system viscoelastoplastic) which works in both quasi-static and quasi-dynamic cases (see, e.g. Needleman, 1988; Wang et al., 1996; Wang et al., 1997). This is many ways similar to what is done in geodynamical models, although it should be noted that inertial forces are typically negligible in models of tectonic deformation and that therefore only the quasi-static results are directly applicable to such setups (unless one wants to include the earthquake cycle). Needleman (1988) performed a one-dimensional analysis of shear localization, and demonstrated that solutions are regularized once a rate-dependent material is added to the setup. He however also demonstrated that the length scale of localization is related to the length scale of the initial heterogeneity, which is unfortunate if applied to geodynamic cases as we can never be sure what the initial heterogeneity length scale was. The results obtained here (Fig. 10) also demonstrate that convergence with respect to strain rate occurs in viscoplastic simulations. Yet, this topic deserves more attention in the future, in particular for viscoelastoplastic rheologies in a geodynamic context. 5. Conclusions Numerical simulations have been performed to address the effects of numerics, geometry and rheology on the orientation of shear bands in geodynamic models. It was demonstrated that it is possible to

obtain the Coulomb orientation. This, however, only occurs if the heterogeneity, which induces this shear band, is resolved with a sufficiently large numerical resolution. Effects of rheology, type of element and iteration strategy were found to play a second order effect. If the heterogeneity is not resolved with sufficiently large numerical resolution (at least 5–10 elements), or if it is situated close to a boundary of the numerical box, resulting shear band orientations are closer to the Arthur or Roscoe angle, independent of the employed resolution. Acknowledgements I thank Hans Muhlhaus, Taras Gerya , and Dave May for helpful discussions and feedback on this work. The constructive reviews of Luc Lavier, Dave May and an anonymous reviewer helped to clarify the text. MILAMIN is available from www.milamin.org. References Arthur, J.R.F., Dunstan, T., Al-Ani, Q.A.J., Assadi, A., 1977. Plastic deformation and failure of granular media. Géotechnique 27, 53–74. Braun, J., Thieulot, C., Fullsack, P., DeKool, M., Beaumont, C., Huismans, R., 2008. DOUAR: A new three-dimensional creeping flow numerical model for the solution of geological problems. Physics of the Earth and Planetary Interiors 171 (1-4), 76–91. Buck, W.R., Lavier, L., Poliakov, A., 2005. Modes of faulting at mid-ocean ridges. Nature 434 (719–723). Buiter, S.J.H., et al., 2006. The numerical sandbox: Comparison of model results for a shortening and an extension experiment. In: Buiter, S.J.H., Schreurs, G. (Eds.), Analogue and Numerical Modelling of Crustal-Scale Processes. Geological Society of London Special Publications. Geological Society of London, pp. 29–64. Byerlee, J.D., 1978. Friction of rocks. Pure and Applied Geophysics 116, 615–626. Coulomb, C.A., 1773. Sur l'application des règles des maximis et minimis à quèlques problèmes de statique relatifs à l'architecture. Mémoires de Mathématique et de Physique. Académie Royale des Sciences 7, 343–382. Cuvelier, C., Segal, A., van Steemhoven, A.A., 1986. Finite element methods and Navier– Stokes equations. Mathematics and its Applications. Reidel Publishing company, Dordrecht. 483 pp. Dabrowski, M., Krotkiewski, M., Schmid, D.W., 2008. MILAMIN: MATLAB-based FEM solver for large problems. Geochemistry Geophysics Geosystems. doi:10.1029/ 2007GC001719. De Borst, R., 1991. Numerical modelling of bifurcation and localisation in cohesivefrictional materials. Pure and Applied Geophysics 137 (4), 367–390. de Souza Neto, E.A., Periae, D., Owen, D.R.J., 2009. Computational methods for plasticity: theory and applications. Wiley. 814 pp. Deborst, R., Sluys, L.J., 1991. Localization in a Cosserat continuum under static and dynamic loading conditions. Computer Methods in Applied Mechanics and Engineering 90 (1–3), 805–827. Fullsack, P., 1995. An arbitrary Lagrangian–Eulerian formulation for creeping flows and its application in tectonic models. Geophysical Journal International 120 (1), 1–23. Gerbault, M., Burov, E.B., Poliakov, A.N.B., Daignieres, M., 1999. Do faults trigger folding in the lithosphere? Geophysical Research Letters 26 (2), 271–274. Gerya, T.V., Yuen, D., 2007. Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems. Physics of the Earth and Planetary Interiors 163 (1–4), 83–105. Hansen, D.L., 2003. A meshless formulation for geodynamic modeling. Journal of Geophysical Research 108 (B11), 2549. Hobbs, B.E., Ord, A., 1989. Numerical simulation of shear band formation in a frictionaldilational material. Ingegnere e l'architetto 59, 209–220. Hughes, T.J.R., 1987. The finite element method. Dover Publications, New York. Kaus, B.J.P., 2005. Modelling approaches to geodynamic processes. PhD-Thesis Thesis, Swiss Federal Institute of Technology, Zurich, 263 pp. Kaus, B.J.P., Podladchikov, Y.Y., 2006. Initiation of localized shear zones in viscoelastoplastic rocks. Journal of Geophysical Research 111 (B04412). doi:10.1029/2005JB003652. Lemiale, V., Muehlhaus, H.B., Stafford, J., Moresi, L., 2008. Shear banding analysis of plastic models for incompressible viscous flows. Physics of the Earth and Planetary Interiors 171 (1–4), 177–186. Lyakhovsky, V., Podladchikov, Y., Poliakov, A., 1993. A rheological model of a fractured solid. Tectonophysics 226 (1–4), 187–198. Mancktelow, N.S., 2006. How ductile are ductile shear zones? Geology 34 (5), 345–348. doi:10.1130/G22260.1. Moresi, L., et al., 2007. Computational approaches to studying non-linear dynamics of the crust and mantle. Physics of the Earth and Planetary Interiors 163, 69–82. Muhlhaus, H.B., Vardoulakis, I., 1987. The thickness of shear bands in granularmaterials. Geotechnique 37 (3), 271–283. Needleman, A., 1988. Material Rate Dependence and Mesh Sensitivity in Localization Problems. Computer Methods in Applied Mechanics and Engineering 67 (1), 69–85. Poliakov, A., Podladchikov, Y., Talbot, C., 1993. Initiation of salt diapirs with frictional overburdens: numerical experiments. Tectonophysics 228, 199–210. Poliakov, A.N.B., Herrmann, H.J., Podladchikov, Y.Y., 1994. Fractal plastic shear bands. Fractals 2, 567–581.

B.J.P. Kaus / Tectonophysics 484 (2010) 36–47 Poliakov, A., Podladchikov, Y., Dawson, J.B., Talbot, C., 1996. Salt diapirism with simultaneous brittle faulting and viscous flow. Geological Society Special Publication 100, 291–302. Popov, A.A., Sobolev, S.V., 2008. SLIM3D: A tool for three-dimensional thermo mechanical modeling of lithospheric deformation with elasto-visco-plastic rheology. Physics of the Earth and Planetary Interiors 171 (1–4), 55–75. Ranalli, G., 1995. Rheology of the earth. Chapman and Hall, London. 413 pp. Roscoe, K.H., 1970. The influence of strains in soil mechanics, 10th Rankine Lecture. Geotechnique 20, 129–170. Rudnicki, J.W., Rice, J.R., 1975. Conditions for localization of deformation in pressuresensitive dilatant materials. Journal of the Mechanics and Physics of Solids 23 (6), 371–394. Scholz, C.H., 1990. The Mechanics of Earthquakes and Faulting. Cambridge University Press. 439 pp.

47

Vardoulakis, I., 1980. Shear band inclination and shear modulus of sand in biaxial tests. International Journal for Numerical and Analytical Methods in Geomechanics 4, 103–119. Vermeer, P.A., 1990. The orientation of shear bands in bi-axial tests. Geotechnique 40 (2), 223–236. Vermeer, P., de Borst, R., 1984. Non-associated plasticity for soils. Concrete and Rock 29, 1–65. Wang, W.M., Sluys, L.J., DeBorst, R., 1996. Interaction between material length scale and imperfection size for localisation phenomena in viscoplastic media. European Journal of Mechanics a-Solids 15 (3), 447–464. Wang, W.M., Sluys, L.J., de Borst, R., 1997. Viscoplastic modelling of localisation problems. Finite Elements in Engineering and Science 463–472. Wolf, H., Konig, D., Triantafyllidis, T., 2003. Experimental investigation of shear band patterns in granular material. Journal of Structural Geology 25 (8), 1229–1240.