Modeling the dynamic rupture propagation on heterogeneous faults with rate-and state-dependent friction

We investigate the effects of non-uniform distribution of constitutive parameters on the dynamic propagation of an earthquake rupture. We use a 2D finite difference numerical method and we assume that the dynamic rupture propagation is governed by a rateand state-dependent constitutive law. We first discuss the results of several numerical experiments performed with different values of the constitutive parameters a (to account for the direct effect of friction), b (controlling the friction evolution) and L (the characteristic length-scale parameter) to simulate the dynamic rupture propagation on homogeneous faults. Spontaneous dynamic ruptures can be simulated on velocity weakening (a < b) fault patches: our results point out the dependence of the traction and slip velocity evolution on the adopted constitutive parameters. We therefore model the dynamic rupture propagation on heterogeneous faults. We use in this study the characterization of different frictional regimes proposed by Boatwright and Cocco (1996) based on different values of the constitutive parameters a, b and L. Our numerical simulations show that the heterogeneities of the L parameter affect the dynamic rupture propagation, control the peak slip velocity and weakly modify the dynamic stress drop and the rupture velocity. Moreover, a barrier can be simulated through a large contrast of L parameter. The heterogeneity of a and b parameters affects the dynamic rupture propagation in a more complex way. A velocity strengthening area (a > b) can arrest a dynamic rupture, but can be driven to an instability if suddenly loaded by the dynamic rupture front. Our simulations provide a picture of the complex interactions between fault patches having different frictional properties and illustrate how the traction and slip velocity evolutions are modified during the propagation on heterogeneous faults. These results involve interesting implications for slip duration and fracture energy. Mailing address: Dr. Elisa Tinti, Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143 Roma, Italy; e-mail: tinti@ingv.it


Introduction
Modeling the dynamic propagation of earthquake ruptures requires the adoption of a constitutive relation, which governs the traction evolution within the cohesive zone (Ida, 1972;(Tse and Rice, 1986;Rice, 1993).Furthermore, they have been applied to model preseismic and postseismic processes, such as earthquake nucleation (Dieterich, 1992;Lapusta and Rice, 2003) and afterslip (Marone et al., 1991).
Another constitutive relation, the Slip-Weakening law (SW), is also currently used in dynamic modeling of earthquake ruptures; according to this constitutive law the traction depends only on the slip (Ida, 1972;Andrews, 1976a,b;Ohnaka and Yamashita, 1989).Although these two constitutive formulations are considered alternative (see Ohnaka, 2003), they both provide a physically correct description of the dynamic propagation of an earthquake rupture (see Bizzarri et al., 2001 and references therein).Recent studies have proposed analytical relations to associate the constitutive parameters of these two formulations (Cocco and Bizzarri, 2002;Bizzarri and Cocco, 2003).The main difference between these two constitutive laws is that SW prescribes the traction evolution with slip, whereas RS do not and traction spontaneously evolves with time and/or slip driven by the state variable evolution (Bizzarri and Cocco, 2003;Cocco et al., 2004).It has to be pointed out that RS laws have been derived from laboratory experiments at low slip velocity (<1 cm/s) and are usually assumed (as in the present study) to be valid also at high slip velocities (∼ 1 m/s), such as those observed during large earthquakes.However, other mechanisms might affect fault friction at high slip velocity such as frictional heating (see Fialko, 2004 and references therein), thermal pressurization (see Andrews, 2002 and references therein) or mechanical lubrication (Brodsky and Kanamori, 2001).
Numerous recent investigations have shown that the rupture history imaged on the fault plane during moderate-to-large magnitude earthquakes is quite complex: slip is non-uniformly distributed on the fault and rupture velocity and slip duration can change dramatically during the dynamic rupture propagation.The heterogeneity of slip and the variations of rupture speed are certainly associated with the complexity of fault geometry, the non-uniform distribution of prestress (Day, 1982;Peyrat et al., 2001) and the variations of frictional properties on the fault plane (Boatwright and Cocco, 1996).The short slip duration and the healing of slip can also be associated with stress heterogeneity (Beroza and Mikumo, 1996;Day et al., 1998).Conversely, in a homogeneous configuration the healing of slip has been modeled by appropriately modifying the constitutive relation (Perrin et al., 1995;Zheng and Rice, 1998;Cocco et al., 2004); this allows a sudden traction re-strengthening.However, it is important to point out that, although the constitutive laws which include self-healing of slip are physically reasonable, they have not been corroborated with laboratory experiments.
The goal of this study is to model the dynamic propagation of an earthquake rupture on a heterogeneous fault using RS constitutive laws.This implies that we have to describe and represent the frictional heterogeneity in terms of non-uniform distributions of RS constitutive parameters along our 2D fault model.Several studies focused on the investigation of the effects of spatial hetero-Table I. Frictional behaviors proposed by Boatwright and Cocco (1996).A= aσn eff  geneities of the constitutive properties (see Tse and Rice, 1986;Rice, 1993;Boatwright and Cocco, 1996 among many others).In particular, Boatwright and Cocco (1996;BC96 in the following) discussed the frictional control of crustal faulting by using a RS law and a single degree of freedom spring-slider dynamic system.They proposed that the fault response to the dynamic stress perturbations can differ depending on the variability of the constitutive parameters.They proposed four frictional fields separating the Velocity Weakening regime (VW) into strong and weak fields and the Velocity Strengthening regime (VS) into compliant and viscous response (see table I).BC96 associated these four frictional regimes with different values of the constitutive parameters and with observed features of seismicity and strain release as depicted in table I. Several authors have suggested that fault frictional parameters can vary with depth (Blanpied et al., 1991(Blanpied et al., , 1995;;Rice, 1993;Lapusta et al., 2000 among many others).BC96 proposed that these frictional parameters can also change along the fault strike, suggesting that frictional properties may control crustal faulting.In this study we used our 2D finite difference algorithm to inves-tigate the rupture propagation on a fault modeled through the different regimes proposed by BC96.We study how the interaction between velocity weakening and velocity strengthening portions of a 2D fault can affect the dynamic rupture propagation.

Methodology
We solve the fundamental elastodynamic equation neglecting body forces for a 2D inplane shear crack for which the displacement and the shear traction depend on time and one spatial coordinate.We solve the fully dynamic, spontaneous problem by using a finite difference method where the fault is simulated as a plane of split nodes in the grid (the Traction-at-Split-Nodes -TSN, numerical technique).The approach is described in Andrews (1973) and in Andrews (1999) and the details of the numerical implementation are in Bizzarri et al. (2001).Our procedure allows the adoption of different constitutive laws.In this study, we use the RS law with slowness evolution equation derived by Dieterich (1986), which has been discussed  in detail by many authors (Gu et al., 1984;Tse and Rice, 1986;Roy and Marone, 1996) .
In (2.1) v is the slip velocity, Ψ is the state variable (which provides a memory of previous slip episodes), µ * and v * are reference values for the friction coefficient and for the slip velocity, respectively.a, b and L are the constitutive parameters: a represents an instantaneous rate sensitivity, that is the direct frictional response to a change in slip velocity, L is the characteristic length that together with b controls the evolution of state variable toward the steady state.σn eff is the effective normal stress; in the literature the parameters A = aσn eff and B = bσn eff are also often used.The length-scale parameter L is different from the characteristic slip-weakening distance Dc (see Cocco and Bizzarri, 2002).We will show the relation between L and Dc in the next section.
All the parameters used to perform the simulations of dynamic rupture propagation presented in this study are summarized in tables II and III: λ and µ are the Lamè constants; vP and vS are P-and S-wave velocities, respectively; vinit is the initial value of slip velocity; ∆ x is the spatial grid size.The adopted discretization and medium parameters guarantee sufficiently accurate numerical experiments: the simulations presented and discussed in this paper have negligible grid dispersion (i.e. the accuracy frequency in space is 50 KHz).Moreover, all the stability and convergence conditions (see Bizzarri and Cocco, 2003) are satisfied.

Dynamic rupture propagation on a homogeneous fault
We present in this section the simulation of the dynamic rupture propagation on a homogeneous fault model with a uniform distribution of constitutive parameters.This fault model can be considered as a reference model for the discussion that will be presented in the following of this paper; this reference model was previously described by Bizzarri and Cocco (2003).The parameters adopted for this simu-lation are listed in table II.The model represents a velocity weakening fault for which the value of b -a is 0.004 (i.e.B -A = 0.4 MPa).The results of this simulation are shown in fig.1a-d and they represent a good example to discuss the typical behavior of a spontaneous dynamic rupture governed by RS.The nucleation strategy is described in details by Bizzarri et al. (2001); the nucleation patch is quite small and the nucleation stage is relatively short.
Figure 1a shows the spatio-temporal evolution of slip velocity, while the bottom panels display the traction evolution for a target point P1 as a function of state variable (fig.1b), slip velocity (fig.1c) and slip (fig.1d).We summarize in the following the major features emerging from this simulation: i) The rupture speed increases during the dynamic propagation and remains for this configuration sub-shear (vcrack ∼ 2200 m/s), although in general it might accelerate to a supershear value, asymptotically approaching Pwave speed (Bizzarri et al., 2001).
ii) Peak slip velocity increases during the dynamic propagation with an increasing spatial distance from the nucleation; the largest increase occurs when the crack accelerates to higher rupture velocities (figs.1a and 2).
iii) Simulations for a homogeneous configuration show no healing of slip: slip velocity does not drop to zero during the rupture propagation (figs.1a and 2).All points have roughly the same final slip velocity value (v 2), which represents the velocity at the new steady state after the dynamic stress release (fig.2).
iv) The dynamic traction evolution shows a characteristic slip-weakening behavior (fig.1d for a target point P1; see also Cocco and Bizzarri, 2002;and Bizzarri and Cocco, 2003).The weakening phase is denoted in this figure by the labels B and D (the latter identifies the point in which dynamic traction reaches the kinetic stress value, τf eq , see fig.1d).The slip required for traction to drop is called «equivalent» slip weakening distance D c eq .v) The phase diagram shown in fig.1c points out that the shear stress reaches the peak value (the maximum stress value τu eq , i.e. the equivalent yield stress, indicated in panel d and characterizing the point B in fig.1b-d) earlier than slip velocity (label C), in agreement with Tinti et al. (2004).
vi) Figure 1b points out that during dynamic propagation the state variable evolves from the assumed initial steady-state value (Ψ init) to a new one (Ψ ss = L / v).This represents the well-known self-damping behavior of the state variable.Bizzarri and Cocco (2003) and Cocco et al. (2004) have shown that it is the state variable evolution (line BC in panel b) that drives the slip acceleration (same line in panel c) and most of the traction drop during the weakening phase (see panel d).
vii) The duration of the weakening phase (line BD in panels c and d) defines the breakdown zone duration, as described by Ohnaka and Yamashita (1989).
viii) Simulations performed with a homogeneous configuration of RS friction with a slowness evolution law (eq.(2.1)) yield a constant weakening rate (see fig. 1d), provided that the spatial and temporal discretization is accurately selected to resolve the state variable evolution and the fast slip acceleration.
The considerations presented above are characteristic of the rupture propagation on a velocity weakening homogeneous fault.It is important to point out that the peak slip velocity, the yield stress, the dynamic stress drop, the rupture velocity, the slip weakening distance are strongly affected by the values of a and b.Moreover, the values of τu eq and τf eq , whose analytical expressions as a function of RS constitutive parameters are computed by Bizzarri and Cocco (2003), do not show relevant variations during the rupture propagation: we observe that, as the crack grows, τu eq slightly increases while τf eq remains nearly constant (see fig. 2).

The cohesive zone
The concept of cohesive zone was originally introduced by Barenblatt (1959a,b) for a tensile crack and subsequently by Ida (1972) for shear cracks in order to remove the physically unrealistic singularity of dynamic stress at the crack tip and to avoid an unbounded energy flux at the rupture front.The cohesive zone is the region of shear stress degradation near the crack tip; it is located just behind the rupture front and it is also named breakdown zone (Ohnaka and Yamashita, 1989).Different physical processes can be responsible of the shear stress degradation within the cohesive zone, and we refer to them as breakdown processes.The most important physical quantity characterizing the breakdown process is the characteristic slip-weakening D c.This parameter represents the amount of slip required to drop the dynamic traction from the upper yield stress (τu eq ) to the kinetic friction level (τ f eq ) and to absorb fracture energy (see fig. 1d).Dc is an input parameter imposed a priori in the well known SW constitutive law (Andrews, 1976a,b).The spatial extension of the cohesive zone scales with Dc and both these parameters are important to assess the required resolution of numerical experiments (see Bizzarri et al., 2001).
As briefly mentioned above, numerical simulations performed by adopting RS constitutive laws display a slip-weakening behavior of dynamic traction (see fig. 1d).Bizzarri and Cocco (2003) refer to this parameter as the equivalent slip-weakening distance (d 0 eq , here indicated as Dc eq ) and suggest that it is related to the charac-teristic length scale parameter (L) of RS formulation (eq.(2.1)) by the following scaling relation: where v 0 is the slip velocity value when the displacement reaches Dc eq (see fig. 1c) and vinit is the initial sliding velocity.This scaling relation was obtained under the assumption that the rupture starts from a steady state value of the state variable, but can be also generalized for an arbitrary initial condition of the state variable (see eq. ( 9) in Bizzarri and Cocco, 2003).
We show in fig. 3 the slip-weakening curves for a set of neighboring points resulting from the simulation illustrated in figs.1a-d and 2 (i.e. the reference model).This figure reveals that, in our simulations, D c eq does not depend on the distance from the nucleating patch.We infer a value of Dc eq larger than the adopted L parameter and the ratio Dc eq /L ranges between 15 and 20, in agreement with Cocco and Bizzarri (2002), Bizzarri and Cocco (2003) and Lapusta and Rice (2003).We observed the slip-weakening curves resulting from several simulations performed by adopting different values of the L parameter (ranging between 8 µm to 15µm): the inferred Dc eq value increases for increasing L in agreement with the scaling relation expressed in eq.(3.1).
All these simulations were performed with parameter values representative of the laboratory scale (L ∼ 10 µm and Dc eq ∼ 0.15 µm).The inferred Dc eq values for actual earthquakes range between several centimeters to meters (see Tinti et al., 2004 and references therein), suggesting that Dc may be a significant fraction of the maximum slip on the fault.However, the scaling of the constitutive parameters from laboratory to real fault dimension is still an open question and different opinions exist about the reliability of these extremely large values of D c for real earthquakes (see Guatteri and Spudich, 2000;Guatteri et al., 2003).

The velocity weakening frictional regime
Because we are interested in modeling the spontaneous dynamic propagation of an earthquake rupture, we first consider the propagation in a velocity weakening frictional regime.This frictional behavior is characterized by (B-A)>0 and relatively small values of the L parameters.We present in this section the results of several simulations performed on a homogeneous fault.The model parameters used for these simulations are listed in table II.The convergence and stability conditions of numerical experiments are discussed in detail in Bizzarri and Cocco (2003).We consider a set of parameters representative of the laboratory scale and we assume a constant normal stress σn eff (σn eff = σn -pfluid; where pfluid is the pore fluid pressure).
The Velocity Weakening (VW) behavior (A< B) represents the unstable regime that causes dynamic slip episodes on the fault plane, dynamic stress drop and the emission of seismic waves.For instance, we can associate VW patches with intermediate depth fault portions, like those located between 3 to 15 km along the San Andreas Fault.

Strong velocity weakening fault
We first model a VW regime characterized by a relatively large difference of constitutive parameters (B-A ≥ 0.6 MPa) and, according to BC96, we refer to it as a strong VW regime.The results of the simulations are shown in fig.4a,b and the input and constitutive parameters are listed in table II.The strength parameter was introduced by Das and Aki (1977a,b) to quantify how a fault area is unstable and ready to fail [S = (τ u -τo)/(τo-τf)].For a 2D fault governed by a linear SW law, Andrews (1976a,b) showed that S expresses also a limit to discriminate if a crack can (S<1.77)or cannot (S > 1.77) propagate with super-shear rupture velocity.For RS constitutive laws it is possible to define a S eq equivalent to S by using the values of parameters τu eq , τo and τ f eq (Bizzarri and Cocco, 2003).The simulation shown in fig.4a,b, performed using the RS constitutive laws for a strong seismic behavior (B -A = 0.75 MPa), is associated with a value of equivalent parameter S eq <1.77.
We can clearly see in fig.4a that the rupture front bifurcates; this was previously obtained by Okubo (1989) and Bizzarri et al. (2001).The spatio-temporal plot shown in this figure has been drawn with a resolution smaller than that adopted for numerical calculations to better illustrate the increase of peak slip velocity and the jump in rupture front speed.Figure 4b shows the temporal evolutions of slip velocity: black lines indicate slip velocity time histories at those fault positions located before the acceleration of the rupture front to super-shear speed; the colored lines show slip velocity at those points where the crack is bifurcated.The increase of peak slip velocity occurs within the cohesive zone and it is more evident for those points located before the rupture front bifurcation.After the bifurcation (x ≥ 4.5 m), peak slip velocity increases during the propagation of the external front (propagating at super-shear rupture speed), but it is nearly constant for the slow internal front.Moreover, these simulations illustrate that when the rupture jumps to a super-shear rupture velocity, the slip velocity drops to zero and thus slip heals before to re-accelerate in the internal rupture front.Thus, a strong VW regime allows the simulation of earthquake ruptures in agreement with theoretical experiments performed with a slip-weakening law (see Bizzarri et al., 2001 and references therein) and with laboratory experiments (Rosakis et al., 2000).
The strong seismic areas have a large breakdown stress drop (defined as the difference between the yield and the frictional stress ∆τb= =τu eq -τ f eq ) and very large coseismic slip: in a strong VW patch, slip and peak slip velocity are three times larger than those obtained for the reference model (fig.1a-d) at the same distances from the nucleation.The rupture velocity is close to 2.6 km/s immediately after the nucleation and the external front moves at about 4.5 km/s after the bifurcation.The separation between the two fronts increases going far from the nucleation, because of the different rupture velocities.

Weak velocity weakening fault
According to BC96 we define a weak fault as a fault characterized by a small difference between A and B (B -A < 0.05 MPa).The results extrapolated by comparing figs.1a-d and 4a,b are representative of simulations characterized by values of B and A parameters selected to have (B -A) > 0.2 MPa.We experienced a severe difficulty in achieving a spontaneous nucleation and then spontaneous dynamic rupture propagation with 0 MPa<B -A < 0.1 MPa performing many numerical experiments.We point out here that decreasing (B -A) we increase the dimension of the nucleation patch, which is defined as (Dieterich, 1992) where η is a geometric constant that depends on the crack type.Therefore, to avoid increasing the size of the nucleation patch we are confined to using a smaller L. We obtained a spontaneous nucleation using (B -A) = 0.1 MPa and L = 0.8 × 10 −6 m (see table II for the whole set of adopted parameters).As expected, under these conditions the slip velocities within the nucleation patch (where the instability is promoted) are much larger than those inferred during the dynamic propagation within the weak fault portion.
Our modeling results suggest that: i) peak slip velocity in weak regions can be even more than 10 times smaller than that simulated for the reference model shown in fig.1a-d; ii) the breakdown stress drop is decreased by 30%, although yield stress and kinetic friction level are both decreased; iii) the equivalent slip weakening distance decreases because both the breakdown stress drop and L decrease; iv) the inferred rupture velocity is smaller than that of the reference model and we never observed a crack bifurcation; v) despite the smaller values of the dynamic parameters we still retrieve an evident slip weakening behavior of the traction versus slip curves.We believe that these results and the inferred behavior of a weak fault portion are physically reasonable because, according to BC96, we expect that weak fault zones undergo dynamic instabilities only during micro-earth-quakes or when they are loaded by a dynamic rupture front propagating in an adjacent strong weakening area.We will model this latter behavior in the next section.

Numerical representation of frictional heterogeneities
The goal of this section is to model the dynamic rupture propagation on a 2D heterogeneous fault.Because in our approach the spontaneous dynamic propagation is governed by a RS friction law, we represent the source heterogeneities in terms of non-uniform distribution of constitutive parameters.As described above, we follow the findings of BC96, who proposed that both the velocity weakening and the velocity strengthening regimes (see Scholz, 1990Scholz, , 1998) ) are separated into two fields depending on the values of constitutive parameters and the response to external loading (see table I).As shown in the previous section an earthquake rupture can spontaneously nucleate only within a Strong Velocity Weakening area (S-VW) characterized by sufficiently large values of the parameters (B -A) and L. However, a dynamic rupture can also propagate in a Weak Velocity Weakening area (W-VW).Nucleation in these regions occurs only for very small earthquakes or when it is forced by a sudden external stress change.We recall here that in their original classification BC96 used the A and B parameters, which in the present study correspond to A = aσn eff and B = bσn eff .
The Velocity Strengthening (VS) behavior (defined by the condition A>B) models a stable sliding, i.e. an aseismic slip.By definition, it is impossible to obtain spontaneous rupture propagation for a homogeneous case.The VS behavior was proposed to describe several portions of the San Andreas Fault characterized by creep events (see Scholz, 1990).A VS area can be characterized by the thickness of unconsolidated sediments (like Southern California) or by unconsolidated gouge within the fault (like in Central California; Marone et al., 1991).According to BC96 the velocity strengthening areas can contribute to the arrest of an earthquake rupture.In particular, compliant areas are velocity strengthening fault regions that slip aseismically but can be driven to instability if they are sufficiently loaded by an abrupt stress increase due to the rupture propagation in an adjacent velocity weakening area (BC96).The VS behavior is characteristic of the RS laws and it cannot be simulated using the SW law.
The two «intermediate» fields proposed by BC96, weak and compliant, have frictional velocity dependencies that are close to velocity neutral: they can modulate both tectonic loading and the dynamic rupture process.Aftershocks can occur on compliant areas around a high slip patch, but most of the stress is diffused through aseismic slip.In the following sections we will present and discuss the results of different simulations of dynamic crack propagation on a heterogeneous fault representing frictional heterogeneity in terms of a non-uniform distribution of either the L or B and A parameters.

Heterogeneous distribution of the parameter L
In this section we will present two numerical experiments in which we consider different values of the characteristic length L along the fault; all models and constitutive parameters are listed in table III.We recall here that the L parameter controls the state variable evolution and affects the size of nucleation zone (eq.(3.2)).Therefore, with the increasing L value, the fault spends a longer time to reach an extension equal to the critical length ,c and to initiate the spontaneous dynamic rupture propagation.In fig.5a,b we show the results of a numerical experiment in which a fault patch with L = 9 µm is surrounded by a region with a greater L (L = 15 µm).In panel a the spatio-temporal evolution of slip velocity is illustrated.In panel b we present the time histories of slip velocity for different points located along the fault.Figure 5a,b points out that the rupture penetrates within the region having a greater L. As previously noted (fig.1ad), during the propagation the peak slip velocity increases as the rupture front moves far away from the nucleation patch and as soon as the crack tip encounters the region with a greater L, peak slip velocity is immediately reduced, but even in this region it starts to grow again.As in the reference model, v2 (i.e. the velocity at the new steady state) is almost the same for all points.We have plotted in fig.6 the slip-weakening curves calculated for different points located along the fault strike for the same simulation shown in fig.5a,b.Because the state variable evolution is different for different values of the L parameter, the dynamic traction also displays a quite diverse behavior.When the rupture enters the region characterized by a greater L, the value of Dc eq and consequently the fracture energy increase.Thus, fig.6

corroborates the linear relation existing between the equivalent slip weakening distance D c
eq and the L parameter (eq.(3.1)).The kinetic friction level is nearly the same in the two regions.We emphasize that the slip-weakening curves maintain the expected trend also during the propagation along a heterogeneous fault but the dynamic physical quantities (such as dynamic stress drop or fracture energy) depend on the constitutive parameters.
The simulation presented above is characterized by a small contrast of the L parameter in the two regions.We show in fig.7a,b the results of a numerical experiment in which the variation of L is more pronounced: we simulate the dynamic propagation along a fault where the inner region has L = 10 µm and the external region has a higher L value (L = 1 mm).In this case the rupture is unable to penetrate the external patch.This configuration was named «barrier model» by Bizzarri et al. (2001).In fact, we observe that during the propagation within the inner region the rupture behavior is identical to the reference configuration of fig.1a.When the rupture front approaches the high-L region (L is 1 mm, thus the contrast is 1000), a back-propagating healing front causes the crack arrest, as the rupture does not have enough energy to break the barrier and propagate inside the high-L area.This process is known as a «barrier-healing».Slip velocity time histories of different fault points are plotted in fig.7b: we can clearly observe the in-crease of the peak as rupture propagates away from the nucleation patch.The back propagating front causes a remarkable decrease of the slip velocity with a consequent healing.In this case our simulation yields a slip velocity time history with a finite duration.

Heterogeneous distribution of (B -A)
In this section we will present and discuss some numerical experiments in which a velocity weakening zone is surrounded by velocity strengthening region.In the first simulation we compute the dynamic propagation of an earthquake rupture along a fault on which a velocity weakening patch is adjacent to a velocity strengthening region.The initial and constitutive parameters are listed in table III. Figure 8a,b shows the results of these calculations: the dynamic rupture nucleates and propagates within the velocity weakening zone, thus it penetrates within the velocity strengthening region for a small distance before being arrested.The rupture arrest is gradual and it generates healing of slip only when rupture is stopped inside the VS region.In panel b we plot the slip velocity time histories, which show the expect-   ates.In this area, the rupture velocity is very low and the crack is almost arrested.Because of the small dimension of the VS region, the crack is able to propagate within the VS patch and it re-accelerates again when it reaches the external VW region.Figure 9b shows the time histories of the slip velocity.It emerges that slip velocity has a finite duration only during the rup-ture propagation within the internal VW region, while in the external VW region it does not return to zero (i.e.no healing of slip).
These simplistic simulations provide a picture of the complex interactions between fault patches having different frictional properties.Our 2D simulations illustrate how the traction and slip velocity evolutions are modified during the propagation on heterogeneous faults.These calculations propose stimulating implications for slip duration and fracture energy.We will discuss these issues in the following sections.

Implications for slip duration
The simulations presented in the previous sections confirm that the dynamic rupture propagation on a homogeneous fault governed by a RS friction law (the slowness formulation defined in eq.(2.1)) does not show the healing of slip (as clearly shown in figs. 2 and 4a,b).This is in agreement with the results of Perrin et al. (1995) and Bizzarri et al. (2001).Although different regularizations or modifications of the constitutive formulations have been proposed to obtain the healing of slip and/or self-healing pulses (see Cocco et al., 2004 and references therein), even source heterogeneities can produce a finite duration of slip velocity and short slip durations (see Beroza and Mikumo, 1996;Day et al., 1998).In the present study we tested two different heterogeneous configurations; the first is based on the variation of the L parameter and the other on the interaction between VW and VS patches represented by changing the values of A and B.
A strong contrast of the L parameter represents a barrier and produces a «barrier healing» (see Zheng and Rice, 1998;Bizzarri et al., 2001) similarly to the arrest on the crack-like rupture propagation.Figure 7b shows that the slip duration depends on the distance from the barrier and it is shorter for closer distances to the barrier.Actually, the slip velocity time histories of the points closest to the barrier have a shorter duration than those located closer to the nucleation region.The interaction between VW and VS still yields finite slip durations (see fig. 8a,b).It is important to emphasize that the time histories of slip velocity simulated for heterogeneous configurations do not resemble a selfhealing pulse, although they have a finite duration.The latter is usually defined as a pulse propagating along the fault with nearly constant slip duration; in this case healing is independent of the crack arrest phase.In other words, while we simulated finite and relatively short slip du-rations, we are unable to generate self-healing pulses.These results are consistent with the findings of Perrin et al. (1995) and corroborate the outcomes of Day et al. (1998), who suggested that source heterogeneity yields healing of slip.An important implication emerging from these results is that healing of slip does not require traction re-strengthening: total dynamic traction remains at the kinetic friction level also when slip velocity tends to zero.

Implications for fracture energy
We computed fracture energy for the different simulations performed in this study.The fracture energy is defined as the amount of energy that the crack spends to advance and increase its length.It has been defined as and is measured as the area below the traction versus slip curve shown in figs.3 and 6 and above the minimum traction τf.For a dynamic rupture propagation on a homogeneous fault governed by a SW law the fracture energy G is constant and known a priori over the whole fault surface [G=1/2(τu−τf)Dc].This is not the case when rupture propagation is governed by other constitutive laws or when the constitutive parameters are not uniform on the fault.
In order to quantify the fracture energy changes during the crack growth, we plot in fig.10a-d the computed G as a function of the position along the fault for four different configurations investigated in this study.Figure 10a shows the resulting values for the homogeneous configuration (the reference model shown in fig.1a-d), which reveals an increasing trend as the rupture propagates far away from the nucleation patch.This is in agreement with the results shown in fig.2, which illustrate that the maximum yield stress increases moving far away form the nucleation patch, while τf eq slightly decreases.Figure 10b-d displays the resulting values computed for three heterogeneous configurations.In particular, fig.10b corresponds to the non-uniform distribution of the parameter L: as the rupture penetrates in the region with higher L, the computed fracture energy increases jumping to a larger constant value; this is also expected looking at the results shown in fig.6, which displays the increase in Dc eq .Figure 10c refers to the rupture arrest within a VS area shown in fig.8a,b: after an initial increase of G with increasing distance from the nucleation patch, consistent with fig.10a, we notice a decrease of G within the VS area.This is due to the decrease of dynamic stress drop and peak slip velocity shown in fig.6.In the VS patch the total amount of work spent to allow the crack advance is lower than that in the VW patch.Finally, fig.10d shows the fracture energy estimated in different fault positions during the rupture propagation in the heterogeneous fault modeled in fig.9a,b: the propagation within the VS area produces a sudden drop in G, which is followed by an increase when it starts to propagate again within the external VW region.
The calculations summarized in fig.10a-d point out that frictional heterogeneity explains the variability of fracture energy on the fault plane, which is associated with both the variations of slip and breakdown stress drop.

Conclusive remarks
One of the main goals of this study is to extend the results of BC96 to investigate the dynamic rupture propagation on a 2D fault with a heterogeneous distribution of constitutive parameters.We used the rate-and state-dependent formulations to characterize fault heterogeneities following the findings of BC96, who proposed to split the velocity weakening and the velocity strengthening regimes into four distinct frictional fields.Our results corroborate the conclusions of BC96 demonstrating that a velocity strengthening area can arrest and can also be driven to a dynamic instability by an earthquake rupture propagating in the adjacent fault patch.We have represented numerically  Our simulations show that the interaction between the propagating dynamic rupture front and the heterogeneous fault patches depends on the values of the constitutive parameters.In particular, we have shown that a variation of the L parameter can modify the peak slip velocity or can arrest the rupture propagation depending on the value of the L contrast.The heterogeneity of the L parameter does not modify the breakdown stress drop nor does it contribute to the variations of velocity if the contrast is smooth.On the contrary, the heterogeneity of the distribution of the difference (B− A) affects the dynamic rupture propagation in a more complex way: dynamic stress drop and strength excess strongly depend on B and A parameters.Moreover, rupture can penetrate within a velocity strengthening area and the heterogeneous distributions of B and A yield complex time histories of slip velocity.
We propose that frictional heterogeneities can explain the observed complexity of slip distribution and the variability of rupture velocity during earthquakes.Our results have important implications for slip durations (i.e.local rise time) and fracture energy.In particular, our results are consistent with the findings of Perrin et al. (1995) and Bizzarri et al. (2001).A large contrast in the L parameter represents a barrier that produces a crack-like solution with variable but finite slip durations.Because the variations of constitutive parameters affect both the critical slip-weakening distance (see Bizzarri and Cocco, 2003 and references therein) and the breakdown stress drop, the inferred fracture energy varies along the fault.We have shown that the increase in the L parameter results in a fracture energy increase and that heterogeneous distribution of (B− A) yields evident variations of fracture energy along the fault.Because we model here a laboratory fault, our estimates (∼10 3 J/m 2 ) of fracture energy cannot be compared with those inferred for real earthquakes.
In this study we mainly focused on the interactions between the propagating dynamic rupture front and the heterogeneous fault patches.According to BC96, the response of fault patches having different frictional properties to a con-stant tectonic load also controls the pattern of seismicity and the behavior of crustal faults.Thus we conclude that fault frictional properties and their variations on the fault plane play an important role in characterizing crustal faulting and the mechanical properties of major fault zones.

Fig
Fig. 1a-d.a) The spatio-temporal evolution of slip velocity during the dynamic rupture propagation on a homogeneous fault.The initial and constitutive parameters used for these calculations are listed in table II.The bottom panels show the traction evolution as a function of state variable (b), slip velocity (c) and slip (d) for a target point P1.Ψinit in (b) represents the initial value of the state variable; v 0 in (c) is the slip velocity when the displacement reaches the critical slip weakening distance Dc eq (d).

Fig. 2 .
Fig. 2. Temporal evolution of dynamic traction and slip velocity for several selected points located at different distances from the nucleation.Dashed lines represent the envelopes of τu eq and τf eq values.The constitutive parameters are the same used for the previous figure and are listed in table II.

Fig. 3 .
Fig. 3. Slip-weakening curves calculated for the reference model shown in the previous figures and listed in table II.The different traction evolutions as a function of slip have been plotted for different points located at different distances from the nucleation patch.

Fig
Fig.4a,b.Dynamic rupture propagation on a strong velocity weakening fault.The model and constitutive parameters used for these calculations are listed in table II.Panel a) shows the spatio-temporal evolution of slip velocity.We have drawn this plot with a reduced resolution in order to emphasize the increase in peak slip velocity.Panel b) shows the time histories of slip velocity in different positions along the fault: black lines identify points located before the point where crack accelerates to a super-shear rupture velocity, while coloured curves identify those points located in the region where the crack has bifurcated.

Fig
Fig. 5a,b.Dynamic rupture propagation along a heterogeneous fault: the adopted constitutive parameters are listed in table III.The fault is represented as two patches having different values of the L parameter (L1 and L3 in table III): the external one has a larger value (1.6 times the inner one).Panel a) shows the spatio-temporal evolution of slip velocity.Panel b) displays the time histories of slip velocity in different positions along the fault: solid lines identify slip velocity computed for points located in the inner region (low L), while dashed lines identify slip velocity computed for those points located in the external region (high L).

Fig. 6 .
Fig. 6.Slip-weakening curves calculated for the simulation shown in the previous figure, whose model and constitutive parameters are listed in table III.The traction evolutions as a function of slip have been plotted for several points located on the two patches having different values of L parameter at different distances from the nucleation patch.

Fig
Fig. 7a,b.Dynamic rupture propagation along a heterogeneous fault: the adopted constitutive parameters are listed in table III.The fault is represented as two patches having different values of the L parameter (L1 and L3 in table III): the external one has a larger value (1000 times the inner one).Panel a) shows the spatio-temporal evolution of slip velocity.Panel b) displays the time histories of slip velocity in different positions along the fault.

Figure
Figure9a,b shows the results of a numerical experiment in which two VW regions are separated by a VS patch located between them.The initial and constitutive parameters are listed in table III.As expected, the rupture initially accelerates within the VW region and it partially penetrates within VS area.Slip velocity is progressively attenuated and the crack tip deceler-

Fig.
Fig. 8a,b.Dynamic rupture propagation along a heterogeneous fault: the adopted constitutive parameters are listed in table III.The fault is represented as two patches having different values of the a and b parameters (a1, b1 and a3, b3 in table III): a velocity strengthening area is adjacent to the velocity weakening patch where the rupture nucleates.Panel a) shows the spatio-temporal evolution of slip velocity.Panel b) displays the time histories of slip velocity in different positions along the fault.

Fig.
Fig. 9a,b.Dynamic rupture propagation along a heterogeneous fault: the adopted constitutive parameters are listed in table III.A narrow velocity strengthening patch separates two velocity weakening areas (identified by the values of the parameters (a1, b1), (a2, b2) and (a3, b3) in table III).The rupture nucleates within the inner velocity weakening patch.Panel a) shows the spatio-temporal evolution of slip velocity.Panel b) displays the time histories of slip velocity in different positions along the fault.

Fig
Fig. 10a-d.Fracture energy values in different positions along the fault for four configurations investigated in this study: a) refers to the reference model shown in fig.1a-d; b) refers to the heterogeneous L distribution shown in figs.5a,b and 6; c) and (d) refer to the heterogeneous distribution of the a and b parameters shown in figs.8a,b and 9a,b, respectively.
by assigning different values of L or (B− A) parameters along the fault line.

Table II .
Model and constitutive parameters used for the simulations on a homogeneous fault.This set of parameters refers to a laboratory-scale experiment.

Table III .
Model and constitutive parameters used for the simulations on a heterogeneous fault.