Weak ductile shear zone beneath a major strike-slip fault: Inferences from earthquake cycle model constrained by geodetic observations of the western North Anatolian Fault Zone

GPS data before and after the 1999 İ zmit/Düzce earthquakes on the North Anatolian Fault Zone (Turkey) reveal a preseismic strain localization within about 25 km of the fault and a rapid postseismic transient. Using 3-D ﬁ nite element calculations of the earthquake cycle in an idealized model of the crust, comprising elastic above Maxwell viscoelastic layers, we show that spatially varying viscosity in the crust can explain these observations. Depth-dependent viscosity without lateral variations can reproduce some of the observations but cannot explain the proximity to the fault of maximum postseismic velocities. A localized weak zone beneath the faulted elastic lid satisfactorily explains the observations if the weak zone extends down to midcrustal depths, and the ratio of relaxation time to earthquake repeat time ranges from ~0.005 to ~0.01 (for weak-zone widths of ~24 and 40 km, respectively) in the weakened domain and greater than ~1.0 elsewhere, corresponding to viscosities of ~10 18±0.3 Pa s and greater than ~10 20 Pa s. Models with sharp weak-zone boundaries ﬁ t the data better than those with a smooth viscosity increase away from the fault, implying that the weak zone may be bounded by a relatively abrupt change in material properties. Such a change might result from lithological contrast, grain size reduction, fabric development, or water content, in addition to any effects from shear heating. Our models also imply that viscosities inferred from postseismic studies primarily re ﬂ ect the rheology of the weak zone and should not be used to infer the mechanical properties of normal crust.


Introduction
Earthquakes occur when the internal elastic shear stress exceeds the frictional strength of an appropriately oriented fault. The internal stress at any time throughout the earthquake cycle is the outcome of a complex balance between far-field loading, the sudden redistribution of stress caused by the earthquake, and viscous creep. The first two factors may be simply parameterized using an elastic model of the crust and are quite accurately constrained by seismological and geodetic measurements. The viscous relaxation, however, depends sensitively on the poorly known viscosity structure within the crust and upper mantle and on the loading history. however, regarding the rheological model that should be used to explain postseismic and/or interseismic crustal displacements in tectonically active regions. Various rheological models have been proposed, including non-Newtonian viscosity [e.g., Freed and Bürgmann, 2004;Fialko, 2012, 2013], the standard linear solid model [e.g., Cohen, 1982;Hetland and Hager, 2005;Ryder et al., 2007], and the biviscous Burger's model [e.g., Pollitz, 2003[e.g., Pollitz, , 2005Hager, 2005, 2006;Pollitz et al., 2006Pollitz et al., , 2008Hearn et al., 2009;Ryder et al., 2011]. Nevertheless, most of these studies (except Hearn et al. [2009] and Takeuchi and Fialko [2012]) have addressed only the postseismic displacements, because there are few cases where both interseismic and postseismic geodetic strain rate measurements with adequate spatial and temporal resolution are available.
Most previous models used to explain earthquake cycle deformation assume that the crust consists of a constant-viscosity layer beneath an elastic lid or consists of constant-viscosity layers. To fit the observations using such models typically requires transient rheologies that represent multiple creep processes. Here we test an alternative hypothesis: spatial variations in the viscosity of the crust can explain the geodetic observations in the context of a Maxwell-type viscoelasticity in which a single viscous creep mechanism responds to the elastic stress field.
In this study, using 3-D finite element calculations, we examine the linear Maxwell viscoelastic response to periodically repeating strike-slip faulting events under a constant far-field loading rate. We consider three basic viscosity models with uniform viscosity (UNV), depth-dependent viscosity (DDV), and a localized weak zone (LWZ) in order to evaluate the impact of viscosity variation on the surface velocity field during the earthquake cycle. Strike-perpendicular profiles of the strike-parallel surface displacement rates before and after an earthquake depend primarily on the ratio of the Maxwell relaxation time to the earthquake cycle period [Savage and Prescott, 1978], but a spectrum of relaxation times may be present because of the variation of the viscosity in the crust. On the basis of these results, we analyze the İzmit and Düzce data sets which Hearn et al. [2009] previously modeled using a transient rheology model in order to test whether a simple spatial variation of Maxwell-type viscosity can explain both postseismic and preseismic creep rates measured across the NAFZ.

Geodetic Observations From the Western North Anatolian Fault Zone
The North Anatolian Fault Zone (NAFZ) in Turkey, a large continental strike-slip fault system that represents the boundary between the Eurasian Plate and Anatolia, is one of the major seismogenic zones in the world, extending more than 1000 km from the Gulf of Saros in the west to the Karhova triple junction in the east (Figure 1) [e.g., Stein et al., 1997;Şengör et al., 2005].
The long-term relative displacement rates along the NAFZ show significant regional differences: 15-25 mm/yr across the eastern NAFZ [Wright et al., 2001b;Walters et al., 2011;Tatar et al., 2012], 20-25 mm/yr across the central NAFZ [Reilinger et al., 2006;Kozacı et al., 2007Kozacı et al., , 2009, and 25-27 mm/yr across the western NAFZ [Reilinger et al., 2000[Reilinger et al., , 2006. The westward increase is part of a regional displacement rate field in which the central and southern Aegean move at~30 mm/yr relative to stable Eurasia [Reilinger et al., 2006]. Instrumental, paleoseismological, and historical seismic data exist to constrain the earthquake cycle to be~150-280 years [Stein et al., 1997;Klinger et al., 2003;Rockwell et al., 2009Rockwell et al., ]. 2001a. Repeated measurements of surface positions using GPS have been used to determine the displacement rate (velocity) fields near the NAFZ before and after the earthquakes [e.g., McClusky et al., 2000;Reilinger et al., 2000Reilinger et al., , 2006Ergintav et al., 2009]. Figure 2 shows the strike-perpendicular profile of strike-parallel surface velocities before and after the 1999 İzmit and Düzce earthquakes. The velocities before the earthquakes were observed during the period 1988-1997[McClusky et al., 2000. The postseismic velocities were estimated by Ergintav et al. [2009] from fits to postseismic GPS time series for the 6 month period following the İzmit earthquake. Because we are at this stage considering the earthquake cycle on infinitely long strike-slip faults, we only use data from the central part of the earthquake ruptures; we project the fault-parallel component of all the GPS sites within 50 km distance of 30.27°E onto the profile. Negative and positive distances from the fault indicate the northern (N) and southern (S) parts of the profile, respectively; positive velocities are to the east.
It is clear that (1) before the earthquakes, displacement rates vary fairly smoothly across a relatively broad region, but the velocity gradient is significantly higher by a factor of~3 within about 25 km of the fault than it is outside of this zone, and (2) 6 months after the earthquakes, the difference in peak velocities either side of the fault (~120-180 mm/yr) is 5-8 times higher than the long-term relative displacement rate. This long-term rate is about 22 mm/yr based on the Anatolia-Eurasia Euler pole calculated by Reilinger et al. [2006] and evaluated at 40°N, 30.25°E, but the grey zones in Figures 2a and 2b show the range of relative displacements 22.5 ± 7.5 mm/yr to facilitate comparison of the two sets of measurements. We ask whether this contrast between preseismic and postseismic relative displacement rates can be explained by any fixed crustal viscosity distribution under conditions of steady far-field loading and periodic fault slip.

The Earthquake Cycle Model
We use a parallelized 3-D finite element code, oregano_ve [Yamasaki and Houseman, 2012a], to solve the linear Maxwell viscoelastic problem in which an instantaneous strike-slip fault occurs at a periodic interval, releasing elastic strain accumulated by regional loading that occurs at a constant rate ( Figure 3). The twolayered model is composed of a faulted elastic layer overlying a viscoelastic layer in which the viscosity parameter may vary with position. The implementation of the fault in the 3-D finite element code by Yamasaki and Houseman [2012a] uses the split node method developed by Melosh and Raefsky [1981].
In most of our calculations, the solution domain has a thickness Z L ′ (= Z L /L 0 ) = 1.0 in vertical direction and covers a horizontal region of X L ′ (= X L /L 0 ) = 25.0 × Y L ′ (= Y L /L 0 ) = 10.0 (Figure 3a), where the reference lengthscale L 0 is assumed to be 40 km, corresponding to the thickness of the crust beneath the İzmit and Düzce region in the western NAFZ [e.g., Vanacore et al., 2013]. We later describe the sensitivity of our solutions to these choices that define the domain extent, comparing experiments with Z L ′ = 1.0 or 2.0 (40 km versus 80 km layer thickness) and those with Y L ′ = 10.0 or 20.0 (400 km versus 800 km across-strike distance).
At the base of the layer, z′ = Z L ′ (= 1 or 2), vertical displacement (U z ′ = U z /S 0 , where S 0 is a reference displacement) and horizontal tractions are set to zero. On the top surface (z′ = 0.0) and at the ends of the box, x′ (= x/L 0 ) = ± 12.5, all components of traction are set to zero. Gravitational forces are omitted from these calculations. The length of the model, X L ′, is great enough that the solutions are close to two-dimensional in spite of the zero traction conditions on the ends, which introduce a small x dependence with X L ′ doubled the maximum postseismic velocity is changed by~0.04%.
Journal of Geophysical Research: Solid Earth 10.1002/2013JB010347 Bürgmann et al., 2002;Feigl et al., 2002], although rupture depths estimated using a layered elastic model exceed 20 km [Hearn et al., 2002;Hearn and Bürgmann, 2005]. H′ is assumed to be 0.3 (H = 12 km), slightly thinner than the depth extent of the fault, so that the elastic layer is entirely within the depth range for which fault slip is constant (± 1).
A constant strike-parallel displacement rate is applied on the boundaries y′ = ± Y L ′/2 of the modeled box. Time (t) is nondimensionalized by the earthquake cycle period Δt: at time t′ the dimensionless displacement on the boundaries y′ = ± Y L ′/2 is U x ′ = ± V xb ′ t′, U y ′ = 0, and U z ′ = 0. Consistent with V xb ′ = 1, the fault slip event is repeated with periodicity Δt′ = 1 (Figure 4).
Beneath the elastic lid, the viscosity is defined by We consider two types of localized weak zone (LWZ) model based on the idea that a long slip history may have induced a strain-dependent weakening in a region that is centered on the fault. In the first LWZ (block)   (2)), with thickness Θ B ′ and width Ω B ′ (Figure 5b). For the Gaussian-type LWZ model the vertical and horizontal variation of dimensionless viscosity are defined by equation (3) for effective thickness Θ G´a nd width Ω G ′ parameters ( Figure 5c). The distribution of weakening factor for LWZ models is shown here for f w = 100.

Journal of Geophysical Research: Solid Earth
10.1002/2013JB010347 model, the viscoelastic layer is subdivided into two domains by a sharp boundary represented by a Heaviside function, with (2) and f = 1, elsewhere (Figures 3a, 3b, and 5b). Thus, Ω B ′ (= Ω B /L 0 ) is the horizontal width of the weak zone, Θ B ′ (= Θ B /L 0 ) is its vertical thickness measured from the base of the elastic layer, and f w is the ratio of viscosities (background: weak zone), which we refer to as the weakening factor.
The second LWZ (Gaussian) model avoids the sharp boundaries shown in Figure 3a by defining where f w is the maximum weakening factor which applies at the top of the viscoelastic layer on the fault (Figures 5a and 5c), i.e., at (y′, z′) = (0, H′). In order to compare results from Gaussian and block-type LWZ models, the Gaussian width parameters σ y and σ z can be related to the thickness Θ G ′ and width Ω G ′ of the zone in which f is greater than f w /2 (based on the half width at half maximum measured from the center point of the Gaussian function) In our tests with the LWZ model, we considered maximum weakening factors f w of 10, 100, 200, or 1000, resulting in minimum viscosity η w ′ = 0.1, 0.01, 0.005, or 0.001, and considered values of Ω B ′ = 0.2, 0.4, 0.6, 1.0, and 2.0, Ω G ′ = 0.1, 0.2, 0.4, 1.0, and 2.0, and Θ B ′ = Θ G ′ = 0.1, 0.2, 0.3, 0.5, and 0.7. The LWZ model is similar to the viscoelastic inclusion models previously investigated by Cohen and Kramer [1984] and Traoré et al. [2014]. Table 1 provides a summary of symbol definitions and relevant parameter values and ranges.

UNV Model
In general, the experiments here evolve through an initial transient phase before establishing a purely periodic behavior. A point on the upper surface moves abruptly when the fault slips at t n = nΔt (for n a nonnegative integer), and between those times its velocity smoothly changes as viscoelastic relaxation occurs. After the transient phase has decayed, a plot of (U x ′ À t′) versus t′ should appear as a simply periodic function.
For example, Figure 6 shows the temporal evolution of (U x ′ À t′) at the surface point (x′, y′, z′) = (0.0, À0.5, 0.0) for UNV models (c = 0.0) with η′ (= η b ′) = 0.01, 0.1, 1, and 10, with Y L ′ = 10, and Z L ′ = 1. The model behavior depends on the dimensionless viscosity η′ = η/(μΔt) = τ/Δt, which is simply interpreted as the ratio of the Maxwell time constant (τ = η/μ) to the interseismic period Δt. For the models with η′ = 0.01 and 0.1 (Figures 6a  and 6b), a cycle invariance is established rapidly, implying that the relaxation of the deviatoric stress induced by fault slip is virtually complete by the end of each interseismic period, and we see the renewed accumulation of elastic strain in the linear change of displacement toward the end of each cycle.
For the models with η′ = 1 and 10 ( Figures 6c and 6d), on the other hand,~15 and~160 earthquake cycles are, respectively, required for cycle invariance to establish. It can be seen in Figures 6c and 6d that the background displacement level increases with time during the transient phase, diagnostic of an increase in the stored elastic stress. For these larger viscosities, only a smaller proportion of the deviatoric stress can be relaxed between seismic events, and the unrelaxed stresses evolve until the cycle invariance is established. The almost constant rate of change of displacement within each earthquake cycle, due primarily to the farfield loading, is now partly offset by viscous relaxation throughout each period. The amplitude of the periodic signal after cycle invariance is achieved depends not only on which point is measured (amplitudes will decay away from the fault) but also on how much of the elastic stress is relaxed during each cycle. These results are consistent with the well-known results of Savage and Prescott [1978].
horizontal length of the model 1000 km, 400 (or 800) km viscosity at the base of the model a factor controlling the viscosity gradient f spatially variable weakening factor f w maximum weakening factor 10-1000 Θ Β thickness of block weak shear zone Figure 7 shows the instantaneous displacement rate (V x ′) profiles along x′ = z′ = 0.0 for UNV models with η′ = 0.01, 0.1, 1, and 10, respectively. The V x ′ profiles immediately before and immediately after an earthquake are compared with the simple linear variation predicted for an elastic (or viscous) sheet under steadily increasing boundary displacement. Immediately after an earthquake, large V x ′ is predicted near the fault  6 8 10 12 14 16 18 20 22 24 26 28 30  0  20 40 60 80 100 120 140 160 180 200 Displacement  (y′ = ±~0.7) for models with η′ = 0.01 and 0.1 (respectively,~16 and~1.8 times higher than the boundary velocity V xb ′) (Figures 7a and 7b). The distance from the fault at which the highest V x ′ is obtained depends primarily on the depth extent of the fault relative to the thickness of the layer. With further distance from the fault V x ′ decays smoothly. The key measures of these numerical experiments later used in comparison with measures from the NAFZ (postseismic peak V x ′, distance to peak V x ′, and preseismic strain rates) are summarized in Table 2.
The high displacement rates decay rapidly for η′ = 0.01 and less rapidly for η′ = 0.1, but, in either case, the displacement rates have returned almost to the linear profile immediately prior to the next earthquake, that is, the relaxation of the stress perturbation caused by the fault slip event is virtually complete (Figures 7a and  7b). For both viscosities, there remains at the end of the cycle no significant indication of strain rate localization related to the fault.
In contrast, the models with η′ = 1 and 10 show only small (negligible for η′ = 10) time variation of the strain rate distribution during the course of the earthquake cycle, that is, little deviatoric stress is relaxed between successive slip events, but, throughout the cycle, they show a strong localization of the strain rate in a narrow region of width a few times the thickness of the elastic layer centered on the fault (Figures 7c and 7d). Outside of this shear zone, the displacement rate of a point on the surface approaches the boundary velocity, and strain rates in this external zone may be smaller by a factor of 10 or more than the strain rates in the near-fault shear zone (Table 2).
The UNV model behavior during the earthquake cycle thus depends primarily on the ratio (τ/Δt) of the Maxwell time constant to the earthquake repeat time as first described by Savage and Prescott [1978]. A high ratio of postseismic to preseismic displacement rates is obtained if τ is less than~0.1Δt, but no localization of strain is found late in the cycle. In contrast, a clear localization of strain occurs throughout the earthquake cycle if τ is comparable to or greater than Δt, but the differences between postseismic and preseismic displacement rates then are small.
In order to examine the impact of domain boundaries, we repeated the UNV model calculations of Figure 7 by separately doubling the cross-strike width of the domain (Figures 8a-8d) and the depth to the lower boundary (Figures 8e-8h). As Figures 8a and 8b show, we see the impact of the lateral boundary conditions: for small Maxwell time constant, postseismic velocity decreases from peak values (that are relatively unchanged) at a slower rate with increasing distance in the far field, and the preseismic velocity gradient is effectively halved as the width of the domain is doubled. Similarly, in the postseismic and preseismic velocity a Postseismic measures are obtained from cycle-invariant solutions immediately following an earthquake and preseismic measures immediately before an earthquake. Model measurements may be compared with dimensionless estimates from the NAFZ (Figure 2) of postseismic peak V x ′ in the range 5 to 8, distance to peak in the range 0.25 to 0.5 (measurements from the 6 month period following the İzmit earthquake), preseismic peak gradient of~1.30 and preseismic far-field gradient in the range 0.34 to 0.45 (measurements from the decade preceding the major earthquakes of 1999 profiles for large Maxwell time constant, the velocity gradient in the near-fault zone is relatively unchanged, but the gradient is significantly decreased in the far-field zone (Figures 8c and 8d). These changes in gradient in the far-field zone all seem consistent with the idea that a larger far-field zone allows lower elastic strain rates to be distributed over a greater area. The impact of domain depth is much less, mainly being evident in a change of velocity gradient in the far-field zone for the UNV models with larger viscosity (Figures 8g and 8h).  Figure 7, but the blue curves are for models in which Y L ′ = 20, with other parameters unchanged (Figures 8a-8d), or Z L ′ = 2, with other parameters unchanged (Figures 8e-8h). The green curves show the analytic solution of Savage and Prescott [1978], which assumes lateral and depth boundaries at infinity beneath a faulted elastic lid of thickness 0.4 (we note here that the Maxwell time constant used by Savage and Prescott [1978] is twice that which we use).

10.1002/2013JB010347
We also compare in Figure 8 our finite element solutions with solutions for an infinite viscoelastic half space beneath a faulted elastic layer obtained using the method of Savage and Prescott [1978]. This comparison validates the accuracy of our calculations in the near-fault zone, showing that peak postseismic velocity and near-fault strain rates are not strongly affected by the distance to lower and lateral boundaries in our calculations. It is also clear that the lateral boundaries can have a significant effect on the far-field strain rates. We note, however, that our calculations with Y L ′ = 10 and Z L ′ = 1 agree well with the half-space model for low viscosities (Figures 8a and 8b) and that present data would be unlikely to discriminate between the small farfield differences predicted for high viscosities (Figures 8c and 8d). In the following sections, we therefore concentrate on models for which Y L ′ = 10 and Z L ′ = 1. Figure 9 shows comparable velocity profiles (preseismic and postseismic) for the DDV model with viscosities at the top of the viscoelastic layer (η h ′) between 10 1 and 10 5 times greater than basal viscosities of η b ′ = 0.001, 0.01, and 0.1, again with Y L ′ = 10 and Z L ′ = 1 (see also Table 2 for summary of key measures). Contrasting the profiles after an earthquake (Figures 9a, 9c, and 9e) with those before an earthquake (Figures 9b, 9d, and 9f), we see again the large contrast between preseismic and postseismic velocities that is developed when the crustal viscosities are relatively low (e.g., η b ′ = 0.001 and 0.01) and the smaller contrast when the crustal viscosities are relatively high (e.g., η b ′ = 0.1). As expected perhaps, we see that the effect of increasing η h ′ with constant η b ′ is to reduce the contrast between preseismic and postseismic profiles. Increasing η b ′ with constant η h ′ also reduces that contrast.

DDV Model
After an earthquake, peak V x ′ is significantly greater than V xb ′ for DDV models with η b ′ = 0.001 and 0.01 and relatively large viscosity contrasts (Figures 9a and 9c), but the location of the peak V x ′ measurement is now significantly further from the fault (y′ = ±~1.0-1.5) than in comparable UNV calculations (y′ = ±~0.5-1.0 in Figures 7a and 7b), indicating that this location is sensitive to the viscosity gradient. In contrast, peak V x ′ is hardly greater than V xb ′ for models with η b ′ = 0.1 (Figure 9e).
Before an earthquake, the strain rate localization (V x ′ gradient) near the fault is significant for any models with η h ′ >~1.0 (see Figures 9b, 9d, and 9f). Consistent with our observation of the UNV calculations, this higher strain rate zone persists because only a small proportion of the deviatoric stress can be relaxed in that upper part of the layer where viscosities are greater than~1. We see that the velocity gradient across the fault before an earthquake increases with η h ′ (compare Figures 9b, 9d, and 9f). Although the DDV model with large viscosity contrast comes closer to predicting the data shown in Figure 2 than did any of the UNV models, difficulties remain: for those parameter sets that predict the required high displacement rates in the postseismic phase, the rates remain too high at distances from the fault greater than 25 km, and the localization of strain near the fault in the preseismic phase is insufficient.

LWZ Model
We now consider the impact of a localized weak zone imposed on a constant background viscosity (η b ′ = 1, c = 0), again with Y L ′ = 10 and Z L ′ = 1. Figures 10a-10d and Table 3 show the impact of a localized weak zone on the preseismic and postseismic displacement rate (V x ′) profiles for the block-type LWZ model, with η w = 0.01 (i.e., f w = 100) for different thicknesses (Θ B ′ = 0.7 and 0.3) and widths (Ω B ′ = 0.2, 0.4, 0.6, 1.0, and 2.0) of the weak zone. These calculations are in some sense intermediate between, and may be compared directly with, those of the UNV models shown in Figure 7a (η′ = 0.01) and 7c (η′ = 1).
If the weak zone extends right through the crustal layer (Θ B ′ = 0.7, Figures 10a and 10c) the widest weak zone that we tested (Ω B ′ = 2.0) produces a postseismic V x ′ profile which is similar to that of the UNV model shown in Figure 7a, but the V x ′ profile clearly is not fully relaxed by the end of the earthquake cycle ( Figure 10c) because of the higher viscosities outside the weak zone. As we decrease the width of the weak zone Ω B ′ to 1.0, 0.6, 0.4, and 0.2, the peak postseismic V x ′ decreases by a factor of about 5 (Figure 10a), while the width of the high-gradient zone in the preseismic V x ′ profile is correspondingly decreased (Figure 10c). A change in strain rate (velocity gradient) is clearly evident above the edge of the weak zone in the postseismic V x ′ profile only for Ω B ′ = 2.0 ( Figure 10a). The width of the high-gradient zone in the preseismic V x ′ profile is influenced by the width of the weak zone, but the change in gradient becomes small as Ω B ′ decreases (Figure 10c).
If the weak zone instead covers only the uppermost part of the viscoelastic layer (H′ < z′ < 0.6; Figures 10b  and 10d), the wider weak zone (Ω B ′ = 2.0) still provides a considerable amplification of the peak postseismic V x ′ (Figure 10b), though less than that of Figure 10a, and the relaxation is much less complete at the end of the earthquake cycle (compare Figures 10c and 10d), causing a significant localization of the strain rate. The edges of the weak zone are evident in the preseismic profiles only for the widest weak zone (Ω B ′ = 2.0); for smaller Ω B ′ the high-gradient zone in the preseismic V x ′ profiles is significantly wider than the weak zone and is effectively a constant shape for Ω B ′ ≤~1.0.
The behavior of the Gaussian-type LWZ model (Figures 10e-10h) is similar to that of the block-type LWZ model (Figures 10a-10d). However, we note the following differences: For comparable maximum weakening factors (with Ω G ′ = Ω B ′ and for Θ G ′ = Θ B ′), the maximum postseismic V x ′ is significantly smaller (Figures 10e and  10f) for the Gaussian-type LWZ because viscosities gradually increase toward the background viscosity. There is a greater dependence of the preseismic shear zone width on the width (Ω G ′) of the LWZ (compare Figures 10c and 10g or 10d and 10h). Significant stress relaxation can occur over a wider region for experiments in which the Gaussian profile climbs gradually to the background viscosity, compared to those in which the viscosity shows a step-like increase (for Ω G ′ = Ω B ′; Figure 5).

10.1002/2013JB010347
We also tested both types of LWZ model using values of η w ′ in the range 0.1-0.001. For given geometrical parameters (Ω′ and Θ′), the peak postseismic V x ′ is approximately inversely proportional to η w ′ but is nearly proportional to Ω′ if Ω′ is sufficiently small (Ω′ < 1) as shown in Table 3. We also see that, for fixed η w ′, the distance to peak postseismic V x ′ increases systematically as Ω′ decreases, approaching values of 1 at Ω′ = 0.2, the smallest Ω′ for which we could obtain reliably resolved solutions.  (c, d, g, h) immediately before an earthquake for block-type LWZ model (Figures 10a-10d) and Gaussian-type LWZ model (Figures 10e-10h). In each case the maximum weakening factor is f w = 100 (η w ′ = 0.01), and the uniform background viscosity is η′ = 1.0. We consider two different weak zone thicknesses (Θ B ′ or Θ G ′) of 0.7 (Figures 10a, 10c, 10e, and 10g) or 0.3 (Figures 10b, 10d, 10f, and 10h), and in each frame we vary the width parameters (Ω B ′ or Ω G ′) as labeled. The grey line is the prediction for a uniform elastic layer subject to continuous loading.

Journal of Geophysical Research: Solid Earth
10.1002/2013JB010347  Figure 5), With Z L ′ = 1 and Y L ′ = 10 a Postseismic measures are obtained from cycle-invariant solutions immediately following an earthquake and preseismic measures immediately before an earthquake. Model measurements may be compared with dimensionless estimates from the NAFZ (Figure 2) of postseismic peak V x ′ in the range 5 to 8, distance to peak in the range 0.25 to 0.5 (measurements from the 6 months period following the İzmit earthquake), preseismic peak gradient of~1.30, and preseismic far-field gradient in the range 0.34 to 0.45 (measurements from the decade preceding the major earthquakes of 1999). Dimensional values are made dimensionless here by distance scale 40 km and velocity scale 11.25 mm/yr.

Application to the North Anatolian Fault Zone
In the context of the UNV model the Maxwell relaxation time τ is central to the response of the viscoelastic system. The high ratio (5-8) of postseismic displacement rates to preseismic displacement rates in the western North Anatolian Fault Zone (Figure 2) clearly indicates that the response is determined by material whose Maxwell relaxation time is small relative to the earthquake cycle period Δt. In the context of a UNV model, such a ratio would imply η′ between 0.01 and 0.1 (Figures 7a and 7b). However, such UNV models are incompatible with the well-defined gradient observed in the preseismic data within about 25 km of the fault, because the relaxation of stress differences occurs too quickly. The localized preseismic strain is produced by material for which η′ > 1.0.
The DDV model, in which viscosity decreases exponentially with depth below the elastic lid, comes closer to explaining both preseismic and postseismic displacement rate distributions (Figure 9). We can predict the high postseismic displacement rates for η b ′ = 0.001 and 0.01, and we can now see the relaxation of stress delayed significantly in the preseismic period for the larger viscosity contrasts (η h ′/η b ′ >~1000), due to the greater Maxwell relaxation time of the upper part of the medium. However, all of the DDV model calculations predict that the highest postseismic velocities occur significantly further from the fault and decay more slowly with increasing distances (e.g., Figure 9a) than is observed (Figure 2a).
The concept of a localized weak zone, whether block-or Gaussian-type, clearly provides the flexibility to reproduce the general features of the geodetic displacement profiles before and after the earthquakes. We use the geodetic profiles to try and constrain the dimensions of a possible weak zone and thereby provide some guidance for the possible spatial viscosity variations which may be needed to explain the data of Figure 2 using a physically grounded rheological model. We compare selected model velocity profiles, dimensionalized for a layer thickness of 40 km and a boundarydisplacement rate of 22.5 ± 7.5 mm/yr with the geodetic observations from the NAFZ for block-type (Figures 11a-11d) and Gaussian-type (Figures 11e and 11f ) LWZ models. We take Δt = 200 years as representative of the earthquake cycle period for the North Anatolian Fault Zone (estimated at~150-280 years by, e.g., Stein et al. [1997], Klinger et al. [2003], and Rockwell et al. [2009]). The viscosity scale then is η 0 = μΔt =~2 × 10 20 Pa s for shear modulus μ = 3 × 10 10 Pa, but the viscosity scale is uncertain by about a factor of 2 due to possible variation of μ and Δt.
The block-type LWZ model with η′ = 1, η w ′ = 0.01, Θ B ′ = 0.3, and Ω B ′ =1.0 (corresponding to Θ B = 12 km and Ω B = 40 km for L 0 = 40 km) almost explains the data, both postseismic ( Figure 11a) and preseismic ( Figure 11b). The most significant remaining misfit is due to the velocity gradient of the high-strain zone of the preseismic data ( Figure 11b) being underestimated. The numerical experiments shown in Figures 10c  and 10d suggest, however, that this type of model does not allow a steeper high-gradient zone in the preseismic data. We tested the idea that reducing the thickness of the upper elastic layer may allow a steeper gradient to develop in the central part of the preseismic profile but found that a decrease from H′ = 0.3 to H′ = 0.2 had negligible effect on this gradient.
Increasing the width (Ω B ′) of the weak zone does not improve the fit to the postseismic rates and significantly degrades the fit to the preseismic velocity gradient across the fault; increasing the thickness (Θ B ′) of the weak zone also degrades the fit to both preseismic and postseismic rates.
On the other hand, a comparably good fit to the data may be obtained for η w ′ < 0.01 by reducing the width (Ω B ′) of the low-viscosity zone, as shown for the block-type LWZ model with η′ = 1, η w ′ = 0.005, Θ B ′ = 0.3, and Ω B ′ = 0.6 in Figures 11c and 11d. This apparent trade-off between viscosity (η w ′) and width (Ω B ′) of the weak zone suggests that an even narrower weaker LWZ may explain the data. Kenner and Segall [2003] also described this trade-off and showed that for shear zone widths less than the thickness of the elastic upper layer, the differences in model surface deformation are minor. We have not tested models with Ω B ′ < 0.2 here, due to restrictions on numerical resolution, but we draw attention to the problem that the distance to peak postseismic V x ′ systematically increases as Ω B ′ decreases. The observation that the peak postseismic V x ′ is found at a distance 15 ± 5 km from the fault (y′ = 0.25-0.5) appears to be at odds with a shear zone narrower than Ω B ′~0.2; our calculations show that distance to peak V x ′ approaches 1.0 for shear zone width~0.2 and is increasing as Ω B ′ decreases further.

10.1002/2013JB010347
Our tests also show that the LWZ model predictions do not change significantly for viscosity values of the background nonweakened domain with η′ ≥ 1.0; we require only that relaxation of the deviatoric stress in that domain is small whether the background viscosity is uniform or depth dependent.
The Gaussian-type LWZ model with similar geometrical parameters (Θ G ′ = Θ B ′ = 0.3, Ω G ′ = Ω B ′ = 1, and η w ′ = 0.01) fits the data less well (Figures 11e and 11f). The postseismic displacement rates are underestimated ( Figure 11e); similarly, the velocity gradient (strain rate) in the preseismic phase is underestimated (Figure 11f). For a given Θ or Ω significant stress relaxation occurs in a broader region for the Gaussian-type LWZ than for the block-type LWZ. Obtaining a better fit to data for this type of model does not seem possible: decreasing the width (Ω) or the thickness (Θ) of the weak zone improves the fit to the preseismic rate but degrades the fit to the postseismic rate.
We also considered the impact of domain width and domain depth on our preferred best fit solution. Doubling the cross-strike width and domain depth increases and decreases the peak postseismic velocity bỹ 10% and~20%, respectively, and changes the near-fault strain rates by a similar factor. The data that we Figure 11. The predictions of the preferred localized weak zone (LWZ) model, compared with the observed GPS velocities (a, c, e) after and (b, d, f) before the 1999 İzmit and Düzce earthquakes; the postseismic velocity profiles are obtained at 6 months after the earthquake. Distances are scaled using L 0 = 40 km, for which the elastic layer thickness (H) is 12 km; the depth extent of the fault (D) is 16 km. Displacement rates are scaled using a far-field relative displacement rate of~22.5 ± 7.5 mm/yr. The preferred model has a weak zone width of Ω B = 40 or 24 km and thickness of Θ B = 12 km for block-type LWZ (Figures 11a-11d) and Ω G = 40 km and thickness of Θ G = 12 km for Gaussian-type (Figures 11e and 11f). Viscosity of weak zone scales to~2.0 × 10 18 Pa s (Figures 11a, 11b, 11e, and 11f) or to~10 18 Pa s (Figures 11c and 11d) and viscosity of background scales to~2.0 × 10 20 Pa s (Figures 11a-11f), based on the estimated earthquake cycle period of~200 years and shear modulus of 3 × 10 10 Pa.
Journal of Geophysical Research: Solid Earth 10.1002/2013JB010347 consider here ( Figure 2) do not sufficiently constrain the preseismic velocity gradient in the far-field zone to enable us to estimate an apparent distance to boundary. We note, however, that Pollitz [2001] finds a better fit to deformation on the San Andreas system using velocity boundary conditions on a block of across-strike width 600 km (we use 400 km) than is possible when the boundaries are at infinite distance.
Finally, we ask what is the impact of assuming an infinite fault length when we know that the 1999 İzmit and Düzce ruptures occurred on fault segments of length~100 km. A limited test of the effect of our infinite faultrupture assumption is obtained by taking our cycle-invariant model, and in place of the next earthquake in the cycle, we substitute an earthquake that ruptures only 100 km of the fault (in the middle of our block). The resulting postseismic velocity profile across the middle of the rupture is predictably decreased in amplitude (by a factor of~1.3) relative to the case of the infinite rupture length. That decrease can be compensated by a relatively minor adjustment of the preferred viscosity of the weak zone in the crust beneath the fault. Figure 12 shows the temporal evolution of (U x ′ À t′) for our preferred LWZ model (η′ = 1, η w ′ = 0.01, Θ B ′ = 0.3, and Ω B ′ =1.0; as Figures 11a and 11b) at several different surface points along x′ = z′ = 0. At each point, the displacement is shown relative to that of the boundary at y′ = Y L ′/2, and we show several earthquake cycles after the initial transient has decayed. In the near field (at scaled distances from the fault of 12, 20, and 40 km), the behavior is similar to that predicted for a UNV model with low viscosity (Figures 6a and 6b), showing again the important role of the weak zone. On the other hand, in the far field (at scaled distances of 80, 120, and 160 km), the displacement rate clearly approaches the background rate of the far-field boundary, reflecting here the importance of the high background viscosity. As was evident in Figure 11, the surface displacement rates still differ significantly from the long-term background rate throughout the earthquake cycle (by~45% at 40 km, by~30% at 80 km, and by~20% at 120 km); the increasing displacement difference is corrected whenever an earthquake occurs. We note that the simplest way to improve the fit to the preseismic data (Figures 11b and 11d) would be to increase the scaling factor for the long-term driving rate to about 30 mm/yr.
Significant asymmetry of the preseismic velocity profiles relative to the North Anatolian Fault [e.g., Le Pichon et al., 2005] is not represented in any of our model calculations. Asymmetry in the data misfit ( Figure 11) has been ascribed to an asymmetric variation of rheological properties (including viscosity [e.g., Malservisi et al., 2001], elasticity [e.g., Lisowski et al., 1991;Le Pichon et al., 2005;Fialko, 2006;Schmalzle et al., 2006;Jolivet et al., 2008;Fulton et al., 2010], and effective elastic thickness [e.g., Schmalzle et al., 2006;Chéry, 2008]) and/or fault geometry [e.g., Fialko, 2006;Jolivet et al., 2008]. Interpreted resistivity sections from this part of the western NAFZ clearly show anomalously low resistivities on the south side of the main active fault strand that passes through Lake Sapanca [Tank et al., 2005]. However, we note that Vaghri and Hearn [2012] have recently concluded that it is difficult to explain such asymmetry by lateral contrasts in viscosity and effective elastic thickness, particularly in cases where strain is localized around the fault late in the interseismic interval. Another source of geometric asymmetry that affects this segment of the NAFZ is the smooth change in relative motion of Anatolia with respect to Eurasia [Reilinger et al., 2006;Floyd et al., 2010] from predominantly strike slip (to the east) to include a larger component of southward displacement (to the west) causing extensional strain that increases westward into the Aegean [Taymaz et al., 1991;Özeren and Holt, 2010]. The impact of these geometrical constraints on the strain field near the NAFZ remains a matter for further quantitative investigations.

Discussion
The compelling feature of the North Anatolian geodetic data set is that the postseismic surface velocities are considerably greater than those observed late in the cycle. This observation requires a model in which a substantial part of the viscoelastic layer has a Maxwell time constant, τ, less than~10% of the interevent time, Δt. On the other hand, the significant localization of preseismic strain rate (velocity gradient) observed near the fault requires that another part of the viscoelastic layer has a time constant τ ≥~Δt. Clearly, a constantviscosity layer cannot satisfy both of these constraints.
Our experiments support the argument that spatial variations in viscosity are required to match the observations. Models with viscosity decreasing exponentially with depth or with a localized weak zone beneath the fault can provide both long-and short-term components in the relaxation of stress during the earthquake cycle. We find that, compared to models with depth-dependent viscosity, the model with a localized weak zone provides a superior fit to the spatial distribution of postseismic velocities. Although we have not explored models which have both of these elements, we anticipate that these may also provide a satisfactory fit to the observations if the stress relaxation in the background nonweakened domain is small.
Our results are consistent with other studies of the impact of a low-viscosity zone beneath the fault by Cohen and Kramer [1984] and Traoré et al. [2014], who also showed that near-field strain rates later in the cycle are relatively high because the stress in the surrounding elastic or high background viscosity is not relaxed during the cycle.
Several physical processes might explain why a low-viscosity zone is present beneath the NAFZ, including non-Newtonian viscosity [e.g., Billen and Houseman, 2004], thermal dissipation [e.g., Yuen et al., 1978;Fleitout and Froidevaux, 1980;Thatcher and England, 1998;Leloup et al., 1999;Savage and Lachenbruch, 2003;Fialko, 2012, 2013], grain size reduction [e.g., Kirby, 1985;Karato and Wu, 1993], water content [e.g., Karato, 1986;Karato et al., 1986;Hirth and Kohlstedt, 1996;Korenaga and Karato, 2008], polymorphic transformation [e.g., Kirby, 1985;Newman et al., 1999], and pore fluid partial pressures [e.g., Raleigh and Paterson, 1965;Brace and Kohlstedt, 1980], or simply structural discontinuity of disparate crustal blocks juxtaposed by accumulated fault displacement. Our models suggest that the viscosity of the weak zone is at least~100 times lower than the background strong material. Anomalous temperatures of up to~300°K are expected in the crust beneath a major strike-slip fault [Yuen et al., 1978;Thatcher and England, 1998;Takeuchi and Fialko, 2012]. Such temperature increases could explain a reduction in effective viscosity by 2 orders of magnitude, depending on the activation energy [e.g., Ranalli, 1995], but the temperature anomaly would decay gradually with distance from the fault [e.g., Thatcher and England, 1998], qualitatively similar to our Gaussian-type localized weak zone model. Our preference for a weak zone bounded by a relatively abrupt change in properties (based on comparison of Figures 11a and 11b, 11c and 11d, and 11e and 11f) perhaps implies that the weakening is more likely to be caused by structural discontinuity or abrupt changes in lithology, grain size, and/or water content.
Several lines of independent evidence point to the existence of shear zones with distinct physical properties beneath the North Anatolian Fault and other major fault zones. A low-resistivity anomaly has been identified beneath the western NAFZ using magnetotelluric data [e.g., Tank et al., 2005;Kaya et al., 2013]. Their measurements show that low resistivities extend to depths of 50 km between and beneath the northern and southern fault branches of the western NAFZ. The hypocenters of the main shock and the aftershocks of the 1999 İzmit earthquake are close to the boundaries of this zone. Wannamaker et al. [2009] found a similar low-resistivity zone beneath the Malborough strike-slip fault zone in New Zealand, and Nakajima et al.
[2010] reported a seismic low-velocity anomaly beneath the Atotsugawa fault in Japan, with a lateral extent of 20-30 km. Bürgmann and Dresen [2008] reviewed conflicting evidence from microstructural and macrostructural studies of exhumed fault zones and suggest that shear zones are fairly ubiquitous, but that their widths may vary from 1 to 2 km, as exposed along the surface trace of the Alpine Fault [Little et al., 2002], to > 100 km based on exhumed mantle xenoliths from the San Andreas Fault system [Titus et al., 2007].

10.1002/2013JB010347
Geodetic observations of postseismic deformation alone have often been used to infer the rheology of the crust and mantle [e.g., Bürgmann and Dresen, 2008;Wright et al., 2013]. Our results suggest that viscosities inferred from postseismic data only should be treated with caution -the magnitude of the postseismic response in our weak zone models is primarily dependent on the viscosity of the weak zone and is relatively independent of the background viscosity of the region. In most cases, earthquake cycle models based on the viscosities inferred from postseismic observations would not reproduce the focused interseismic strain observed at most faults [Wright et al., 2013] later in the deformation cycle.
The screw dislocation model popularized by Savage and Burford [1973] is commonly used to model interseismic deformation [Wright et al., 2013]. This model contains no time dependence-it consists of an elastic half space with uniform slip beneath a locked elastic lid-but the slip rates inferred are often considered to be the long-term steady state slip rates for the observed faults. In our preferred weak zone model for the North Anatolian Fault, the postseismic displacement rates decay to near-steady state values after 10-15% of the earthquake cycle (20-30 years for the NAFZ; Figure 13). Once this postseismic period has passed, the screw dislocation model provides a reasonable fit to the velocity versus distance profiles ( Figure 13), with best fit relative displacement rates within 10% of the imposed background displacement rate. Geodetic estimates of slip rate determined using the steady state assumption of Savage and Burford [1973] therefore likely reflect actual long-term relative displacement rates, as shown by the general agreement between geologic and geodetic estimates [Thatcher, 2009;Meade et al., 2012].

Conclusions
As has previously been argued, producing rapid postseismic transients early in the earthquake cycle and focused interseismic strain late in the cycle is not possible in a material with uniform viscosity. In the case of the North Anatolian Fault, matching the observed velocities in the period following the 1999 earthquakes requires relaxation times in the viscoelastic region that are <~10% of the interevent time, whereas focusing strain near the fault late in the cycle requires relaxation times that are equal to or longer than the interevent time. We have shown here that this apparent paradox can be reconciled if the viscosities of the substrate are spatially variable. We argue that a weak shear zone in the midcrust beneath the fault is the simplest and most likely the way to explain the observations. For the North Anatolian Fault, our best fit to the observed postseismic transients is obtained for a Maxwell relaxation time in the weak zone~100 times less than the interevent time, implying viscosity~2 × 10 18 Pa s. The inferred weakening factor trades off against the width of the shear zone to some extent as narrower shear zones require weaker material to match the observations. To obtain focused strain late in the cycle, the relaxation time of the upper crust away from the fault zone must be equal to or longer than the interevent time, implying viscosities ≥ 2 × 10 20 Pa s. Shear heating can provide some of the weakening required beneath the fault, but additional mechanisms appear also necessary-these are likely to include lithological discontinuity, grain size reduction, elevated water content, and/or an apparent weakening due to a power law rheology.
The high deformation rates in the postseismic period in our models depend primarily on the relatively low viscosities of the weak zone. Conclusions about the viscosity of the crust should therefore not be drawn from models that fit postseismic deformation alone. Counterintuitively, focused strain in the interseismic period does not result primarily from a weak zone beneath the fault but instead arises because the background viscosity is high. Faster viscous relaxation in the weak shear zone reduces the preseismic strain rate in the near-fault zone. Although a significant weak zone is needed to explain the postseismic strain rates, it must be restricted in both lateral and depth extent. High viscosities below and adjacent to the LWZ are needed to explain the localization of strain rate in the interseismic period.
Our best fit model, in which the background material is strong, relaxes to a near-steady state velocity profile after 10-15% of the earthquake cycle, after which the localization of strain near the fault is broadly consistent with the classic screw dislocation solution frequently used to estimate strain rates on strike-slip fault systems. This near-fault localization of strain implies a Maxwell time constant of the crust (away from the inferred weak zone) that is generally large relative to the length of the earthquake cycle.