Pinch and swell structures: evidence for strain localisation by brittle–viscous behaviour in the middle crust

The flow properties of middle crustal rocks are commonly represented by viscous flow. Examples of pinch and swell structures found in a high strain zone at St. Anne Point (Fiordland, New Zealand) and Wongwibinda (N.S.W., Australia) suggest pinch and swell structures may be initiated by brittle failure of the more competent layer in conjunction with subsequent material softening. On this basis we develop a numerical model where Mohr–Coulomb constitutive strain localising behaviour is utilised to initiate pinch and swell structure development. Results show that pinch and swell structures develop in a competent layer in both Newtonian and non-Newtonian flow, provided the competent layer has sufficient viscosity contrast and can localise strain to form shear bands. The flow regime and strain localising characteristics of the surrounding country rock appear not to impact pinch and swell structure formation. The degree of material softening after the initial strain localising behaviour is shown to impact pinch and swell characteristics, while extensive material softening causes the formation of thick necks between swells by limiting the focused localisation of strain into shear bands. To aid analysis of the structures and help derive the flow properties of rocks in the field, we define three stages of pinch and swell development and offer suggestions for measurements to be made in the field. Our study suggests that Mohr–Coulomb strain localising behaviour combined with viscous flow is a viable alternative representation of the heterogeneous rheological behaviour of rocks seen in the middle crust. This type of mid-crustal rheological behaviour can have significant influence on the localisation of strain at all scales. For example, inclusion of Mohr– Coulomb strain localising behaviour with viscous flow in just some mid-crustal layers within a crustal-scale model can result in significant strain localisation, extending from the upper crust into the middle crust. This localisation also influences the development of near-surface structures.


Introduction
Rock deformation in the Earth's crust is known to be highly heterogeneous (Rudnick and Gao, 2003;Taylor and McLennan, 1985) with localised domains of higher strain seen in metamorphic zones from the upper to lower continental crust.Crustal-scale numerical models of tectonic processes using rheological flow properties have provided important insights into scenarios such as continental break-up and collision, and their surface expression in the current landscape.While there is ongoing discussion of appropriate rheological models for the upper crust (e.g.Moresi and Mühlhaus, 2006), it is generally accepted that the most suitable rheological model for the middle to lower crust is thermally activated viscous flow (e.g.Molnar and Houseman, 2004).However, the assumption of viscous flow for the middle continental crust is now increasingly being questioned (Brantut et al., 2013;Bürgmann and Dresen, 2008;Karrech et al., 2011;Mancktelow, 2006Mancktelow, , 2008Mancktelow, , 2009;;Mancktelow and Pennacchioni, 2005;Pennacchioni and Cesare, 1997;Pennacchioni and Mancktelow, 2007;Regenauer-Lieb and Yuen, 2008).
For decades geologists have used field evidence of different structures to derive the flow characteristics of crustal rocks.For example, Treagus and Treagus (2002) used deformed pebbles, Ramsay (1980) used shear zones, whereas boudins and mullions have been used in several studies (e.g.Kenis et al., 2004;Urai et al., 2001).Boudins (from French R. L. Gardner et al.: Strain localisation in pinch and swell structures boudin, "sausage") are chains of rock fragments sometimes joined by thin necks giving them the appearance of a string of sausages (Lohest et al., 1908;Ramsay, 1866).They form when a more competent layer undergoes layer parallel extension in a weaker matrix, and can range in size from microscopic to outcrop to regional scales (Goldstein, 1988;Klepeis et al., 1999;St-Onge et al., 2009).Boudins have been used to suggest that rock units are deformed by a mixture of brittle and/or semi-brittle to ductile flow (Abe and Urai, 2012;Arslan et al., 2008;Komoróczi et al., 2013;Marques et al., 2012;van der Molen, 1985).Pinch and swell structures, also called drawn or necked boudins, are a subset of boudins that retain continuity of the drawn layers but show variable thinning of the original layer thickness.They have been reported from a large variety of tectonic settings and crustal levels (Arslan et al., 2008;Goscombe et al., 2004;Hobbs et al., 2010;Kenis et al., 2004;Ponce et al., 2013;van der Molen, 1985).
Investigations into the necessary flow properties that allow pinch and swell formation have encompassed three general areas: (i) theoretical analysis (Emerman and Turcotte, 1984;Smith, 1975Smith, , 1977) ) using pure viscous flow, which suggests that non-Newtonian flow and material softening is required for pinch and swell formation; (ii) analogue modelling, which generally confirmed the earlier theoretical concepts (Kobberger and Zulauf, 1995;Mengong and Zulauf, 2006;Neurath and Smith, 1982;Zulauf and Zulauf, 2005); and (iii) numerical modelling, which has commonly implemented solutions where a geometric instability initiates continuous necking of a power law material undergoing viscous flow based on the linear stability analysis of Fletcher (1974).This has been used to investigate folding scenarios (e.g.Llorens et al., 2013;Schmid et al., 2004) and has also been extended to pinch and swell structures, where an initial irregularity on the competent layer to matrix interface (here, referred to as the competent layer edge) or within the layer is allowed to grow and dominate any other irregularities to allow strain to localise, thereby causing the pinch and swell structure to initiate (e.g.Passchier et al., 2005;Schmalholz and Maeder, 2012;Schmalholz et al., 2008;Smith, 1977).Abe and Urai (2012) and Komoróczi et al. (2013) create boudins using a model where initial brittle deformation caused by simulated high fluid pressure is used to create 90 • fractures (mode I; Koehn et al., 2005;Pollard and Segall, 1987;Scholz, 2002), which initiate the boudin structure, followed by viscous flow to allow the structures to develop.This model readily creates torn boudins, a subset of boudins where the separated ends of the competent layer fragments are rectangular.
The numerical scheme presented here aims to represent an appropriate alternative to the commonly used viscous flow model of middle crustal rocks.We use a rheological model combining Mohr-Coulomb strain localisation and subsequent viscous flow allowing simulation of the described wide range of strain localisation and subsequent softening behaviours observed in nature.Once initiated through strain localisation, pinch and swell structure evolution is influenced by the flow properties of all the layers (e.g.Schmalholz and Maeder, 2012;Smith, 1977).
Viscous flow is assumed for post-failure flow and in its simplest form is described by where ˙ is the strain rate, σ is differential stress, n is the stress exponent and A is a constant which reflects the material properties at a specific temperature and pressure.We utilise field examples of several well-exposed chains of pinch and swell structures to determine the prerequisites for a numerical model that allows realistic representation of the initiation and developing characteristics of the structures observed.
We systematically test the impact of (I) the stress exponent and Mohr-Coulomb strain localising behaviour in the competent layer, (II) initial viscosity ratio of the competent layer to matrix, (III) competent and surrounding layerdependent stress exponent and Mohr-Coulomb strain localising behaviour, and (IV) material softening following Mohr-Coulomb strain localisation, on the formation and evolution of pinch and swell structures.We use our models to derive a scheme of simple field measurements on a pinch and swell outcrop which can provide insights into relative rheological properties of the competent layer and the surrounding layers.
The model is then extended to a crustal-scale model to test the effect of the presence of strain localising followed by viscous behaviour in a subset of the middle crust layers.In these exploratory tests of crustal extension, we investigate the strain localisation throughout the model and the characteristics of the developing surface structures.
Our results highlight the importance of including rock units with contrasting rheological properties in terms of strain localising behaviour, material softening and viscosity differences in numerical models on the thin-section (mm) to map-scale (km), as all these factors impact the dynamics of the developing structures.

Pinch and swell structures: field and thin section observations with initial interpretation
Brittle deformation resulting in the formation of pinch and swell structures is readily visible in the plagioclase dominated pegmatites found in the Wongwibinda Metamorphic Complex, New England Orogen, N.S.W., Australia (30 • 18 41.6 S, 152 • 8 35.3 E).On the outcrop scale, large plagioclase grains up to 3 cm have brittle fractures of variable angles of 30 to 55 • from the principal stress direction (Fig. 1a).On the thin-section scale, brittle failure can also be seen as distinct fractures and zones of fine-grained recrystallised grains, with a dominance of angles of 45 to 55 • to the external foliation (Fig. 1b).The Wongwibinda pinch and swell structures show little impact from fluid influx and are formed at temperatures of ∼ 600 • C and shallow crustal depths of ≤ 200 MPa (Craven et al., 2013).
The rock platform at St. Anne Point, Fiordland, New Zealand (44 • 34 19.9 S, 167 • 46 55.3 E) displays an extensive variety of mafic pinch and swell structures, many in long chains with irregular spacing developed in a metasedimentary country rock (Fig. 2b).Parallel, 2 mm to 50 cm thick mafic layers run SW-NE across the rock platform, generally parallel to the overall steeply dipping foliation of the surrounding metasedimentary layers.Different mafic layers exhibit different structures (Fig. 2).Some layers form welldefined pinch and swell chains (Fig. 2b), while others display only minor undulation of the mafic to felsic edge.Very narrow (Fig. 2a) to wider (Fig. 2b) neck thicknesses as well as short to long swell separations develop (Fig. 2a and b) in different layers.Numerous examples can be seen of apparently brittle failure within the mafic layers, where pinch and swell structures typically have residual evidence of brittle failures at 30 to 40 • to the foliation (Fig. 2c, thick dashed line) and general layer orientation.Within the swell structures, the internal foliation has been variably deflected, deviating up to 15 • from the orientation of the general matrix foliation (Fig. 2c).
Sample pairs were taken from the swell centre and neck of the mafic competent layer in multiple pinch and swell chains.A summary of microstructural features is provided here (see Supplement for detailed thin section descriptions and a description of the surrounding metasedimentary layers with additional photomicrographs).The main difference between swell centre and neck samples is the relative abundance of phases.The swell centre samples have more amphibole + garnet but less biotite + quartz + plagioclase compared to neck samples.In all samples, amphibole grains are fractured.However, fractures are more common in neck samples.In all cases, the fractures are filled with quartz, biotite and/or plagioclase (Fig. 2d, e; white arrows).Strain shadows around garnet grains contain the same minerals and indicate brittle failure and ductile fabrics formed under similar amphibolite facies conditions (Klepeis et al., 1999).In addition, neck samples have smaller grain size and exhibit welldeveloped shear bands.
Brittle fractures in the centre of the swell (Fig. 2d) are at right angles to the general foliation and are not associated with any shear bands.By contrast, similar to outcropscale brittle features, shear bands and brittle deformation in neck samples are at 30-40 • angles to the general foliation (Fig. 2e).The spatial association of increased brittle deformation, the presence of hydrous minerals and shear band formation at the neck suggests that the localisation of strain in the mafic layer occurred where the necks have formed and is associated with mode II failure, fluid influx and a transition from brittle to ductile deformation, in agreement with Brander et al. (2012), Fusseis et al. (2006) and Mancktelow (2006).Fluid influx was aided by local porosity and permeability increase due to the brittle failure.The presence of numerous grains displaying brittle failure in the neck suggests that once initiated, these narrow, intensely sheared regions could develop into a through-going shear band (Fig. 2e).These shear bands, in turn, are inferred to potentially develop into outcrop-scale failures in the pinch and swell chain, where enough strain is localised into a neck.For the example shown in Fig. 2e, not enough strain was localised to form the neck.The angles of the residual brittle failures at 30 to 40 • to the foliation suggests these may have rotated as part of the continuous development of the pinch and swell structure.
The observed assemblage as well as the steeply dipping, SW-NE trending foliation indicate that the deformation occurred during the regional D 3 dextral transpression event (Klepeis et al., 1999).This event occurred at upper greenschist to amphibolite facies conditions at ∼ 8.7 kbar and ∼ 585 • C (Klepeis et al., 1999).These conditions are consistent with amphibole and garnet undergoing brittle deformation (Passchier and Trouw, 2005), while quartz and plagioclase deform in a generally ductile manner (Passchier and Trouw, 2005;Rybacki and Dresen, 2004).Quartz and biotite (Fig. 2e) did not deform brittlely, but instead underwent dislocation creep, as evidenced by the presence of undulose extinction, subgrains and irregular boundaries (Stipp et al., 2002).Hence, these phases would have deformed by non-Newtonian flow (e.g.Ranalli, 1997).

Conceptual model based on field analysis
Ideally, in any model, the pinch and swell structures should initialise and develop with minimum manual specification of the initial location and orientation where strain concentrates to form the necks.The Wongwibinda and St. Anne Point samples both suggest that Mohr-Coulomb mode II brittle failure (Koehn et al., 2005;Pollard and Segall, 1987;Scholz, 2002) is characteristic of the initiation mechanism producing the failure planes seen in the field examples.
Field examples also suggest that material softening occurs after the initial brittle failure and this softening process needs to be included in the model.This material softening should incorporate both lower viscosity, simulating a higher percentage of "soft" phases and lower cohesion, simulating a higher degree of fracturing (grey dashed lines, Fig. 3a).The St. Anne Point samples display non-Newtonian flow, but, Newtonian flow, that is, grain-size-sensitive flow is commonly in-ferred for natural rocks (Svahnberg and Piazolo, 2010), so both Newtonian and non-Newtonian flow should be tested in the model.
Based on these observations, conceptually, a numerical model to describe the growth of pinch and swell structures should encompass: (1) no specified geometric perturbation (e.g.sine waves or notches) in the competent layer; (2) the ability to deform by brittle failure, resulting in pairs of fault planes that occur dynamically; (3) the ability to soften the material to variable degrees after brittle failure has occurred and (4) the ability to implement both Newtonian and non-Newtonian viscous flow.
Using the numerical model, we investigate the effect of the rheological properties of stress exponent, strain localising behaviour, initial viscosity ratio and material softening on the formation of pinch and swell structures.

General model set-up and initial testing
We have used the numerical modelling platform Underworld (Moresi et al., 2003(Moresi et al., , 2007)), a 2-D and 3-D simulation package that uses Lagrangian particle-in-cell finite element methods for modelling large-scale geodynamic processes on geologic timescales to investigate the formation and subsequent evolution of pinch and swell structures.Theoretically, as stress is increased, the layers should at first undergo purely elastic deformation at very low strain (cross-hatched and dotted areas in Fig. 3c), followed by viscous deformation (Griggs and Handin, 1960).While elastic behaviour is only important at very low strain (< 1.03 stretch, e.g.Ranalli, 1997), the model has been simplified by assuming a combi-nation of localising and viscous behaviour without explicit elastic behaviour, as we are interested in behaviour postyielding.Thus, stress accumulation, elasticity and/or elastically stored energy are ignored in the model (Fig. 3c).
Strain localising behaviour was included to simulate the brittle failure inferred for the St Anne Point samples.We have used the Underworld implementation of Mohr-Coulomb constitutive behaviour by Moresi and Mühlhaus (2006) (discussed in Ord, 1991;Rudnicki and Rice, 1975;Vermeer and De Borst, 1984), which uses an intrinsically unstable incompressible material with no dilation and no elasticity.Other material models, for example, Drucker-Prager or von Mises materials could have been used, but we considered the Mohr-Coulomb constitutive behaviour, which readily localises strain at unspecified locations, to best represent the brittle failure inferred for the Wongwibinda and St Anne Point pinch and swell chains.
The Moresi and Mühlhaus (2006) implementation also optionally incorporates material softening for the particles that have undergone yielding.This provides the capability of a layer localising strain by plastic processes, possibly including brittle fracturing, then flowing in a viscous manner.The material fails immediately on initiation of the model run, and stress gradually decreases as the run continues (Fig. S2 in the Supplement).Failure is critically dependent on the initial cohesion of the competent layer (Fig. 3a; c i ).In these models, it has been set just below the maximum stress in the system.It is also necessary to define the cohesion of the layer after failure and softening (c s ).We define a cohesion ratio (R Co ): (2) R Co values are dimensionless and thus provide relative rather than absolute information (Fig. 3a).This ratio is kept constant in the models except where the impact of c s is specifically tested in analysis IV.
We implement 2-D models with two different general geometries: (i) a three-layer model which is the simplest possible geometry, where the central layer (Fig. 3b, layer B) is more competent, and (ii) a multilayer model with a brittle upper layer (layer C), i.e. upper crust and variably competent layers in the lower parts (layer B'), i.e. middle crust of the model (Fig. 3d).Passive black marker lines are included in each layer to show the movement of pre-existing planar features such as bedding or foliation as the simulation proceeds (Fig. 3b and d).All simulations are based on an incompressible Uzawa Stokes Flow calculation, starting with an initial 1 × 1 square that is extended horizontally using a stretch of −0.5 % on the left-hand side and +0.5 % on the right-hand side at each step to simulate extension.Top and bottom boundaries are non-periodic with free slip.In all experiments, the models were run over 200 iterations or steps, to a stretch (s) of 3.2, with each step increasing the stretch as calculated from the elongation of the original 1 × 1 square based on where L is the length at that step and L i is the initial length (Twiss and Moores, 1992).Non-Newtonian flow simulations use a stress exponent of three (n = 3, Eq. 1) as we consider this to closely represent the flow of naturally occurring rocks (Table S1 in the Supplement).For the three-layer models, we use a 1 × 1 unit square, allowing for comparison of relative properties between the layers and between models, while for the kilometre-scale models, we use SI units.Details of parameter values used in the three-layer model are provided in Table 1.Healing was not included in the models.Gravity has been included for the kilometre-scale models, but omitted in the outcrop-scale models.
Each layer material has characteristics that combine to give the material constant (A, see Eq. 1) for the viscous flow behaviour in that layer.The material constant and viscosity are interrelated, with viscosity known to be dependent on grain size, temperature and pressure, activation energy and volume, and the presence or absence of melt and/or water (Hirth and Kohlstedt, 2003).In our models, these variables have been incorporated into a single viscosity value (η); layer B is more competent than the surrounding layers A (Fig. 3c).An initial viscosity ratio (R V ) is calculated to be where η Bi is the initial viscosity of the competent layer B and η Ai is the initial viscosity of the matrix layers A. The viscosities for both layers will change, depending on the strain rate as the non-Newtonian flow simulations progress.By contrast, in the linear (n = 1) Newtonian flow simulations, the viscosities and R V are constant throughout the simulations (per Eq. 1).
In our model, the angle of failure on the particle scale (from Moresi and Mühlhaus, 2006;Eq. 14) is determined by using θ = ±1/2 arctan(1/tanµ). (5) A history of angle failure is recorded on each particle in the model and is used to determine the location and angle of the shear bands formed.The macroscopic shear bands are an emergent feature that take a number of iterations to form (for details see Moresi and Mühlhaus, 2006).Typically, the shear bands have formed by a 4 % stretch.All of our models use a friction coefficient of 0.6, as this is recommended for crustal rocks by a number of researchers (e.g.Byerlee, 1978;Ord, 1991).Mesh tests undertaken on the reference model (cf.Figs.5a,  S3) show that the tortuosity and width ratio values used in this paper are generally mesh-independent, while the numbers and spacings of swells are mesh-dependent.In the following, we do not use the mesh-dependent parameters (number and spacing of swells), other than to indicate the spacing is irregular and based on where strain localises.

Model analysis
Numerical results were extracted from the model at regular intervals throughout the simulations.The results are plotted as spatial distributions of different properties, for example, different particle layers, strain rate, stress and viscosity.Results are extracted as stretch increases and are used to represent the model evolution.To monitor the behaviour of the system on a local scale, "tracer" particles on a 9 × 5 grid across the model were marked and their evolution was followed throughout each numerical simulation.Characteristics were collected at each step and have been used to produce the graphs presented.In addition, model results were analysed using canny edge detection methods in the Image Processing Toolbox TM in Matlab ® to determine swell and neck widths and tortuosity.For the three-layer models, average tortuosity is the average of the ratios of the two competent layer edge lengths to the straight line length.An average tortuosity of > 1.1 was taken to indicate probable formation of pinch and swell structures.An additional ratio of the minimum neck width (W N ) to the maximum swell width (W S ) was measured using An R W value < 0.4 was taken to indicate successful pinch and swell initiation, a value < 0.2 indicates successful formation of the pinch and swell structures and a value < 0.05 indicates attenuation of the neck.For the landscape models, tortuosity of the surface (T S ) and R W for the upper crust layer were measured.

Three-layer model
In order to evaluate the model behaviour and investigate the impact of material property variations of the surrounding country rock and competent layers on pinch and swell structure characteristics, we systematically test their development and evolution.

Analysis I: Impact of strain localising behaviour and stress exponent
To test if pinch and swell structures can be initiated by Mohr-Coulomb strain localising behaviour and form irrespective of the stress exponent, simulations were completed where layer B (Fig. 3b) had Mohr-Coulomb strain localising behaviour, and all layers had either Newtonian flow (n = 1) or non-Newtonian flow (n = 3).
Strain softening was included and R V was set at 20 for this analysis.Other parameters in the models were as specified in Table 1.R W was used to indicate successful pinch and swell initiation (R W < 0.4) and formation (R W < 0.2).

Analysis II: Impact of relative initial viscosity
To test the impact of the initial viscosity ratio (R V ) of the competent layer to the surrounding layers on the formation of pinch and swell structures, we undertook additional simulations in non-Newtonian flow.R V was increased until pinch and swell structures were successfully formed.Variation of initial viscosity ratios up to 2 orders of magnitude have been commonly used in numerical modelling (e.g.Llorens et al., 2013;Passchier et al., 2005;Takeda and Griera, 2006), so a wide range of R V values, from 10 to 160 for n = 3 were included in this analysis.Again, R was used to indicate successful vs. unsuccessful pinch and swell initiation and formation.

Analysis III: Impact of layer-dependent flow properties
To test the impact of Mohr-Coulomb strain localising behaviour and flow regime across all layers in the model, four separate series of simulations with R V set at 10 were undertaken (Table 2).We systematically varied the stress exponent and Mohr-Coulomb strain localising behaviour: (1) all layers were given Newtonian flow; (2) matrix layers (A) were given Newtonian flow while layer B was given non-Newtonian flow (n = 3); (3) layer B was given Newtonian flow while the matrix layers (A) were given non-Newtonian flow (n = 3); (4) all layers were given non-Newtonian flow (n = 3).Mohr-Coulomb strain localising behaviour was varied within each of these series (

Analysis IV: Impact of degree of material softening
To test the impact of material softening on the evolution of pinch and swell structures a series of models with Newtonian flow (n = 1), an R V of 20 and Mohr-Coulomb strain localising behaviour in only layer B was used.Material softening was varied using the initial cohesion (c i ) from the previous analyses.R Co values (Eq.2) ranged from 40, where there is a high degree of material softening on those particles that yielded, through 20, 10, 4, and 2, to 1, where there is no softening.

Multilayer model
To explore the impact of including rock layers with strain localising behaviour in the middle crust, on the overall pattern of strain localisation and near-surface structural development, a multilayer model was developed to represent a crustal-scale extensional scenario.In this model, a compressible layer at the top (Fig. 3d layer D) allows topography to develop (Hetényi et al., 2011) and a layer (Fig. 3d layer C) with Mohr-Coulomb strain localising behaviour represents the upper continental crust reaching the surface.Five layers (Fig. 2d, A' and B' layers) of alternating material properties, similar to our three-layer model, represent the middle continental crust.The relative widths and viscosity ratio (R V = 20) of these A' and B' layers are the same as for the three-layer reference model.We present results from three sets of simulation pairs at 10 × 10 km, 20 × 20 km and 40 × 40 km, where the rheological properties of the lower five layers (i.e.A' and B' layers) are varied.Within each set, (i) layers A' and B' have purely viscous rheological properties with no strain localising behaviour, and (ii) the B' layers were given Mohr-Coulomb strain localising behaviour.The model was run with non-Newtonian (n = 3) flow and gravity has been included.

General model behaviour
Flow property testing of the Underworld stress exponent and yielding parameters on a single material 1 × 1 model ensured flow behaviour according to theoretical predictions (Fig. S2).
Particle layer model plots (Fig. 4a) show the development of the pinch and swell structures and the rotation of the planar bedding features.The neck and swell development can be seen most clearly in the plot of the second invariant of strain rate (Fig. 4b).Neck and swell development is less obvious in the corresponding stress invariant plot (Fig. 4c) as there is less variation of stress across the whole model.Therefore, strain rate invariant plots have been used in the subsequent figures.The viscosity, strain rate and differential stress values were plotted (Fig. 4d, e and f) for selected particles within the competent layer (Fig. 4a, b and c: red and green stars from the swell centre and black star from the neck).An inverse relationship between strain rate and stress (per Eq. 1) can be seen in Fig. 4b and c.Where strain rate is higher (Fig. 4b dark bands), stresses are relatively low (Fig. 4c light areas) and the planar features (Fig. 4a) become distorted.The black parti-  cle has relatively higher strain rate and lower stress as it is in an area that subsequently develops into a neck.By contrast, the red and green particles are in relatively low strain, higher stress areas that develop into swells.The example shown has a stress exponent of three (n = 3), so, per Eq. ( 1), there is an inverse relationship between viscosity and strain rate (Fig. 4d  and e).Graphs show the model is working as theoretically expected for a yielding material (compare Fig. 4f with Fig. 3c).

Analysis I: Impact of strain localising behaviour and viscous flow characteristics
The progressive development of the pinch and swell structure is shown in Fig. 5 at R V = 20 for both Newtonian (Figs. 5a, 6a dashed black line) and non-Newtonian (Figs.5b, 6a solid black line) flow.The positions for the formation of the necks are defined by the strain localisation surfaces initiated at low stretch value, early in the simulations.These failure surfaces, where the strain localises to form local shear bands at ∼ 30 • to the principal stress direction, can be seen in the strain rate invariant results at each stretch value (Fig. 5).As the stretch increases (Fig. 5), all layers are elongated and thinned.The necks thin more quickly than the swells, causing the marker beds to distort near the localisation of strain in the neck (Figs.4a, 5).This analysis allows three stages of pinch and swell structure to be distinguished (Fig. 6a) as the structures develop with progressive deformation, that is, with decreasing R W and increasing stretch: 1. Stage 1 is where zones of localised strain are establishing and the structures have not started to form: 0.4 < R W < 0.75 (e.g.Fig. 5a stretch 1.2 to 1.6; Fig. 6a).
2. Stage 2 is where the strain has localised and pinch and swell structures have successfully been initiated and are starting to separate: 0.2 < R W < 0.4 (e.g.Fig. 5a stretch 1.6 to 2.5; Fig. 6a).
3. Stage 3 is where the pinch and swell structures have successfully formed and are well-defined and separated: R W < 0.2 (e.g.Fig. 5a stretch 2.5 to 3.2; Fig. 6a).Attenuation of the neck is developed where R W < 0.05.
The pinch and swell structures are more obvious at higher stretch values (Figs.5a and 6a, n = 1 dashed black line) where the necks have become attenuated into thin threads and R W < 0.4 (Fig. 6a).The bands of localisation have initial angles to the vertical (that is to the model principal shortening direction), of approximately ±28 • , with some variation discernable in the strain localisation plots (for example, Fig. 5a stretch 1.6, strain rate invariant plot).As the model is elongated, the shear bands rotate towards the horizontal (Fig. 6b,  dashed line).Shear bands displaying strongest localisation have the largest degree of rotation.Some swell structures display rotation of the passive marker beds as the strain localisation is established (Figs. 5  and 6b, to stretch 1.6).Where the strain localisation is regularly distributed (± ∼ 30 • to vertical) this results in swells with the passive beds rotating up to ±5 • (Fig. 5b).At higher stretch values, these become curved parallel to the edge of the swell to matrix interface.The strain localisation not being regularly distributed either end of the swell, that is, one direction dominating (e.g.Fig. 6b swell 1 and 2), results in swells where the passive beds rotate up to ±15 • .
Figure 5 shows that the localisation of strain on just a few failure planes is more effective in Newtonian (Fig. 5a) than non-Newtonian (Figs.5b and 6a) flow regimes in this model.The increased number of strain localisation planes for our non-Newtonian flow model results in weaker isolation of Where the competent layer had no strain localising behaviour (Fig. 6a black dotted line), the competent layer was evenly stretched and thinned and no pinch and swell structures were formed, that is, R W remained at 1.

Analysis II: Impact of relative initial viscosity ratio
The viscosity of the matrix layers (A) was varied to see the impact of initial viscosity ratio (R V ) on the formation of pinch and swell structures (Fig. 6a).The ratio of the neck width to swell width (R W ) for the series (Fig. 6a) shows that R W decreases as R V increases, indicating that pinch and swell structures are initiated (R W < 0.4) and form (R W < 0.2) at lower stretch when R V values are higher.For Newtonian flow, pinch and swell structures were successfully initiated (i.e.R W < 0.4) from R V = 10, while for non-Newtonian flow this was not the case until R V = 125 (Fig. 6a) at stretch ∼ 2.4.

Analysis III: Impact of layer-dependent flow property behaviour
A summary of results where the Newtonian and non-Newtonian flow and Mohr-Coulomb strain localising behaviour of the surrounding layers (A) vs. competent layer (B layer) was systematically varied is provided in Table 2. Figure S4 shows numerical particle layer results at stretch 2.3.Pinch and swell structures form if layer B has Mohr-Coulomb strain localising behaviour (Table 2, rows i and iii) and do not form if it does not (Table 2 rows ii and iv).The results show that the strain localising behaviour of the surrounding layers (A) does not influence the formation of pinch and swell structures in the competent layer (B).The results also show pinch and swell structures can form in both Newtonian and non-Newtonian flow in all flow regime combinations of matrix vs. competent layer (Table 2, rows i and iii) if at least layer B can localise strain by Mohr-Coulomb constitutive strain localising behaviour.

Analysis IV: Impact of degree of material softening
Pinch and swell structures successfully initiate and develop (R W < 0.2) for all models where R Co > 1 (Fig. 7a to e and h).
The average tortuosity of the two competent layer edges in each model (plotted in Fig. 7h dashed line) shows a gradual increase with R Co to R Co = 10 where it plateaus.This corresponds to an increasing number of interacting strain localisation planes causing more complicated pinch and swell structure shapes in the R Co > 10 models (Fig. 7a, b and c).
A plot of differential stress vs. stretch for a representative point at x, y co-ordinates (0.4, 0.5) for the different cohesion ratios is graphed (Fig. 7g) and shows a sharp drop in stress after the material has failed.Where the material is significantly softer after fracturing (higher R Co values) the stress is lower, though high levels of softening (R Co = 20 and 40) show similar differential stress levels.The softened material shows some variability immediately after yielding, perhaps a numerical artifact, which then stabilises to constant stress at larger stretch values.High levels of material softening (R Co = 40) also appear to impact strain localisation, as an area of thick neck forms, where no single pair of strain localisation planes has dominated to form a pinch and swell structure (Fig. 7a, black arrow).

Multilayer model: Impact of strain localising behaviour on surface topology
Model results at stretch 2.3 for scales of 10, 20 and 40 km (Fig. 8) show the development of near-surface structure, topography, layer shape and the distribution of strain rate throughout the crustal-scale model.In addition, the development of surface layer tortuosity and R W is plotted.At low levels of stretch there is little difference seen in the surface tortuosity or neck to swell width ratios between all models presented.Where the lowest five layers, i.e. the middle crust, have no strain localising behaviour (Fig. 8a(i), b(i) and c(i)), strain is concentrated in the upper crustal layer, leading to fewer areas of focused localisation.This causes very simple near-surface strain localisation and a relatively low R W value, allowing the lower layers to be closer to the surface at all scales of the model.R W and T S values show little variation with the changing scale of the model.Where some of the middle crustal layers (i.e.B' layers), were given strain localising behaviour (Fig. 8a(ii), b(ii) and c(ii)) representing different rock units with contrasting rheological behaviour in the middle crust, strain localisation is strongest in the upper crust, but also occurs across all layers, leading to areas of relatively dispersed strain localisation and areas of homogeneous strain in the middle crust.The models show the upper crust has additional areas of interacting nearsurface strain localisation and R W values that increase as the scale increases as a result of the lower layers being further from the surface.Surface tortuosity (T S ) for the 10 × 10 km model (Fig. 8a(iii) dashed black line) is similar to that seen in the model where no localising behaviour was included in the middle crust (Fig. 8a(iii) dashed grey line).T S is decreased in the 20 × 20 km and 40 × 40 km scale models (Fig. 8b(iii) and show surface tortuosity (dashed lines) and R W (solid lines) for the upper crustal (blue) layer from the (i) and (ii) models.R W for models with no M-CB in the lower layers (i) increases as the model scale increases, while remaining similar for all scales where M-CB has been included in the lower layers.Surface tortuosity for models with no M-CB in the lower layers (i) remains the same for all model scales, while decreasing as the scales increase where M-CB has been included in the lower layers.Greater upper crust and surface variation is seen where Mohr-Coulomb strain localising behaviour is included in the middle to lower crust (ii models).c(iii)) where strain localising behaviour has been included in the B' layers.

Pinch and swell initiation by brittle deformation
Our observations from the two field examples of pinch and swell chains suggest brittle failures with variable amounts of material softening have initiated strain localisation in these structures.Field examples show that brittle failure is followed by viscous flow which has subsequently caused the formation of the pinch and swell structures.The numerical model used successfully initiates pinch and swell structures using Mohr-Coulomb constituent strain localising be-haviour and allows the structures to develop in subsequent viscous flow, simulating the pinch and swell structures seen in the St Anne Point and Wongwibinda field examples.The combination of field evidence and numerical model results confirm that combining Mohr-Coulomb strain localising behaviour (which includes brittle deformation) and subsequent viscous behaviour is a plausible rheological model for producing pinch and swell structures.
Our results support Mancktelow (2006) who concluded that ductile shear zones form after Mohr-Coulomb brittle failure at middle to lower continental crust depths.Mancktelow (2006) suggests that ductile flow with fluid influx requires ongoing pressure-dependent brittle failures, including micro-cracking for fluid to continue to flow into the shear zone.The thin-section-scale brittle failures observed in the St. Anne Point samples (Fig. 2d, e) are likely examples of the micro-cracking invoked.This micro-cracking develops into through-going micro-shear bands (Fig. 2e, yellow lines), which are then inferred to develop into outcrop-scale shear bands if enough strain is localised.
The discrete element numerical model used by Abe and Urai (2012) and Komoróczi et al. (2013) simulates mode I brittle failures, initiating pinch and swell structures by varying the competent layer properties from being semi-ductile, producing pinch and swell structures to being brittle, producing torn boudins.Our results are not directly comparable to these results, due to the contrasting numerical methods; our pinch and swell structures are initiated by strain localisation processes which do not readily form torn boudins.

Effect of stress exponent, strain localising
behaviour, R V and R Co for the competent layer and matrix Strain localisation of at least the competent layer is required for pinch and swell structures to form and these form readily in Newtonian flow (Fig. 5a) from viscosity ratios of ∼ 10, but in non-Newtonian flow (Figs.5b and 6), they form only at much larger initial viscosity ratios of ∼ 125.This result contrasts with other researchers who suggest that pinch and swell structures form where the stress exponent (n) is larger than 1 (Schmalholz and Maeder, 2012;Smith, 1977).It should be noted that our results are not directly comparable to these studies as their layers are purely viscous and the competent layer required an initial geometric irregularity, or instability, to grow to initiate the pinch and swell structure.By contrast, we use Mohr-Coulomb strain localising behaviour for pinch and swell structure initiation, causing irregular swell intervals and viscous flow only for post-failure structure development.Our results are also not directly comparable to the experiments by Abe and Urai (2012) and Komoróczi et al. (2013), as they used mode I brittle fractures to initiate the structures, while we have modelled strain localising behaviour to simulate the St Anne Point mode II brittle failures.Komoróczi et al. (2013), who simulated mode I brittle frac-ture for pinch and swell structure initiation, indicated that the matrix viscosity impacts the number of swell structures formed.Our model suggests this is possibly due to the impact of the viscosity ratio on the localisation of strain that initiates the pinch and swell structure formation.
In agreement with our results, many authors have previously suggested material softening enhances pinch and swell formation (Hobbs et al., 2009;Neurath and Smith, 1982;Smith, 1977;van der Molen, 1985).Our results are in general agreement.In addition we show that the degree of softening after failure has a marked effect on the nature of the developing pinch and swell structure (Fig. 7).At low levels of material softening (low R Co ), simple symmetric pinch and swell structures are formed, initiated by simple pairs of brittle fractures that develop into shear bands, while at higher levels of material softening, more interacting shear bands form, causing variable pinch and swell structure shapes and distributions to develop (Fig. 7).At even higher material softening levels, strain localisation is impaired as softening in the competent layer is very effective and no more strain localisation is necessary, causing a wide zone of weak material to govern the material behaviour and shape development.These differences in behaviour can be utilised in the interpretation of pinch and swell structures.For example, the presence of thick necks in the pinch and swell chain indicates that simple strain localisation planes were not initiated, possibly due to material softening within the competent layer accommodating all of the strain, thereby limiting strain localisation processes.Asymmetry, or rotation of planar features within a single swell, can be due to the uneven direction and strength of strain localisation at either end, forming the necks.This too can be used to indicate uneven distribution of strain localisation across an outcrop.

Application of the numerical model to field interpretation
The numerical model can provide some additional insights into the processes forming the pinch and swell structures: (i) rotation > 5 • of once planar structures in the swell structures suggests uneven strain localisation on the strain localising plane pairs; (ii) tortuosity of the competent layer edges > 1.2 and/or complex pinch and swell formations suggest significant material softening after strain localisation in the competent layer; (iii) the presence of a thick neck suggests strong material softening within the competent layer has reduced the ability of the material to localise strain; (iv) the ratio of neck width to swell width can be compared with the three stages of pinch and swell structure development (Fig. 6a) and a general idea of the relative layer viscosities can be deduced; and (v) shear band angles and rotation can provide an idea of the extent of stretch undergone.However, for a more accurate determination, both Newtonian vs. non-Newtonian flow and stretch should be independently validated by complementary methods (e.g.microstructural analysis).
The applicability of our model results is supported by comparison with the natural examples of pinch and swell structures observed on the St. Anne Point rock platform.Our model suggests that brittle failure of at least the competent layer is required for the St. Anne Point pinch and swell structures to form.The natural pinch and swell structures show both stage 2 development, with R W < 0.4 (Fig. 2c) and stage 3 development, with R W < 0.2 (Fig. 2a and b).Rotation of the bedding greater than 5 • can be seen in Fig. 2c, similar to that seen in the models where fracturing was dominated by one direction (Figs.5a and 6b).This, together with the presence of both stage 2 and 3 development, suggests that strain localisation was spread unevenly across the rock platform.The inferred positions of previous brittle failure planes (Fig. 2 red dashed lines) appear to have rotated away from the principal stress direction in a similar manner to the shear bands in the modelling (Fig. 6b, dashed line).
Thin section analysis showed the St. Anne Point samples display non-Newtonian flow.Therefore, comparison of model results (Fig. 6a) and nature (Fig. 2a and b) predicts a high viscosity ratio (R V > 125) between the mafic pinch and swell layer and the surrounding metasediments.At the same time, the St. Anne Point pinch and swell chains (Fig. 2a  and b) have well-separated swell structures, some with attenuated necks (Fig. 2a and b black arrows) and locally show thick necks (Fig. 2b black arrow).Accordingly, a high degree of softening i.e. larger R Co value (Fig. 7a, black arrow) has occurred variably across these rocks.

Implications of the presence of brittle-viscous behaviour in the middle continental crust on crustal-scale structures and topography
Pinch and swell structures are commonly observed in amphibolite facies terranes (e.g.Goscombe et al., 2004;Hobbs et al., 2010;van der Molen, 1985), supporting the interpretation that middle crustal rocks can be heterogeneous, with different rock units having contrasting rheological behaviour in terms of both viscous flow properties and strain localising behaviour.This is consistent with the pinch and swell structures seen at St Anne Point and the growing evidence that the flow behaviour of the middle continental crust should not be considered as exclusively viscous (e.g.Mancktelow, 2006).Our crustal-scale numerical models where two layers in the middle crust have heterogeneous rheological properties (Fig. 8a(ii), b(ii) and c(ii)) suggest that these layers impact both surface topography and strain localisation across the middle to upper crust.In our 20 and 40 km scale models, areas of heterogeneous and homogeneous strain are seen in the middle crust (Fig. 8b(ii) and c(ii)), possibly similar to the variations seen in seismic reflectors of the middle crust.Our numerical models create a viable scenario at the kilometre scale, suggesting that surface features, fluid influx pathways and, therefore, ore genesis, may be highly dependent on the heterogeneous nature of the middle and lower crust in terms of strain localising and viscous behaviour.If this is the case, then strain localising behaviour should be included in both small-and crustal-scale models of deformation within the Earth's crust.

Conclusion
Our numerical models simulate the brittle deformation we see in the Wongwibinda and St Anne Point pinch and swell chains by using an inherently unstable Mohr-Coulomb constituent material to initiate strain localisation followed by viscous flow to form the pinch and swell structures.The model results indicate: (i) these structures can form in either Newtonian or non-Newtonian flow if the competent layer can localise enough strain, contrary to previous conclusions, where viscous flow in a power law material had been assumed; (ii) higher viscosity ratios are required for pinch and swell formation in non-Newtonian regimes; (iii) low levels of material softening in the competent layer enhance strain localisation and pinch and swell development; (iv) higher levels of material softening can limit strain localisation, resulting in thick necks and (v) the country rock stress exponent and strain localising behaviour have little influence on pinch and swell formation.
We suggest that measuring planar bedding/foliation rotation, neck to swell thickness ratio (R W ), tortuosity and shear band rotation on a pinch and swell chain can provide important insights into the localisation of strain and the flow properties (e.g.stress exponent, relative viscosity, extent of stretching and degree of softening) of both the competent and country rock layers.
In addition, in crustal-scale numerical simulations, the incorporation of strain localising behaviour in some rock units at middle crustal levels can have a significant impact on the pattern of strain localisation throughout the continental crust and the development of near-surface structures.Consequently, we suggest some strain localising behaviour should be included with viscous flow when utilising such numerical models.
The Supplement related to this article is available online at doi:10.5194/se-6-1045-2015-supplement.

Figure 1 .
Figure 1.Pinch and swell structure characteristics, examples from Wongwibinda Metamorphic Complex, New England Orogen, N.S.W., Australia (30 • 18 41.6 S, 152 • 8 35.3 E).(a) Outcrop-scale examples of pinch and swell structures with positions of brittle failures shown as red, green and yellow lines, dependent on the angle from the primary stress direction, inferred to be perpendicular to the foliation.(b) Photomicrographs of a swell (sample WJ1655D2); brittle fractures and areas of fine grain recrystallisation, indicated by yellow lines, show extensive brittle failure in the competent plagioclase-rich layer.Brittle failure zones have likely rotated towards the horizontal during extensional deformation.

Figure 2 .
Figure 2. Pinch and swell structure characteristics, examples from the St. Anne Point rock platform, Fiordland, New Zealand (44 • 34 19.9 S, 167 • 46 55.3 E). (a-c) Outcrop-scale examples of pinch and swell structures; inferred positions of previous brittle failure shown as red dashed lines; (a) mafic layer showing pinch and swell structure with attenuated neck; (b) pinch and swell chain with attenuated and thick necks; red arrows show location of thin sections depicted in (d) and (e); (c) pinch and swell structure in mafic layer showing brittle deformation and rotation of planar features; (d) Photomicrographs of the centre of swell (sample AS1331I); two sets of stars (black and white) indicate brittlely fractured amphibole grains; foliation and stress direction indicated.Softer minerals biotite (Bt) and quartz (Qz) have filled brittle deformation localities.(e) Photomicrographs of the swell near the neck (sample AS1331G); four sets of stars (white, red and black) indicate brittlely fractured amphibole grains; yellow lines show residual fractures with direction of movement (white arrows).Thicknesses of grains either side of the shear band can vary slightly due to the rotation of the grain fragments in 3-D.Similar brittle fractures are inferred to have localised enough strain to form a through-going fracture, thereby allowing the formation of a neck.Two phases of biotite (Bt A and Bt B ) and quartz (Qz A and Qz B ) are indicated.The later-phase biotite (Bt B ) and quartz (Qz B ) are concentrated in the fracture zones.Minerals are labelled per Whitney and Evans (2010); modal percentages of the soft minerals decrease towards the swell centre.

Figure 3 .
Figure 3. (a) Examples of classic Mohr diagrams showing the theoretical faulting angle (θ) for two different materials, shown as black solid and grey dashed lines; red X signifies mode II failure; orange lines (solid and dashed) show failure planes used in numerical modelling with c i and c s denoting cohesion for the material before and after yielding.(b) Numerical model set up for three-layer model used in analysis tests I to IV; layer B more competent than matrix layers A, with passive marker beds (black lines) and stretch of −0.5 % on left and +0.5 % on the right at each model step.(c) Theoretical strain vs. stress diagram for elastic (cross-hatching and dots) and viscous flow (after Griggs and Handin, 1960); behaviour of layers A and B (with respect to three-layer model (b)) are shown: competent layer B -solid line, less competent matrix layers Adotted line; dashed line shows a layer with brittle deformation and material softening.Orange cross-hatching shows area of the theoretical model covered by the numerical model.Only yielding and post-yielding behaviour is investigated.(d) Multilayer crustal-scale model used to explore the effect of brittle behaviour on strain localisation and surface topography.The model includes an air layer (D) to allow realistic topographic feature development, an upper crust brittle layer (C) and five layers with two different material properties (layers A' and B') representing the middle to lower crust.B' layers are more competent than A' layers.Note: the models depicted in (b) and (d) show the particles overlaying the model framework; they do not represent a geometrical perturbation.

Figure 4 .
Figure 4. General model behaviour: snapshots (a-c) and graphs (df) showing characteristics of three sample particles from the n = 3, R V = 160, three-layer model.Numerical plotted results are shown for stretch 1.3 (grey dashed line on the graphs).Red and green stars indicate the location of two "tracer" particles in a low strain-higher stress area which develops into a swell, while the black star indicates the location of a particle in the higher stain-low stress area of the fracture which develops into a neck of the pinch and swell chain.Viscosity (d), strain rate (e) and differential stress (f) for these three particles are plotted against stretch.Note that the viscosity and strain rate graphs (d) and (e) show an inverse relationship and the stress invariant plot (c) shows less variation of stress across the whole model.See text for further explanation.

Figure 5 .
Figure 5. Analysis I results: series showing pinch and swell structure formation using the three-layer model, where layer B has Mohr-Coulomb strain localising behaviour, R V = 20; particle layer plots and strain invariant plots are shown at the stretch specified in the left hand column.In (a) all layers have Newtonian flow (n = 1) and in (b) all layers have non-Newtonian flow (n = 3).Pinch and swell structures are successfully formed in Newtonian flow, where strain is localised into a limited number of failure planes.Range of values is as specified while range of colour is the same for all model results showing strain rate invariant (˙ ); colour bars are included at stretch 1.0.Note: model results from stretch 1.6 have been zoomed to fit in the figure.

Figure 6 .
Figure 6.Analysis I and II results; (a) the effect of viscosity ratio (R V ) in both Newtonian (n = 1) and non-Newtonian (n = 3) flow on pinch and swell structure formation in the three-layer model.R W is minimum neck width/maximum swell width at that stretch.Pinch and swell structures are deemed to have initiated when R W < 0.4 and formed when R W < 0.2.R W decreases as R V increases indicating pinch and swell structures form more readily at higher R V values.Three stages of pinch and swell growth are shown.(b) Rotation of passive marker beds up to 15 • can be either positive (swell 1, black solid line) or negative (swell 2, grey solid line) as the model develops.Shear bands rotate from ∼ 55 • towards horizontal as stretch increases (black dashed line).Plotted data are from Newtonian (n = 1) series in Fig. 5a.

Figure 7 .
Figure 7. Analysis IV results: (a-f) model results at stretch 2.3 showing the effect of varying cohesion after softening in Newtonian flow with R V = 10.R Co is the ratio of initial cohesion to cohesion after softening (Eq.2); high R Co (a, b) indicates strong material softening while R Co = 1 (f) has no material softening.The black arrow in (a) shows formation of a thick neck.Strain rate invariant (˙ ) ranges of values are as specified while range of colour is the same for all images; colour bars are included on (a) and (e).Panel (g) shows stretch vs. differential stress graph showing reduction in differential stress as material softening (R Co ) is increased.Stress data are taken from a point towards the centre of the model at x, y co-ordinates of 0.4, 0.5.Differential stress drops after fracturing and fluctuates until reaching a steady state depending on the R Co .Panel (h) shows impact of R Co on competent layer edge tortuosity and width ratio (R W ; Eq. 6).Tortuosity increases as R Co increases to R Co = 10, reflecting more complex pinch and swell structures where R Co > 10.Pinch and swell structures are formed where R Co > 2 (i.e.R W < 0.2).

Figure 8 .
Figure 8. Multilayer (a) 10 × 10 km, (b) 20 × 20 km and (c) 40 × 40km model results: non-Newtonian multilayer models of the continental crust at stretch 1.5, showing differences in topography development (layer affinity results) and strain rate localisation (strain rate invariant results) where the middle to lower crust has a difference in viscosity (layer B' : layer A' is R V = 20); (a(i), b(i), c(i)) B' layers (pink) have no Mohr-Coulomb strain localising behaviour (M-CB); and (a(ii), b(ii), c(ii)) B' layers (pink) have initial Mohr-Coulomb strain localising behaviour, causing strain localisation through the whole model and more complex topography to develop.Panels (a(iii)), (b(iii)) and (c(iii))show surface tortuosity (dashed lines) and R W (solid lines) for the upper crustal (blue) layer from the (i) and (ii) models.R W for models with no M-CB in the lower layers (i) increases as the model scale increases, while remaining similar for all scales where M-CB has been included in the lower layers.Surface tortuosity for models with no M-CB in the lower layers (i) remains the same for all model scales, while decreasing as the scales increase where M-CB has been included in the lower layers.Greater upper crust and surface variation is seen where Mohr-Coulomb strain localising behaviour is included in the middle to lower crust (ii models).

Table 1 .
Moresi and Mühlhaus, 2006ameters used for Mohr-Coulomb strain localising behaviour in Analysis I (seeMoresi and Mühlhaus, 2006and text for details).In Analysis II R V values of 20, 40, 80, 100, 125 and 160 were used.In Analysis IV R Co values of 1, 2, 4, 10, 20 and 100 were used.defined as R V = η Bi /η Ai , where η Bi is the initial viscosity of the competent layer B and η Ai is the initial viscosity of the matrix layers A b defined as R Co = c i /c s , where c i is the initial cohesion of the competent layer and c s the cohesion after failure and softening a

Table 2 .
Analysis III results showing effect of systematically varying Newtonian/non-Newtonian flow and Mohr-Coulomb strain localising behaviour (M-CB) characteristics of the more competent layer B, with respect to matrix layers (A) on pinch and swell development (detailed model results are provided Fig. S2); an x indicates Mohr-Coulomb behaviour; ticks indicate pinch and swell structures formed in layer B for that simulation.