Scaling of viscous shear zones with depth-dependent viscosity and power-law stress–strain-rate dependence

SUMMARY One of the unresolved questions concerning fault deformation is the degree and cause of localization of shear at depth beneath a fault. Geologic observations of exhumed shear zones indicate that while the motion is no longer planar, it can still be localized near the down-dip extension of the fault; however, the degree of localization is uncertain. We employ simple analytic and numerical models to investigate the structural form of distributed shear beneath a strike-slip fault, and the relative importance of the physical mechanisms that have the potential to localize a shear zone. For a purely depth dependent viscosity, η = η 0 exp ( − z / z 0 ), we ﬁnd that a shear zone develops with a half-width δ w ∼ √ z 0 for small z 0 at the base of the layer, where lengths are non-dimensionalized by the layer thickness ( d km). Including a non-linear stress–strain-rate relation (˙ (cid:4) ∝ σ n ) scales δ w by 1 / √ n , comparable to deformation length scales in thin viscous sheet calculations. We ﬁnd that the primary control on the shear-zone width is the depth dependence of viscosity that arises from the temperature dependence of viscosity and the increase in temperature with depth. As this relationship is exponential, scaling relations give a dimensional half-width that scales approximately as where T 12 (K) is the temperature at the midpoint of the layer, R (J mol − 1 K − 1 ) the gas constant, Q (J mol − 1 ) the activation energy and β (K km − 1 ) the geothermal gradient. This relation predicts the numerical results for the parameter range consistent with continental rheologies to within 2–5 per cent and shear-zone half-widths from 2 to 6 km. The inclusion of shear-stress heating reduces δ w by only an additional 5–25 per cent, depending on the initial width of the shear zone. While the width of the shear zone may not decrease signiﬁcantly, local temperature increases from shear-stress heating range from 50 to 300 ◦ C resulting in a reduction in viscosities beneath the fault of several orders of magnitude and a concomitant reduction in the stresses needed to drive the motion.


I N T RO D U C T I O N
We observe faults as quasi-planar structures at the Earth's surface, but what governs the distribution of shear at depth? Geologic observations of exhumed shear zones (Poirier 1980;White et al. 1980;Rutter 1999;Norris & Cooper 2003) indicate that while the motion is no longer planar, it can still be localized near the downdip extension of the fault; however, the degree of localization is still uncertain (Wilson et al. 2004;Bürgmann & Dresen 2008). In this paper, we are interested in the structural form of distributed shear beneath a strike-slip fault and the physical mechanisms that have the potential to localize a shear zone. Previous work has fo-cussed on the effect of thermomechanical coupling. Yuen et al. (1978) used a 1-D model and found that once a shear zone forms it will remain localized due to shear-stress heating. Thatcher & England (1998) expanded upon this by looking at the effect of shear heating on a 2-D idealized strike-slip fault in an elastic lid over a linear, Newtonian, viscous layer. They found that, while a broad range of behaviour could be simulated, for reasonable parameter choices shear zones are narrow. The key conclusion was that shear localization is driven by dissipative heating near the axis of the shear zone, causing a localized reduction in viscosity. Takeuchi & Fialko (2012) used a 3-D model with periodic boundary conditions to simulate the 2-D case of a rigid lid over a viscoelastic layer and concluded that a thermally coupled non-linear power-law rheology could also generate such localization. They found a maximum thermal anomaly of 240 • C, in agreement with Thatcher & England (1998), capable of reducing the viscosity by up to 2.5 orders of magnitude in the region immediately beneath the fault tip. More recently, Duretz et al. (2014) modelled the development of shear zones in a 2-D viscous sheet undergoing compression. They concluded that the degree of shear-zone localization depends both on thermomechanical coupling and a non-linear rheology, and they predict half-widths of 1-7 km with typical temperature increases of 50-200 • C due to thermomechanical coupling for crustal shear zones.
The degree of localization of shear at depth and the effective viscosity structure beneath the fault have important implications for the interpretation of geodetic observations of surface deformation throughout the earthquake cycle. If shear is localized at depth, then the region of applicability for elastic half-space models of fault creep at depth (Savage & Burford 1973) is greater than previously thought. A localized shear zone would allow interseismic strain accumulation resembling the well-known arctangent distribution of surface displacements. The results we present on the localization of shear at depth complement previous work by Hetland & Hager (2005, 2006, Barbot et al. (2008) and others, where they develop analytical solutions for viscoelastic earthquake cycle models by allowing for the development of a weak zone beneath a fault. The development of a weak zone may also explain the relatively low apparent viscosities in the vicinity of a fault required to fit postseismic motion by Yamasaki et al. (2014).
In this paper, we investigate viscous shear zones and examine the interplay between three controls on the localization of shear in a viscous shear zone: the depth dependence of the viscosity, shearstress heating and a power-law stress-strain-rate relation. We start from the model used to investigate ductile shear zones in Thatcher & England (1998). Scaling laws are derived which can be used to infer shear-zone widths for a given viscosity structure and we show that shear-zone localization can be produced by a linear rheology with a depth-dependent viscosity alone.

MODEL
The model construction is that of an idealized strike-slip fault, infinite in length, embedded in a rigid lid over a viscous layer, as presented in Thatcher & England (1998). The lid physically represents the seismogenic brittle upper crust and the viscous layer the mid to lower crust undergoing ductile deformation. We treat this as a 2-D antiplane problem (˙ xx =˙ yy =˙ zz =˙ xz = 0, where˙ i j is the strain-rate tensor), so the region of interest becomes a vertical slice taken perpendicular to the fault ( Fig. 1) with all velocities assumed to be orthogonal to the plane and a function of position, u(x, z). As the problem is antisymmetric about x = 0, we need only analyse one half of the domain. Unless otherwise stated, all length scales are non-dimensionalized by the thickness of the viscous layer (d) and velocities by the driving velocity (v 0 ). The direction z is defined positive downwards, so the domain of u is (0 ≤ x ≤ ω, 0 ≤ z ≤ 1), with range 0 ≤ u ≤ 1 and aspect ratio ω = x max /d.
The main assumption of this model is that the system is fully relaxed beneath the rigid lid so we can neglect the effect of elasticity in this region, modelling it as a purely viscous process. This assumption is valid as long as the effective viscosities are low enough in the viscous layer that the relaxation timescale for a viscoelastic material, η/μ, is less than the earthquake recurrence interval. Motion is driven by fixed far-field velocities and the overlying lid is assumed to have block motion in a time-averaged earthquake cycle. The base of the layer is assumed to be shear-stress free and motion in the viscous layer is governed by the force balance equation with a viscous rheology (viscosity η) which for a 2-D antiplane problem reduces to at Bodleian Library on July 14, 2015 http://gji.oxfordjournals.org/ Downloaded from (∂ x u) 2 + (∂ z u) 2 . Rheological parameters used in the calculations are taken from Hirth & Kohlstedt (2004) and Hirth et al. (2001).

Uniform properties
Depth-dependent viscosity Arrhenius law Constant viscosity analytical solution, eq. (5). Contours are plotted at 10 per cent intervals of the driving velocity, with a dashed contour for 50 per cent, the aspect ratio ω = 4 and the shear-zone half-width δ w = 0.561.
At this stage eq. (3) is still general, but may be analysed further once a viscous rheology is selected (see Table 1). The boundary conditions are We also investigated the second model from Thatcher & England (1998) in which the lid is allowed to be passively displaced (shearstress free on the z = 0 surface), but it may be shown that in the absence of shear-stress heating the velocity field is a linear ramp (u(x) = x/ω) for any choice of rheology. However, even with shear-stress heating, for a wide range of parameters the shear zones widen to form a broad region of distributed deformation, as noted in Thatcher & England (1998). Thus, we did not develop that model further and instead concentrated on the model presented above. Finally, to aid direct comparison between the results of our models, a definition of shear-zone width (δ w ) is required; we use the halfwidth at the base of the layer, u(δ w , z = 1) = 1/2. We selected the half-width at the base of the layer as this aids comparison of the width of the shear zone between the different mechanisms and avoids second-order structural considerations.
In the following sections we analyse a constant viscosity layer, exponentially depth dependent viscosity and an Arrhenius law with a linear stress-strain-rate relationship. We then repeat the analysis for a non-linear stress-strain-rate relationship and finally examine the impact of shear-stress heating in both linear and non-linear regimes.

Uniform properties
For reference, we first investigate the case of a constant viscosity (η = η 0 ). As there are no gradients in viscosity, eq. (3) reduces to solving Laplace's equation, ∇ 2 u = 0, on a rectangular domain. This has a Fourier series solution, and is plotted for ω = 4 in Fig. 2. This solution is pointwise convergent, exhibiting the Gibbs phenomenon at the point x = 0, z = 0 due to the step discontinuity there. The shear-zone half-width for a constant viscosity layer, δ w = 0.5611, can be used to directly compare results for different rheologies and provides a maximum width of the shear zone. For ω 1, the shear-zone width should be independent of ω, governed only by the thickness of the viscous layer, d. As the solution is a Fourier series, it cannot be directly inverted to calculate δ w . However, we can investigate the behaviour of the limiting case of ω → ∞ by numerically calculating δ w as a function of ω. The Fourier series was evaluated with 50 000 terms and δ w (ω) calculated for ω ∈ [1, 4000] using a numerical minimization routine in Mathematica (Wolfram Research 2012). The limiting behaviour is illustrated in Fig. 3, and for comparison with the numerical models we select ω = 4 to give δ w = 0.5611 with precision of four significant figures (4sf).  . Contours are plotted at 10 per cent intervals of the driving velocity, with a dashed contour for 50 per cent, the aspect ratio ω = 4, and plots are labelled with z 0 and the shear zone half-widths, δ w .

Exponentially depth dependent viscosity
The next rheology we investigate is that of an exponentially depth dependent viscosity; as we show later, this is approximately what arises from the geothermal gradient with a temperature-dependent viscosity. Substituting the depth-dependent linear viscosity structure (η D ) from Table 1 into eq. (3) gives This equation also admits an analytical solution, in the form of a Fourier series: with the full derivation given in Appendix A. Visualizing these solutions (Fig. 4) indicates a clear dependence of the shear-zone half-width (δ w ) on the e-folding length scale (z 0 ). Velocity profiles at z = 1 are shown in Fig. 5, with the black profile representing the constant viscosity solution. Taking the limit as z 0 → ∞ in eq. (7) recovers eq. (5), as expected. The strong dependence of the shear-zone width on the vertical e-folding length scale prompted us to perform a simple scaling analysis to investigate their relationship. In the model, the shearstress free bottom boundary condition is such that it forces all Figure 5. Velocity profiles at z = 1 for exponentially depth dependent viscosities, eq. (7), with constant viscosity solution in black. Vertical dashed grid lines correspond to shear-zone half-widths (δ w ) for each z 0 . velocity contours to be perpendicular to the base of the layer, where the differential operators can be approximated under the assumption that u varies over the length scale d in the vertical direction. If we assume that z 0 is small then ∂ zz u min{∂ xx u, (1/z 0 )∂ z u}, where ∂ i /∂ ii indicates a first/second-order partial derivative operator in the ith direction. The ∂ zz u term may then be neglected in eq. (6), so that Replacing the derivatives by order of magnitude estimates gives at Bodleian Library on July 14, 2015 http://gji.oxfordjournals.org/ Downloaded from Figure 6. Shear-zone half-width for depth-dependent viscosity and Arrhenius law as a function of z 0 and effective z 0 (eq. 14), respectively. Scaling analysis of the depth-dependent viscosity equation shows that for small z 0 the fit should be proportional to √ z 0 . The red line is c √ z 0 , where the coefficient c = 0.937 was found empirically using a least-squares minimization for z 0 < 0.3. The limiting case of a constant viscosity is included as the dashed grey line. and the resulting scaling relation This scaling relation is plotted in Fig. 6 versus the shear-zone halfwidth obtained by numerically solving eq. (7) for x when z = 1 and u = 1/2.

Arrhenius temperature dependence
The next increase in complexity is to incorporate an Arrhenius law, η A from Table 1, which allows us to model the temperature dependence of viscosity. We choose to neglect the effect of pressure on the Arrhenius law as it is a minor contribution at crustal depths. At this stage, we are also neglecting the feedback mechanism due to shear-stress heating; for a temperature-dependent viscosity, η(T), eq. (3) becomes and substituting for η = η A gives, We assume a linear geotherm, where d e is the dimensionless lid thickness (t e /d), T is in K, and eq. (12) becomes We assume a linear geotherm for simplicity and because we are primarily interested in upper to mid crustal deformation where the geotherm is approximately linear; the same analysis would apply to more complex geotherms. As the depth, z, appears in the denominator of the coefficient of the ∂ z u term, eq. (13) requires a numerical method of solution. The numerical codes to solve eq. (13) utilize a successive over-relaxation finite-difference method written in Mathematica 9.0 (Wolfram Research 2012) and were benchmarked using the analytical solutions, eqs (5) and (7). Rheological parameters were selected to cover a wide range of lithologies and continental geothermal gradients (Hirth et al. 2001;Hirth & Figure 7. Log plot of viscosities for Arrhenius law, with zeroth-(equivalent to depth-dependent viscosity) and first-order Taylor approximations. Parameters are for wet quartzite (Q = 247 kJ mol −1 , B = 2.6 × 10 −4 K m kg −1 ) with layer thickness d = 30 km, elastic lid t e = 15 km and effective z 0 = 0.045 from eq. (14). Dimensional values are used to illustrate the viscosity range for a wet quartzite rheology in 30 km of crust.
Kohlstedt 2004). Resolution tests for our numerical models were run with grid spacings from 1/10 to 1/200 and they showed no more than 5 per cent variation in δ w , with that solely due to interpolation errors between grid points, where quadratic interpolation was used. Thus, the numerical solution fields are essentially independent of the grid resolution. Tests on non-linear geotherms show that the velocity structure is sensitive to the thermal gradient, with higher gradients producing more localization of shear.
Upon close examination of the numerical solutions it is apparent that the shear-zone half-width (δ w ) varied in a similar fashion, for certain parameter choices, to that of an exponentially depth dependent viscosity. Motivated by that similarity, we attempted to identify an e-folding length scale equivalent to z 0 . The zeroth-order Taylor approximation of eq. (13) is equivalent to eq. (6). The equation is Taylor expanded about z = 1/2 as this gives the closest approximation to the gradients of the viscosities for 1/2 < z < 1 (Fig. 7). This was tested by computing ∂ z ln η for the approximation and the Arrhenius law, and calculating the root mean square deviation as a function of expansion location (0 ≤ z ≤ 1), which showed a minima at z ≈ 0.5. The effective length scale is then If we use eq. (14) to define an effective z 0 , we find a velocity profile at z = 1 that is visually indistinguishable between the numerical solution and the analytical solution for the exponentially depth dependent viscosity, eq. (7), with that value of z 0 . Additionally, the shear-zone half-widths, δ w , for the numerical solutions are found to be equivalent, to four significant figures, with that of the corresponding analytical solution. The velocity fields from these numerical solutions are included in the first column of Fig. 8, with a dimensional shear-zone half-width,δ w = δ w d. A comparison of δ w from the numerical Arrhenius solutions and the analytical depthdependent viscosity solution, eq. (7), is illustrated in Fig. 6, where eq. (14) is used to define an effective z 0 for the Arrhenius model. The scaling of the shear-zone half-width then follows as eq. (10) with z 0 given by eq. (14). This empirically confirms our choice of expansion location, as a plot of δ w against effective z 0 only produces close agreement with the depth-dependent viscosity results when the expansion is taken about z = 1/2.  (Hirth et al. 2001;Hirth & Kohlstedt 2004). Each square plot has height and width d, while the rectangular plots have width 2d. Each plot is labelled with a dimensional shear-zone half-width,δ w , for d = 30 km. The shear-zone half-widths of the Arrhenius and the non-linear Arrhenius models agree withδ w predicted by eq. (33) to the number of digits stated. Contours are plotted at 10 per cent intervals of the driving velocity, with a dashed contour for 50 per cent, the aspect ratio ω = 4 and plots are labelled with the dimensional shear-zone half-widths,δ w .
If we compare the velocity field for the Arrhenius model to the equivalent exponentially depth dependent model (Fig. 9), we note that the Arrhenius model exhibits significant narrowing for z < 1/2. The equivalent exponentially depth dependent solution accurately captures the behaviour of a shear zone governed by an Arrhenius law for z > 1/2 with an absolute deviation between the velocity fields of less than 6 per cent, illustrated in Fig. 10. The narrower shear zone at shallow depths for the Arrhenius model is due to the much larger viscosity gradients at shallow depth (Fig. 7). To explore the behaviour for z < 1/2, a first-order Taylor expansion of eq. (13) about z = 1/2 is solved in Appendix B, which quantitatively matches the numerical results and is also included in Fig. 9.

N O N -L I N E A R D U C T I L E S H E A R Z O N E S
The shear-zone model may be extended to examine the effect of a non-linear power-law stress-strain relation˙ ∝ τ n . If we define the second invariant of the strain-rate tensor to be J 2 = 2 − 1 2 (∂ x u) 2 + (∂ z u) 2 , then non-linear equivalents to the rheologies of Section can be constructed with effective viscosities given by Table 1.

Uniform physical properties
Taking the effective viscosity, η n , from Table 1 and substituting it into eq. (3) gives the governing equation for the non-linear rheology with uniform physical properties, This equation is solved numerically using a successive overrelaxation finite-difference method written in Mathematica 9.0 (Wolfram Research 2012). The velocity fields are illustrated in Fig. 11 with the shear-zone half-width at the base of the layer (δ w ). Inspecting the results of the variation of δ w with n on a loglog plot (Fig. 12) leads us to infer a fit of the form δ w ∝ 1/ √ n. Using the numerical solutions as a guide to the analysis of eq. (15), it is apparent that for the majority of the shear zone the horizontal gradients in velocity are significantly greater than the vertical. In Appendix C, we use the assumption that this is generally true to derive the approximate governing equation, which has the analytical solution, Eq. (17) produces a velocity field that is in agreement within the accuracy of the numerical experiments (3sf) on all grid points for 1 ≤ n ≤ 100, where the Fourier series is truncated at an upper Nyquist frequency appropriate for the grid spacing. It is apparent from eq. (17) that the length scale in the Fourier coefficient is modified by a factor of 1/ √ n, matching that found empirically from the numerical solutions. It also exactly reduces to the analytical solution for the linear case when n = 1. The other limit of interest is that of purely plastic behaviour, or n → ∞, where it can be shown that eq. (17) reduces to Physically, as we approach plastic behaviour, all motion is localized onto a plane on the down-dip extension of the fault, with δ w = 0. Numerical solutions also show shear collapsing onto a plane in the limiting case of n → ∞ for eq. (15). A scaling analysis of the approximate governing eq. (16) produces a relation which is commensurate with both the numerical results (Fig. 12) and the analytical limit, eq. (18). The result of eq. (19) is similar to that found in both England et al. (1985) and Harig et al. (2010). England et al. (1985) investigated the equation for creeping flow in a thin viscous sheet with a non-Newtonian viscosity with imposed velocity boundary conditions, and their main assumption was that the velocity components normal to the imposed velocity boundary condition could be neglected. Though they were modelling deformation in plan view, and their starting equations were somewhat different, their equation could at Bodleian Library on July 14, 2015 http://gji.oxfordjournals.org/ Downloaded from Figure 11. Uniform parameter non-linear velocity field, eq. (17). Contours are plotted at 10 per cent intervals of the driving velocity, with a dashed contour for 50 per cent, the aspect ratio ω = 4 and the shear-zone half-widths, δ w , for each n are included. also be approximated as a quasi-Laplace equation modified by n, resulting in deformation decaying away from the driving boundary on a length scale dependent on 1/ √ n. Harig et al. (2010) investigated thinning of mantle lithosphere due to a Rayleigh-Taylor instability and found that the degree of localization of the downwelling and upwelling also varied as 1/ √ n.

Exponentially depth dependent viscosity
The modifications required to include depth dependence in the problem are relatively straightforward. Substituting the viscosity structure η Dn from Table 1 into eq. (3) produces which is solved numerically. The half-width of the shear zone at the base of the layer is then plotted in Fig. 13 as a func- As expected, this simplifies to give the solution for the depthdependent viscosity in the limit of n → 1, eq. (7), and the approximate solution for a non-linear ductile shear zone with uniform properties in the limit of z 0 → ∞, eq. (17). Eq. (22) also produces a velocity field that is in agreement within the resolution of the numerical experiments (3sf) on all grid points for 1 ≤ n ≤ 100. The scaling of the shear-zone half-width with both exponentially depth dependent viscosity and non-linear rheology follows as:

Arrhenius temperature dependence
The governing equation for this rheology is obtained by substituting η An from Table 1 into eq. (3), and for the temperature structure T(z) = T s + βd(d e + z) this becomes which is also solved numerically. Velocity fields for a non-linear rheology with Arrhenius temperature dependence are plotted in Fig. 8, with dimensional shear-zone half-widths,δ w = δ w d, for d = 30 km ranging from 2.34 to 3.52 km. From examining the numerical solutions it is apparent that the scaling relations already presented, eqs (14) and (23), predictδ w with an accuracy greater than 98 per cent. The approximate analytical solution for an exponentially depth dependent viscosity, eq. (22), combined with the effective e-folding length scale provided by eq. (14) also provides a qualitatively similar velocity field to the numerical solutions. However, the same caveat as illustrated in Fig. 9 applies in that the approximate analytical solution produces a wider shear zone than the numerical simulation for z < 1/2. Again, a first-order Taylor approximation improves this considerably, at the expense of significantly increased complexity, with these equations given in Appendix D.

S H E A R -S T R E S S H E AT I N G A N D T H E R M O M E C H A N I C A L C O U P L I N G
The last physical mechanism for shear-zone localization we investigate is thermomechanical coupling. The equations to be solved are based upon the linear and non-linear eqs (11) and (24) where the temperature field is now allowed to vary laterally. The solution for the temperature field arises from shear-stress heat production and heat diffusion governed by the heat equation, using the rate of shear heating at each time step, to modify the temperature field and hence the viscosities. The dissipative heating term, eq. (27), is dependent on both the velocity field, u(x, z), and the temperature field through the definition of viscosity, η. We must also include the rigid lid (Fig. 1) as part of the domain to be solved because heat will be lost via conduction through this layer, and we neglect the contribution of shear heating in the brittle fault zone. Frictional heating on the fault plane does not significantly affect the temperature field in the viscous layer as it is dominated by the heat generated by the stress singularity at the base of the fault. For the solution of eq. (26) we choose z = 0 to lie on the interface between the rigid lid and the viscous layer, so the domain of T is where d e is the non-dimensionalized elastic layer thickness, t e /d. A linear geothermal gradient, T(t = 0) = T s + βd(d e + z), is used as the initial condition on temperature for the heat equation, and dissipative heating, eq. (27), allowed in the region 0 ≤ z ≤ 1 only. The boundary conditions used for the heat eq. (26) are

Arrhenius temperature dependence
For the Arrhenius temperature dependence with a linear stressstrain-rate relationship we replace η in eq. (27) with η A from Table 1. The initial velocity condition for the shear heating term, u(t = 0), is then taken to be that of the Arrhenius solution with a linear geothermal gradient. Eq. (13) is solved first, and then the heat eq. (26). The steady state solution is obtained by iteration between solving eqs (12) and (26). The numerical codes for solving the heat equation utilize the method of lines with a pseudospectral Chebyshev grid for the spatial discretization, as implemented in the NDSolve routine from Mathematica 9.0 (Wolfram Research 2012). A Chebyshev grid is used to ensure that small thermal transients near the axis of heating are accurately captured. The model is run until it reaches an approximately steady state for both temperature and velocity fields. We define steady state as less than a 0.1 per cent change in the thermal and velocity fields over 1 Ma. The time evolution is controlled by the thermal diffusion timescale, which for our model ranges from 0.5 to 1 Ma. The derivation of this timescale can be found in Appendix E, where the region of shear heating is approximated by an infinite half-plane source with depth-dependent strength.
A comparison of the velocity fields with and without shear-stress heating are given in Fig. 8, while Fig. 14 illustrates the temperature increase due to shear-stress heating at steady state (t ∼ 5τ ). The variation in shear heating is strongly dependent on the initial width of the shear zone which is itself controlled by the rheology. If the strain rates are not large enough, or equivalently the effective z 0 not small enough (z 0 < 0.2), then no significant shear heating will be generated on the thermal diffusion timescale, and the temperature field remains relatively unperturbed. Thus, for thermomechanical coupling to be a major contributing factor in shear-zone localization, at Bodleian Library on July 14, 2015 http://gji.oxfordjournals.org/ Downloaded from it requires rheological parameters that produce a narrow shear zone initially, with effective z 0 < 0.2. The rheological parameters that produce a shear zone this narrow and generate the most shear-stress heating are those associated with dry and wet olivine, lithologies that we would not expect to find in the upper crust. It should also be noted that even when z 0 = 0.1 and significant shear heating is generated, such as for dry olivine, the change in the shear-zone width is only ∼20 per cent.

Non-linear rheology with Arrhenius temperature dependence
Thermomechanical coupling for a non-linear Arrhenius law is modelled in the same way as described in Section 5.1, except that the non-linear effective viscosity η An is used in place of the linear η A (Table 1). If we look at the shear heating term, eq. (27), the dependence on the gradients in velocity is actually weaker for a non-linear rheology than a linear one. Substituting η An into eq. (27) gives, for 0 < z < 1, The exponent for the last term in eq. (30), (n + 1)/2n, is a maximum for n = 1 and decreases monotonically, approximating 1/2 in the limit n → ∞; so the strongest dependence of the shear heating production on the velocity gradients is for a linear rheology, though the initial width of the shear zone decreases with increasing n. The numerical solvers iterate between eq. (24) and the heat eq. (26) with the shear heating term given by eq. (30). Velocity profiles for a non-linear Arrhenius rheology with shear heating are plotted in Fig. 8, next to the equivalent model without shear heating. Fig. 14 illustrates the temperature increase due to shear-stress heating, in comparison to a linear rheology. One consequence of our model is that despite starting the iteration for the thermomechanically coupled model with a narrower shear zone due to the non-linear rheology, which in a linear rheology would result in greater shear heating, less shear heating is generated overall because of strain-rate softening. However, because of the non-linearity of eq. (24) significant localization of shear can still occur, in some cases more than for a linear rheology despite having a lower maximum temperature increase, T max .

D I S C U S S I O N
We have investigated the interplay of three physical mechanisms responsible for the localization of shear at depth beneath a strike-slip fault: depth-dependent viscosity, a non-linear stress-strain relation and thermomechanical coupling driven by shear-stress heating. Velocity fields illustrating the effect of each physical mechanism are summarized in Fig. 8. Rheological parameters were selected to simulate a wide range of possible behaviour and are not necessarily indicative of the lithology that would be expected at those depths. For example, while we would not expect a predominantly dry olivine composition at 15 km, it elucidates the behaviour of any material with a similar activation energy.
The strongest control on the shear-zone width is the effective length scale of vertical variation of viscosity, as a result of an Arrhenius law with a linear geothermal gradient. Simplifying the force balance eq. (1) under the assumption of antiplane conditions, and substituting in eq. (2) produces From eq. (32) it can be seen that when the right-hand side is large due to rapid variations in viscosity, to balance it on the left-hand side requires a reduction in the horizontal length scale. This leads to the scaling relation where the horizontal length scale is inversely proportional to the square root of the vertical gradients in viscosity.
Relative to the constant viscosity solution, the depth dependence of viscosity alone can produce a reduction in width of 60-75 per cent, while a non-linear rheology with n = 3 produces a further reduction in width of (1 − 1/ √ 3) ≈ 42 per cent. Shear-stress heating can add an additional 20-25 per cent to these effects. For a range of rheological properties and geothermal gradients, shear zones with a half width of 1.7-6 km can develop in continental crust 30-40 km thick, compared to 17 km for a constant viscosity layer or 10 km for a layer with uniform non-linear properties and n = 3. These widths scale linearly with the thickness of the crust, so thinner crust would result in correspondingly narrower shear zones.
The half-width of the shear zone with an Arrhenius rheology is illustrated in Fig. 15. To look at the scaling of the shear-zone width for a non-linear Arrhenius rheology we can combine eqs (14) and (23) using the empirical scaling coefficient 0.937 from Fig. 6. This is further simplified by noting Q RT1 2 , where T1 2 is the temperature in the middle of the layer, T s + (d e + 1/2)βd, producing a dimensional shear-zone half-width This relationship predicts the calculated half-widths from Fig. 8 to within 2 per cent, and is accurate to within 5 per cent for the full range of parameters consistent with earth rheologies. All linear and non-linear Arrhenius calculations with effective z 0 < 0.4 (equivalent to Q > 100 kJ mol −1 for β > 10 • C km −1 ) are plotted against this predicted width in Fig. 16, which shows the departure from eq. (33) at large z 0 (outside the range of reasonable parameters). For most of the parameter regime appropriate to the continental crust, z 0 < 0.2. The introduction of shear-stress heating and thermomechanical coupling is important in the cases where the shear-zone width is already narrow due to the depth dependence of viscosity and a nonlinear stress-strain-rate relation. The addition of both a non-linear rheology and shear-stress heating to the depth dependence due to the geothermal gradient is capable of producing a shear zone of only 1.7 km. Depending on the choice of rheological parameters, the viscosities may already be low, in which case the change in temperature due to shear-stress heating is also relatively low and the reduction in viscosity correspondingly smaller. Fig. 17 illustrates the effective viscosities for each of the parameter sets and physical mechanisms from Fig. 8 relative to an Arrhenius temperature dependence. It should be noted that the greatest reduction in effective viscosity requires a dry olivine composition, which we would not expect to find at these depths, though it provides an upper bound to the viscosity reductions that may be generated by shear heating.
The combined power law and Arrhenius dependence show a large reduction in viscosity in the immediate vicinity of the shear zone, which returns to a background uniform depth dependence outside of the shear zone. The reductions here are significantly more localized than those generated by shear-stress heating, which produce greater changes in viscosity but over a broader region. The gradients in viscosity, which are the primary control on the shear-zone width, are not as large. For both the non-linear rheology and the non-linear rheology with shear-stress heating, the viscosities are reduced in a relatively narrow vertical region approximately parallel to the downdip extension of the fault (Fig. 17). This reduction in viscosities directly underneath the fault should be resolvable when examining post-seismic displacements. To construct an earthquake cycle model using a Maxwell viscoelastic rheology consistent with both interseismic and post-seismic observations of displacement around the western North Anatolian Fault Zone, Yamasaki et al. (2014) required the addition of a rectangular region of low viscosity directly beneath the fault, with relatively sharp boundaries, to explain the amplitude of the post-seismic transient in the near-field (<30 km from the fault). Our results suggest that a non-linear thermomechanically coupled rheology could generate such a region, with the effect becoming more pronounced as the shear zone approaches steady state. The non-linear rheology would ensure relatively sharp boundaries, while the thermomechanical coupling would localize the low-viscosity zone to a region surrounding the fault tip (Fig. 17). A wet quartzite rheology produces a viscosity change of 10 −2 to 10 −4 , in agreement with the viscosity change they require of 10 −2 . A non-linear rheology on its own produces a vertically extended zone of reduced viscosity, which would not necessarily fit these data.
The effect of strain localization due to shear-stress heating can only be included in the case of the Arrhenius rheology as temperature dependence is required for the feedback mechanism. Practically, this means that the shear zone has already been localized significantly due to the depth dependence of the viscosity before thermomechanical coupling is considered. The additional reduction in shear-zone width is dependent upon this initial width and in most cases only provides between 5 and 20 per cent of additional narrowing, producing shear-zone half-widths of 3-6 km for linear or 1.7-3.5 km for non-linear rheologies (n = 3). However, this reduction in shear-zone width would be less pronounced if heat were advected away from the fault, such as in the case where the fault is allowed to migrate laterally, though the depth dependence of viscosity will still generate a narrow shear zone.
Our shear-zone half-widths are consistent with those published in Thatcher & England (1998), where their non-dimensional constant is equivalent to the reciprocal of our length scale (Q/(RT Max ) ≈ 1/z 0 ). After conversion their results scale as δ w ≈ 1.6 √ z 0 , with the larger coefficient a direct result of their definition of shear-zone half-width u(δ w , 1) = 0.8. Their work attributed the localization of shear to thermomechanical coupling, while Takeuchi & Fialko (2012) attributed the localization of shear to both thermomechanical coupling and a non-linear rheology. In considering a linear Maxwell rheology, Takeuchi & Fialko (2012) restricted their models to the case of constant viscosity only, where they found no significant localization. However, in order to exclude a linear Maxwell rheology, it is first necessary to consider the effect of the depth dependence of viscosity. We include viscosity variations with depth in our linear model and suggest that the primary control on the localization of shear beneath a strike-slip fault is the dependence of the viscosity on the geothermal gradient, with secondary contributions from non-linear creep mechanisms and shear-stress heating.
The main assumption of our model is that the earthquake cycle can be approximated, over geologic timescales, by the block motion of a rigid lid over a viscous layer. For the case of a pure Maxwell rheology this approximation is valid. However, it is not entirely clear in the case of a non-linear rheology that we can neglect the contribution of the build-up of elastic strains and treat the solution as a fully relaxed viscoelastic material behaving as a viscous media. Takeuchi & Fialko (2012) found that steady state for a nonlinear rheology took considerably longer to achieve than for any linear rheology, and that a non-linear rheology with uniform properties would lead to an unfeasible build-up of stresses (many GPa) at Bodleian Library on July 14, 2015 http://gji.oxfordjournals.org/ Downloaded from Figure 17. Comparison of viscosity fields for given material parameters (Hirth et al. 2001;Hirth & Kohlstedt 2004). Arrhenius temperature dependence is shown as absolute viscosities, while additional mechanisms are plotted as viscosity change compared to background Arrhenius temperature dependence. Each plot has height and width d. Contours are plotted at intervals of two orders of magnitude. due to the elastic strains, which would not be the case for a temperature-dependent rheology. Once their models are re-run with a depth-dependent Maxwell rheology, we can directly compare their time-dependent earthquake cycle model with our approximation.
When discussing creep mechanisms, the dependence of diffusion creep upon grain size should not be neglected. Cataclasis and dynamic recrystallization (Goodwin & Wenk 1995) are both mechanisms that can cause a reduction in grain size in mylonites (Sammis et al. 1986;Biegel et al. 1989;Marone & Scholz 1989;Blenkinsop 1991;Ray 1999). The effective viscosity is proportional to the grain size, η ∝ d m , with the exponent m between 0.7 and 3 (Gueydan & Précigout 2014). The equilibrium grain size depends on both the shear-strain rate (d ∝ 1/˙ for nonnegligible strain rates) and temperature (d ∝ T). We would expect the flow profile and effective viscosities to demonstrate further localization of shear, and reduced effective viscosities, close to the down-dip extension of the fault where the strain rates are high-est. This should act in a similar fashion to including a non-linear stress-strain relation.
Without further calculations, we cannot rule out the possibility of a grain size dependence resulting in a widening of the shear zone at depth. This is possible if the grain size temperature dependence is stronger than the temperature dependence of the Arrhenius law and the strain rates are low. An increase in grain size with depth would result in an increasing viscosity with depth, producing a negative effective z 0 . Evaluating eqs (6) and (7) with a negative z 0 produces δ w greater than that for a constant viscosity (δ w = 0.5611).

C O N C L U S I O N S
The localization of shear at depth arises directly from the dependence of the viscosity on the geothermal gradient, with secondary contributions from non-linear rheologies and thermomechanical at Bodleian Library on July 14, 2015 http://gji.oxfordjournals.org/ Downloaded from coupling. For a purely depth dependent rheology the range of parameters consistent with continental lithologies all result in shear zones of less than 6 km for continental crust less than 40 km thick. A non-linear rheology with uniform properties can generate shear zones of ∼11 km; narrower than that of a constant viscosity layer (18 km). Simple scaling relationships are derived which can be used to predict shear-zone widths for given parameter combinations, these relationships agree with the numerical results for continental rheologies to within 5 per cent. The results of Yamasaki et al. (2014) indicate that a region of low viscosity beneath the rupturing fault is required to fit both pre-seismic and post-seismic rates of deformation on the North Anatolian Fault, and our results suggest that a non-linear rheology plus thermomechanical coupling could produce such a reduction in viscosity in this region. The narrowest shear zones from our model include contributions from thermomechanical coupling, a non-linear rheology and the depth dependence of viscosity to produce shear zones of 1.7-3.5 km for a wide range of parameter choices, significantly narrower than the 17 km shear zone produced by a constant viscosity layer. Nonetheless, the primary control on shear-zone width remains that of the dependence of the viscosity on the geothermal gradient, especially as that is a pre-requisite for significant shear-stress heating.

A C K N O W L E D G E M E N T S
This work was supported by the Natural Environment Research Council through a studentship to James Moore, the Looking into the Continents from Space project (NE/K011006/1) and the Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET).
We thank Philip England for helpful discussions during the course of this work.