N Traoré

L Le Pourhiet

J Frelat

F Rolandone

B Meyer

Does interseismic strain localization near strike slip faults result from boundary conditions or rheological structure?

Keywords: Numerical approximations and analysis, mechanics, theory and modelling, continental tectonics: strike slip and transform, fault zone rheology, creep and deformation, dynamics and mechanics of faulting

Geodetic observations across most of the major strike-slip faults show an interseismic strain rate which presents a sharp localization of elastic shear strain in the fault vicinity (20-60 km).

The screw dislocation model of Savage & Burford is commonly used to fit these geodetic data and to retrieve the far field velocity and the locking depth. This model is very popular because of its inherent simplicity to derive fault slip rates, and mostly because it predicts locking depths which are of the same order of magnitude as the base of the seismogenic zone (5-20 km). A first issue with the screw dislocation model is that localization is paradoxically introduced by imposing a step function in the velocity field at a depth where the crust is otherwise recognized to behave following a viscous rheology. A second issue with this model is that it is not consistent with the rheological model of the crust that is valid for both postseismic (1 -10 yr), interseismic (1-several kyrs) and long term geodynamic (several kyrsseveral Myrs) time-scales. Here we use numerical models to study how alternative and more geologically realistic boundary conditions and rheological structures can lead to the localization of elastic

N.Traoré et al.

strain at the earth surface during the interseismic period. We find that simple elastic models resembling the Savage & Burford model but driven by far field plate velocity are inefficient at localizing strain unless this driving velocity is transmitted by a rigid indentor. We also find that models including a weak viscous heterogeneity beneath the fault zone are able to produce appropriate localization of the deformation near the fault. This alternative class of model is shown to be pertinent in regards to the boundary conditions and geological observations along exhumed ductile strike slip shear zone.

INTRODUCTION

The study of geodetic data across strike-slip fault zones is used for unraveling the geological and mechanical characteristics of these zones, and plays a key role in our understanding of the lithosphere behavior. Most of the major strike slip faults, such as the San Andreas fault [START_REF] Lyons | Fault creep along the southern San Andreas from interferometric synthetic aperture radar, permanent scatterers, and stacking[END_REF], the North Anatolian fault [START_REF] Mcclusky | Global Positioning System constraints on plate kinematics and dynamics in the eastern Mediterranean and Caucasus[END_REF][START_REF] Provost | 3d mechanical modeling of the GPS velocity field along the North Anatolian fault[END_REF][START_REF] Flérit | The mechanical interaction between the propagating north anatolian fault and the back-arc extension in the aegean[END_REF] or the Sagaing fault [START_REF] Vigny | Present-day crustal deformation around Sagaing fault, Myanmar[END_REF] show a surface deformation that presents a strong localization along the fault which corresponds to a steep gradient in fault parallel velocity across the fault.

The [START_REF] Savage | Geodetic determination of relative plate motion in Central California[END_REF] model is the simplest approach to compute synthetic interseismic velocity field across strike slip faults that approximates the interseismic velocity, obtained by GPS or InSAR data, across the major faults [START_REF] Lisowski | The velocity field along the San Andreas Fault in Central and Southern California[END_REF][START_REF] Smith-Konter | Locking depths estimated from geodesy and seismology along the San Andreas Fault System: Implications for seismic moment release[END_REF]). This model relies on two geological hypothesis. Firstly, the seismogenic part of the fault of thickness d does not slip during the interseismic period. Secondly, the fault itself extends at depth beneath this so-called locked part of the fault, and accommodates the plate relative motion s in a localized but aseismic manner.

Mechanically this system is analogous to a screw dislocation in an elastic solid. Considering this screw dislocation occurs in a semi-infinite half-space, [START_REF] Weertman | Elementary dislocation theory[END_REF] have shown that the velocity v at the surface of this infinite half-space can be described as a function of the distance to the fault x following

v = - s π arctan x d . (1)
[Figure 1 about here.]

The geodetic observations provide easily two parameters: the far field velocity and the width of the zone where the deformation is mainly localized. In the screw dislocation model, the localization width does not appear explicitly. However, using the method of the tangent described in Appendix A and Fig. 1, we find that this width, W c is directly proportional to the locking depth d

W c = πd, (2)
and therefore depends neither on the rheological properties nor on the far field velocity. On Earth, the observations infer a width of the localization zone ranging between 20 to 60 km, so that the retrieved locking depth d always remains between 5 and 20 km in good agreement with depth distribution of earthquakes along strike slip fault zones [START_REF] Miller | Thermal-mechanical controls on seismicity depth distributions in the San Andreas fault zone[END_REF]. Moreover, the blocks outside a fault zone are usually assumed to behave rigidly, that is also a characteristic of the screw dislocation model (Appendix A). Together with its simplicity, these two properties of the [START_REF] Savage | Geodetic determination of relative plate motion in Central California[END_REF] model have largely contributed to its popularity.

Following this approach, other elastic models originating from crystalline physics have been transferred to geological applications. They fall in two categories, the dislocation methods and the back-slip methods. On the one hand, the dislocation methods consist in considering one or several dislocations embedded in an elastic material, as [START_REF] Savage | Displacement field for an edge dislocation in a layered half-space[END_REF] and Le Pichon et al. ( 2005) using edge dislocation models, or [START_REF] Savage | Dislocation pileup as a representation of strain accumulation on a strike-slip fault[END_REF] via the introduction of dislocation pileups. Multiple screw dislocations [START_REF] Fay | Fault slip rates, effects of elastic heterogeneity on geodetic data, and the strength of the lower crust in the Salton Trough region, southern California[END_REF][START_REF] Brink | Three-dimensional models of deformation near strike-slip faults[END_REF] have also been used to produce more complex faults system and study the interactions between fault segments. On the other hand, the back-slip method consists in subtracting the coseismic slip rate, given by the elastic analytical solution of [START_REF] Okada | internal deformation due to shear and tensile faults in a half-space[END_REF] for finite dislocation, from the long term geologic stepwise velocities [START_REF] Bennett | Global Positioning System constraints on fault slip rates in southern California and northern Baja, Mexico[END_REF][START_REF] Meade | Estimates of seismic potential in the marmara sea region from block models of secular deformation constrained by global positioning system measurements[END_REF][START_REF] Becker | Constraints on fault slip rates of the southern California plate boundary from GPS velocity and stress inversions[END_REF]. All these models assume a disconti-nuity in the horizontal velocity field at the base of the brittle crust during the interseismic period.

The mechanical appropriateness of this particular basal boundary condition is yet questionable.

Imposing the velocity underneath the locked fault geologically implies that the aseismic parts of the crust and the lithosphere impose a basal drag on the brittle crust during the interseismic period. This requires the viscosity of the ductile crust to be large enough to promote a strong viscoelastic coupling between the brittle and ductile crust. Paradoxically, the velocity profile in this part of the viscous crust is assumed to be a step function, requiring the viscosity to be extremely small bellow the locked fault.

In order to study transients in the velocity gradients across a strike slip fault during the seismic cycle, the Savage & Burford solution was extended for a viscoelastic space by [START_REF] Savage | Asthenosphere readjustment and the earthquake cycle[END_REF], and later by many others [START_REF] Johnson | Influence of lithosphere viscosity structure on estimates of fault slip rate in the mojave region of the san andreas fault system[END_REF][START_REF] Hetland | Interseismic strain accumulation: Spin-up, cycle invariance, and irregular rupture sequences[END_REF][START_REF] Pollitz | Viscoelastic shear zone model of a strike-slip earthquake cycle[END_REF][START_REF] Pollitz | A physical model for strain accumulation in the san francisco bay region[END_REF] in order to account for more complex vertical layering of the viscoelastic lower layer.

However, all these viscoelastic coupling models (e.g. [START_REF] Savage | Viscoelastic coupling model of the San Andreas Fault along the big bend, southern california[END_REF][START_REF] Savage | Viscoelastic-coupling model for earthquake cycle driven from below[END_REF] share a lot of similarities with the original Savage & Burford model when it comes to the treatment of the elastic upper layer during the interseismic period. This layer is indeed always considered as an horizontally homogeneous elastic layer containing a screw dislocation. The main difference between the original Savage & Burford model and the viscoelastic coupling models is that the amplitude of the displacement discontinuity imposed by a Burgers vector below the elastic layer is depth-and time-dependent and accounts for postseismic relaxation.

The present work focuses on alternative mechanical models that permit strain localization at the surface without imposing it through a Burgers vector located underneath a locked fault. We are not focusing on postseismic deformations and related transients in the velocity field, instead we aim at finding a model satisfying both the secular interseismic loading and the long term observations. We therefore assume that our models represent the part of the seismic cycle that is recorded late after the last large earthquake and explore how the choice of the boundary conditions and the rheological structure of the model impacts strain localization at the surface. We have designed alternative models in which both the boundary conditions and rheological structure are geologically and mechanically justified. The resulting surface velocity fields from these alternative models are obtained using the finite element method.

NUMERICAL METHODS

Constitutive equations

The numerical experiments were run with Cast3M (C.E.A / D.M.T / L.A.M.S), a Finite Element software developed by the CEA, the French Alternative Energies and Atomic Energy Commission.

Cast3M permits both the definition of the models (geometry, meshing, rheology) and the solution of the numerical problem.

We consider the continuous problem to find the stress tensor σ that solves the equilibrium equations, associated with boundary conditions. In linear elasticity, the stress tensor is defined by the generalized Hooke's law:

σ = C : , (3)
where C is the fourth order rigidity tensor and the Cauchy's strain tensor for small deformations defined by = 1 2 (∇u + ∇u t ), where u is the displacement field. In viscoplasticity, the stress tensor is given by the constitutive equation:

σ = C :( -v ) , (4)
the time derivative of the viscous strain being a function of the stress tensor v = f (σ).

The problem has been then discretized using the classical Finite Element method that leads to solve the linear system

Ku = F , (5)
or K u = Ḟ in the viscoplastic case, where K the stiffness matrix and F the force vector that takes into account the displacement loadings as well as the contribution of the viscous deformation term in the viscoplastic case. The plastic yield value is equal to zero that corresponds to the effective viscoelastic case.

The transient state is out the scope of this paper and we consider the part of the interseismic cycle that is far beyond the coseismic and the postseismic periods for which one can posit the solution corresponds to a steady state in the viscoelastic case. In order to compare the model with geodetic data, we study the discrete elastic problem in velocity rather than in displacement, assuming a constant strain rate which is also consistent with our steady state approximation.

General description of the numerical models

The geometry chosen is a rectangular box. The Young's modulus of the medium is equal to 10 10

Pa and the Poisson's ratio is 0.25. All the numerical models presented here share some common boundary conditions (Fig. 2), i.e. a zero velocity is imposed on the left and right lateral boundaries in the x direction and for the basal surface in the z direction (free slip conditions).

[Figure 2 about here.]

Periodic models have the velocity of the vertical plane at the back of the model that is equal to the velocity of the vertical plane at the front. This type of model can be viewed as a slice of an infinitely long domain in the y direction, where therefore the fault is infinite.

The domain is meshed with elements that are triangular prisms (15 nodes per elements). The mesh is refined in the vicinity of the fault plane, the typical element size near the fault is 0.5 km, and 1 km near the domain boundary. The number of elements in total, for instance for our benchmark, is 65640, that represents 224507 nodes. We tested the sensitivity of the models to changes in the number of elements (half the initial number) with no significant results variation.

Benchmarks

First we compare the analytical solution given by Eq. ( 1) with a benchmark of models built with Cast3M.

The geometry of the benchmark model is given in Fig. 3(a). The benchmark model is 100 km wide in y but boundary conditions are periodic in this direction so that the fault can be considered as infinitely long along strike. The thickness of the locked segment is h 1 =1 2km for all the models of the benchmark. The boundary conditions consist of imposing the sliding velocity on the vertical surfaces of the unlocked portion of the fault, with a value of 1 cm/yr on the left side and -1 cm/yr on the right side. Between the locked part of the fault where sliding velocity is zero and the lower part of the fault where velocity is ±1 cm/yr, there is a point of singularity. In order to obtain a reliable numerical solution, it was necessary to add an intermediate domain of thickness 0.1 km (which is not represented in the figure) where free slip conditions are imposed on the fault plane.

The benchmark consists of 24 simulations for which we systematically varied the thickness of the sliding part of the fault, (h 2 =10-20 -30 -40 -50 -100 -200 -1200 km), and the width of the model (L 0 =400-800 -1200 km). The range for the sliding thickness of the fault, h 2 , was chosen to cover geologically acceptable values of 10 -50 km to a value of 1200 km intended to simulate an infinite half-space. The range for the width of the model, L 0 , was chosen to correspond to the extent of typical geodetic data set or to simulate numerically an infinite half-space.

For all the models, the deformation is observed on a reference line that lies horizontally on the upper free surface, perpendicularly to the fault. The parameter h 1 corresponds to the parameter d in the analytical solution given by Eq. ( 1), and the total imposed velocity (2 cm/yr) corresponds to the parameter s. Matching the deformed reference line with the arctangent functions defined by Eq. ( 1), we retrieved the values of the parameters d and s which provided the best fit, following the method of least square regression detailed in Appendix B. Then we computed the relative errors (Fig. 3 b) between the numerical input parameters (12 km and 2 cm/yr) and the parameters retrieved from the analytical solution (d et s) following the formulas |12 -d|/12 for the depth error and |2 -s|/2 for the slip error.

[Figure 3 about here.]

The results of the benchmark show clearly that the relative error obtained with Cast3M diminishes as the model approaches the dimension of a semi-infinite half-space. The break in the convergence curve occurs as the steady slipping part of the fault approaches 200 km in thickness. The effect of the width of the model produces a rather small error, which corresponds to the residual error for an unlocked slipping fault greater than 200 km in thickness. The model with L 0 =1 2 0 0km width and h 2 =1 2 0 0km reaches an error of less than 3%, which allows us to conclude that our model is sufficiently well resolved for studying the effects of other geometries and/or different boundary conditions. However, the relative differences between the analytically retrieved locking depth and the locking depth retrieved from the synthetics produced by a more realistic unlocked fault thickness (h 2 ≤ 100 km) are very large and deserve a discussion.

The finite space approximation

The differences between the numerical locking depth h 1 , being equal to 12 km, and the values of d, retrieved from the analytical solution are given in Tables 1 and represented on Fig. 3(b). They display about 2% to 57.5% of relative error depending on the thickness on which the localized slip is imposed.

For a 400 km long model, the thickness of the free sliding fault zone, h 2 , must be higher than 100 km to get a depth error inferior to 15% and a slip error inferior to 10%. The best fitting arctangent gives a locking depth of 5.1 km for h 2 =1 0km and 9 km for h 2 =5 0km, instead of the 12 km input in the model, which corresponds to more than 55% of relative error for the former, and to 25% of relative error for the latter. The retrieved slip is 0.94 cm/yr for the h 2 =1 0km model and 1.7 cm/yr for the h 2 =50km model, instead of the imposed 2 cm/yr, which represents a slip relative error of 53 to 15%. For the 800 km long and 1200 km long models, the results are roughly the same, which implies that the numerical models are rather independent of the choice of the width.

In conclusion, it requires to account for a thick zone of localized slip (≥ 200 km) for the retrieved value of the locking depth to correspond to the actual locking depth imposed in the numerical model. When this critical thickness is not reached, the locking depth retrieved using [Table 1 about here.]

Implication for locking depth assessment

A mechanical definition of the locking depth is the depth extent within the crust which behaves elastically during the interseismic period. As this portion is loaded, it accumulates elastic strain and therefore potential energy, which is then released during earthquakes. Part of the seismic risk assessment is to evaluate this depth. Based on the rate of loading and the elapsed time since the last earthquake on a segment, it gives an estimate of the seismic moment accumulation M 0 and therefore a raw estimate of the magnitude M W of the next earthquake.

There are two main ways in evaluating the locking depth. The locking depth may be estimated using the seismogenic zone, i.e. the distribution of micro-seismicity, assuming that this corresponds to the brittle part of the fault. This depth inferred from seismologic data is usually called the seismogenic depth. Alternatively, it is possible to retrieve a geodetic locking depth along each segment of the fault using surface displacement measurements obtained from geodetic data, either GPS or InSAR, across the major faults together with the [START_REF] Savage | Geodetic determination of relative plate motion in Central California[END_REF] model. The standard approach is to fix the slip rate to the long terms velocities deduced from various geological studies, and to use the Savage & Burford model to retrieve the locking depth by a least square fitting approach.

These two approaches in estimating the locking depth are fundamentally very different in their assumptions and therefore should generally reach different results. However, this is not always the case and many studies actually give consistent values using both of these methods to determine the locking depth [START_REF] Smith-Konter | Locking depths estimated from geodesy and seismology along the San Andreas Fault System: Implications for seismic moment release[END_REF]. This consistency explains why the seismogenic depth and the geodetic locking depth are sometimes considered as providing the same information.

However, taking a closer look at the outliers in the data, Table 2, shows that, when the estimates clearly disagree, the locking depth inferred from geodetic data systematically underestimates the seismogenic depth as it was predicted by the finite space model. As these underestimations of the depth of the locked zone have been shown to be up to a factor of two in error (Sec. 2.4), the retrieved geodetic depth should be used carefully, specifically when considering this depth to be representative of that of the seismogenic zone.

[Table 2 about here.] This limitation of the original Savage & Burford model was somewhat accounted in the viscoelastic coupling model of [START_REF] Savage | Asthenosphere readjustment and the earthquake cycle[END_REF]. In this model, right after the earthquake, the sudden stress change causes the viscoelastic substrate to respond elastically to infinite depth and the retrieved velocity field at the surface, when averaged over the interseismic period, corresponds to the exact [START_REF] Savage | Geodetic determination of relative plate motion in Central California[END_REF] solution. However, as time goes on during the interseismic period, part of the elastic strain energy is dissipated into permanent plastic deformation at shallower and shallower depth. As a result, long time after the last large earthquake, the elastic part of the viscoelastic lower layer approaches zero and the strain localization at the surface is limited to that of a finite plate as described in here.

ALTERNATIVE MODELS

Validation criteria

The alternative set of boundary conditions and initial geometry presented in this study form the basis of idealized models which serve to help our understanding of the origin of the localization zone. These models were designed to minimize the number of parameters and to reproduce the general shape of geodetic data across strike slip fault zones. In order to evaluate the first order pertinency of our numerical models, we wish to find criteria derived from general observations that are simple enough to enable an immediate validation.

[Figure 4 The characteristic velocity ratio indicates the ratio of the velocity reached in the central deformation area to the far field velocity. A high characteristic ratio (higher than 80% of the total velocity) reflects then that the velocity reaches an important value in the vicinity of the fault, there-fore the velocity outside this area is not increasing significantly to reach the far field velocity and the velocity field shows a flat repartition far from the fault, of type (a) or (b) on Fig. 4. In this case, the medium can be assumed to behave almost rigidly outside the fault area which is concordant with the assumption of plate rigidity within the context of plate tectonics. Low characteristic ratios of type (c) and (d) break this assumption.

Geologically, only a sharp and flat repartition (type (a)) reflects that the deformation is localized around the fault and that the blocks outside the fault zone behave rigidly. The three others types represent distributed deformations that are not observed near strike-slip fault. In the following, a model is considered as valid upon its ability to meet three criteria, (i) the characteristic width W c must range from 20 to 60 km, (ii) the characteristic velocity ratio R c must be higher than 80%, and (iii) the boundary conditions and the rheological structure must be mechanically consistent with the basic knowledge of long and short term rheology of the crust and lithosphere.

Basal model

Numerous previous studies used basal boundary conditions [START_REF] Fay | Fault slip rates, effects of elastic heterogeneity on geodetic data, and the strength of the lower crust in the Salton Trough region, southern California[END_REF], or combined them with lateral boundary conditions [START_REF] Brink | Three-dimensional models of deformation near strike-slip faults[END_REF][START_REF] Roy | Crustal rheology and faulting at strike-slip plate boundaries: 2. effects of lower crustal flow[END_REF] to model the interseismic velocity of a strike slip fault. We perform numerical simulations in order to evaluate how accounting for simple basal boundary conditions impacts the surface deformation in an elastic model. As during the benchmark, we extract the horizontal surface velocity from the simulations, but here we test different locking depths.

The basal model (Model 1) is a rectangular box with two elastic layers. The fault is locked through the whole upper layer, whose depth, h 1 , varies from 5 km to 20 km. The lower layer is 20 km deep, and the fault can slide freely through this layer. Model 1 is periodic, and the boundary conditions are applied underneath the box, with an opposite velocity field imposed on the left and the right planes (Fig. 5a) of the basal surface.

[Figure 5 about here.]

The shape of the horizontal velocity field at the surface of the model approximates well the arctangent shape, as shown in Fig. 5(b), for h 1 =1 0km. The results presented in Table 3 show the values of the localization width evaluated by following the method of tangents described in Appendix A. For h 1 =10km, the localization zone is about 41.1 km width, which is a reasonable width according to our first criterion. In addition, the velocity becomes almost uniform far from the fault which is shown by a characteristic velocity ratio higher than 98%, the material can be considered then as a rigid block outside the local deformation zone. The model based on basal boundary conditions gives a velocity field that matches two of our validation criteria.

[Table 3

Lateral model

In order to overcome the previous issue, we choose to implement a simple elastic model that is driven laterally (Model 2). The model is periodic, with opposite velocities imposed on the lateral sides of the box (Fig. 6a). This model is intuitive since it considers that the large scale motion of the plates is a driving force inducing a motion on both side of the domain. Because of the nature of these boundary conditions, this kind of model is sensitive to the choice of the width of the model, but not to its thickness. An approach that overcomes this limitation is to impose directly the driving forces on the lateral boundaries, but it would require to invert the values of such forces from the far field velocity field [START_REF] Pollitz | A physical model for strain accumulation in the san francisco bay region[END_REF]. Here we chose to simply apply the far field velocities that are easily deduced from geodetic data, and to ensure that these velocities are applied sufficiently far from the fault (here at a distance of 100 km to the fault).

[Figure 6 about here.]

The geometry of Model 2 is a rectangular box with two distinct layers. The depth of the upper layer, h 1 , varies from 5 km to 20 km, while the thickness of the lower layer is fixed and equal to 150 km. The fault is locked inside the upper layer but is free to slip inside the lower layer. Model 2 is periodic.

The retrieved surface velocity fields are presented in Fig. 6(b). The deformation is distributed almost linearly across the domain and is not well localized along the fault. Although the localization width is reasonable (W c is about 52 km for h 1 =1 0km), as shown by Fig. 6(c) and Table 4, the characteristic velocity ratio is never higher than 64%, that is compatible with neither a rigid plate-tectonic type behavior, nor the geodetic data available on Earth. Therefore we also tested the case where the fault frictional strength is almost zero, considering a very shallow locked fault of 0.5 km. Even in this case, the distribution is not localized but distributed (the velocity ratio being only 71%).

[Table 4 about here.]

Contrarily to the conclusion of [START_REF] Savage | Viscoelastic-coupling model for earthquake cycle driven from below[END_REF], we find that the solution is radically different whether loading is applied as a displacement at the lateral boundary or from below. This strikingly different conclusion is the result of different modeling approaches. [START_REF] Savage | Viscoelastic-coupling model for earthquake cycle driven from below[END_REF] model is indeed built on the [START_REF] Savage | Asthenosphere readjustment and the earthquake cycle[END_REF] model, i.e. using a Burgers vector to load the fault, and then post processing the stress in the far field or at the base of the lithosphere in the models a posteriori. With this semi-analytical approach, strain localization indeed does not depend on where the loading is integrated, it is only the result of applying this load at the tip of the dislocation.

Here, we build a screw dislocation model and apply displacement on the boundary of the model (lateral model, Fig. 6) or at the base (basal model, Fig. 5) and find, like [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF] that when driving stresses are transmitted laterally along the elastic upper layer, the strain hardly localizes self consistently around a screw dislocation.

Indentor model

As in the previous models, the indentor model (Model 3) mimics a strike slip fault, locked down to depth h 1 and free to slip below, embedded inside a box whose base is a square of 200 km long.

In this model however, the motion is driven by the extrusion of an indenting block located at the front of the model on the left side. The indentor is 100 km long normal to the fault and 20 km wide in the direction of the fault. A normal velocity of 2 cm/yr is applied on the indentor front surface, which is perpendicular to the fault plane. In order to ensure the symmetry of the surface velocity field, free slip boundary conditions are applied on the lateral faces of the model (parallel to the fault) and rigid boundary condition are applied at the front of an elastic block that has been added at the front of the model symmetrically to the indentor (transparent gray block on Fig. 7a).

[Figure 7 about here.]

Contrary to the previous models, the back side of the model is a free surface. As a result, this model is not periodic and the surface velocity depends on the distance to the indentor. To allow quantitative comparison with the previous models, the results presented in Fig. 7(b) show the relative velocities with respect to the axis of the fault.

The results in Fig. 7(b) show an appropriate localization of the deformation in the vicinity of the fault and the blocks behave rigidly. For instance, the localization width is about 54 km and the velocity ratio is about 83% (Table 5) at a distance of 40 km from the indenting block, when the depth of the locked fault is 5 km. For a deeper locking depth, the deformation is accommodated over a wider zone (from 65 to 79 km), but is still quite localized (about 80% of the total velocity is reached within the localization zone). The far field velocity decreases with increasing distance from the fault tip (Fig. 7c, from 40 to 160 km) as the deformation is being absorbed elastically by the model along the strike of the fault.

[Table 5 about here.]

Viscoelastic notch model

Model 4 is fully viscoelastic, with a rigid part, in which the viscosity is 10 25 Pa.s, surrounding a notch of low viscosity (10 19 Pa.s) in the lower layer. The model is periodic, i.e. it corresponds to an infinite fault, and is formed by a rectangular box that is 200 km long by 100 km wide. As for Model 2, the velocity imposed laterally amounts to 2 cm/yr and simulates the general plate motion that drags the box on both sides with opposite directions. Contrary to the previous models, the displacement at depth does not take place along a discrete fault. Instead, it is distributed within a broader weak zone which corresponds to the notch.

[Figure 8 about here.] Fig. 9 shows the horizontal velocities obtained for several depths of locked fault, with a notch width L of 50 km and a viscosity µ of 10 19 Pa.s. The results, for the locking depths from 5 to 15 km, show a good localization of the deformation near the fault. The localization width is about 63 km for h 1 =5km and a notch width of 50 km, and the velocity ratio is 87% (Table 6). A shallow notch is therefore more efficient than a deep one to help the localization of the deformation at the surface.

[Figure 9 about here.]

[Table 6 about here.]

Then, we tested different values of the notch's width L, with a locked fault of thickness h 1 =5

km and a viscosity of 10 19 Pa.s. Fig. 10 shows the surface velocity fields for L =1 0 , 50 and 100 km. Narrow notches (L<20 km) induce a localized deformation within a 30 km wide zone, but the velocity ratio does not reach 80% (Table 7), which means that the surrounding of the fault zone does not behave rigidly. Wide notches (L>50 km) do not produce localized deformation (W c > 60 km) and tend to result in a linear velocity field when the notch's width tends to the total width of the model.

For L in the range [20, 50] km, the deformation is well localized near the fault in a 30 to 60 km wide zone, and the velocity inside the localization zone reaches 80% of the far field velocity.

Thus, it is the lateral contrast of viscosity in the fault vicinity, with a lower viscosity beneath the fault, that favors the localization.

[Figure 10 about here.]

[Table 7 about here.]

The models obtained by choosing the parameters h 1 in the range [5; 15] km, L in the range This set up is similar to the Zatman (2000) model because it accounts for lateral reduction of the viscosity in the lower crust beneath the fault (Fig. 8), but, contrary to the [START_REF] Zatman | On steady rate coupling between an elastic upper crust and a viscous interior[END_REF] model, we impose that the viscosity outside the notch is high enough for the elastic stress not to be relaxed at the time scale of interseismic deformation. As a result, the elastic stresses are not relaxed at depth on the lateral part of the model and drive a long term localized ductile shear flow beneath the fault zone. The pattern is in the end extremely similar to [START_REF] Bourne | The motion of crustal blocks driven by flow of the lower lithosphere and implications for slip rates of continental strike-slip faults[END_REF] model, except that the viscous shear flow is driven laterally by the strong part of the lithosphere rather than from below. Like [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF], we find that accounting for a local drop of the elastic thickness of the lithosphere, here through a localized viscous notch, permits strain localization at the surface to be self-consistently driven by stress transmitted along the elastic upper layer.

DISCUSSION

Localization and boundary conditions

To model strain localization around strike slip faults during the interseismic period, most of the models impose the slip on the fault itself, at depth. This type of loading is inconsistent with the coseismic and long-term effective rheology of the lithosphere. In this contribution, we have tested by the mean of numerical experiments alternative modes of loading and studied specifically how they affects shear localization in the vicinity of strike slip faults.

We have found that with screw dislocation type models, strain localization only occurs when the velocity applied on one boundary includes a step function. This step function can be imposed either on the basal plane (model 1) or by an indentor on the lateral boundary (model 3). However, when the shear is applied by a velocity in the far field in order to mimic intuitively a plate tectonic like loading (model 2), the elastic deformation is so important over the whole domain that very little localization occurs in the vicinity of the fault zone. With this type of boundary conditions, the deformation is always distributed unless the frictional strength of the fault vanishes, i.e. when the locking depth is lower than 1 km. The fact that strongly localized deformation cannot be achieved within a screw dislocation model by far field loading unless the fault is not locked at interseismic time scale raises the question whether these models have a geologically relevant rheological structure.

The indentor models (Model 3) might be meaningful in numbers of geological context, such as the North Anatolian fault or some of the Tibetan faults where the old and rigid Arabian and Indian plates respectively might induce a step function as velocity field on the boundary [START_REF] Flérit | The mechanical interaction between the propagating north anatolian fault and the back-arc extension in the aegean[END_REF][START_REF] Peltzer | formation and evolution of strike-slip faults, rifts, and basins during the India-Asia collision: an experimental approach[END_REF][START_REF] Peltzer | formation and evolution of strike-slip faults, rifts, and basins during the India-Asia collision: an experimental approach[END_REF][START_REF] Socquet | India and Sunda plates motion and deformation along their boundary in Myanmar determined by GPS[END_REF]). However, the elastic screw dislocation models with a step function applied at the base of the upper crust (Model 1) suffer from the same limitation as the original Savage & Burford model. They indeed localize the deformation with a similar pattern as what is observed in nature but they present mechanical incompatibility between the strong viscoelastic coupling that is assumed in order to impose velocity boundary conditions at depth and the infinite velocity gradient assumed to occur on the aseismic part of the fault. This mechanical inconsistency can be eliminated if one assumes that the deepest part of the fault is a thin vertical layer made of either a frictional material with low effective friction, or a low viscosity material, which is embedded between blocks of stiffer material. In other words, lateral rheological variations must exist at depth beneath the fault.

In order to account for this type of lateral variations, it is possible to design a model (model 4), which includes a viscous notch beneath the fault. This model is found to be quite efficient to localize the deformation when the notch is shallower than 10 km. The simulations have shown that, for a narrow notch (< 50 km), most of the deformation is accommodated in the vicinity of the fault. However, when the notch is embedded deeper than 10 km, the elastic deformation becomes predominant over the whole domain and the assumption of rigidity of plates is lost.

One should also note again that in the limit of a very thin notch with an extremely low viscosity, both the notch and the screw dislocation model become equivalent and none of them can actually localize the deformation in the vicinity of the fault zone when subjected to far field loading. Therefore, if one posit that the loading of strike fault is imposed in the far field by plate tectonics, one must admit that this loading is always imposed by an indenting block or that the rheological structure of the crust around strike slip faults includes a local zone of weakness, which must be of the approximate same width as the zone on which the steep velocity gradient is observed during the interseismic period. This weakness could correspond either to elastic damage in the brittle part of the crust [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF] or to a local 5 to 50 km wide low viscosity zone in the ductile part of the crust.

Geological evidence for viscous notch

The only class of models, for which laterally driven boundary conditions lead under some conditions, to a surface velocity field that approaches the geodetic observations, are the models in which we introduce a viscous notch in the brittle crust. As it was already noted by [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF][START_REF] Platt | The mechanics of continental transforms: An alternative approach with applications to the San Andreas system and the tectonics of California[END_REF] and [START_REF] Platt | Where is the real transform boundary in California?[END_REF], the high interseismic strain rate can be explained by the low rigidity of the fault zone. By introducing a notch in the model, we locally lower the integrated elastic strength in the vicinity of the fault. In this paper, we have assumed the presence of this low viscosity zone and we have provided quantitative limits on the viscosity and the size of such a weak zone in order to reproduce the general shape of geodetic data.

Yet, one may wonder what is the geological meaning or existence of such a weak zone. This discussion focuses on the case of the San Andreas Fault because it is one of the best-documented strike slip faults. Using numerical modelling, previous studies conclude that significant vertical [START_REF] Rolandone | The distributions of slip rate and ductile deformation in a strike-slip shear zone[END_REF][START_REF] Fay | Fault slip rates, effects of elastic heterogeneity on geodetic data, and the strength of the lower crust in the Salton Trough region, southern California[END_REF] and horizontal [START_REF] Chéry | An integrated mechanical model of the San Andreas fault in central and northern California[END_REF] variations in viscosity provide the strain localization that is observed in the San Andreas fault vicinity. Strain localization of ductile deformation could be due to the effect of shear heating on the local viscosity [START_REF] Rolandone | The distributions of slip rate and ductile deformation in a strike-slip shear zone[END_REF], as also modeled by [START_REF] Dayem | Localization of shear along a lithospheric strength discontinuity: application of a continuous deformation model to the boundary between Tibet and the Tarim Basin[END_REF] for the Altyn Tagh fault zone.

Moreover, magnetotelluric sounding infers the presence of a 40 km wide low electrical resistivity zone at around 20 km depth beneath the San Andreas Fault, which has been interpreted as marking the presence of hot saline fluid-filled porosity network at depth. [START_REF] Becken | Correlation between deep fluids, tremor and creep along the central san andreas fault[END_REF].

Similarly, shear wave splitting measurements shows lateral variation in the vicinity of the San Andreas Fault [START_REF] Bonnin | Uppermantle deformation beneath the north american-pacific plate boundary in california from sks splitting[END_REF][START_REF] Becker | Mantle flow under the western United States from shear wave splitting[END_REF][START_REF] Baldock | Width of mantle deformation across a continental transform: Evidence from upper mantle (Pn) seismic anisotropy measurements[END_REF][START_REF] Savage | Strain modelling, seismic anisotropy and coupling at strike-slip boundaries: Applications in New Zealand and the San Andreas fault[END_REF],

where the deformation zone was again marking a 40 km wide high shear strain zone in the mantle.

In these studies, the data seem to infer a relatively localized low viscosity zone located beneath the active fault zone.

In the ductile layer, at high stress, grain boundary sliding creep and dynamic grain size reduction [START_REF] Precigout | Mantle weakening and strain localization: Implications for the longterm strength of the continental lithosphere[END_REF][START_REF] Gueydan | Grain-size-sensitive flow and shear-stress enhancement at the brittle-ductile transition of the continental crust[END_REF] provide two valid mechanisms to localize strain in zones as narrow as 40 km. In the crust, it is useful to study field analogues of large exhumed strike slip fault zone to understand how localization operates below the seismogenic crust.

Examples of large exhumed strike slip faults are found in the remnant of the variscan orogeny and in the vicinity of the Hymalayas. These strike slip ductile shear zones all display narrow and elongated domal shape [START_REF] Leloup | The Ailao Shan-Red River shear zone (Yunnan, China), Tertiary transform boundary of Indochina[END_REF], for the Red River fault zone) generally attributed to a transpressional regime. Reconstruction of the thermal history of the rocks exhumed in these fault zones as well as structural geology analysis display that the shear strain is localized within elongated domes of exhumed lower crustal material and gneisses [START_REF] Denèle | Formation des dômes gneissiques hercyniens dans les Pyrénées : exemple du massif de l'Aston-Hospitalet[END_REF] for the Pyrenées;

Corsini & Rolland 2009 for the Maures Massif; West & Hubbard 1997 for the Norumbega fault zone). Similar elongated domes are found in Norway along the exhumed strike slip segments of the main detachment [START_REF] Labrousse | Pressure-temperature-time deformation history of the exhumation of ultra-high pressure rocks in the Western Gneiss Region, Norway, in Gneiss Domes in Orogeny[END_REF]. Recently, 3D thermo-mechanical models have found that these domes can form self consistently over a Myr time scale within pure strike slip boundary conditions within the transtensional relay of strike slip faults (Le [START_REF] Pourhiet | Kinematic interpretation of the 3d shapes of metamorphic core complexes[END_REF].

An example of application to the San Andreas Fault System

Here, we focus on the interseismic velocity distribution across the Carrizo segment of the San Andreas Fault in central California. Different viscoelastic models [START_REF] Segall | Integrating geologic and geodetic estimates of slip rate on the San Andreas fault system[END_REF][START_REF] Fulton | Reconciling patterns of interseismic strain accumulation with thermal observations across the carrizo segment of the san andreas fault[END_REF] also chose this segment of the San Andreas Fault and were able to fit the present-day distribution of geodetic data with an appropriate choice of parameters. This segment is chosen because it belongs to the locked portion of the San Andreas fault system, which last ruptured in the 1857 Fort Tejon earthquake. Therefore, it corresponds to the main assumptions made in our purely interseismic model where the transient effects related to postseismic relaxation are negligible compared to the interseismic signal.

We use the geodetic data of the Carrizo segment (USGS public data) to infer the best parameters for our viscoelastic notch model. Fig. 11 presents the GPS fault parallel velocity components given by the stations placed in a 75 km wide band southeast of Cholame, perpendicularly to the fault direction. The best prediction of our notch model provides a 10 km locking depth and a 50 km wide viscoelastic notch, with a reasonable RMS misfit (0.105 cm/yr).

The geodetic data are also overlaid by two solutions of the standard and widely used Savage & Burford model. One, obtained using a locking depth of 12 km and a slip rate of 3.4 cm/yr, leads to a RMS misfit of 0.274 cm/yr. The second curve is the best fitting arctangent in the sense of least square regression and is obtained for an unrealistic locking depth of 32 km, that might be underestimated according to our previous conclusions of Sec. 2.4. If one assumes the locking depth reflects the seismogenic depth, this underestimation leads to depths that have no geological meaning.

[Figure 11 about here.]

Finally, the viscoelastic notch model is shown to reproduce the signal around locked faults system with a similar RMS misfit than the best fitting arctangent of the Savage & Burford model, except that it predicts a more realistic locking depth.

CONCLUSIONS

This contribution questions the geological relevance of the boundary conditions imposed in the [START_REF] Savage | Geodetic determination of relative plate motion in Central California[END_REF] model to reproduce the general form of the velocities observed in the geodetic data. This model and other more complex numerical models, which account for multiple faults (e.g., [START_REF] Brink | Three-dimensional models of deformation near strike-slip faults[END_REF], always impose the localization of the deformation at the surface by imposing a step function in the velocity in the ductile crust. Different boundary conditions for which the loading is applied on the resistant part of the crust have therefore been tested by the mean of numerical modeling.

Indentor type models are found to allow strain to be localized around strike slip faults because, as in the [START_REF] Savage | Geodetic determination of relative plate motion in Central California[END_REF] model, a step function is imposed in the boundary conditions.

However, these provide a limited understanding of the origin of the strain localization.

In order to localize the strain using far field laterally driven boundary conditions, it is necessary to introduce a local mechanical heterogeneity in the crust. [START_REF] Chéry | Geodetic strain across the San Andreas fault reflects elastic plate thickness variations (rather than fault slip rate)[END_REF] However, as was pointed out in the discussion, the screw dislocation model is similar to a notch model in the limit where the notch becomes extremely thin and weak. Therefore, the screw dislocation model remains a valid first approximation to retrieve far field velocity from GPS data. Nevertheless, using the [START_REF] Savage | Geodetic determination of relative plate motion in Central California[END_REF] model to estimate the locking depth from geodetic data has been shown to be problematic. The locking depth retrieved with this approach systematically underestimates the effective depth of locking with up to 50% error, because the finite thickness of the sliding part of the fault is neglected. As a result, seismic moment rates and earthquake magnitudes estimated from geodetic locking depths are probably systematically underestimating the risk assessment.

We regard numerical models to be appropriate, when W c ranges from 20 to 60 km, and when R c ranges from 80% to 100%. The combination of these two criteria ensures that the deformation is accommodated locally around the fault zone, and that the blocks outside this zone can be considered as rigid bodies. The arctangent function fulfills the first criterion for d in the range 6 to 19 km, the second criterion being always true.

APPENDIX B: THE LEAST SQUARE REGRESSION

The least square regression method, used in order to choose the best fitting arctangent, is detailed here. Considering a reference line, lying at the surface of the numerical models perpendicularly to the fault plane, this line is discretized into a set of N points x 1 ,... x N . The deformation of this reference line, in the direction Oy that is parallel to the fault, is represented by a vector, The goal of the method is to find the parameters s atan and d atan that have to be input in the N equations above so as to minimize the quadratic error between the vector V num and the vector V atan given by

e q = 1 N N i=1 (v num i -v atan i ) 2 . (B.2)
We perform the calculation over a rectangular range of values for the couple (s atan ,d atan ) and retrieve the best possible couple that minimizes the error e q .

11 Fault parallel velocity field of the Carrizo segment of the SAF: GPS velocities gathered from USGS public data (solid dots), velocity fields given by the Savage and Burford model with a locking depth d = 32 km and a total slip s = 3.5 cm/yr for the best fit (dashed line) and with a locking depth d = 12 km and a total slip s = 3.4 cm/yr for the average central California parameters (solid thin line) and velocity field given by the viscoelastic notch model when the notch is embedded below the depth h 1 = 10 km, its viscosity is µ= 10 19 Pa.s and its width is L = 50 km, the embedding block is of viscosity 10 25 Pa.s (solid thick line). . The three other cases represent distributed deformations: b) a high characteristic width and a high velocity ratio reflect a large deformation zone with a flat field far from the fault; c) a small characteristic width and a small velocity ratio indicate that the field is sharp in the fault vicinity but the deformation is also accommodated far from the fault; d) the combination of the two previous defaults (high characteristic width and small velocity ratio) leads to a deformation that is only accommodated far from the fault (geologically unrealistic). In our study, we consider that the deformation is localized when the characteristic widths are between 20 to 60 km and the velocity ratios are higher than 80% of the total far field velocity. . Fault parallel velocity field of the Carrizo segment of the SAF: GPS velocities gathered from USGS public data (solid dots), velocity fields given by the Savage and Burford model with a locking depth d = 32 km and a total slip s = 3.5 cm/yr for the best fit (dashed line) and with a locking depth d = 12 km and a total slip s = 3.4 cm/yr for the average central California parameters (solid thin line) and velocity field given by the viscoelastic notch model when the notch is embedded below the depth h 1 = 10 km, its viscosity is µ= 10 19 Pa.s and its width is L = 50 km, the embedding block is of viscosity 10 25 Pa.s (solid thick line).

Table 1. Results of the benchmark: the locking depth (d) and far field velocity (s) retrieved from the best fitting arctangent approximation, the relative errors with respect to the imposed locking depth (12 km) and imposed slip (2 cm/yr) for different thicknesses h 2 of the steadily sliding fault at depth and different widths L 0 of the model.

  least square fitting of the Savage & Burford model are systematically underestimating the locking depth. It is very unlikely that a 200 km discrete plan of localized slip exists under any strike-slip fault on Earth and we therefore conclude that the locking depths retrieved by fitting geodetic data with the original Savage & Burford model systematically underestimate the actual locking depth of the fault.

  about here.] In Appendix A, we define a characteristic width W c and a characteristic velocity ratio R c that are used for validation. The characteristic width indicates the width of the central deformation area. A small characteristic width will then give a sharp repartition in the fault vicinity, of type (a) or (c) on Fig. 4. In order to match the geodetic observations, we accept characteristic widths between 20 to 60 km.

[

  20; 50] km and µ =1 0 19 Pa.s for 10 25 Pa.s outside meet the criterion of a velocity field coherent with the observations, with a sharp gradient near the fault and a rigid blocks behavior around. It results that the notch width L must be large enough to make the blocks behave rigidly, but small enough to keep a sharp velocity gradient in the fault vicinity. The lateral contrast of viscosity plays a key role by weakening the fault zone locally and inducing the deformation of the lithosphere around.

  approach consisted in reducing the effective elastic thickness of the crust by dropping the elastic coefficients, we instead introduce a viscous notch at the base of the fault. The notch model constitutes a possible alternative model to the screw dislocation model because (i) it provides an acceptable fit to the geodetic data, (ii) it is compatible with long term deformation along strike slip faults and (iii) its geometry and boundary conditions are justified by geophysical and geological evidences. The notch model has the advantage over the screw dislocation models and the indentor models, because it does not produce unrealistic displacement profiles (step function) at the brittle ductile transition.

V

  num , whose components are the velocities of the N points, v num i , i ∈{1,N}, extracted from the numerical simulations. Let's call V atan =( v atan i ) i∈{1,N } the interpolation of the function defined by Eq. (1) on the discrete set {x 1 ,... x N }, i.e.

Figure 1 .Figure 2 .

 12 Figure1. Representation of an arctangent function of parameters d (the locking depth) and s (the total slip rate), and its characteristic lines (the two asymptotes and the tangent at the origin). The localization width of an arctangent function is πd, its characteristic velocity ratio is 100%.

Figure 3 .Figure 4 .

 34 Figure 3. Benchmark models: a) geometry of the numerical models; b) relative errors between the numerical input parameters (12 km and 2 cm/yr) and the parameters retrieved from the analytical solution (d et s) for models varying in width (L 0 = 400 -800 -1200 km) and in sliding fault thickness (h 2 = 10 -20 -30 -40 -50 -100 -200 -1200 km).

Figure 5 .

 5 Figure 5. Elastic Model 1: basal boundary conditions. a) Geometry of model; b) horizontal velocity of the reference line for the locking depths h 1 =5, 10, 15 and 20 km; c) horizontal velocity and tangent lines for the locking depth h 1 = 10 km, in this case, the localization width is about 41 km and the velocity ratio is 98.5%.

Figure 6 .

 6 Figure 6. Elastic Model 2: laterally driven model. a) Geometry of the model; b) horizontal velocity obtainedfor the locking depths h 1 =0.5, 5, 10, 15 and 20 km; c) horizontal velocity and tangent lines for the locking depth h 1 = 10 km, in this case, the localization width is about 52 km and the velocity ratio is 60.5%.

Figure 7 .

 7 Figure 7. Elastic Model 3: indentor model. a) Geometry of the model; b) horizontal velocity obtained for the locking depths h 1 =5 , 10, 15 and 20 km, at the distance y = 80 km from the fault tip; c) horizontal velocity obtained for the locking depth h 1 =5km at the distances y = 40, 80, 120 and 160 km from the fault tip.

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. Viscoelastic model: a notch of lower viscosity µ and of width L is embedded inside a viscoelastic box below the depth h 1 down to the depth h 1 + 150 km. The boundary conditions are applied laterally.

  Table 2. locking depth inferred from geodetic data versus locking depth inferred from seismologic data for the North Anatolian fault, the San Andreas fault and the Dead Sea fault. NAF SAF DSF Geodetic depth (km) 2-10 a 9 d ; 2-6 e ; 14 f ;6 g 8-15 i Seismogenic depth (km) 17 b ; 20 c 10-20 h 20-32 j a From Le Pichon et al. (2003) b From Bulut et al. (2009) c From Ito (2002) d From Savage & Burford (1973) e From Chinnery (1961) f From Smith-Konter et al. (2011) g From Smith-Konter et al. (2011), average of the outliers h From Nazareth & Hauksson (2004) and Lin et al. (2007) i From Al Tarazi et al. (2011)

	Steadily sliding fault thickness h 2 (km) 10	20	30	40	50	100 200 1200
	L 0 = 400 km
	d (km)	5.1	6.7	7.6	8.4	9.0 10.4 11.3 11.4
	s (cm/yr)	0.94 1.31 1.50 1.62 1.70 1.87 1.95 1.97
	depth error (%)	57.5 44.2 36.7 30.0 25.0 13.3 5.8	5.0
	slip error (%)	53.0 34.5 25.0 19.0 15.0 6.5	2.5	1.5
	L 0 = 800 km
	d (km)	5.5	6.8	7.7	8.4	9.0 10.4 11.5 11.7
	s (cm/yr)	0.95 1.31 1.50 1.62 1.70 1.87 1.96 1.99
	depth error (%)	54.2 43.3 35.8 30.0 25.0 13.3 4.2	2.5
	slip error (%)	52.5 34.5 25.0 19.0 15.0 6.5	2.0	0.5
	L 0 = 1200 km
	d (km)	5.4	6.7	7.7	8.4	9.0 10.7 11.5 11.7
	s (cm/yr)	0.95 1.31 1.50 1.62 1.70 1.88 1.96 1.99
	depth error (%)	55.0 44.2 35.8 30.0 25.0 10.8 4.2	2.5
	slip error (%)	52.5 34.5 25.0 19.0 15.0 6.0	2.0	0.5

j From

[START_REF] Aldersons | Lower-crustal strength under the dead sea basin from local earthquake data and rheological modeling[END_REF]

Table 3 .

 3 Characteristic width W c and velocity ratio R c estimated from the method of tangents for the basal model (Model 1), for different locking depths h 1 .

	h 1 (km)	5	10	15	20
	W c (km) 35.4 41.1 46.6 51.9
	R c (%)	99 98.5 98	97

Table 4 .

 4 Characteristic width W c and velocity ratio R c estimated from the method of tangents for the laterally driven model (Model 2), for different locking depths h 1 .

	h 1 (km)	0.5	5	10	15	20
	W c (km) 36.7 43.5 51.7 58.7 64.2
	R c (%)	71	64 60.5 58.5 56.5

Table 5 .

 5 Characteristic width W c and velocity ratio R c estimated from the method of tangents for the indentor model (Model 3), for different locking depths h 1 .

	h 1 (km)	5	10	15	20
	W c (km) 54.2 64.6 72.6 78.5
	R c (%)	83	81 79.5 78

Table 6 .

 6 Characteristic width W c and velocity ratio R c estimated from the method of tangents for the notch model (Model 4) for a notch width L = 50 km and when the notch is embedded below different depths h 1 .

	h 1 (km)	5	10	15	20
	W c (km) 62.7 69.5 74.0 77.6
	R c (%)	87	80 74.5 70.5

Table 7 .

 7 Characteristic width W c and velocity ratio R c estimated from the method of tangents for the notch model (Model 4) for different notch widths L and when the notch is embedded below the depth h 1 =5km.

	L (km)	10	20	30	40	50	75	100	150
	W c (km) 30.2 34.5 43.7 53.1 62.7 86.6 110.2 155.5
	R c (%) 72.5 78	82	85	87 90.5 92.5	94.5

ACKNOWLEDGMENTS

We acknowledge funding from Émergence-UPMC-2010 research program, and two anonymous reviewers for their constructive comments and suggestions. We thank the USGS for giving access to their GPS data and the CEA for providing the software Cast3M.

APPENDIX A: THE CHARACTERISTIC WIDTH AND THE VELOCITY RATIO

Some elementary mathematics can be used to retrieve the characteristic width and the characteristic velocity ratio. On the fault parallel velocity field, three lines of interest can be identified: two asymptotes when |x|→∞and the tangent line at the origin (Fig. 1). The tangent line crosses the two asymptotes at two points P 1 and P 2 . If the equation describing the asymptotes are of the form y = a 1 x + b 1 and y = a 2 x + b 2 respectively, and the tangent line is of the form y = a 0 x + b 0 , the components of P 1 and P 2 are

and the characteristic velocity ratio R c by

where v f is the far field velocity.

Considering the case of an arctangent function, of the form y = -s/π arctan(x/d), the tangent line at the origin is of the form y = -sx/πd and crosses the two horizontal asymptotes at y = s/2 and y = -s/2 respectively. Inserting a 0 = -s/πd, b 1 = s/2 and b 2 = -s/2 and b 0 = a 1 = a 2 =0 into Eq. (A.1) and Eq. (A.2), we can conclude that the characteristic width is given by

and that the characteristic ratio is given by

Representation of an arctangent function of parameters d (the locking depth) and s (the total slip rate), and its characteristic lines (the two asymptotes and the tangent at the origin). The localization width of an arctangent function is πd, its characteristic velocity ratio is 100%. 2

Common boundary conditions of all the numerical models: free slip boundary conditions. The deeper plane is locked in the vertical direction, the lateral sides that are parallel to the fault plane are locked in the perpendicular direction of the fault plane. 3

Benchmark models: a) geometry of the numerical models; b) relative errors between the numerical input parameters (12 km and 2 cm/yr) and the parameters retrieved from the analytical solution (d et s) for models varying in width (L 0 =400-800 -1200 km) and in sliding fault thickness (h 2 =10-20 -30 -40 -50 -100 -200 -1200 km). 4

Validation criteria: a) a localized deformation is characterized by a small characteristic width W c (the curve is sharp in the fault vicinity) and a high velocity ratio R c (the curve is flat outside the fault vicinity). The three other cases represent distributed deformations: b) a high characteristic width and a high velocity ratio reflect a large deformation zone with a flat field far from the fault; c) a small characteristic width and a small velocity ratio indicate that the field is sharp in the fault vicinity but the deformation is also accommodated far from the fault; d) the combination of the two previous defaults (high characteristic width and small velocity ratio) leads to a deformation that is only accommodated far from the fault (geologically unrealistic). In our study, we consider that the deformation is localized when the characteristic widths are between 20 to 60 km and the velocity ratios are higher than 80% of the total far field velocity. 5

Elastic Model 1: basal boundary conditions. a) Geometry of model; b) horizontal velocity of the reference line for the locking depths h 1 =5 , 10, 15 and 20 km; c) horizontal velocity and tangent lines for the locking depth h 1 =10km, in this case, the localization width is about 41 km and the velocity ratio is 98.5%. 6

Elastic Model 2: laterally driven model. a) Geometry of the model; b) horizontal velocity obtained for the locking depths h 1 =0 .5, 5, 10, 15 and 20 km; c) horizontal velocity and tangent lines for the locking depth h 1 =10km, in this case, the localization width is about 52 km and the velocity ratio is 60.5%. 7

Elastic Model 3: indentor model. a) Geometry of the model; b) horizontal velocity obtained for the locking depths h 1 =5 , 10, 15 and 20 km, at the distance y =8 0km from the fault tip; c) horizontal velocity obtained for the locking depth h 1 =5km at the distances y =40, 80, 120 and 160 km from the fault tip. 8

Viscoelastic model: a notch of lower viscosity µ and of width L is embedded inside a viscoelastic box below the depth h 1 down to the depth h 1 +150km. The boundary conditions are applied laterally. 9

Viscoelastic model (Model 4): horizontal velocity obtained with a notch embedded below the depths h 1 = 5 to 20 km, the notch width is L = 50 km and its viscosity is µ = 10 19 Pa.s. The embedding material is of viscosity 10 25 Pa.s. 10 Viscoelastic model: horizontal velocities obtained with the notch widths L = 10, 50 and 100 km, the notch is embedded below the depth h 1 = 5 km and its viscosity is µ = 10 19 Pa.s. The embedding material is of viscosity 10 25 Pa.s.

LIST OF TABLES

1

Results of the benchmark: the locking depth (d) and far field velocity (s) retrieved from the best fitting arctangent approximation, the relative errors with respect to the imposed locking depth (12 km) and imposed slip (2 cm/yr) for different thicknesses h 2 of the steadily sliding fault at depth and different widths L 0 of the model. 2 locking depth inferred from geodetic data versus locking depth inferred from seismologic data for the North Anatolian fault, the San Andreas fault and the Dead Sea fault.

3

Characteristic width W c and velocity ratio R c estimated from the method of tangents for the basal model (Model 1), for different locking depths h 1 . 4

Characteristic width W c and velocity ratio R c estimated from the method of tangents for the laterally driven model (Model 2), for different locking depths h 1 . 5

Characteristic width W c and velocity ratio R c estimated from the method of tangents for the indentor model (Model 3), for different locking depths h 1 . 6

Characteristic width W c and velocity ratio R c estimated from the method of tangents for the notch model (Model 4) for a notch width L =5 0km and when the notch is embedded below different depths h 1 . 7

Characteristic width W c and velocity ratio R c estimated from the method of tangents for the notch model (Model 4) for different notch widths L and when the notch is embedded below the depth h 1 =5km.