We study the formation of localized shear zones during the layer-parallel extension of viscous multi-layers using two-dimensional numerical simulations based on the finite-difference method. For power-law viscous layers and a linear viscous embedding medium, the extended multi-layer develops boudins due to necking. For power-law viscous layers embedded in a power-law viscous medium, the extended multi-layer develops first distributed necks, and subsequently a localized shear zone with a vertical offset (with a size of several layer thicknesses) along the multi-layer. During the extension, the deformation style switches from distributed and symmetric necking to localized and asymmetric shearing. A localized shear zone develops in the viscous multi-layer although the rheology is everywhere strain-rate-hardening (power-law stress exponent >1) and no material softening and/or energy feedback mechanism (e.g., shear heating combined with a temperature-dependent viscosity) is applied. The shear localization is caused by structural softening because the formation of a localized shear zone decreases the bulk resistance and hence the work required to deform the multi-layer. A localized shear zone forms in the multi-layer when the distance between the stiff layers is approximately equal to or less than the layer thickness. The shear localization was observed in multi-layers with nine and with only three stiff layers.


During the layer-parallel extension of a stiff ductile layer, the strain can become localized due to a mechanical instability termed necking (e.g., Hutchinson and Neale, 1977; Smith, 1977). In geology, necking is often called boudinage, and the resulting structures are termed boudins or pinch-and-swell structures (e.g., Ramberg, 1955; Goscombe et al., 2004; Fig. 1A). In a viscous layer, a necessary condition for necking is that the rheology of the layer is nonlinear, e.g., power-law viscous with a power-law stress exponent n > 1 (e.g., Smith, 1977). Strain localization by necking can, hence, occur in a strain-rate-hardening material (stress increases with increasing strain rate; n > 1), and a softening behavior is not necessary (stress decreases with increasing strain rate; n < 0). Necking generates symmetric zones of reduced layer thickness (Fig. 1A), can occur on the outcrop scale (Fig. 1A), and is also an important mechanism during lithospheric rifting (e.g., Fletcher and Hallet, 1983; Zuber and Parmentier, 1986) or slab detachment (e.g., Duretz et al., 2012). While necking in viscous layers is relatively well understood, the formation of localized and asymmetric ductile shear zones is less well understood (Fig. 1B). Several mechanisms have been proposed for the generation of ductile shear zones, e.g., shear heating (e.g., Regenauer-Lieb and Yuen, 1998), weakening due to metamorphic reactions (e.g., Brodie and Rutter, 1987), or weakening by grain-size reduction (e.g., Platt and Behr, 2011). All of these mechanisms require material softening and/or energy feedbacks during the extension. Here, we investigate shear zone formation in viscous multi-layers without any material softening or energy feedback mechanism.

In both laboratory and numerical modeling studies, multi-layers can represent either any type of lithological layering (e.g., Fig. 1), the mechanical layering of the lithosphere due to the vertical variation of the yield strength, or simply the overall anisotropy of rock volumes. The behavior of multi-layers under extension has been studied experimentally (e.g., Cobbold et al., 1971; Allemand and Brun, 1991; Kidan and Cosgrove, 1996; Cosgrove, 2007; Gomez-Rivas et al., 2015), analytically (e.g., Bassi and Bonnin, 1988), and numerically (e.g., Schmalholz and Maeder, 2012). Schmalholz and Maeder (2012) studied multi-layer extension in power-law viscoplastic layers with a von Mises–type yield strength. The experimental and numerical studies reported the development of localized shear zones associated with boudinage, but the conditions for shear localization remain so far unclear. Here, we study the formation of localized asymmetric shear zones in power-law viscous multi-layers under extension with two-dimensional (2-D) numerical simulations based on the finite-difference method.


We solve the 2-D Stokes equations for incompressible flow without gravity using a finite difference–marker-in-cell code (Duretz et al., 2014). The governing partial differential equations are discretized on an Eulerian staggered grid, and the material properties are advected using Lagrangian markers. The Eulerian grid resolution is 801 × 1001 nodes (∼2.4 × 106 degrees of freedom). Such resolution is necessary in order to accurately resolve the shear zone formation (see the Appendix). Each cell initially contains 16 markers, which are regularly distributed throughout the mesh. Approximately 3000 time steps are required for each numerical simulation.

We consider a multi-layer composed of nine competent layers each of thickness H (as in Schmalholz and Maeder [2012]). The layer spacing (i.e., thickness of matrix inter-layers) D is set to 0.5H in the reference simulation. The impact of D on the model results is discussed in the Appendix. The model has an initial width of 12H and a height of 25H (Fig. 2). The distance separating the top competent layer from the model top and the bottom competent layer from the model bottom is 6H. The model is extended with a dimensionless constant background bulk strain rate forumla = 1, and free slip is imposed on all boundaries. For homogenous material, these boundary conditions provide a pure shear deformation. Both the layer and matrix rheology are described with a viscous power-law flow law with an effective viscosity, η, given by (Schmalholz and Schmid, 2012): 
where η0 is the reference viscosity and forumla is the second invariant of the strain-rate tensor. We apply a ratio of reference viscosities between layers and matrix of 50 (η0 = 1 in the matrix); the value of n for the layers is 10 and for the matrix is either 1 or 3. Two types of initial perturbations (see the Appendix) are considered in order to initiate the necking instability.


During the initial stage of extension of multi-layers with a pseudo-random geometrical perturbation of the layer interfaces, the bulk deformation is dominated by pure shear. During this stage, necking preferentially develops in the outer layers of the multi-layer. This first appearance of necking in the outer layers is expected from the folding and necking stability analysis, which predicts faster folding and necking growth rates in layers surrounded by a thicker matrix (e.g., Schmalholz et al., 2002). After further extension, each layer is independently necking and develops pinch-and-swell structures. For a linear viscous matrix (n = 1; model 1, Fig. 3A), the multi-layer is initially thinning. Multi-layer boudinage only starts developing for large extensional strain (>200%). Although pinches and swells have a tendency to align with progressive bulk extension, no prominent shear localization develops, and the necks remain symmetric. For the same model configuration, simulations with a power-law viscous matrix (n = 3; model 2, Fig. 3B) also first exhibit homogenous thinning. However, after a limited period of symmetric necking, localized shear zones develop across the multi-layer. The deformation around the shear zone is dominated by simple shear, and hence switches to an asymmetric mode. The intensity of layer necking is thus moderate, as most deformation is accommodated along a prominent shear zone. The shear zone causes a vertical offset along the layers (Fig. 3B), while for a linear viscous matrix, the average height of the layers remains the same (Fig. 3A). The switch from pure shear– to simple shear–dominated deformation, and the localization of strain into shear zones, is solely due to the power-law viscous rheology of the matrix. The effective viscosity decrease in the matrix allows the alignment of necks across the multi-layer and the strain localization within narrow regions of the matrix.

Additionally, we performed a simulation with a linear viscous matrix and a strong rheological contrast between layer and matrix (effective viscosity ratio = 100, n of layer = 20, and initial inter-layer spacing = 2H; model 3, Fig. 3C) as well as a model with power-law viscous matrix (n = 3) but with only three layers (effective viscosity ratio = 30, n of layer = 10, and initial inter-layer spacing = 0.5H; model 4, Fig. 3D). Model 3 shows that a strong rheological contrast enhances necking of the competent layers but does not promote shear localization. Model 4 exhibits shear localization into a single shear zone despite the restricted number of layers. In general, the bulk extension required to initiate the shear zones depends on the rheological contrast between layers and matrix and on the amplitude of the initial geometrical perturbation. For highly nonlinear layers and/or a large initial viscosity contrast, necks and shear zones develop after less bulk extension (<100%, Fig. 3C; Fig. DR2 in the GSA Data Repository1).

Strain localization for a power-law viscous matrix is characterized by an increase of the strain rate with respect to the bulk extension rate within the developing shear zones. The maximum value of forumla is ∼500× larger than the bulk extension rate forumla (Fig. 4A). In contrast, models with a linear viscous matrix do not localize considerably, and provide negligible amplification of the maximum (Fig. 4A). Monitoring the apparent area-averaged normal viscosity, 
(with V being the model area, and τxx the horizontal deviatoric stress), and the viscous dissipation (representing the incremental mechanical work), 
during extension highlights the intensity of the so-called structural softening. During structural softening, the bulk resistance (quantified by ηNA) of a rock volume decreases due to the formation of structures, here by the formation of boudins, but structural softening also occurs due to the formation of folds (e.g., Schmalholz et al., 2005; Schmalholz and Schmid, 2012). The material properties η0 and n remain constant during structural softening; i.e., the material remains strain-rate-hardening at each point. A significant decrease of ηNA and Φ (Figs. 4B and 4C) is observed for models with a linear matrix and a strong rheological contrast (model 3) or in models with a power-law viscous matrix (models 2 and 4). The strong layers get effectively mechanically disconnected by either necking or shear localization, which causes the structural softening. The decrease in ϕ shows that the development of boudins and shear zones decreases the amount of work needed to deform the multi-layer.


Strain localization can occur in homogeneous linear viscous material if (1) the applied boundary conditions (e.g., a velocity discontinuity; Brun, 2002) or (2) the geometrical configuration (e.g., flow around a rigid obstacle) impose a priori such localization independent of the material properties and rheology. This type of localization is termed kinematic strain localization (e.g., Brun, 2002; Braun et al., 2010). In our simulations, neither the pure shear boundary conditions nor the geometrical configuration impose a priori strain localization, because for a linear viscous rheology or a viscosity ratio between layers and matrix that equals one, no (dynamic) strain localization will occur. The presented simulations show the spontaneous development of localized shear zones in the absence of material softening or energy feedbacks. The shear zone formation arises naturally in a power-law viscous multi-layer under extension, because shear zone formation decreases the work required to deform the multi-layer. The multi-layer power-law behavior may be the simplest mechanism for dynamic shear localization in viscous material under extension, which requires the fewest model assumptions and fewest material properties (only two constant parameters, η0 and n). The applied material properties are relatively well constrained from laboratory rock deformation experiments. The applied effective viscosity ratio of 50–100 is reasonable for natural rock, and the applied values of n for the matrix (3) and the layer (10 and 20) are representative for dislocation creep and low-temperature plasticity (e.g., dislocation glide), respectively (e.g., Schmalholz and Fletcher, 2011). Low-temperature plasticity has been documented in laboratory rock deformation experiments for various lithologies such as calcite (e.g., Renner and Evans, 2002), dolomite (e.g., Davis et al., 2008), and olivine (e.g., Mei et al., 2010). Flow laws for low-temperature plasticity are usually described with an exponential flow law (e.g., Karato, 2008) which can be transformed into a power-law flow law (Schmalholz and Fletcher, 2011). The effective stress exponent corresponding to low-temperature plasticity can be considerably larger (>10) than the one for dislocation creep (∼3–5; e.g., Renner and Evans, 2002; Davis et al., 2008; Schmalholz and Fletcher, 2011).

The presented modeling results could explain the structures shown in Figure 1. Whereas Figure 1A depicts multi-layer necking with relatively large layer spacing (D > H), Figure 1B shows shear localization across a restricted number of layers with a tight spacing (D < H). Because both natural examples exhibit boudinage, the results suggest that the rheology of the boudinaged (necked) layers was nonlinear viscous (e.g., low-temperature plasticity). The overall structure in Figure 1A is similar to the one of model 3 (Fig. 3B), which encompassed large layer spacing (D = 2H) and linear viscous matrix. The results therefore suggest that the rheology of the natural matrix (fine-grained marble, Fig. 1A) was close to linear viscous (e.g., diffusion creep), which could explain the absence of a localized shear zone. The overall structure of Figure 1B is similar to the that of model 4 (Fig. 3D). which is characterized by a tight layer spacing (D = 0.5H) and a power-law viscous matrix. The results therefore suggest that the rheology of the natural matrix (carbonate) was likely nonlinear viscous (e.g., dislocation creep) allowing for the development of a localized shear zone across a few competent calcite layers. A thin calcite vein above the calcite multi-layer including the shear zone has not been offset by the shear zone, which indicates that the shear zone is restricted to the multi-layer (Fig. 1B). Microstructural analysis (Schmalholz and Maeder, 2012) was performed for the competent calcite layers shown in Figure 1B. This analysis indicates twinning of the calcite grains located in the swells, suggesting low-temperature plastic deformation in agreement with stress exponents >>1.

Structural softening due to the formation of necks and/or shear zones can occur in strain-rate-hardening, layered power-law viscous materials. Localization of deformation into necks can occur in power-law viscous single (Schmalholz et al., 2008) and multi-layers. However, localization of deformation into shear zones only occurs in multi-layers with power-law viscous layers and matrix.

We thank J.-P. Brun, M. Jessell, an anonymous reviewer, and editor R. Holdsworth for their insightful comments and advice. This research was supported by the University of Lausanne, Switzerland.

We consider two types of initial geometrical perturbations on the layer interfaces. The first type is a pseudo-random perturbation of similar amplitude (H/6) where each layer interface has a different perturbation. This type of perturbation exerts less conditioning on the final model result and allows for variable periodicity of structures (cf. Fig. 1). The second type of perturbation is sinusoidal. Such perturbation allows for reproducibility of model results and is useful to test model grid convergence (cf. Fig. DR1). In this case, the perturbation takes the form: 
with A = H/6, λ = 2H, Δx = tan(π/4). A and λ correspond to the amplitude and wavelength, respectively, of the perturbation; x and z represent the spatial coordinates The index k is 1 for the upper layer interface and 2 for the lower interface, and the index i is the layer number, which varies from 0 to 8 in a case with nine competent layers.
1GSA Data Repository item 2015245, Figure DR1 (the effect of layer spacing on the style of extension) and Figure DR2 (the effect of numerical resolution), is available online at www.geosociety.org/pubs/ft2015.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.