Discrete Fault Models

Fault surfaces are characterized by an inhomogeneous friction distribution, that can be represented with asperity models. Fault mechanics is dominated by asperities, so that a fruitful approach is to use discrete models, where asperities are the basic elements and the state of the fault is described by the average values of stress, friction and slip on each asperity. Under reasonable assumptions, the equations of motion can be solved analytically, with a deeper understanding of the behavior of the system. Fault dynamics has a sticking mode, where asperities are stationary, and a number of slipping modes, corresponding to the separate or simultaneous motion of asperities. Any seismic event is a sequence of slipping modes and a large variety of source functions is possible. Many large earthquakes are observed to be the consequence of the failure of two asperities: a discrete two-asperity model shows a rich dynamics and allows a detailed study of interaction between asperities. In this framework, fault evolution during coseismic and interseismic intervals can be calculated in terms of fault slip, stress state, energy release and seismic spectrum, including viscoelastic relaxation, fault creep and stress perturbations from other faults. Discrete models may include interaction between neighboring faults, allowing to assess conditions for the occurrence of seismic sequences in a fault system. A review of recent work on this subject is presented with applications to real earthquakes.


Introduction
Fault surfaces usually present a strongly nonuniform friction distribution, producing a remarkable complication in fault mechanics. An important achievement of seismic source theory is that friction distribution can be represented by asperity models (Lay et al., 1982;Ruff, 1983;Ruff & Kanamori 1983;Scholz, 1990). Such models assume that earthquakes result from slip of a small number of fault patches, characterized by high static friction and velocityweakening dynamic friction.
An inhomogeneous friction distribution entails a nonuniform distribution of coseismic slip. Seismic source functions retrieved from seismometric data show that fault slip usually takes place irregularly and suggest a complex interaction between different fault patches, resulting from a continuous stress transfer between them. One of the first events for which this behavior was observed is the 1992 Landers, California, earthquake (Kanamori et al., 1992;Wald & Heaton, 1994).
Heterogeneous fault models have been proposed by many authors. Models with nonuniform friction distributions under a uniform shear stress were considered by Miyatake (1978, 1995) and Beroza and Mikumo (1996). Somerville et al. (1999) showed that in most cases the slip distribution can be represented in terms of a small number of asperities. Pisarenko (2002) suggested that velocity weakening friction is a consequence of the dynamical failure of asperities. A stochastic model of fault slip complexity was developed by Mai and Beroza (2002). The case of nonuniform initial stress was considered by Cochard and Madariaga (1994), Ripperger et al. (2007) and Bizzarri et al. (2010). A multiple-asperity model was considered by Johnson (2010). Zielke et al. (2017) considered friction deriving from geometrical roughness of the fault surface. Numerical solutions were obtained in most of these studies.
A different approach is suggested by the evidence that asperities play a crucial role in fault mechanics. Processes controlling fault mechanics are stress accumulation on asperities, asperity strength and slip, as well as stress transfer between asperities. Therefore, the fault can be considered as a discrete dynamical system where asperities are the basic elements (Ruff, 1992;Rice 1993;Turcotte, 1997). Fault dynamics coincides with asperity dynamics and the state of the fault can be described by a finite number of variables characterizing asperities.
In fact, fault mechanics is controlled by a small number of asperities, so that models with a small number of degrees of freedom can be employed. Discrete fault models reproduce the essential features of the seismic source, but avoid the more complicated description based on continuum mechanics. The small number of degrees of freedom allows visualization of the state of the fault and of its evolution by inspection of orbits of the representative point in the phase space.
It should be noticed that discrete models considered here are different from models made of a large number of blocks discretizing the fault, such as those originally proposed by Burridge and Knopoff (1967). Those models approximate a continuous system with a system made of a finite, but very large number of degrees of freedom. On the contrary, models considered here are devised to study those aspects of fault mechanics resulting from a structure made of few strong patches, i.e. asperities. Average values of stress, friction and slip on each asperity are considered.
Originally, these models were conceived as loworder analogues of real faults, with asperities described as blocks sliding on a rough plane. However, equations governing these systems can be viewed as approximations of continuum mechanical equations, with a correspondence between model parameters and continuum mechanical quantities (Dragoni & Santini, 2012). In this formulation, an asperity is assumed to be a compact and simply connected subset of the fault surface and asperity failure can reproduce any kind of source mechanism.
In this framework, a fault with n asperities is a nonlinear, dissipative, piecewise smooth dynamical system with n degrees of freedom (di Bernardo et al., 2008). The evolution of such systems in the phase space is characterized by smooth intervals separated by rapid transitions, that can be considered as instantaneous events. Smooth intervals are the dynamic modes of the system and each mode is associated with a different set of differential equations. Such systems are named Filippov systems (Filippov, 1988).
A fault with n asperities has 2 n dynamic modes: a sticking mode, corresponding to stationary asperities, and 2 n À 1 slipping modes, corresponding to individual or simultaneous asperity motion. If asperities are labeled with integers from 1 to n, each mode can be indicated by a string of n digits i 1 i 2 . . .i n , that may assume values 0 or 1, where 0 denotes that the asperity is stationary and 1 that the asperity is slipping. A seismic event generated by the fault is a sequence of slipping modes.
Seismological data show that several large earthquakes occurred in the last decades can be ascribed to the failure of two asperities, from the 1964 Great Alaska Earthquake (Christensen & Beck, 1994) to the 2010 Maule earthquake (Delouis et al., 2010). A discrete fault model with two asperities has two degrees of freedom and four dynamic modes, and exhibits a very rich dynamics. Its evolution can be studied in detail during coseismic slip and interseismic intervals. Complex seismic events, that are sequences of several slipping modes, can be generated and initial conditions producing such events can be determined. Slip rates and shear stresses on asperities can be calculated as functions of time and their relationships can be investigated. Energy release, seismic moments and seismic wave spectra can be also calculated as functions of initial conditions.
The present paper summarizes the results obtained in a number of papers published on this subject in the last decade. Section 2 presents the model: a plane fault with two asperities subject to uniform strain rate from the motion of tectonic plates. The system variables and parameters are defined and the dynamic modes of the system are introduced. The general formulae for slip rates, moment rates, shear stress, seismic efficiency, seismic moments and spectra are given.
Section 3 illustrates the dynamics of the twoasperity fault, starting from the simplest case of identical asperities. The properties of long-term behavior and the conditions for the occurrence of large earthquakes are investigated. The case of two asperities with equal areas, but different strengths is then considered and the variety of possible source functions that it can generate is illustrated. Effects of wave radiation are highlighted by calculation of seismic energy and spectrum as functions of model parameters. The more general case of asperities with different areas and strengths is then considered and the role of stress perturbations in controlling fault evolution is examined. Finally, the effect of viscoelastic relaxation of the lithosphere during interseismic intervals is treated. Section 4 considers particular aspects of fault mechanics, such as the effect of variable strain rate, fault creep and afterslip, fault interaction and the occurrence of seismic sequences in a fault system. Section 5 gives a summary of the main results.

The Asperity Model
A plane fault with two asperities is considered. The asperities are named 1 and 2 respectively: they are disjoint subsets of the fault surface, with areas A 1 and A 2 (Fig. 1). It is assumed that the fault lies in a shear zone that is a homogeneous and isotropic Poisson solid with rigidity l. The shear zone has thickness d and is subject to a uniform strain rate _ e imposed by the motion of two tectonic plates at constant relative velocity v.

The Dynamical System
The state of an asperity at any time t is described by its slip deficit, defined as the slip that the asperity should undergo in order to recover the relative plate displacement occurred up to that time (Dragoni & Santini, 2015). The reason for using slip deficit is that it changes continuously even during interseismic intervals, when asperities are stationary. Accordingly, the state of the fault is described by slip deficits x(t) and y(t) of asperities 1 and 2, respectively. The model assumes the no-overshooting conditions implying that asperity slip does not exceed the plate displacement accumulated up to that time.
Since it is assumed that asperities move as rigid surfaces, their dynamics can be described by forces instead of tractions. Asperity slip is induced by forces exerted by the surrounding medium and is controlled by friction. Let f 1 and f 2 be the tangential forces applied to asperities in the slip direction. They can be written as (Lorenzano & Dragoni, 2018a) where terms ÀK 1 x and ÀK 2 y are the effect of tectonic loading, with Terms AEK c ðx À yÞ are contributions of stress transfer between asperities, with a coupling constant where s is the shear traction (per unit seismic moment) that the slip of one asperity imposes to the other, calculated at the asperity centroid (Lorenzano & Dragoni, 2018a). Let a be the distance between the asperity centroids. At distances a [ 1:5 ffiffiffi A p , the traction produced by a finite dislocation source of area A is indistinguishable from that of a point-like doublecouple source (Dragoni & Lorenzano, 2016). Accordingly, shear tractions are The two-asperity fault model. A rectangular fault surface with two asperities of areas A 1 and A 2 is shown. Asperities may have different strengths. The fault is subject to uniform strain rate, but shear stress is typically nonuniform on the fault. Adapted from Lorenzano and Dragoni (2018a) Vol. 179, (2022) Discrete Fault Models 3099 for a strike-slip mechanism and for a dip-slip mechanism (Love, 1944). Finally, terms Ài 1 _ x and Ài 2 _ y in (2) and (3) are due to radiation damping, where i 1 and i 2 are impedances, and contribute only during fault slip (Rice, 1993).
It is assumed that asperities obey a velocityweakening frictional law, with static frictions f s1 and f s2 and dynamic frictions f d1 and f d2 (Dragoni & Santini, 2012), that can be considered a simplified version of the more general rate and state-dependent law (Ruina, 1983;Dieterich, 1994;Scholz, 1998). Accordingly, conditions for failure of asperities 1 and 2 are, respectively The fault has four dynamic modes: a sticking mode (indicated by the string 00), where asperities are stationary, and three slipping modes, corresponding to slip of asperity 1 (10), slip of asperity 2 (01) and simultaneous slip (11). Each mode is described by a system of two differential equations. Equations for the sticking mode are and equations for the slipping modes are where m 1 and m 2 are masses associated with the asperities (Lorenzano & Dragoni, 2018a). Hence, the evolution of the system is described by (9) and (10) in mode 00, by (11) and (10) in mode 10, by (9) and (12) in mode 01 and by (11) and (12) in mode 11.

Nondimensional Formulation
It is useful to express the model variables and parameters in nondimensional form. To this aim, nondimensional variables and time are introduced (Lorenzano & Dragoni, 2018a): Accordingly, the state of the fault is described by nondimensional slip deficits X(T) and Y(T). The following nondimensional parameters are also introduced: where a is a coupling parameter for asperities; b is the ratio between the frictional stresses of asperity 2 and the same quantities of asperity 1; c is a measure of impedance associated with asperity 1; is the ratio between dynamic and static friction for both asperities; n is the ratio between asperity areas; V is the nondimensional velocity of tectonic plates. It is assumed that masses m 1 and m 2 are proportional to the respective asperity areas and that the impedance per unit area is the same for both asperities, so that Nondimensional forces are defined as Thanks to (2), (3) and (14), they can be written in terms of slip deficits as where dots indicate differentiation with respect to time T. Therefore, the evolution equations for the sticking mode become and those for the slipping modes become 3100 M. Dragoni Pure Appl. Geophys.
where a 0 ¼ a n ð23Þ Since the system has two degrees of freedom, its phase space is a 4-manifold S. The evolution of the system can be described by the orbit of its representative point (X, Y, _ X, _ Y) in S. The effect of radiation on fault dynamics is expressed by parameters (Dragoni & Santini, 2015;Lorenzano & Dragoni, 2018a) where If m(t) is the seismic moment produced in a seismic event as a function of time, the nondimensional moment is Let w(t) be the energy of the system at time t and r(t) be the radiation produced by asperity slip after a time t from the beginning of the event. The corresponding nondimensional quantities are Finally, if x is the angular frequency of seismic waves, the nondimensional frequency is

Initial Conditions
For most time, the system is in the sticking mode 00, where _ X and _ Y are equal to the velocity V of tectonic plates, hence extremely small with respect to the values they assume in slipping modes. Therefore, the representative point of the system in mode 00 can be assumed as belonging to the plane XY. The onset of seismic events is controlled by forces that are applied to asperities in mode 00. The noovershooting condition (1) implies that both F 1 and F 2 are in the direction of velocity v, or From (8), the conditions for failure of asperities 1 and 2 are respectively Therefore, the state of the system in mode 00 is constrained to a subset of the plane XY, whose borders are defined by conditions In terms of slip deficits, Eq. (32) can be written as that are the equations of four lines defining a The sticking region Q of the system in the plane XY. The sides of Q are defined by conditions on the forces F 1 and F 2 acting on asperities. The dashed line is Y ¼ X þ p 0 , separating initial conditions leading to failure of asperity 1 or 2 respectively. The region is drawn for a ¼ 1, b ¼ 0:5, c ¼ 0:5, ¼ 0:7, n ¼ 1.
Adapted from Dragoni and Tallarico (2016) Vol. 179, (2022) Discrete Fault Models parallelogram Q (Fig. 2). The subset Q, corresponding to stationary asperities, is called the sticking region (di Bernardo et al., 2008). Vertices of Q are the origin O and points P a , P b and P with coordinates (Lorenzano & Dragoni, 2018a) If we consider an initial point P 0 ¼ ðX 0 ; Y 0 Þ 2 Q, the evolution of the system from P 0 is given by solutions of (19) and (20), or These are the parametric equations of a line If one calculates the difference F 1 À F 2 at P 0 from (29), one obtains so that p is proportional to the difference between forces acting on asperities: hence it expresses the degree of stress inhomogeneity on the fault (Dragoni & Santini, 2010). A seismic event occurs when the representative point of the system, moving along line (40), intersects one of lines (34). The subsequent evolution depends on the value of p, expressing the initial condition for the seismic event. Values of p range in the interval Different initial conditions may produce very different seismic events. In particular, an event begins with mode 10 when p\p 0 and with mode 01 when p [ p 0 , where In many cases, a 1-mode seismic event, involving motion of a single asperity, is produced. More complex events are generated if p belongs to the narrower interval ½p 1 ; p 2 , where (Lorenzano & Dragoni, 2018a) where In this case, the motion of an asperity triggers the motion of the other one, originating a stress interchange between them: events that are sequences of three or more slipping modes are produced. In the particular case p ¼ p 0 , a 2-mode event 11-01 takes place: this is the largest event produced by the 2-asperity fault.

Seismic Events
During a seismic event, a continuous change in slip deficits X and Y, as well as in slip rates _ X and _ Y, takes place. Therefore, the event is represented by an orbit in the 4-manifold S. However, a seismic event can be more easily characterized by drawing the projection of its orbit on the plane XY, where the system dwells for most of its lifetime.
Each event is a sequence of m slipping modes and a segment of curve in the plane XY is associated with each mode. The functions X(T) and Y(T) relative to that mode are the parametric equations of the segment. If one introduces normal coordinates the solution for a generic mode starting at T ¼ 0 can be written as (Dragoni & Tallarico, 2016) where h 1 , h 2 , a 1 , a 2 , b 1 and b 2 are constants depending on the particular mode and on initial conditions, while X 1 and X 2 depend on model parameters. These are the parametric equations of a damped Lissajous curve (Lawrence, 1972). The orbit representing an m-mode event is the union of m segments of such curves, separated by singular points. The sequence of modes in the event and the associated slip amplitudes are evident from inspection of orbits in the plane XY. Let T i (i ¼ 1; 2; . . .m) be the instant of time when the system enters the i-th mode and T mþ1 be the end time of the event. The duration of an m-mode seismic event is then The duration of each mode, as well as the total duration of the event, may vary sensibly from one event to the other (Dragoni & Tallarico, 2016). Let X i and Y i be the slip deficits at time T i . The slip amplitudes in the i-th mode are then and the total slip amplitudes are respectively The dynamics of a seismic event is evident if we consider the slip rates of asperities and the moment rate of the event. In each slipping mode, slip rates are given by different functions for each asperity. Let X i ðTÞ and Y i ðTÞ be the slip deficits in the i-th mode. The slip rates of asperities are The slip rates in an m-mode seismic event are then where H(T) is the Heaviside function. For each event, the moment rate can be calculated from (55) and (56). Let M 1 be the seismic moment of a 1-mode event 10 in the absence of radiation (c ¼ 0), corresponding to the slip amplitude U defined in (46). The moment rate of a generic m-mode event is then Very different moment rates result depending on initial conditions (Dragoni & Tallarico, 2016). The total seismic moment of an m-mode event is The moment rate spectrum S of a seismic event is defined as the magnitude of the Fourier transform of _ MðTÞ: The energy of the system at a point ðX; YÞ 2 Q is (Dragoni & Santini, 2015) EðX Let E i be the energy at the beginning of the ith mode. The energy change in the ith mode is then and the total energy change is or, thanks to (52), In the ith mode, seismic energy is released at a rate (Dragoni & Santini, 2015) and the energy released in the ith mode is The total seismic energy released in the event is Finally, the seismic efficiency can be defined as

Dynamics of a Two-Asperity Fault
The dynamics of a two-asperity fault can be highlighted by considering a number of cases of increasing complexity. In all cases, analytical solutions for the dynamic modes of the system can be obtained.
According to Turcotte (1997), chaotic behavior is attained for higher values of the coupling parameter a. From (14), (4), (5) and (6) or (7), it results that a is proportional to the ratio A=a 3 , where A ¼ A 1 or A 2 . Since asperities are disjoint subsets of the fault surface, larger areas imply larger values of distance a. Therefore the high values of a required for chaotic behavior are hardly attained and this kind of behavior is not considered here. Analysis of seismological data is not conclusive in this regard (e.g. Marzocchi et al., 1997).

Identical Asperities
The simplest case is a fault with two identical asperities (Dragoni & Santini, 2010). In this case Then the parameters b and n are equal to 1. For the sake of simplicity, wave radiation is neglected, a reasonable assumption in view of the low seismic efficiency of faults (Kanamori, 2001): hence c ¼ 0. Therefore, Eqs. (21) and (22) The system has a rich variety of behaviors even in this simpler case. An analysis of orbits in the phase space shows that the system has an infinite number of limit cycles, each one characterized by the variable p defined in (41). An initial value p þ U gives rise to the same cycle, but run in the opposite way (Fig. 3). An interesting aspect is that the recurrence pattern of seismic events depends on the initial stress distribution on the fault (Dragoni & Santini, 2010). Each cycle is made of two 1-mode events 10 and 01, separated by interseismic periods of variable duration. The slip amplitude U given in (46) is associated with each event. From (59) In the particular case when one interseismic period is equal to zero, a single 2-mode event 10-01 or 01-10, with seismic moment 2M 1 takes place in the cycle.
There is no simultaneous slip of asperities.
Orbits reach a limit cycle only after entering a particular subset L of the sticking region Q. This means that only a subset of initial stress distributions on the fault allows the system to enter a limit cycle.

Figure 3
A limit cycle of the two-asperity fault, projected on the plane XY. Segments P 1 P 2 and P 3 P 4 represent 1-mode seismic events 10 and 01 respectively. Segments P 2 P 3 and P 4 P 1 represent interseismic intervals of different durations. The cycle can be run in both senses depending on the initial value of the variable p (À0:17 or 0.13), expressing stress distribution on the fault. Model parameters are Dragoni and Santini (2010)  In this case, the behavior is periodic, with alternate asperity slip. When the representative point of the system is outside L, the motion of one asperity can trigger the motion of the other one and simultaneous motion takes place (mode 11). This happens because, while an asperity is slipping, the condition for slip of the other asperity is attained. In this case, (69) and (70) must be solved simultaneously (Dragoni & Santini 2011). When simultaneous slip takes place, asperity 1 may trigger the failure of asperity 2 or vice versa: the two cases produce earthquakes with similar seismic moments, but different epicentres. The seismic moment of events including mode 11 is always larger than the maximum value 2M 1 produced in a limit cycle. Events including simultaneous slip of asperities are the largest earthquakes that can be generated by the fault.
In the long term, the representative point reaches in any case the subset L and enters a limit cycle, with periodic behavior. Of course, this contrasts with observation, showing that the seismic activity of a fault is aperiodic and produces earthquakes of different magnitudes. This behavior can easily result from the model if one assumes that the fault is not isolated, but may receive stress transfers from the seismic activity of neighboring faults.
The system is sensitive to small perturbations, because the sticking region Q can be divided into narrow stripes leading to different evolutions (Dragoni & Santini, 2011). If the system is in L, a small stress perturbation can shift the system from one limit cycle to another having a different recurrence pattern. The system may also be shifted outside L: in this case, one or more larger events, with simultaneous asperity motion, will be produced, until a stress distribution leading to periodic behavior is restored. Therefore periodicity could not be observed.

Asperities with Different Strengths
In the case of a fault with two asperities with equal areas, but different strengths, we have still n ¼ 1 but b 6 ¼ 1. The equations for slipping modes are Solution of these equations shows that the system has a greater range of possible behaviors (Dragoni & Santini 2012). As in the case of equal strengths, the fault can produce earthquakes from the failure of one asperity or may involve failure of both asperities, but the symmetry of section 3.1 is lost. When one asperity is involved in a seismic event, the following event is originated in most cases by the failure of the other asperity, but in some cases failure of the same asperity may occur. In particular cases (Fig. 4), an earthquake may be due to slip of one asperity followed by slip of the other one and again by slip of the first one (Dragoni & Santini, 2012). At any instant of time, the evolution of the fault depends on the variable p, that is proportional to the difference between forces F 1 and F 2 according to (42). For the largest part of the interval ½p a ; p b , the difference is remarkable. It is concluded that the initial stress distribution on the fault is not uniform in most cases.
As an example, Dragoni and Santini (2012) considered the 1964 Great Alaska Earthquake, one of the largest events in the last century, with An example of orbit for a 3-mode seismic event 01-11-01, projected on the plane XY. Initial condition is p ¼ À0:15. In this case p 2 ½p 1 ; p 2 , so that the event includes mode 11. Each segment of the orbit is labeled with the string indicating the slipping mode. The dashed line is Y ¼ X þ p 0 , where p 0 ' À0:17. Model parameters are a ¼ 1, b ¼ 0:5, c ¼ 0:5, ¼ 0:7, n ¼ 1. Adapted from Dragoni and Tallarico (2016) Vol. 179, (2022) Discrete Fault Models magnitude 9.2. Seismological, geodetic and tsunami data show that the earthquake was due to the failure of two large asperities (Christensen & Beck, 1994;Holdahl & Sauber 1994;Santini et al. 2003;Zweck et al., 2002). On the basis of data, a coupling parameter a ¼ 0:1 was assumed. Friction distribution was estimated from the ratio of maximum slip amplitudes for the two asperities, yielding b ¼ 0:75. With an appropriate value of p, the event can be modelled as a mode sequence 10-01. If one assumes that after an earthquake the state of stress is controlled only by tectonic loading, the model predicts that next earthquakes will involve only one asperity: hence they will be smaller than the 1964 one, until conditions for a new large event are attained. This prediction leans on simplifying assumptions and approximate data, but it shows that discrete fault models, focusing on large-scale properties, may disclose mechanisms controlling the longterm evolution of faults (Dragoni & Santini, 2012).
The source function of a seismic event is usually expressed by the moment rate _ M as a function of time. A study of the moment rate of events produced by a two-asperity fault shows how the moment rate changes as a function of initial conditions and model parameters (Dragoni & Santini, 2014). For given values of parameters, the moment rate function depends on the sequence of slipping modes in the event, which is controlled by the value of p. The graph of moment rate appears as a sequence of humps, each one corresponding to a slipping mode, that may partially overlap in connection with simultaneous slip (Fig. 5).
A seismic event including simultaneous asperity slip was the 2010 Maule earthquake, a magnitude 8.8 thrust event occurred in central Chile (Delouis et al., 2010;Lay et al., 2010a;Vigny et al., 2010). Slip concentrated on two main asperities located south and north of the epicentre. According to observations, the event can be modeled as a mode sequence 10-11-01, involving asperities with similar strengths (b ¼ 1). The model provides a good fit of the observed moment rate and an estimate of stress distribution on the fault before and after the earthquake (Dragoni & Santini, 2014).
A consequence of friction heterogeneity on faults is that complex events can be produced (Dragoni & Tallarico, 2016). With the exception of the case p ¼ p 0 , all events start with the motion of one asperity. The initial mode is 10 or 01 according to whether p is smaller or greater than the value p 0 given in (44). After some time, the asperity motion triggers the motion of the other asperity and the system passes to mode 11. When asperities slip simultaneously, they never stop at the same time, so that mode 11 is followed by mode 10 or 01.
When p has the smallest values in the interval ½p a ; p b , the mode sequence is 10-11-01. As p increases, mode 11 becomes prevailing over modes 10 and 01, until an inversion in the order of arrest takes place in mode 11: asperity 2 stops earlier than asperity 1, so that the sequence is 10-11-10. When p ¼ p 0 , the initial mode 10 disappears and asperities slip simultaneously from the onset of the event. Asperity 1 stops earlier, so that a 2-mode sequence 11-01 results. For p [ p 0 , the event starts with mode 01 and the sequence is 01-11-01. A greater complication results for the largest values of p, such as a 5-mode sequence 01-11-10-11-01. Different sequences are possible depending on the values of model parameters (Dragoni & Tallarico 2016).
The evolution of stress on the fault is described by forces F 1 and F 2 on asperities, given in (18). They increase with time on both asperities in the sticking mode, but oscillate in the slipping modes, owing to two processes: stress drops on asperities and stress transfers between them (Fig. 6). These processes determine the state of stress at the end of the event. Depending on the initial stress distribution and on the mode sequence, the final stress can be more or less heterogeneous than the initial one.
An example of complex event is 1992 Landers, California, earthquake, a magnitude 7.3 event originated by strike-slip faulting (Kanamori et al., 1992;Olsen et al., 1997;Peyrat et al., 2001;Wald & Heaton, 1994). The observed slip distribution is very heterogeneous, but can be modelled in terms of a friction distribution made of two asperities with different strengths (b ¼ 0:5). The model shows that the initial stress distribution was strongly heterogeneous, with a reduced heterogeneity at the end. With an appropriate choice of model parameters, the moment rate for a 2-mode event fits the observed function (Dragoni & Tallarico, 2016).

Seismic Energy and Spectrum
Consideration of wave radiation introduces further effects. Radiation during slipping modes can be taken into account by addition of a term proportional to slip rate (Rice, 1993). If asperities are assumed to have the same area (n ¼ 1), equations for slipping modes are In the presence of radiation, fault dynamics changes both in the sticking and in the slipping modes, in a measure that depends on the impedance c (Dragoni & Santini, 2015). Emission of radiation changes the evolution of the system from a given state, because it shifts the boundaries between different subsets of the sticking region Q. In particular, the presence of radiation changes the width of the interval ½p 1 ; p 2 leading to simultaneous asperity slip. In the presence of radiation, slip amplitude is smaller, while slip duration is longer. The moment rate depends on seismic efficiency g and has its maximum in the first half of source duration DT.
The smaller value of slip implies a smaller seismic moment M 0 , that decreases with increasing g, at constant radiated energy. It is found that (Dragoni & Tallarico, 2016) for a 1-mode event 10 and for a 1-mode event 01, where j 1 is given in (24). Events including mode 11 have a moment where f(p) must be evaluated numerically (Fig. 7).
The greatest values of M 0 are achieved by these events, with a maximum at p ¼ p 0 . The seismic spectrum SðXÞ defined in (60) exhibits the classical Brune (1970) shape (Fig. 8). It has an infinite number of relative minima, their positions and values depending on b. They are zeroes only in the case of homogeneous friction (b ¼ 1). The corner frequency is (Dragoni & Santini, 2015) where x 1 is defined in (25) and DT is the source duration. It depends on parameters a, b and c and can sensibly change as a function of initial conditions. If one considers again the 1964 Great Alaska Earthquake, mode durations, slip distribution, moment rate and seismic moment can be calculated in the presence of radiation (Dragoni & Santini 2015): they are consistent with observed values (Ichinose et al., 2007) provided an appropriate value of p is chosen. The state of the fault after the earthquake is different from that obtained in the absence of radiation, implying that the long-term evolution is different.
During fault slip, elastic strain energy is partly dissipated into heat and partly radiated in the surrounding medium (Kanamori & Heaton 2000;Kanamori & Rivera 2006;Rudnicki & Freund, 1981). The model shows that the energy release DE decreases, while the radiated energy DR increases with increasing impedance c. Accordingly, seismic efficiency g given by (68) increases with increasing c, its maximum value depending only on the ratio between dynamic and static friction (Dragoni & Santini, 2015).
The effect of friction heterogeneity on seismic energy and spectrum can be studied by considering events produced by the failure of asperities with different strengths (Dragoni & Santini 2017). As underlined above, the stress distribution on the fault, expressed by b, is in general strongly heterogeneous at the onset of a seismic event. Seismic energy DR decreases with increasing b, while seismic efficiency g is constant. An equation relating g to the parameters of the friction law can be obtained, showing that g is maximum for smaller values of . The model provides a relation between DR and the seismic moment M 0 , that is consistent with the empirical relation between the two quantities (Kanamori, 2001). It results that heterogeneity introduces a correction to the value of energy radiated by a homogeneous fault (Dragoni & Santini, 2017).
As an example, the 1965 Rat Islands, Alaska, earthquake, was considered, a magnitude 8.7 event occurred at the Alaska-Aleutinian Trench (Beck & Christensen 1991;Wu & Kanamori, 1973). The event can be modeled as due to the separate failure of two asperities with different strengths, with b ¼ 0:67. For this event, moment rate, stress evolution, radiated energy and seismic spectrum can be calculated on the basis of the model, showing how friction heterogeneity controls the characteristics of a seismic event (Dragoni & Santini, 2017).

Asperities with Different Areas and Strengths
In the case of a fault with two asperities of different areas and strengths, the equations for slipping modes have the general forms given in (21) and (22). There are remarkable differences with the case of asperities with equal areas (Lorenzano & Dragoni, 2018a).
Force rates on asperities are not equal to each other and their difference is not constant during interseismic intervals. The sticking region Q grows with the total asperity area, while the probability that the fault produces events involving simultaneous asperity slip decreases. If one considers events produced by the failure of a single asperity, slip duration and amplitude increase with asperity size, while corner frequency decreases. In 2-mode events involving consecutive failure of asperities, seismic efficiency g depends on the total asperity area.
As an application of the model, Lorenzano and Dragoni (2018a) considered the 2007 Pisco, Peru, earthquake, a magnitude 8.0 event at the border between the Nazca and South American plates (Lay et al., 2010b). Seismological and geodetic data indicate the presence of two asperities (Sladen et al., 2010). By an appropriate choice of initial conditions, the event can be modelled as a 2-mode sequence 01-10, starting with slip of the weaker asperity and followed by slip of the stronger one, with b ¼ 0:5 and n ¼ 0:6. The moment rate fits reasonably well the observed function.
The case of simultaneous asperity motion was considered by Santini and Dragoni (2020). The aim was to model the observed moment rate and seismic moment of the 2018 Gulf of Alaska earthquake, a magnitude 7.9 event originated by strike-slip faulting (Krabbenhoeft et al., 2018;Lay et al., 2018;Ruppert et al., 2018). Observations indicate that the earthquake was due to the failure of two main asperities and suggest an intermediate phase of simultaneous slip. Accordingly, the earthquake can be ascribed to a 3-mode sequence 10-11-01. Asperity sizes are very different from each other, with n ¼ 0:25. An estimate of frictional strengths yields b ¼ 0:45. The moment rate calculated from the model fits observations reasonably well.

Effects of Stress Perturbations
In a fault system, any fault is subject to stress perturbations due to earthquakes generated by neighboring faults. The stress redistribution produced by each earthquake affects occurrence times and magnitudes of following earthquakes. Therefore, stress transfer between faults plays an important role in fault systems (Stein et al., 1992;Harris, 1998;Stein, 1999;Steacy et al., 2005;Tallarico et al., 2005). Belardinelli et al. (2003) discussed the effects of stress perturbations in the case of a homogeneous fault. Dragoni and Piombo (2015) considered the case of a fault containing two asperities. Fault heterogeneity produces effects that are not present in the case of a homogeneous fault. This occurs because the stress field of a dislocation is inhomogeneous, so that asperities belonging to a fault are subject to different stress changes. As a consequence, a stress perturbation may not only advance or delay the next earthquake, but it may change the sequence of dynamic modes in the event, thus changing the hypocentre position, source duration and seismic moment.
The proximity of a fault to failure can be expressed by its Coulomb stress. If r t is shear stress in the direction of fault slip and s s is static friction on the fault surface, the Coulomb stress is (Gomberg et al., 2000;Stein, 1999) Accordingly, r C is negative during interseismic intervals and seismic events occur when r C vanishes.
Since there are two asperities, a value of Coulomb stress must be assigned to each of them. If forces are used instead of stresses, nondimensional Coulomb forces on the two asperities are where F 1 and F 2 are given by (2) and (3) and n ¼ 1 is assumed for the sake of simplicity. While F C 1 and F C 2 measure the proximity of asperities to the respective failure conditions, the difference F C 2 À F C 1 indicates the proximity to the condition of simultaneous slip.
When the fault is subject to a stress perturbation, the changes DF C 1 and DF C 2 in Coulomb forces are different from each other, so that the evolution of the fault is controlled by the difference expressing the change in the asymmetry degree of stress and friction on the fault. These results explain why earthquakes generated by a fault are not only aperiodic, but are different from each other as to hypocenter position, slip amplitude and involved area. Differences are connected with the kind of stress perturbations occurring during interseismic intervals (Dragoni & Piombo, 2015). Vol. 179, (2022) Discrete Fault Models 3109 The model has been applied to the fault of the 2010 Maule earthquake, that was considered in Sect. 3.2, in order to investigate the effect of the perturbation produced by the 1960 Great Chilean Earthquake (Dragoni & Piombo, 2015). As noted above, the stresses imposed on the asperities of the Maule fault were necessarily different from each other: this altered the stress distribution on the fault, with important consequences for its evolution.
The model predicts that, in the absence of the 1960 earthquake, the Maule earthquake would have occurred several decades later, with a different sequence of dynamic modes and a different seismic moment. The 1960 earthquake increased stress inhomogeneity on the Maule fault, preparing the conditions for the 2010 earthquake.
The effects of viscoelastic relaxation on a twoasperity fault were considered by Amendola and Dragoni (2013) in the case of identical asperities and by Dragoni and Lorenzano (2015) in the case of asperities with different strengths. As reported in Sect. 3.1, in the case of purely elastic coupling the long-term behavior of the system is a limit cycle with a recurrence pattern of earthquakes depending on the degree of stress inhomogeneity on the fault. In the presence of viscoelastic relaxation, this simple behavior is modified, because stress transferred from one asperity to the other partially relaxes during interseismic intervals.
In order to take into account viscoelasticity, a third variable z, expressing viscoelastic deformation, is introduced in the model. In nondimensional form Hence the system has three degrees of freedom and phase space S is 6-dimensional. In the sticking mode, tangential forces on asperities are where terms AEaZ are contributions of stress transfer between asperities in the presence of viscoelastic deformation. If asperities with equal areas are considered for the sake of simplicity (n ¼ 1), conditions for the onset of slip are from (31) or, thanks to (84), that are the equations of two planes in the space XYZ.
If one considers the no-overshooting condition (30), the sticking region is defined from the six inequations that define a convex tetrahedron T (Fig. 9). A seismic event takes place when the representative point of the system reaches face ACD or BCD, entering mode 10 Figure 9 The sticking region T of the two-asperity fault in the presence of viscoelastic relaxation. T is a tetrahedron ABCD in the space XYZ. Seismic events take place when the representative point of the system reaches face ACD or BCD. Model parameters are a ¼ 1, b ¼ 1, c ¼ 0, n ¼ 1. Adapted from Dragoni and Lorenzano (2015) 3110 M. Dragoni Pure Appl. Geophys. or 01. If it reaches the edge CD, the system enters mode 11. The simplest viscoelastic model describing longterm behavior of lithospheric rocks is the Maxwell body (e.g. Ranalli, 1995), characterized by a relaxation time s, that is expressed in nondimensional form as The evolution equations in mode 00 are then where Z is determined by the Maxwell constitutive equation. Hence where X 0 , Y 0 and Z 0 are the coordinates of the system at T ¼ 0. From (86) and (90), the time T 1 required for the fault to reach the condition for slip of asperity 1 is where W is the Lambert function with argument An analogous formula holds for asperity 2, with where From (91) and (93), it is easy to find that initial states belonging to T lead to mode 10 or to mode 01 according to whether they are below or above the surface From (84) and (90), forces applied to asperities in mode 00 are A comparison with the purely elastic case shows that there are additional nonlinear terms, so that forces change non-monotonically during the interseismic interval. Additionally, the difference F 1 À F 2 changes in time, entailing a change in stress distribution on asperities. In summary, seismic events are anticipated or delayed with respect to the elastic case, the importance of viscoelastic relaxation being controlled by the product VH (Dragoni & Lorenzano, 2015). In order to illustrate this effect, the fault that originated the 1964 Great Alaska Earthquake can be considered again. Being a large-size event, the earthquake was followed by a remarkable postseismic deformation (Zweck et al., 2002). As reported in Sect. 3.2, the earthquake source can be represented as a 2-mode event 10-01. Consideration of viscoelastic effects shows that stress relaxation controls the occurrence times of earthquakes produced by that fault (Amendola & Dragoni, 2013).
Fault evolution depends on the state from which the 1964 event originated, that is inferred from the observed moment rate. This constrains the evolution of the system to a subset of phase space. Knowledge of the moment rate of the next earthquake would further constrain the orbit, and so on (Dragoni & Lorenzano, 2015).
There is a complex interplay between stress perturbations and viscoelastic relaxation. Following a stress perturbation, the change in Coulomb stress on a given asperity determines the anticipation or delay of slip of that asperity, if a purely elastic behavior is assumed. This property no longer holds in the presence of viscoelastic relaxation (Lorenzano & Dragoni 2018b). In the latter case, from (81) and (84), Coulomb forces on asperities are Following a perturbation in normal stress, static frictions on asperities are modified. If f 0 s1 and f 0 s2 are the new static frictions on asperity 1 and 2, respectively, the following parameters are defined Vol. 179, (2022) Discrete Fault Models 3111 Changes in static friction imply that conditions (85) for the onset of asperity slip become meaning that changes in normal stress change the sticking region. Then the changes in frictions can be written as If DF 1 and DF 2 are the changes in tangential forces due to the perturbation, the changes in Coulomb forces are The signs of DF C 1 and DF C 1 determine whether the perturbation takes an asperity closer to or farther from failure conditions. However, there is not a direct connection between the signs of the two quantities and the anticipation or delay. The effect depends on the state of the fault immediately before and after the stress perturbation (Lorenzano & Dragoni 2018b).
An interesting case is the stress perturbation produced by the 1999 Hector Mine, California, earthquake on the fault originating the 1992 Landers earthquake, the latter being followed by a remarkable viscoelastic relaxation (Lorenzano & Dragoni 2018b). The 1999 Hector Mine earthquake was a magnitude 7.1 event generated by strike-slip faulting located about 20 km from the Landers fault (Jónsson et al. 2002;Salichon et al., 2004). As reported in Sect. 3.2, the Landers earthquake source can be modelled as a 2-mode 01-10 event. The stress transfer produced by the 1999 Hector Mine earthquake can be calculated and the complex effects of the stress perturbation on the future activity of the Landers fault can be explored.

Other Aspects of Fault Mechanics
The model developed in previous sections can be adjusted in order to investigate other aspects of fault mechanics, such as the effects of variable strain rate, the presence of fault creep and the occurrence of seismic sequences in fault systems.

Variable Strain Rate
In the study of fault behavior it is generally assumed that faults are subject to a constant strain rate. However, velocities of tectonic plates may change in the very long term, probably due to changes in the rate of mantle convection (e.g. Iaffaldano & Bunge, 2009;King et al., 2002). Therefore strain rates controlling seismic activity are also bound to change in the long term.
Studies of long-term correlations of earthquakes have been based on different approaches: addition of Brownian perturbations to steady tectonic motion (Matthews et al., 2002;, use of the concept of self-organized criticality (Abaimov et al., 2007;Baiesi, 2009) or study of stress evolution in discrete fault models (Ben-Zion et al., 2003;. Dragoni and Piombo (2011) considered the effect of a slowly variable strain rate on the activity of a fault with a single asperity. The state of the fault is described by slip deficit X(T) and can assume two dynamic modes, denoted by 0 and 1 respectively. Two cases were considered: a sinusoidal oscillation in velocity of tectonic plates and a slow transition between two velocity values. It was assumed that such oscillations or transitions have small amplitudes with respect to the average velocity and longer periods and durations than recurrence times of earthquakes.
In the case of constant plate velocity V, equations governing fault dynamics in the two modes are respectively (19) and (69) with a ¼ 0 (radiation is neglected): In this case, the fault produces seismic events with a recurrence period where In the case of variable plate velocity, slip deficit changes with a variable rate _ XðTÞ in mode 0. Since interest is focused on the interseismic intervals, seismic events can be considered as instantaneous.
In the case of sinusoidal oscillation with frequency X, fault activity is made of cycles including several seismic events and repeating periodically with a period that is a multiple of 2p=X (Dragoni & Piombo, 2011). Within each cycle, recurrence times of events oscillate about an average value equal to DT 0 and the oscillation amplitude is proportional to that of strain rate oscillations. The number of events in a cycle depends on the ratio between X and 2p=DT 0 . In the case of monotonic transition between different values of strain rate, recurrence times change gradually from an initial to a final value.
If the fault is subject to stress perturbations, the subsequent earthquake is anticipated or delayed. In the case of sinusoidal oscillations in strain rate, the perturbation will interrupt the current seismic cycle and will start a new one. Hence the pattern of seismic cycles determined by strain rate oscillations is destroyed by frequent stress perturbations, as may occur in systems made of several faults. In the case of a monotonic transition in strain rate, perturbations similarly alter the pattern of recurrence times and the number of events occurring during the transition (Dragoni & Piombo, 2011).
Therefore, even though slow variations in strain rate are difficult to observe in seismicity records, they significantly contribute to the aperiodicity of seismic events in the long term.

Fault Creep and Afterslip
It is often observed that fault slip continues for some time after an earthquake, although at a decreasing rate, a phenomenon called afterslip. Afterslip is interpreted as aseismic slip of a velocity-strengthening region of the fault (Belardinelli & Bonafede, 1995;Marone et al., 1991;Scholz, 1990). This is confirmed by seismic and geodetic observations, indicating that faults can accommodate tectonic motion with stable, quasi-static slip or with fast slip and production of seismic waves. Such observations can be accounted for if one considers two kinds of regions on the fault surface: stable regions, which mostly creep, and unstable regions, producing earthquakes (Johnson, 2010).
The time dependence of aseismic slip has been mostly described by empirical relationships. An exponential function approaching a constant value was proposed by Nason and Weertman (1973). Later, observations and theoretical considerations suggested a logarithmic function (Marone et al., 1991). A review of different time functions was given by Barbot et al. (2004).
The interaction between two fault segments due to aseismic slip was studied by Dragoni and Tallarico (1992) and Tallarico et al. (2002) in the framework of continuum mechanics. A discrete fault model with two mechanically different regions was considered in Dragoni and Lorenzano (2017). The two regions are an asperity and a weak region and the state of the fault is described by their slip deficits x and y.
In this section, quantities referring to the asperity and to the weak region are denoted by indices 1 and 2, respectively. Let f 0 1 and f 0 2 be the frictional resistances of the two regions. For the asperity, a velocity-weakening law is assumed as above, characterized by a static friction f s and a dynamic friction f d . For the weak region, a velocity-strengthening law is assumed where f 0 is the steady-state dynamic friction and K is a constant. Nondimensional slip deficits are where K 1 is defined as in (24). Nondimensional parameters are defined as where a new parameter k has been intoduced and the others have the same meaning as in two-asperity models. The evolution equations in the interseismic intervals are Vol. 179, (2022) Discrete Fault Models 3113 When the two regions slip, the evolution equations are respectively The equations are solved analytically for interseismic intervals, asperity slip and afterslip in the weak region (Fig. 10). During interseismic intervals, the asperity is stationary, while the weak region is creeping. Stress accumulates on the asperity, until it is released when frictional threshold f s is exceeded. Asperity slip transfers stress to the weak region and afterslip takes place. In turn, afterslip transfers stress back to the asperity, determining the conditions for the next earthquake (Dragoni & Lorenzano, 2017). According to the model, afterslip is a dynamic mode of the fault and its time dependence can be approximated by a function where k is a function of a, k and n. In the long term, afterslip approaches an asymptotic value U a and may have any duration, depending on the intensity of velocity strengthening. The amount of afterslip U a is proportional to seismic slip of the asperity, in agreement with observations (Dragoni & Lorenzano, 2017). The model was applied to the fault of the 2011 Tohoku-Oki earthquake (Ide et al., 2011;Maercklin et al., 2012;Simons et al., 2011), that was followed by a prolonged afterslip episode (Ozawa et al., 2011). According to observations, seismic slip concentrated at shallow depth, while afterslip took place downdip (Lay et al., 2012;Silverii et al., 2014;Wei et al., 2012).
On these grounds, a fault with a single shallow asperity and a downdip weak region was considered (Dragoni & Lorenzano, 2017). With a suitable choice of model parameters, a good fit of the seismic moment rate is obtained. The model suggests that afterslip dominated the first months after the event, while later postseismic deformation was due to bulk viscoelastic relaxation (Sun et al., 2014;Yamagiwa et al., 2015).

Seismic Sequences
A typical manifestation of seismic activity is the occurrence of seismic sequences. This term is used to mean a series of earthquakes produced by sources located in a relatively small region and occurring in a time interval much shorter than the intervals between sequences (Dragoni & Lorenzano, 2016). Sequences take place in fault systems producing earthquakes with similar mechanisms and magnitudes. A sequence is typically made of a small number of larger events having medium magnitudes, accompanied by a greater number of smaller events.
Since the faults of the system are close to each other, fault interaction plays a major role, concentrating the events in a shorter time interval. The distribution in time of seismic sequences has been investigated by multifractal analysis (e.g. Telesca et al., 2004;Telesca & Lapenna, 2006).
If we consider a system made of n faults, the state of the system can be described by n variables that are

Figure 10
Representation of the cycle made of seismic slip, afterslip and interseismic creep in the space XY. The dashed and dotted lines represent conditions for asperity failure and fault creep, respectively. Points P 1 , P 2 and P 3 represent the state of the fault at the beginning of the seismic event, at the end of the event and at the end of afterslip, respectively. Segments P 1 P 2 and P 2 P 3 give the amplitudes U s and U a of seismic slip and afterslip, respectively. From P 3 to P 1 , the fault is subject to creep. Adapted from Dragoni and Lorenzano (2017)  the Coulomb stresses of the faults (Dragoni & Lorenzano, 2016). If faults are ordered according to the magnitude of their Coulomb stresses, each state of the system can be expressed by a permutation of the first n integers, such as where fault i 1 has the maximum value of Coulomb stress and fault i n has the minimum. This permutation changes whenever a fault produces a seismic event, so that the evolution of the system can be described as a sequence of permutations. When a seismic sequence of n events takes place, the order of fault activation is a consequence of initial stress state of the fault system and of fault interaction. One may conceive an order implicit in the initial state, that is modified due to changes in Coulomb stresses occurring whenever an event takes place. There are n! possible sequences, differing for the order of fault activation, and a sequence itself can be expressed as permutation where faults i 1 , j 1 ...k 1 are in the order of activation (Dragoni & Lorenzano, 2016). Thanks to the model, the stress state of a fault system as a function of time can be retrieved from observation of the order of fault activation in a seismic sequence. It results that consecutive seismic sequences originated by a fault system are necessarily different from each other, because the state of the system at the end of a sequence is always different from the initial one.
The model has been applied to the 2012 Emilia (Italy) seismic sequence, that was made of seven events with magnitudes between 5 and 6 (Castro et al., 2013;Pezzo et al., 2013;Scognamiglio et al., 2012). According to the model, stress evolution in the fault system was conditioned by the first and fourth events, that were greater than the others, and produced a greater stress heterogeneity at the end of the sequence. The model predicts that, in the absence of external perturbations, the next sequence will take place after a few centuries and will be different from the 2012 one (Dragoni & Lorenzano, 2016).

Conclusions
Discrete fault models based on the concept of asperity have proven to be a useful tool in the study of fault mechanics. In such models, the state of the fault is described by few variables expressing the state of asperities, so that the dynamical system is characterized by few degrees of freedom.
The present review has focused on the dynamics of a fault surface with two asperities. This is a common case, because many large earthquakes can be interpreted as due to the failure of two asperities. A two-asperity model is apparently simple, but reveals an unexpected richness of complicated behaviors.
Due to the presence of friction, the system is nonlinear and evolves as a sequence of dynamic modes, where each mode is described by a different set of differential equations. Analytical solutions can be obtained and the evolution of the system can be visualized in a low-dimensional phase space.
The main conclusions are: 1. The fault can produce very different seismic events according to initial conditions, expressed by the stress distribution on the fault. Events are sequences of slipping modes, involving the separate or simultaneous slip of asperities, with durations ranging over a wide time interval. 2. The system is sensitive to small changes in initial conditions and to small stress perturbations. Hence its evolution is unpredictable in the long term, even if there is no chaotic behavior. 3. Most initial conditions of seismic events correspond to strongly nonuniform stress distributions on the fault. This should be considered the typical case in fault activity. 4. The moment rate of seismic events may have very different shapes, according to initial conditions. Observation of the moment rate of a real event allows to constrain initial conditions to a subset of the larger set of possible conditions. 5. Interaction between asperities plays a crucial role during fault slip: continuous stress transfer between asperities controls the sequence of dynamic modes in a seismic event. 6. The stress distribution on the fault after a seismic event is always different from that before the Vol. 179, (2022) Discrete Fault Models 3115 event: hence the next event will be different from the previous one. 7. Radiation of seismic waves has an influence on fault dynamics and the frequency spectrum of emitted waves exhibits remarkable variations as a function of initial conditions. 8. Viscoelastic stress relaxation occurring during interseismic intervals plays an important role, because it changes the stress distribution on the fault and may lead to a completely different evolution. 9. If one of the asperities is replaced with a velocity strengthening region, the latter is subject to afterslip as a consequence of coseismic slip of the asperity. Afterslip is a possible dynamic mode of the fault. 10. Generalization of the model to a fault system shows that the order of events in a seismic sequence is a consequence of initial stress state and of interaction between faults. The stress states of the system can be retrieved from observation of the order of fault activation in the sequence.
Discrete models are of course a simplification of real faults. Their strength is to catch the essential features of fault dynamics and to provide analytical solutions, enlightening the details of processes occurring both in coseismic and interseismic intervals. The possibility of visualizing the evolution in the phase space adds a long-term overview that is precluded to models based on continuum mechanics.

Acknowledgments
The author is grateful to the editor Carla F. Braitenberg and to anonymous reviewers for acceptance of the paper.

Funding
Open access funding provided by Alma Mater Studiorum -Università di Bologna within the CRUI-CARE Agreement.

Declarations
Conflict of interest No funds, grants, or other support were received during the preparation of this manuscript. The author has no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.