Rake rotation introduces ambiguity in the formulation of slip-dependent constitutive models : slip modulus or slip path ?

The linear slip–weakening (SW) law, predicting that the traction decreases for increasing fault slip, is one of the most widely adopted governing models to describe the traction evolution and the stress release processes occurring during coseismic slip failures. We will show that, contrary to other constitutive models, the SW law inherently poses the problem of considering the Euclidean norm of the slip vector or its cumulative value along its path. In other words, it has the intrinsic problem of its analytical formulation, which does not have a solution a priori. By considering a fully dynamic, spontaneous, 3–D rupture problem, with rake rotation allowed, in this paper we explore whether these two formulations can lead to different results. We prove that, for homogeneous configurations, the two formulations give the same results, with a normalized difference less than 1%, which is comparable to the numerical error due to grid dispersion. In particular, we show that the total slip, the resulting seismic moment, the fracture energy density, the slip–weakening curve and the energy flux at the rupture front are practically identical in the two formulations. These findings contribute to reconcile the results presented in previous papers, where the two formulations have been differently employed. However, this coincidence is not the rule. Indeed, by considering models with a highly heterogeneous initial shear stress distribution, where the rake variation is significant, we have also demonstrated that the overall rupture history is quite different by assuming the two formulations, as well as the fault striations, the traction evolution and the scalar seismic moment. In this case the choice of the analytical formulation of the governing law does really matter.


Introduction
It is well known that numerical experiments represent, in addition to theory and to laboratory experiments, a suitable way to explore natural phenomena.In seismology, due to the obvious impossibility to plan direct experiments at the real world scale and due to the intrinsic, technical limitations of the laboratory machines, the numerical computations are extremely powerful and useful tools to explore, reproduce and model, in realistic conditions, the energy-dissipating processes occurring during faulting [e.g., Bizzarri 2009].
In pioneering, self-similar, elliptical expanding crack models [e.g., Burridge and Willis 1969] the slip is assumed to be everywhere and always parallel to the direction of the initial shear stress (the pre-stress).On the other hand, Madariaga [1976] shows that for a finite circular crack the rupture process introduces a component perpendicular to the direction of the shear prestress, which in general is quite small, but non zero.
One of the outcomes of realistic numerical experiments is that the direction of the fault slip (namely, the azimuth of the fault slip vector) can vary during the dynamic propagation of an earthquake rupture, even in homogeneous conditions.This phenomenon, also named rake rotation, appears both in the time domain (in a given fault node the slip direction changes as long as the rupture dynamically propagates) and in the spatial domain (at a given time the slip direction is not the same over the slipping portions of the fault).This is a wellknown result after the 2-D mixed-mode models by Andrews [1994] and the 3-D fault models by Bizzarri and Cocco ([2005]; readers can refer to Section 2 of Bizzarri [2011] for a compendious summary of the naming conventions).As pointed out also by Guatteri and Spudich [1998], the amount of the rake rotation is inversely proportional to the overall shear stress acting on the fault surface, making it possible to discriminate between the so-called high-stress and low-stress faults.Indeed, Bizzarri and Cocco [2005] show that, even if the strength parameter S is the same (namely, if the ratio between the strength excess and the dynamic stress drop is kept constant), the rake rotation is smaller if the absolute values of the stress levels increase (see their Figure 13).
We emphasize here that the rake rotation is, by definition, explicitly neglected in fault models where the slip is not allowed in the direction perpendicular to the direction of the initial shear stress (far of being ex-haustive, see for instance Aochi et al. [2000aAochi et al. [ , 2000b]], Fukuyama and Madariaga [2000], Madariaga et al. [1998], Nielsen and Olsen [2000]).The same occurs in fault models where the governing law is written in a vectorial form (i.e., it is written independently for each component of the physical observables), but only one component is non null [e.g., Olsen et al. 1997, Fukuyama and Madariaga 1998, Fukuyama et al. 2003].
Rake rotation is often connected to the phenomenon of fault striations [e.g., Etchecopar et al. 1981, Cashman and Ellis 1994, Otsuki et al. 1997, Haeussler et al. 2004] and it has been documented from the analysis of geodetic and ground motions data [Ide et al. 1996, Wald 1996, Yoshida et al. 1996].Moreover, it can be regarded as a possible way to estimate the level of shear stress prior to a slip instability event [Spudich 1992].
When rake rotation occurs it is clear that there is a difference between the slip modulus and the slip path.At a given time, the first quantity simply represents the actual magnitude of the slip vector, while the second one represents the fault slip cumulated during its evo-lution up to the considered time.This situation is clearly described in Figure 1 (see also Figure 2 of Spudich et al. [1998]) and can be written in analytical form as it follows: (1) and (2) respectively (in Equation (2) v is the time derivative of u).Equations ( 1) and ( 2) are written for a 3-D problem, in which the solutions of the elastodynamic equation depend explicitly on two spatial coordinates (in the special case of a vertical, planar fault normal to the x 2 , as such considered here, the axis x 1 is along the strike and the x 3 axis is down the depth) and on the time t.
Among the large class of governing models proposed in the literature to describe the chemical and physical processes potentially taking place during faulting (see Bizzarri [2011] for a review), only the slip-de-, , , , In the example displayed, at time level 3, u (mod) is identified by the magnitude by the thin blue vector u (3) , while u (path) is identified by the length of the thick green line.In the sketch it is also emphasized that, in general, u and v are not collinear; this is apparent from the difference of the two azimuths { u and { v .

#
pendent constitutive models suffer from the ambiguity stated by Equations ( 1) and (2).In other words, the ambiguity exists because the magnitude of the shear stress, T, can be function of u (mod) or of u (path) .On the contrary, when T does not depend on the slip, but on the slip rate (the fault slip velocity, v) the ambiguity automatically disappears.
Within the class of the slip-dependent laws, the linear slip-weakening law (SW thereinafter; see Equation (25) in Bizzarri [2011]) is probably one of the most popular (and widely employed) constitutive equation applied in dynamic model of faults, starting from the pioneering works by Ida [1972] and Andrews [1976aAndrews [ , 1976b] ] where the concept of slip-weakening mechanism has been introduced and from the laboratory experiments by Ohnaka and co-workers (see Ohnaka [2003] for a comprehensive review).The wide use of the SW governing model is because of its inherent simplicity from a modeling point of view (all the stress levels are a priori defined) and from a computational point of view (its numerical implementation is straightforward).The body of the literature here is immense, but the most interesting point is that some authors adopt the slip modulus (and hence the cumulated fault slip u is given by Equation ( 1)), while others use the slip path (Equation ( 2)).Just for example, Olsen et al. [1997], Bizzarri and Cocco [2005, 2006a, 2006b], Bizzarri and Belardinelli [2008] and Bizzarri and Das [2012], fall in the first category, while Day [1982], Day et al. [2005], Dalguer and Day [2006], Bizzarri et al. [2010] and Zhang et al. [2014] fall in the second category.
The goal of the present study is to explore whether the two formulations (which are different from a mathematical, conceptual and geometrical point of view) lead to different results in terms of rupture history, stress release, final (total) slip, seismic moment, etc.To properly handle this question we will consider both homogeneous and heterogeneous configurations (see Sections 3 and 4, respectively).

Statement of the problem and methodology
In the present paper we consider planar, vertical, strike slip faults of finite width.The ruptures start from the imposed hypocenter ((x 1 H , x 3 H )) and then spread spontaneously and obey the linear SW law.In models with homogeneous initial shear stress during their early stages the ruptures are forced to expand at a constant speed in order to obtain the earthquake nucleation, as described in Bizzarri [2010].In the model discussed in Section 4.2 the initial shear stress distribution contains an asperity located in the imposed hypocenter which guarantees the rupture propagation, as the shear stress exceeds the maximum yield stress.All the parameters are tabulated in Table 1.The borders of the computational domain have been selected in order to not interfere with the supershear transition, which occurs in the homogeneous cases.We have deliberately chosen supershear ruptures (Models A and B) because in the supershear regime the rake rotation is greater than for subshear ruptures [Bizzarri and Cocco 2005; see also Bizzarri and Das 2012] and therefore we maximize the possible differences between the two formulations (1) and (2).
As stated above, we consider a 3-D problem, which is spontaneous after the nucleation and it is fully dynamic (the inertial term is always considered); the elastodynamic problem, neglecting body forces, is solved numerically, by adopting the second order-accurate, conventional grid, finite difference approach described in Bizzarri and Cocco [2005].

Reference case
In the present section we will consider the most idealized configuration, in which the whole fault surface is characterized by the parameters listed in Table 1, Model A. More elaborated models, including spatial heterogeneities of the frictional parameters and initial shear stress, will be discussed in Section 4. The slip at t = 4.2 s resulting from the spontaneous rupture propagation in the case of formulation ( 1) is reported in Fig- ure 2a.The result from the formulation (2) is visually identical; in Figure 2b we report the normalized difference (misfit) between the developed slip from a numerical experiments adopting Equation ( 1) with respect to that from Equation (2).We can conclude that the variations are negligible, except near the bottom of the fault.The misfit between the developed fault slip reported in Figure 2b is the same over the whole duration of the rupture and not peculiar to the selected time.Small differences (of the order of a few percent) only emerge near the bottom of the fault; in this region the end of the traction-at-slip-node area is located and therefore here the arrest mechanisms occur.The latter cause the abrupt cessation of fault slip (from a non null value down to zero) and this cause in turn rake rotation.Incidentally, we note that the numerical method does not influence these differences, since independingly on the adopted numerical scheme and type of fault implementation, at the end of the computational domain an unbreakable zone is usually imposed and therefore arrest process will take place.Indeed, the scalar seismic moment and the corresponding magnitude of the event are practically the same for the two cases; in the case of formulation (1) we have: M 0 = 5.014 × 10 19 Nm (or equivalently M w = 7.067) and in the case of formulation (2) we have: M 0 = 5.015 × 10 19 Nm (the corresponding moment-magnitude is the same as in the previous case).
In Figure 3 we report the fault striation (i.e., the time evolution of u 3 vs.u 1 ) in a target fault receiver.We select a fault node near the supershear rupture front, where we know that the rake rotation (and the possible differences between formulations (1) and ( 2)) is maximum [Bizzarri andCocco 2005, Bizzarri andDas 2012].Nevertheless, we can clearly see that the differences of u between the two formulations are negligible (see inset of Figure 3a), in that the interval from 0 to d 0 (recall that after d 0 the magnitude of the shear traction is exactly at the residual level x f which is independent on the fault slip and thus the distinction between formulations ( 1) and (2) naturally disappears) the variations of u 3 are so small that, at each time, Equations ( 1) and (2) give the same results.Note that u 3 remains one order of magnitude below u 1 ; see inset in Figure 3a.Correspondingly, the slip-weakening curves (Figure 3b) are practically indistinguishable, so that the so-called local fracture energy density E G is the same for the two numerical simulations (E G = 10.45MJ/m 2 ).
An useful indication of the energy balance of the SLIP MODULUS OR SLIP PATH? two models arises from the computation of the energy flux F at the rupture tip (see Equation (8) of Bizzarri [2013] and references cited therein).This quantity is able to concurrently capture the speed at which the rupture front advances and the stress release occurring in the process zone.Contrary to the E G which is a local estimate (because it is computed at a given fault node), F takes into account all the slipping fault nodes at a given time.The time evolution of the energy flux is reported in Figure 4a and the misfits is plotted in Figure 4b.It is apparent that as long as the rupture remain subshear (t < 1.57 s) the solutions are identical.After the supershear transition the differences only emerge in a late stage of the rupture (nearly from t = 2.1 s; see Figure 4b).The maximum normalized variation of F is less than 1%, which is comparable to the numerical error due to grid dispersion.

Low initial stress case
In the model A presented in Section 3.1, the initial shear stress is relatively high, compared to the dynamic stress drop (we have: x 0 /Dx b ≈ 4).We then consider another homogeneous model (Model B; see Table 1), in which the frictional parameters have been chosen in order to reduce that ratio; in this case we have: x 0 /Dx b = 2 (in model B x 0 = 20 MPa, instead of 73.8 MPa, as in Model A).As discussed above (Introduction), we know that a high initial shear stress can reduce the rake rotation during the propagation of the rupture (and thus can mask possible differences between formulations (1) and ( 2)).The results of this numerical experiment are reported in Figure 5. Also in this case the variations between the two models are small, as small is the variation of u 3 with respect to its initial value of 0. The fault striations are very similar in the two formulations (1) and (2) (see Figure 5a), as well as the traction evolution (see Figure 5b).The fracture energy density is slightly different in this Model B (E G = 4.208 MJ/m 2 in case of ( 1) and E G = 4.203 MJ/m 2 in case of ( 2)).This makes sense, in that in the case of slip path the traction reaches earlier the residual level x f , as also illustrated in the inset of Figure 5b, so that the area under the slip-weakening curve is smaller in this case.As in the previous Model A, also in this case the total seismic moment in the case of slip path is slightly larger than that pertaining to the slip modulus (M 0 = 1.859 × 10 19 Nm and M 0 = 1.858 × 10 19 Nm, respectively; this corresponds to a difference ~5%).

Results for a heterogeneous rheology
The selected configurations discussed in Sections 3.1 and 3.2 clearly indicate that in the canonical, homogeneous configurations, even under supershear conditions, even with a relatively large value of characteristic SW distance, the two formulations (1) and (2) produce practically the same results.This is basically due to the fact that the rake rotation is not so pronounced to induce variations in the computation of the frictional resistance specified by the linear SW law.
The obvious problem is then: what happens if the rake rotation is indeed significant?In order to answer to this problem in the present section we explore more interesting configurations, including heterogeneities in the characteristic SW distance (Section 4.1) and in the initial shear stress (Section 4.2).High frequency (namely, greater than 1 Hz) patterns emerging in the radiated waves can be caused by frictional inhomogeneities (asperities and barriers), geometrical irregularities (such as bending and branching) or both.The characterization of the heterogeneities in dynamic model has received much attention of many authors [e.g., Lavallée et al. 2006, Bizzarri et al. 2010, Andrews and Barall 2011, Song and Dalguer 2013].Here we will consider heterogeneous distributions of frictional parameters (Section 4.2) and of the initial spatial distribution of the shear stress (Section 4.2).

Heterogeneities of d 0
In this section we consider a heterogeneous distribution of the characteristic SW distance d 0 , in the same way as Tinti et al. [2005]; in particular, for | we set d 0 = 1 m.All the other parameters are the same as those of Model B, except for a region of 1 km below the free surface, having a high strength (unbreakable barrier); in such a way we prevent the local supershear acceleration, which has been observed in numerical models when the rupture front interacts with the free surface [Olsen et al. 1997, Chen and Zhang 2006, Zhang and Chen 2006, Day et al. 2008, Bizzarri 2010, Kaneko and Lapusta 2010, Bizzarri and Das 2012, Zhang et al. 2014].Once the rupture penetrates the region with a high d 0 it decelerates and the propagation of the rupture is disturbed (see Figure 6a), so that the slip distributions are different compared to the homogeneous configuration of Model B. The misfit between formulations ( 1) and ( 2) (reported in Figure 6b) indicates that in the present case the differences are a little bit more significant with respect to the homogeneous configuration, especially in the correspondence of the primary (subshear) and secondary (supershear) rupture fronts.By comparing Figure 2b (homogeneous case) against Figure 6b we have that the misfit of the total fault slip predicted by the two formulations is <0.5% in the homogeneous case and ~2% in the present heterogeneous one.The resulting seismic moment and magnitude are: M 0 = 1.629 × 10 19 Nm (or equivalently M w = 6.741) for Equation ( 1) and M 0 = 1.631 × 10 19 Nm (or equivalently M w = 6.742) for Equation ( 2).The time histories of the solutions in a selected fault node (Figure 7) basically confirm the previous findings obtained for the homogeneous configurations (Models A and B); although d 0 in the selected receiver is now larger, the difference between the two formulations ( 1) and ( 2) is only in the slip-weakening curve (see the inset in Figure 7b).As a net results the fracture energy density is slightly different: E G = 7.007 MJ/m 2 in Figure 6.The same as in Figure 2, but now for a model having d 0 = 1 m at distances from the hypocenter lager than 5 km.The other parameters are the same as in Model B, with the exception of an unbreakable barrier extending 1 km below the free surface.In panel (b) the scale is clipped at 2% just to emphasize the (positive) differences near the supershear front.
Figure 7.The same as in Figure 3, but now for the heterogeneous models of Figure 6.  1) and E G = 7.000 MJ/m 2 in the case of formulation (2).This difference is comparable with that of previous Model B.

Heterogeneous distribution of the initial shear stress
In all the numerical experiments discussed so far the pre-stress was uniform and aligned along the strike direction (x 1 ).We build now a model in which the initial shear stress is highly heterogeneous and both its two components are non null; the spatial distributions are reported in Figures 8a and 8b.In this case the fault is larger than in the previous case and the rupture propagation is no longer symmetric with respect to the hypocenter.
The resulting distribution of the fault slip at t = 3.8 s is reported in Figures 9a and 9b for the formulations (1) and ( 2) respectively.The obvious difference between the two maps is quantified in Figure 9c, which reports the misfit .In general, we can see that, at a given time, the slip obtained with the formulation (2) (Figure 9b) is greater than that pertaining to the formulation (1) (Figure 9a).This is confirmed by the different values of the scalar seismic moment: (M 0 = 6.634 × 10 19 Nm and M 0 = 5.518 × 10 19 Nm, respectively) or momentmagnitude (M w = 7.15 and M w = 7.09, respectively).Note the difference between these values, compared to the small differences found in Sections 3.1, 3.2 and 4.1.
A first outcome of Figure 9c is that, contrarily to the previous cases (Sections 3.1, 3.2 and 4.1) in which the | m | is of the order of few percents, the variations can now exceed 75%.Moreover, we can note that the differences between the two models depend on the specific location of the fault.This is not surprising, giving the fact that the rake rotation is different on the whole fault plane, owing to the heterogeneous distribution of the initial shear stress.Finally, we observe that the differences tend to enhance as long as the rupture propagates.Indeed, near the external boundaries of the cracked region reported of Figure 9c m can even reach the maximum value 100%, indicating that, in a given fault node, u (path) is non zero, whereas u (mod) = 0.In the sliding logic of the linear SW law, this means that in the formulation (2) this specific fault node is slipping (and releasing stress), while in the formulation (1) that node is still at rest (i.e., not slipping).This translates into a time shift of the rupture fronts (or analogously into a shift of the rupture onsets), which increases through time, because of the effects of the dynamic stress redistribution on the fault plane.Overall, we can conclude that the whole rupture history is different in the two cases.
Figure 10a reports the fault striations in the receiver R indicated in Figure 9c.Contrarily to previous cases (Sections 3.1, 3.2 and 4.1), now the two components of the fault slip have roughly the same magnitude (this is an obvious consequence of the fact that both the components of the initial shear stress are non zero; see Figure 8).Remarkably, the two curves are now significantly different one from the other also below d 0 (see inset in Figure 10a).Moreover, they are not straight lines, suggesting that relevant rake rotation occurs in this fault node.The rake variation, reported in Figure 10b for both the formulation, is expressed as D{ = { − { 0 , where { 0 is the initial value of rake and the actual value { is computed as it follows: (3) (In such a way we capture also the rake rotation before the onset of the rupture, which corresponds to the first (negative) peak in each curve.)The shift between the two curves reflects the shift of the rupture onset we observe at t = 3.8 s (see Figure 9c).Remarkably, the most significant rake rotation (occurring in the correspondence of the traction degradation from x u to x f , Bizzarri and Cocco [2005]) can reach a maximum absolute value of the order of 20°.A relevant rake rotation is the cause of a different temporal evolution of the frictional resistance in the two cases.This in turn translates into a globally different behavior of the dynamic rupture, leading to the different total slips and seismic moments.From this simple example we can therefore conclude that the formulations (1) and ( 2) can lead to different results in case of highly heterogeneous simulations.
SLIP MODULUS OR SLIP PATH?

Conclusions
Contrary to rate-and state-dependent friction laws [e.g., Ruina 1983], the slip-dependent laws inherently pose a formal, geometrical problem in the definition of the independent physical observable on which the magnitude of the shear traction depends, the fault slip u.Indeed, we can define the amount of u at a given time t and in a given node of the fault surface (x 1 , x 3 ) in terms of its Euclidean norm (Equation ( 1)) or in terms of slip cumulated along its path (Equation ( 2)).
Since a priori we do not know -and we can not predict -whether these two formulations can lead to different results in temm of the resulting dynamic response of the fault system, in the present study we have quantified the differences.First of all, we remark that the ambiguity between slip modulus, u (mod) , and slip path, u (path) , only exists in 3-D or in 2-D mixed-mode problems (see Section 2.1 in Bizzarri [2011] for details), where we can have two non null components of slip, slip velocity and shear traction.In the case of 2-D problems (pure in plane or pure anti plane problems), originally used in the first applications of dynamic ruptures obeying the linear SW friction law [e.g., Andrews 1976aAndrews , 1976b]], u (mod) and u (path) are identical, if motion reversal (or back slip) is not allowed.Moreover, by definition of the linear SW law, the possible differences appear only during the breakdown phase of the rupture, i.e., when the fault traction evolves from the maximum level x u down to the residual level x f .Indeed, once the traction is at x f it does not depend on u anymore, so that the distinction between formulations (1) and (2) naturally disappears.It is obvious that the potential differences are affected by the value of the characteristic SW distance d 0 ; for small values of d 0 we expect that Equations ( 1) and (2) will provide practically the same results.This is important, since we know that additional mechanisms, such as chemical reactions (or thermochemical pressurization) can significantly reduce the characteristic SW distance in case of non linear slip-dependent friction [Chen et al. 2013].
When no rake rotation occurs, formulations (1) and (2) coincide by definition.We already know that homogeneous, subshear ruptures tend to have very small, or negligible, rake rotation.On the other hand, in homogeneous conditions, supershear earthquakes exhibit variations in the rake only near the supershear front [Bizzarri andCocco 2005, Bizzarri andDas 2012].The numerical results presented here show that, in homogeneous conditions and with a relatively large d 0 the misfit between the two solutions is less than 0.5% (see Figure 2b), the resulting total slip is the same, and so is the scalar seismic moment.Moreover, the fault striations are the same (see Figure 3a), as well as the fracture energy density (see Figure 3b) and the energy flux at the rupture front (see Figure 4).These results are also confirmed for configurations having a small initial stress value compared to the dynamic stress drop (Model B; see Figure 5).The general conclusion here is that, for homogeneous configurations, formulations (1) and (2) predict practically the same behavior of the dynamic rupture and thus we expect to have practically identical resulting ground motions.
These results are confirmed also in heterogeneous conditions, in which a high value of d 0 (1 m) is assumed in a portion of the fault.Compared to the homogeneous simulations more differences emerge near the  3a, but now in the case of Model C, having heterogeneous initial shear stress (see Table 1).The selected fault receiver is located at (26.25,7) km.(b) Time evolution of the rake variation (with respect to the initial value) for the two formulations.
(a) (b) rupture fronts (compare Figures 2b and 6b).In this case the differences are of the order of few percents, which is comparable to the grid dispersion error.
It is important to remark that the coincidence between the results pertaining to the two formulations (1) and ( 2), however, is not the rule.As a counterexample, in a case with a highly heterogeneous initial shear stress (Model C; see Figures 9 and 10) the rake rotation (see Figure 10b) is relevant enough to cause significant differences.The whole rupture history of the fault is quite different, indicating a more unstable behavior in the case of formulation (2).In other words, the slip modulus case tends to underestimate the stress release, the cumulated fault slip and thus the scalar seismic moment.Moreover, given the differences of the on-fault processes, we can expect to have some variations also in the synthetic seismograms on the free surface, especially at closer distances from the fault trace.
As an overall conclusion, we have demonstrated here that, in homogeneous cases, the ambiguity between slip modulus and slip path in the framework of the SW law is only formal and that the two different formulations provide the same results (with a tolerance of the errors of the numerical code).This contributes to reconcile previous papers, where the two formulations have been differently employed, depending on the preferences of the authors.From another point of view, this can be regarded as a future perspective that in homogeneous cases the analytical formulation of the linear SW law does not really matter.
On the other hand, in highly heterogeneous simulations, the normalized difference can even reach 100% (see Figure 9c), indicating that the choice of the analytical formulation of the SW plays some role and it affects the whole rupture history.It is difficult, in this case, try to highlight a general guidance.First, it should be noticed that kinematic inversions can provide, at maximum, with the time evolution of the shear traction change in some target fault points, but they are unable to clarify whether the slip path or slip modulus formulation is preferred.However and at a more fundamental level, if we consider the physical processes at the micro-asperities level and if we interpret their evolution only in terms of the cumulated fault slip (i.e., in terms of a macroscopic, averaged, slip-dependent law), then the ambiguity between u (mod) and u (path) would virtually disappear, in that the slip path would be more preferable than the slip modulus from a theoretical point of view and in agreement with the rate-and state-dependent friction governing models.

Figure 1 .
Figure 1.Schematic representation of the difference between slip modulus and slip path stated by Equations (1) and (2).Thin vectors represent the fault slip vector u at 3 subsequent, hypothetical time levels (denoted by superscripts).Thick vectors indicate the vectors vDt.Following the standard explicit Euler scheme, the 1-component of the fault slip is updated as u 1 (m) = u 1 (m−1) + y 1 (m−1) Dt, as reported in the sketch.(For the 3-component the situation is the same.)We have that u (mod) (m) = | |u (m) | |, while u (path) (m) = | |u (1) | | | |v (n) | |Dt.In the example displayed, at time level 3, u(mod) is identified by the magnitude by the thin blue vector u(3) , while u(path) is identified by the length of the thick green line.In the sketch it is also emphasized that, in general, u and v are not collinear; this is apparent from the difference of the two azimuths { u and { v .

Figure 2 .
Figure 2. (a) Distribution of the fault slip at time t = 4.2 s for a fault with homogeneous properties (the parameters are listed in Table 1) and in the case of the formulation (1).(b) Corresponding normalized difference between formulation (1) and formulation (2).In each fault node (x 1 ,x 3 ) we plot the misfit (*).The locations of the hypocenter (H) and of the fault received (R) are also indicated.Note that due to the symmetry of the problem only one half of the fault is reported.(*) m u u u u 100 ( ) ( ) ( ) ( ) mod mod

Figure 4 .Figure 3 .
Figure 4. (a) Time evolution of the energy flux F for the two models (F is computed from Equation (8) of Bizzarri [2013]).As for Figure 3 the thick black lines refer to the formulation (1), while the thin gray lines refer to the formulation (2).(b) Normalized difference between F from the two models; namely, in the ordinate axis we plot, for each time t, the misfit (*).In both the panels the vertical dashed gray line indicates the first time when the supershear transition occurs.(*) m F F F F 100 ( ) ( ) ( ) ( ) mod mod

Figure 5 .
Figure 5.The same as Figure 3, but now in the case of Model B. The fault receiver is the same as in Figure 3.

Figure 8 .
Figure 8. Spatial distribution of the initial shear stress adopted in Model C. (a) Along strike component.(b) Along depth component.

Figure 10 .
Figure 10.(a) The same as Figure3a, but now in the case of Model C, having heterogeneous initial shear stress (see Table1).The selected fault receiver is located at (26.25,7) km.(b) Time evolution of the rake variation (with respect to the initial value) for the two formulations.

Table 1 .
Parameters adopted in the present study.
* For models A and B only the x 1 -component of the initial shear traction vector T is non zero.