Strain localisation in mechanically layered rocks beneath detachment zones: insights from numerical modelling.

. We have designed a series of fully dynamic numerical simulations aimed at assessing how the orientation of mechanical layering in rocks controls the orientation of shear bands and the depth of penetration of strain in the footwall of detachment zones. Two parametric studies are presented. In the ﬁrst one, the inﬂuence of stratiﬁcation orientation on the occurrence and mode of strain localisation is tested by varying initial dip of inherited layering in the footwall with regard to the orientation of simple shear applied at the rigid boundary simulating a rigid hanging wall, all scaling and rheological


Introduction
Structural analysis of ductile deformation in crustal rocks reveals drastically different structures developing in isotropic (mylonite fabric in granites, Berthé et al., 1979) or layered material (extensional crenulation cleavage in micaschists, Platt and Vissers, 1980). In the latter case, relative geometry of inherited rheological layering and shear bands show systematic patterns, implying a mechanical impact of layering orientation on the development of shear bands. Nevertheless, exhumed field examples exhibiting the finite result of cumulated strain through time and space do not always allow the deciphering of initial layering impact in the localisation of strain. The present studies concentrate on the impact of such a pre-existing layering on occurrence, dynamics and preferred orientation of newly formed shear bands and associated finite strain in medium-(10 m) to large-scale (1 km (Handy et al., 2005), extension: e.g. Betics (Platt and Vissers, 1980) and Aegean (Jolivet et al., 2010), strike-slip: e.g. Red River Fault (Leloup et al., 1995) and Cap de Creus (Fusseis et al., 2006)) and are known to exist at larger scale from seismic reflection data (Torvela et al., 2013). Focus has been put here on detachment zones because their footwalls exhibit exhumed rocks deformed at the brittle-ductile transition and also because we can use active analogues to calibrate strain rates in the models. Mechanical layering is present at all scales in rocks as a result of sedimentation, diagenesis as well as inherited structures. Shear localisation in layered media results in strain partitioning between high strain domains, where initial layering is completely overprinted (shear zones) and low strain domains where the initial layering in the host rock is deflected but preserved. Cobbold et al. (1971) models show that strain localisation occurs along kink bands or shear bands in mechanically stratified media where strong and weak isotropic viscous layers alternate. At low strain, and for simplified boundary conditions, orientation of these bands can be successfully reproduced by mathematical models based on linear stability analysis established for anisotropic media (Biot, 1964;Cosgrove, 1976). Large-strain numerical models of folding and convection in anisotropic media (Muhlhaus et al., 2002(Muhlhaus et al., , 2004 have shown that anisotropy exacerbates shear instabilities. One may wonder then, under which conditions these local shear instabilities favour strain localisation at larger scale. Following previous studies (Mancktelow, 2006;Frehner et al., 2011;Schmalholz and Maeder, 2012), we make use of a numerical code that solves the momentum equation for stress and strain to perform dynamic numerical modelling of layered rocks. This approach allows measuring the variation of effective strength of the bulk shear zone as the initial layering deforms and monitoring the orientation of neoforming structures as well as stress distribution in the different layers.
After describing the model itself, we systematically describe the effect of initial orientation of the layering on strain localisation and orientation of new localised shear bands, when they occur. The modelling results are analysed in terms of finite strain and evolution of pressure and deviatoric stress in strong and weak layers during shear localisation. In a second part, the orientation of layering versus shear direction is fixed and chosen to favour shear localisation. Boundary conditions and constitutive mechanical properties are then changed to focus on the impact of length scale, effective viscosity contrast, plasticity, and strain rate on the degree of localisation. Strain localisation in the footwall of a detachment fault can eventually be considered as resulting from a combination of inheritance, boundary conditions and competition between ductile and brittle deformation mechanisms.

Geological motivations
We aim at modelling strain localisation within a shear zone forming in a mechanically layered medium located at the top of a detachment footwall (Fig. 1a). We focus on the inner structure of the shear zone and particularly on the relationships between initial layering, orientation and depth of penetration of newly formed shear bands as a function of the initial orientation of the mechanical layering. The analysis encompasses a geometrical description as well as a focus on the system dynamics at small scale (pressure jumps and microstructures) and large scale (structural softening and hardening within the shear zone). Our model is precisely designed to reproduce the western Gulf of Corinth detachment zone, among the best monitored in the world. This active shallow-dipping crustal-scale shear zone accommodates extension at a shear rate of 1 cm yr −1 , distributed at 10-15 km over a 1 km thick zone (Avallone et al., 2004). Crustal-scale thermo-mechanical models of this detachment have shown that kinematics of high angle normal faults in its hanging wall as well as the observed formation of out-of-sequence faults can only be reproduced for low effective viscosity of the order of 10 19 to 10 20 Pa s within the shear zone (Le Pourhiet et al., 2004) . However, with such a large strain rate, brittle-ductile interactions are still expected to occur at small scale, as attested by the repeated occurrence of small-magnitude earthquakes (Rigo et al., 1996). Lecomte et al. (2012) have shown that these events can be interpreted as slip on neoforming shear bands within a kilometer-scale shear zone. This active example can easily be compared to exhumed examples from the Aegean, which display similar rates of shearing close to the brittle-ductile transition (Gueydan et al., 2005).
Within the shear zone, the initial dip δ 0 of the layers (Fig. 1b) is counted positive when dipping towards the right of the model. The grey and white layers in Fig. 1a, b represent the alternating low and high competence layers. At such scale, each layer is 30 m thick. In these conditions, the background strain rate approaches 3 × 10 −13 s −1 . Confining pressure, set to 400 MPa, is intended to simulate a kilometer-scale shear zone buried at 15 km depth. To keep the results general and independent from the orientation of the shear direction versus gravity and thermal gradients, the problem has been deliberately simplified by removing gravity and temperature dependence of viscosity.

Choice of the rheological model
Plastic yielding is an important localisation parameter close to the brittle-ductile transition (Mancktelow, 2006;Gueydan et al., 2003). This is the motivation for introducing yield strength in the models. Since the plastic strain equation is solved numerically, we prefer introducing a realis-tic pressure-dependent finite yield strength calculated with the Mohr-Coulomb criterion (Vermeer and De Borst, 1984;Le Pourhiet, 2013), rather than an approximation. Indeed, a power law viscous rheology or a von Mises criterion would neglect the effects of pressure on strength, even if they facilitate the establishing of analytical solutions (Schmalholz and Fletcher, 2011). For similar reasons, elastic effects are not neglected in the models. All layers follow the same non-linear rheology formulation including elasticity, linear viscosity and pressure-dependent maximum yield strength. Elastic and plastic parameters are, however, kept constant with internal friction angle (φ) set to 30 • , cohesion (C 0 ) set to 20 MPa, while elastic bulk (K) and shear (G) moduli are set to 50 and 30 GPa, respectively. The competence contrast between weak and strong layers is achieved by varying only one parameter, namely their viscosity, fixed to a contrast of 100 in all experiments with a weak layer at η ref w = 10 19 Pa s and a strong layer at η ref s = 10 21 Pa s.

Parametric approach
In the first part of the paper, we explore the effect of the original dip of the layer on the shear bands that form in the detachment vicinity. We therefore describe in detail the results of twelve simulations in which the only varying parameter is the angle between the layering and the detachment zone. These models are named after their initial layering dip δ 0 , which has been varied taking values of 10, 20, 30, 40, 45 and 50 • in the direction of shearing and −10, −20, −30, −40, −45 and −50 • in the direction opposite to shearing. Experiments with higher angles were deliberately excluded from the present study due to numerical resolution issues.
In the second part of the paper, the model sensitivity to layering length scale, viscosity contrast or strain rate is explored. The results are presented using seven supplemental models in which layering dip is kept constant at 40 • . One run without the implementation of Mohr-Coulomb plasticity is presented. These extra runs are denominated by letters, run f being the reference model with δ 0 = 40 • . Detailed parameters are listed in Table 1.

Boundary conditions
To simulate a rigid hanging wall, shear is imposed by fixing a constant horizontal velocity along the top of the model, as in Gueydan et al. (2004). At the bottom of the models a constant velocity of opposite sign is applied in the horizontal direction but vertical velocity is not fixed. Instead, a normal stress corresponding to the confining pressure is imposed onto a thin elastic layer to allow the shear zone to thicken or thin vertically as it is sheared (Fig. 1b- keeping the model size constant is not straightforward. We therefore chose to perform the simulation on a 4 × 1 aspect ratio modelling box and use the results located only at the center of the model, far from the anomalous rotation effect at the lateral boundary conditions (Frehner et al., 2011). Similarly, we only use the upper half of the models close to the rigid boundary condition rather than the lower half, which is located close to the deformable elastic basement which has no geological meaning and was only introduced to impose confining pressure. The presented part of the models is represented by a grey box in Fig. 1c. Given these simplified but realistic boundary conditions, no analytical solution can be proposed for that problem, which may only be solved using analogue material or numerical modelling. The later solution permits to monitor both bulk and local stress as well as variations of these parameters with time. To do so, we solve the conservation of momentum equation for displacements using a slightly compressible formulation, which relates volume change to pressure via the media elastic compressibility (See Appendix A).

Numerical solution
The equations are solved on a 400 × 100 elements mesh so that each layer is 3 elements thick. The numerical code Flamar, based on Paravoz, uses an explicit time marching algorithm, which limits time step to half of the smallest Maxwell relaxation time in the model (Poliakov et al., 1993a). Although each increment of strain is small, Jaumann corotational corrections are applied to ensure the objectivity of the stress/strain update at each time step. The code is implemented in a fully Lagrangian formulation. Therefore, largestrain simulations require re-meshing to ensure the mesh quality is sufficient to obtain an accurate solution at each time step. However, re-meshing introduces two kinds of errors during the simulation: stress oscillations and numerical diffusion of material phases. Correction of the stress tensor components after remeshing remains diificult because they cannot be interpolated. However, the oscillations are easily recognized from the longer-term signal in the results. Throughout this study we made no attempt to interpret these oscillations as being associated with any physical process. Numerical diffusion is reduced using the re-meshing algorithm described in Yamato et al. (2007), which is based on the implementation of passive markers. However, we acknowledge that this methodology, together with the regular quadrilateral mesh, is not sufficient to accurately capture the initial rate of growth of viscous shear instability, such as folding and boudinage. Therefore, we do not attempt to do so in these studies and instead concentrate our interpretation on shear banding and the orientation of shear bands, which are well numerically resolved.

Representation of the results
Passive markers implemented by Yamato et al. (2007) and a geometric method derived from Ramsay and Huber (1983) are used to compute the finite strain field (See Appendix B and Fig. 2). This method has the advantage over the Frehner and Schmalholz (2006) method that it is not a punctal method but actually considers the deformation within a volume. This method avoids dropping some non-linear terms at large strain, it is general and it does not require the usage of interpolants of strain rate within the element.
Results of models are represented by superimposing lithologies in grey with the finite strain pattern in colour. Black and grey areas represent the weak and strong materials, respectively, and define the deformed initial layering. For a specific region of the domain (white box in Fig. 1c), we also represent the local maximum stretching direction (MSD, white lines in Fig. 2). Since models presented here are 2-D plane strain, the MSD is equivalent to the trace of the incipient schistosity in rocks deformed in ductile regimes.
The colour scale denotes the amount of finite strain, with blue and red indicating low-and large-strain domains, respectively. For all models, we also represent the total horizontal displacement along an initially vertical line as a function of depth beneath the detachment. This allows for qualitatively discussing the degree of localisation at the scale of the shear zone (large scale, Fig. 1d).

Fig. 3.
Results of the parametric study at γ = 1, accompanied by sketches of results. Black and grey indicate, respectively, the location of weak and strong layers; colour scale indicates strain intensity. Red circles refer to specific structures described in the text. P and T are the instantaneous shortening and stretching axis at large scale, respectively. Maximum large-scale shear directions are indicated in magenta.

Geometric and mechanical analysis of models
At γ =1, different large-and small-scale structures form. At large scale, a shear zone forms below the rigid top boundary. At small scale, shear bands, folds and necks in the strong layers are observed within this shear zone (Fig. 3).
Four different evolutions have been distinguished, primarily controlled by the initial layering dip δ 0 . Types I and II correspond to positive values of δ 0 and include models in which layers are initially stretched along layering, whereas Types III and IV correspond to negative values of δ 0 and include models in which layers are initially shortened (Fig. 3). For low obliquities (|δ 0 | < 30 • , Types I and III), shear strain primarily occurs in weak layers, and small-scale structures accommodating the strain are discontinuous. For large obliquities (|δ 0 | > 30 • , Types II and IV), small-scale structures are continuous and little strain occurs in weak layers.
In this section, we describe and analyze the structures and the strain partitioning on the one hand and the link between penetration of deformation, strain localisation and evolution of shear stress on the other hand.

Structures and strain partitioning revealed by finite strain field
In Type I models ( Fig. 4a), strain is accommodated at small scale by shearing along weak layers and by necking of strong layers (δ 0 = 10 • and 20 • ). In case δ 0 = 30 • , top-tothe-left shear bands crosscutting a few strong layers also form. Movement along these shear bands and sliding along weak layers define a conjugate network. Finite strain reflects a change in strain partitioning with increasing δ 0 . This partitioning is outlined both by the large amount of finite shear strain in weak layers ( Fig. 4a) and by the oblique orientation of MSD versus layering (Fig. 4b). This is particularly visible at δ 0 = 10 • , where MSD is systematically parallel to layering in strong layers (indicating a local large pure shear component) and oblique to layering in weak layers (indicating a local large simple shear component). For a larger initial dip, finite strain intensity is more homogeneous and the MSD orientation trends parallel to layering. In Type II models (Fig. 5a), discrete, localised and continuous top-to-the-left shear bands crosscutting layering develop. They induce the formation of small isolated sigmoidal  bodies of strong material. Dip of these shear bands depends on initial dip of the layering: shallowly dipping to the left for δ 0 = 40 • , horizontal for δ 0 = 45 • and shallowly dipping to the right for δ 0 = 50 • . Shear bands form at an angle approximately 45 • to initial layering. Within one shear band, finite deformation intensity is large and MSD tends to become parallel to shearing. Between two successive shear bands, finite strain intensity is reduced and the MSD shows a sigmoidal pattern consistent with the top-to-the-left kinematics. In Type III models (Fig. 6a), deformation is partitioned between shearing along weak layers and thickening of strong layers (δ 0 = −10 • and −20 • ). This partitioning is outlined by large finite strain intensity in weak layers and MSD roughly normal to strong layers. Shortening is also accommodated by small-scale folding with axial planes dipping against bulk sense of shear (i.e. to the right). Such instabilities develop sooner for large obliquities (δ 0 = −30 • ), indicating that growth rate of folds rises with increasing obliquity. Hinge zone of these folds being very localised, these can be described as kinks. Finally, shear bands do not form in Type III models. , 4, 135-152, 2013 www.solid-earth.net/4/135/2013/ In Type IV models (Fig. 7a), shearing is first accommodated by rotation of layering and thickening of layers within the shear zone. This rotation produces a large-scale fold, of which the upper limb corresponds to the sheared portion of rocks. Once the layers are rotated enough, the upper limb of the fold is stretched, after which shear bands form with characteristics similar to those of Type II models. In the three cases presented, at γ = 1, shear banding is only incipient for δ 0 = −40 • and δ 0 = −50 • . For δ 0 = −45 • , shear banding occurs along a horizontal direction only after γ = 1.3 (not shown). The incipient shear bands (dashed lines in Fig. 7a) dip to the right in δ 0 = −40 • model and slightly to the left in δ 0 = −50 • case. The dipping direction of shear bands obtained is opposite to those obtained for the same obliquity in Type II models with a roughly similar absolute dip value.

Solid Earth
At low obliquity (|δ 0 | < 30 • ), the weak layers are almost parallel to the upper boundary condition. Shearing along them is therefore favoured. In these cases, only small-scale discontinuous structures develop to accommodate layerparallel stretching (necking in the Type I models) or layer parallel shortening (folding in the Type III models). For larger obliquity, shearing along weak layers is not favoured. Continuous shear bands crosscutting the whole layering develop instead. If layering dips against bulk sense of shear, these shear bands form rapidly. In the opposite case, they form after the layering is overturned in the shear zone. Consistently, these shear bands form at δ 0 = 45 • to layering and have the same top-to-the-left kinematics as boundary conditions.

Penetration depth of deformation, strain localisation and evolution of the shear stress
All models show strain localisation directly below the upper limit of the model. The rigid upper boundary imposes the same strain rate in both the weak and strong layers close to it. This creates strong shear stress differences that trigger shear instabilities and bring stress states close to yield strength in strong layers. In contrast, far from the upper boundary, layers are free to deform at different strain rates, which results in a more diffuse strain pattern. The degree of localisation within the shear zone can evolve with time, as shown by the depth/displacement profiles (vertical black lines in Figs. 4a, 5a, 6a and 7a). It is closely related to the evolution of shear stress in strong layers, which can be used as a proxy for measuring the strength of the shear zone. If it increases with increasing bulk strain, the shear zone displays a strain hardening behaviour. In the opposite case, it displays a strain softening behaviour.
All the Type I models show very little softening with bulk strain. This evolution is consistent with the distributed pattern of finite strain in these models. In detail, the increasing obliquity also correlates with an increasing depth of penetration of deformation of strong layers: they are deformed in the whole investigated region for model δ 0 = 30 • while only their top 10 % to 20 % are deformed for δ 0 = 10 • (Fig. 4a). Depth of penetration of deformation is also correlated with the rise of average shear stress in strong layers as initial obliquity increases (Fig. 4c). Top-to-the-left displacement is distributed over all depths from the earliest deformation stages (curves in Figs. 4a and 4d for γ = 0.3) and the displacement rate at a given depth is constant with increasing bulk strain. Even if displacement rates may either be constant with depth (case δ 0 = 10 • ), resulting in a linear finite displacement profile at γ = 1, or concentrated at the top of the column (case δ 0 = 30 • ), resulting in a non-linear depth/displacement profile, there is no proper localisation of strain with time, i.e. no abandoned shear zones and no increase of convexity in displacement profiles with time. This steady-state behaviour is consistent with the absence of strain weakening deduced from the constant shear stress measured in strong layers.
In contrast, progressive strain localisation highlighted in the depth/displacement profiles from Type II models is compatible with the slight strain weakening observed in the evolution of the deviatoric stress with rising strain (Fig. 5c). This weakening is maximum for δ 0 = 45 • . At γ = 1, the penetration of deformation with depth is relatively similar in all models but the number of the shear bands and their spacing change with the dip of shear bands. Increasing convexity of depth/displacement profiles with increasing strain in all models illustrates this localisation process. For δ 0 = 45 • , profiles do not change in the lower half of the investigated area after γ bulk is higher than 80 %, meaning the structures at this level are not active anymore.
Within Type III, all depth/displacement profiles indicate ditribution of deformation from the earliest phase of deformation. Those profiles have a tendency to linearise with increasing bulk strain (especially δ 0 = −20 • ), which would imply that the shear rate increases downward with increasing bulk deformation to become perfectly distributed through the whole medium. This indicates that a natural shear zone forming in such a configuration would widen with increasing strain. This evolution is consistent with the strain hardening behaviour of the small-obliquity models (δ 0 = −10 • and δ 0 = −20 • ). Surprisingly, for an obliquity of δ 0 = −30 • a slight strain softening occurs (Fig. 6c). This stress decrease does not correspond to strain localisation but to the rotation of the soft layers towards directions that are more compatible with the bulk strain direction.
Shear bands form late in Type IV models. They appear as layering is rotated and are marked by inflections in the depth/displacement profiles. In these cases, younger shear bands are formed downward with migration of the layering fold hinge. Resulting shear band geometry is not as continuous as in Type II models. This lack of microstructures at δ 0 = −45 • prevents the initial structural softening that affects the strong layers in the two other cases (Schmalholz et al., 2005). As a result, strong layers remain stronger and rotate passively in a simple shear regime, as outlined by the   obliquity of MSD to layering in strong material. In weak layers, shearing is so intense that MSD is parallel to layering (Fig. 7b).
The numerical models show that the development and evolution through time and space of shear bands in sheared layered rocks is highly controlled by the initial orientation of layering regard to shearing direction. Shear bands develop and progressively localise strain in the vicinity of a highstrain zone only for layering orientations close to 45 • to shear sense. If the angle between shear and layering is low, shear is accommodated by conjugated activation of discontinuous shear bands and sliding in weak layers, resulting in boudi-nage of strong layers. Layering will passively overturn before being crosscut by shear bands if layering is initially dipping at high angles opposite to the sense of shear.

Pressure jumps
The Mohr-Coulomb yield criterion used in this study is pressure-dependent. To understand how the plastic deformation localises, the pressure pattern in the layers has to be  described first. The weak and strong layers in our models sustain very different pressure, which deviates by one fourth to one half from the confining pressure applied at the boundaries. The pressure jump amplitude depends on the effective strength of strong layers and the gradient sense depends on layering orientation with regard to shear direction with similar trends as Schmid (2005) found around rigid inclusions.
Each type of model is characterized by a pressure pattern and evolution. In Type I models (Fig. 4c), pressure in weak layers is slightly higher than the confining pressure (400 MPa) and significantly lower in strong layers (300 MPa). This results in a constant pressure jump of 100-150 MPa which equates approximately the amplitude of maximum shear stress in strong layers. Type II models (Fig. 5c) display similar characteristics between strong and weak layers except that pressure is sensitively higher than confining pressure in weak layers and that it tends to decrease with bulk strain with a maximum slope for an initial dip of 45 • . The pressure jump and maximum shear stress in strong layers are again in the same range.
In Type III models, pressure is very close to confining pressure in weak layers and 100 to 150 MPa higher in strong layers. This difference is also close to shear stress in strong layers. In the hardening cases (δ 0 =-10  pressure tends to increase in strong layers with increasing bulk strain (Fig. 6c) and layers do not deform significantly by folding. In the δ 0 = −30 • case, pressure decreases during folding with the same trend as the shear stress in the strong layer. This effective structural softening of the shear zone occurs as large-amplitude folding develops (Schmalholz et al., 2005). We do not discuss the rate of this case of instability in the present paper because Schmalholz and Schmid (2012) did so with a better numerical method and a setup that is very similar to δ 0 = −30 • . Type IV models show a reversal of pressure difference between strong and weak layers during the early deformation stages and then mimic a Type II case. This reversal corresponds to the onset of stretching in strong layers.
To provide a theoretical background to the impact of pressure and shear stress distribution and the localisation of shear bands in the models, it is usefull to think of the problem in terms of continuity of stress using Mohr circle constructions (Fig. 8). , 4, 135-152, 2013 www.solid-earth.net/4/135/2013/

Mohr Circle construction
The local stress state within strong layers is represented by a circle in the Mohr space. Strength of weak layers being two orders of magnitude lower than strength of strong layers, their Mohr circle has a very small radius and can be approximated by a point on the horizontal σ n axis (magenta square in Fig. 8). For reasons of stress continuity across the interface between strong and weak layers, shear stress and normal stress acting on the layer interface must be the same in both layers. This implies that, at the layer scale, the local layering direction (L) must belong to the circle that describes the stress state in the strong layers. Taking these two conditions into account, one of the principal stress axes within the strong layers must be aligned with the weak layer. For the Type I and II models, the layer initially lengthens and one expects σ 1 to be normal to the layers. For the Types III and IV, the layer initially shortens and one expects σ 3 to be normal to the layers. Under these conditions and looking at the circle construction, it is obvious that continuity of stress implies that the maximum shear stress in the strong layer must be equal to the gap of pressure between the strong and weak layer, and that the sign of the pressure jump defines if the layering (L) is perpendicular either to σ 1 or to σ 3 .
Accounting for this local orientation of the principal stress axis, we may draw the maximum shear-stress directions within the strong layers (lines a and b in Fig. 8) and predict the kinematics of the corresponding shear zones according to the local principal stress orientation.

Effect of Mohr-Coulomb plasticity
In the models presented here, the Mohr-Coulomb yield criterion was used. This choice can be discussed since it is pressure dependent while localised shear bands tend to form throughout the entire lithosphere in nature i.e. independent from pressure. This apparent contradiction is solved when we consider that the low permeability of the lower crust and upper mantle favours undrained conditions and that nonhydrostatic fluid pressure can limit the depth dependence of the frictional yield criterion (Sibson, 1994). Non-newtonian stress-dependent rheology or fixed yield-strength rheology like Peierls creep (Goetze and Evans, 1979) at depth could also allow deformation to localise faster than in viscous media (Schmalholz and Fletcher, 2011). These mechanisms are not explored in this study focused on shear band development at the brittle-ductile transition where Mohr-Coulomb plasticity is more appropriate to describe the brittle behaviour.
The description of the previous experiments showed that the localisation of shear bands only occurs when pressure in the strong layers is lower than the confining pressure, i.e. when these layers are being stretched (Types I and II) while they deform with only little strain localisation during initial shortening (Types III and IV). In Type IV, strain localisation occurs as layers start stretching near the rigid upper bound- ary. This asymmetry between stretching and shortening is due to the pressure dependent term of the Mohr-Coulomb yield criterion in the model. At the brittle-ductile transition, whether shear stress in strong layers exceeds the plastic yield criteria or not depends merely on whether strong layers are stretched or shortened along their dip, and on the sign of the pressure jump.
To confirm this hypothesis, we have taken the best localising experiment, i.e. δ 0 = 40 • , and ran a simulation which is identical except that plastic yielding was disabled. The results (model c, Fig. 9) lacks the formation of shear bands and the depth/displacement profile shows a low degree of localisation. Similarly, models in which the strain rate is lowered by a factor of 10 (model b, Fig. 9) or in which the viscosity of the strong layer is reduced by a factor of 10 (model g, Fig. 9) do not produce shear bands or localisation. In both cases, the maximum viscous shear stress of the strong layer is 30 MPa, hence one order of magnitude lower than the confining pressure. Oppositely, modifying the maximum viscous shear stress of the strong layer too largely exceeds the yielding criteria, does not significantly change the results (models d and e, Fig. 9). In all the cases, brittle-plastic deformation is a necessary condition for the localisation of shear bands for γ < 1.

Maximum shear-stress directions, and the orientation and occurrence of shear bands
The maximum shear-stress directions a and b are reported on the representation of shear strain at γ = 0.1. For Types I and II, strain starts to localise along shear bands parallel to the maximum shear-stress direction from the very early stages of the deformation onward. Type IV shows poorly localised strain at small scale for γ = 0.1, although the b direction dipping into the direction of shear is aligned with local concentration of shear strain (Fig. 8). Contrasting with Types I and II, this local concentration is not delineating a shear band. However, the b direction does correspond to the orientation of the shear band that will eventually develop at larger strain, as shown in Fig. 3 local shear strain concentration only accommodates the passive rotation of layers in the first stages. However, in the later stages the shear bands develop following the initial alignment of fold hinges. For the case of Type III, however, none of the a and b axes are activated in the early stages of the deformation and Fig. 3 outlines that no shear bands form within this orientation range at larger strain. In this case, the pressure is always higher than confining pressure in the strong layer. This configuration does not allow for plastic yielding to occur and therefore favours folding at least up to γ = 1.
Comparing with the preferred shear orientations deduced from Mohr circle analysis with the intensity of strain at γ = 1 (Fig. 3), it is clear that the shear band orientations, when they form, correspond only to one of the two conjugate local maximum shear stress orientations. The orientation chosen is always the one compatible with the larger scale shear direction.
The present numerical models with large-strain computation confirm that shear localisation directions are acquired at very small strain and therefore that the small strain formulation of Biot (1964) allows retrieving the orientation of large strain features (Cosgrove, 1976). As a result, although simple shear was applied at the boundary, the shear bands dipping towards the sense of shear display a normal kinematics and those with an opposite dipping sense display a reverse kinematics.
In conclusion, while boundary conditions control the overall kinematics, the direction of initial mechanical layering controls the orientation of shear bands forming in the media as well as the effective strength of the sheared media. Accounting for Mohr-Coulomb plasticity instead of Peierls creep (Schmalholz and Fletcher, 2011) Table 1. to occur, the models show that the layering passively rotates by folding towards more favourable orientations at the scale of the shear zone. Therefore, at larger strain, the layering always becomes favourably oriented for shear banding to occur. If layering is present at different scales in the medium, the rocks may go through cycles of hardening (folding) and weakening (shear banding), giving rise to anastomosing shear zone as was already pointed out by Ghosh and Sengupta (1987).

Scaling of the model and strain localisation close to natural detachments
Given the realistic rigid boundary conditions used in our study and the usage of elasticity and Mohr-Coulomb plasticity, as well as the locally quite large strain displayed, it is impossible to conduct a complete scaling analysis. However, we did perform some models to address the limit of applicability of our results at different length scales, background strain rates, and viscosity contrasts. It is beyond the scope of this study to include the effect of varying the thickness ratio between the weak and strong layers and we refer the reader to Schmid and Podladchikov (2006) for detailed insight into the effective behaviour of layered media as a function of spacing. However, within that frame, our analysis holds as long as the shear stress in the weak layer is negligible enough to ensure that the local principal stresses are aligned with the layer (cf. Sect. 4.2). This implies that if the thickness of the weak layers is decreased as compared to the thickness of the strong layer, their viscosity has to be dropped to obtain the same results. In this part, we concentrate on the model δ = 40 • because, as we have shown earlier, shear banding only occurs when the layering is being stretched, and when it does not occur, large amplitude folding tends to rotate part of the layering towards orientations that are favourable for shear banding.

Scale invariance
To test whether the results are independent of the length scale, we have run an experiment in which the dimensions were divided by a factor 100, yielding a total model dimension of 10 by 40 m and initial layer thickness of 30 cm (Fig. 9a). We have kept the viscosities and background strain rate similar as in the reference experiment (Fig. 9f). The results obtained are similar enough to conclude that our models are valid at least from the kilometers scale to the outcrop scales.
The experiment of Fig. 9h was produced by increasing the viscosity of both the weak and strong layers by a factor of 100 and decreasing the background strain rate by a factor of 100 compared to experiment i. The results obtained are very similar and lead to the conclusion that our results are applicable to a wide range of viscosities and strain rates as long as the viscous stress in the strong layer exceeds the plastic yield strength.
In all the models where shear bands were forming, the same alignement of large-strain structures is observed with a dip of approximately 10-15 • towards the left, which indicates that the dip of the shear bands is controlled by the dip of the initial layering rather than the scale, the strain rate, or the viscosity contrast. At constant strain rate, our models are length-scale independent. At large strain rate, shear bands are favoured but the strain is more distributed in the model. As the strain becomes more penetrative, flow pattern becomes more linear at large scale.

Bounds on viscosity
Our models also show that localised deformation crosscuting the initial foliation, i.e. shear bands, only occurs for a restricted range of initial orientation of the layering (Types II and IV). On the one hand, the occurrence shear bands requires that the strong layers yield plastically, implying that where η s , τ y andγ stand for the viscosity of the strong layers, the plastic yield strength and the background strain rate, respectively. On the other hand, large-scale strain localisation requires that the weak layer does not yield plastically, and hence where η w is the viscosity of the weak layer. Therefore, providing that both the plastic yield strength and the largescale kinematics are known, the microstructure and the observation of strain localisation at larger scale in turn provide bounds on the viscosity of the strong layers. Conversely, one could use the laboratory data on creep parameters and plastic yield strengths to estimate the large-scale strain rate.

Penetrative deformation versus strain localisation
The series of experiments framed in the red box in Fig. 9d, e, f and g aims to characterise the impact of the strong-layer viscosity on the outcome of the models. In the experiment in which the strong-layer viscosity is decreased by a factor 10 ( Fig. 9g), the displacement profile indicates that strain localisation is weak at the scale of the experiment and the finite strain does not feature shear bands. In that case the degree of localisation is small. For larger viscosities, shear bands always form. Their orientations and the spacing do not change significantly although there is one more shear band in the experiments d and e than in the reference model f. The displacement profiles outline a larger degree of localisation for viscosity contrast greater or equal to 1000, but the depth of penetration of the localised deformation is not affected by the viscosity of the strong layer. The series of experiments framed in the blue box in Fig. 9b, f and i) show that strain localisation at the scale of the shear zone is highly senstive to the background strain rate. At lower strain rate (Fig. 9b), shear banding does not occur and the degree of localisation is small, as attested by the straight depth/displacement profile. The results obtained for low strain rate are very similar to those obtained turning off plasticity in the model (Fig. 9c). At high strain rate, shear banding instability does not drop the strength of the shear zone sufficiently and the deformation becomes more penetrative. For these higher strain rates (Fig. 9i), the shear stress in the weak layers is on the order of 30 MPa and is therefore no longer negligible compared to the plastic yield strength of the strong layers. As a result, the effective viscosity contrast between the layers drops and much less strain occurs in the weak layers, probably causing the lack of localisation at large scale. However, the orientation of the shear bands is not affected by strain rate.

Application to strain localisation near detachment faults
Under geological conditions, the maximum shear velocity is limited. Strain rates are therefore more sensitive to the thickness of the shear zone than to the change in velocity. Comparing the results of low, moderate and high strain rate experiments is interesting when trying to upscale the rheology of the shear zone and considering a plate boundary velocity of ∼1 cm yr −1 . The low strain rate experiment corresponds to a 10 km thick shear zone, the moderate strain rate experiment corresponds to a 1 km thick shear zone and the high strain rate experiments corresponds to a 100 m thick shear zone. In the first case, the deformation distributes itself within the model with a low degree of localisation; in the second case, the degree of localisation is maximum and the depth/displacement profile shows that at γ = 1 most of the deformation occurs in 150-200 m layer located close to the hanging wall along 20 m thick shear bands. In the case with higher strain rate ( Fig. 9h and i), the deformation pattern obtained in the high strain rate actually also resembles more the pattern observed in mylonites and cataclasites, and although plastic yielding  Fig. 10. Sketch model for shear localisation beneath a detachment fault. Layered rocks are initially at high temperature and large confining pressure corresponding to model c. The rapid cooling during exhumation causes the strong layers to start yielding. According to model f, this results in shear banding and rapid localisation of strain (γ < 1). Within the shear bands, the strain rate is approximately increased by a factor of 10, corresponding to model i and resulting in more penetrative deformation and onset of mylonitisation. occurs and shear bands form, the depth/displacement profile shows the lowest degree of localisation of all the models and is indeed more linear than in the purely viscous model (Fig. 9c).
In our numerical experiments, where the finite strain is relatively small (γ = 1), strain localisation only occurs in a specific regime in which brittle layers coexist with effectively viscous ones. We do not dismiss that strain localisation may occur without plasticity at larger γ because the depth/displacement profile of the low strain rate and viscous experiments are not perfectly linear. There is indeed a tendency towards strain localisation but the rate of strain localisation is much smaller than in our models. At higher strain rate, the strain is very penetrative and does not lead to localisation because all the layers are yielding and the shear zone behaves at constant shear stress, i.e. with an effective perfectly plastic behaviour.
Indeed, when warm layered material is exhumed from depth along a detachment fault, the deformation first localises at small rate due to shear instabilities such as folding and boudinage. However, once the material has reached a smaller confining pressure and a lower temperature, the effective rheology of the layered media becomes extremely non-linear due to the onset of plastic yielding in the stronger layers that are kinematically stretched by shearing. As a result, the effective shear stress drops in the shear zone and this mechanical instability results in a rapid strain localisation. In response to strain localisation, the strain rate increases, and the deformation becomes more penetrative. Penetrative deformation being effectively linear, it cannot cause further strain localisation, stabilising the system again (Fig. 10).

Conclusions
We have run numerical experiments to simulate the deformation of a mechanically stratified medium in simple shear to provide quantitative insights on the orientation of shear bands forming beneath a detachment fault as a function of the initial orientation of the layers relative to sense of shear imposed by the rigid hanging wall of the detachment. The simulations were separated into four groups according to the absolute value of the initial dip, δ 0 , and its sign. The absolute value of δ 0 controls strain partitioning between the weak and strong layers and its sign controls the effective shortening or stretching of the layers as they rotate within the shear zone. For Type I models, the resulting pattern is relatively similar to the extensional crenulation cleavage described by Platt and Vissers (1980). In Type II and IV cases, the chosen orientation corresponds to one of a neoformed shear band with S-C structures. Type III models do not lead to shear banding.
Near the brittle-ductile transition, pressure jumps influence strain localisation in rocks by slowing it down and favouring folding over shear banding for orientation of Types III and IV or by leading to rapid localisation of the deformation for Types I and II. To our knowlegde the deformation and the mechanical layering close to detachment faults mostly fall within Type I and II geometries and constitutes an argument in favour of the existence and the importance of dynamic pressure in deformation processes at the brittle-ductile transition. This local pressure jump contributes to enhance chemical potential differences and increases mineral solubilities favouring local pressure-solution (Fletcher, 1977) over advective mass transfer, which requires seemingly unrealistic quantities of external water (Ferry and Dipple, 1991) to forms veins. Following this hypothesis, our model would thus predict that elements migrate toward the strong layers in Type I and II regimes and that quartz veins may form in the weaker layers in Type III and in the early deformation stages for Type IV orientations. In both cases, the synkinematic paragenesis will form in the part of the model where pressure is nearly lithostatic, providing an explanation for the scarcity of petrologically proven examples of tectonic over/under pressure in rocks (Vrijmoed et al., 2009).
It was already shown in Le Pourhiet et al. (2004) that at crustal scale, dipping mechanical layering in presence of a thermal gradient alone can cause the localisation of normal faults in the upper crust at the location where the weakest layer crosses the brittle ductile transition. Huet et al. (2011) have shown that these results could be upscaled and that mechanical layering due to nappe stacking facilitates the formation of metamorphic core complexes exhumed along detachment faults. Here, we have explored the effect of mechanical layering at smaller scale and we show that it causes www.solid-earth.net/4/135/2013/ Solid Earth, 4, 135-152, 2013 This direction is represented with a white line on the details snapshot of the models. The intensity of stretching in that direction is and the intensity of maximum shortening (normal to the direction of θ ′ ) is with = a 2 − b 2 + c 2 − d 2 2 + 4 (ab + cd) 2 . The ellipticity is a measure of the intensity of the finite strain, which varies between 1 (no strain) and +∞ (infinite strain). For visualisation purposes, we prefer to use which renormalises the finite strain on a hyperbolic scale that varies between 0 (no strain, blue in the figures) and 1 (infinite strain, red in the figures).