Conditions for the occurrence of seismic sequences in a fault system

We consider a fault system producing a sequence of seismic events of similar magnitudes. If the system is made up of n faults, there are n! possible sequences, differing from each other for the order of fault activation. Therefore the order of events in a sequence can be expressed as a permutation of the first n integers. We investigate the conditions for the occurrence of a seismic sequence and how the order of events is related to the initial stress state of the fault system. To this aim, we consider n coplanar faults placed in an elastic half-space and subject to a constant and uniform strain rate by tectonic motions. We describe the state of the system by n variables that are the Coulomb stresses of the faults. If we order the faults according to the magnitude of their Coulomb stresses, a permutation of the first n integers can be associated with each state of the system. 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. A crucial role is played by the differences between Coulomb stresses of the faults. The order of events implicit in the initial state is modified due to changes in the differences between Coulomb stresses and to different stress drops of the events. We find that the order of events is determined by the initial stress state, the stress drops and the stress transfers associated with each event. Therefore the model allows the retrieval of the stress states of a fault system from the observation of the order of fault activation in a seismic sequence. As an example, the model is applied to the 2012 Emilia (Italy) seismic sequence and enlightens the complex interplay between the fault dislocations that produced the observed order of events.


Introduction
Seismic sequences are a characteristic aspect of seismic phenomenology.Recent examples in Italy are the 1997-1998 Umbria-Marche sequence (Morelli et al., 2000;Salvi et al., 2000;Santini et al., 2004) and the 2012 Emilia sequence (Scognamiglio et al., 2012;Castro et al., 2013;Pezzo et al., 2013).We call "seismic sequence" a series of earthquakes generated by sources located in a relatively small region (of the order of 100 km) and occurring in a time interval (of the order of a few months) much shorter than the intervals during which the system is at rest.The interval elapsing between two seismic sequences in the same region is of the order of several decades at least.We do not include in this definition aftershock sequences following a greater event that may have similar features but are strongly conditioned by the main shock.
Sequences are originated by fault systems that produce similar earthquakes as to mechanism and magnitude.A sequence is typically made up of a small number (< 10) of larger events having a medium magnitude, in general between 5 and 6, plus a greater number of smaller events.We take into account the larger events, neglecting the smaller ones.
Since the faults of the system are close to each other, it is believed that fault interaction plays a major role in seismic sequences.Fault interaction and its role in earthquake triggering have been widely studied (e.g.Stein et al., 1992;Harris, 1998;Stein, 1999;Gomberg et al., 2000;Belardinelli et al., 2003;Steacy et al., 2005).Whenever a fault slips, it transfers stress to neighbouring faults, thus anticipating the instant of time when they slip in their turn.Therefore the interaction has the effect of concentrating the events in a shorter time M. Dragoni and E. Lorenzano: Seismic sequences in a fault system interval, hence reducing the duration of the sequence (e.g.Tallarico et al., 2005).
In the present paper, we consider a model of a fault system and investigate the conditions under which the system may originate seismic sequences with the characteristics described above.In particular, we ask the following questions.
1. Which are the stress conditions under which a sequence may take place?
2. What determines the order of seismic events in the sequence?
3. What makes the order of events change from one sequence to the following one?
4. Is the observed order of events informative about the state of stress before and after the sequence?
The model is applied to the 2012 Emilia (Italy) seismic sequence that was made up of seven major events with similar focal mechanisms and with magnitudes between 5 and 6, occurring in a time interval of 15 days.

The model
We consider a system made up of n plane faults that we assume to be coplanar and lined up, with the same strike and dip angles (Fig. 1).The fault system is placed in an elastic half-space with Lamé constants λ and µ.We introduce a coordinate system (x, y, z) such that the x axis coincides with fault strike, y is the horizontal direction perpendicular to strike and z is depth.Let δ be the dip angle of the faults.
We number the faults from 1 to n, starting from one end of the system.Let A i be the area of the ith fault and r ij be the distance between the centres of the ith and j th faults.We introduce the following assumptions: 1. the fault system is subject to a strain rate ė that is constant in time and uniform in space; 2. the onset of seismic events is controlled by the average values of tangential traction and static friction on fault surfaces; 3. fault slip is a step function of time and does not produce overshooting; 4. each fault slips once and only once during a sequence; 5. there is no simultaneous slip of two or more faults and a finite time interval elapses between the failures of any two faults; 6. the duration of a sequence is much shorter than the interval between two consecutive sequences; and 7. the system is not subject to external stress perturbations.
Most of these assumptions are commonly introduced in seismic source models.Assumption 1 is reasonable, since by definition the n faults belong to the same seismogenic region, for which the same tectonic mechanism is observed.Assumptions 2 and 3 are based on the fact that we are not interested in the details of each event, which has a much shorter duration than the duration of the sequence, but rather in the relationship between the n events.Assumptions 4, 5 and 6 are suggested by the features of the sequences we are describing.Sequences are made up of distinct events, each one associated with the failure of a distinct fault in the system, and there is no observation of reactivation of the same fault during a sequence.The duration of a sequence is typically several weeks or a few months, while the interval between two consecutive sequences may even be centuries long (Rovida et al., 2011).
As to assumption 7, it is a fact that the evolution of a fault system can be altered by external perturbations.Any fault system is not isolated, but is surrounded by other faults that may transfer stress to it whenever they slip (e.g.Dragoni and Piombo, 2015).Generally, contributions from external faults may be numerous during an interseismic interval, but they are smaller than contributions from faults belonging to the system, due to greater distances and to different orientations of fault surfaces.Such contributions may also partially cancel each other.
In the case of normal and reverse faults, we assume plane strain, according to the Anderson model (Anderson, 1951;Sibson, 1974;Turcotte and Schubert, 2002).The nonvanishing strain components are e yy = ėt, where ė is positive for tensile strain and negative for compressive strain.The stress components are where ν is the Poisson modulus.We introduce the stress rate The rates of normal and tangential traction on the faults are then where the upper sign in σ t is for normal faults and the lower sign is for reverse faults.In the case of transcurrent faults, we consider simple shear, with strain and stress components e xy = ėt, In this case, we define and the rates of normal and tangential traction on the faults are Let σ i be the average tangential traction applied to the ith fault in the slip direction and τ i be the average static friction of the ith fault.We define the Coulomb stress (Stein, 1999) of the ith fault as Since the σ i are always positive or zero, the x i range between −τ i and zero.When x i = 0, an earthquake is generated by the ith fault.The rates of σ i and τ i are, respectively, where κ is the coefficient of static friction.Then the rate of Coulomb stress is where for normal and reverse faults and for transcurrent faults.Then, in the absence of earthquakes, the Coulomb stress of the ith fault changes in time as where x 0i is the Coulomb stress at an arbitrary time t = 0. Due to the presence of friction, the set of n faults is a nonlinear dynamical system.The Coulomb stresses x i can be considered as the components of an n-dimensional vector x(t) describing the state of the system as a function of time.
The possible states of the system belong to an n-dimensional parallelepiped S, defined by the n disequalities According to assumption 5, all the components of x are different from each other.Therefore one (and only one) component will vanish first, generating the first event in the sequence.Whenever an earthquake occurs, the fault dislocation produces a static stress field that is transferred to the system and modifies the Coulomb stress of all faults, producing a sudden change in x.
In general, the change in Coulomb stress on the j th fault due to the failure of the ith fault can be written as where σ ij is the coseismic change in tangential traction, H is the Heaviside function and σ ij is the change in traction due to time-dependent processes, including pore fluid diffusion, after-slip and viscoelastic relaxation.Since we have assumed that faults are coplanar, there are no changes in normal stress on the fault plane.
The traction σ ij could be calculated from the formulae for a rectangular dislocation in an elastic half-space collected by Okada (1992).However, if r ij > 1.5 √ A i , the traction of a finite dislocation source is virtually indistinguishable from that of a point-like double-couple source in an unbounded medium (Love, 1944) and the latter simpler formula can be used (Appendix A).
A Poisson solid (λ = µ) is considered.Accordingly, if m i is the seismic moment of the ith fault, we have in the case of a strike-slip mechanism and in the case of a dip-slip mechanism.As to the stress change of the ith fault, it is where σ i is the static stress drop that can be estimated from the average slip u i and the fault area A i as where C is a nondimensional constant of the order of unity determined by the geometry of the fault (Kanamori, 2001).
According to discrete fault models (e.g.Dragoni and Piombo, 2015), the stress drop is a fraction of static friction, where is the ratio between the average values of dynamic and static frictions that we assume to be the same for all faults.
If the medium is porous and saturated with fluids, the coseismic stress field induces a fluid flow that changes the stress field in turn (e.g.Wang, 2000;Piombo et al., 2005).As shown in Appendix B, the effect of fluid diffusion is at least 1 order of magnitude smaller than coseismic stress transfer: for the sake of simplicity, we do not consider it in the following.However, in some cases pore fluid diffusion may have a role in the evolution of a seismic sequence (e.g.Convertito et al., 2013).
As to after-slip and viscoelastic relaxation, the events we are considering are relatively small, and it is assumed that they do not produce appreciable after-slip or impose considerable stress on deeper ductile regions that may relax it afterwards.Viscoelastic relaxation of lithospheric rocks may change the stress distribution in the long term as a consequence of larger earthquakes (e.g.Dragoni and Lorenzano, 2015).
3 Evolution of the system Let t k be the occurrence times of the events in the sequence (k = 1, 2, . . .n), so that the durations of the interseismic intervals are Then the initial state is x(t 1 −).If the first event is due to the failure of the i 1 th fault, x has a sudden change and its k component becomes Afterwards, x changes continuously in time, as a consequence of tectonic loading, according to At t = t 2 , the second event takes place, due to the failure of the i 2 th fault, so that x has another sudden change, and so on.
At the end of the sequence, the state is that can be written as where is the duration of the sequence that can be written as In Eq. ( 25) the difference between the final and initial states is made up of two terms: the first one is tectonic loading during the time interval t; the second one is the effect of earthquakes.The latter term has the effect of concentrating in a shorter time interval a series of events that otherwise would be farther in time.The shortening in duration is obtained by calculating how much the instant t n of the last event is anticipated.The decrease in t n is due to the sum of the stresses that are transferred to the i n th fault from the other n − 1 faults.From Eq. ( 27), the duration of the sequence in the absence of interaction is The interseismic intervals in Eq. ( 21) can be calculated as During the interseismic intervals, the representative point x moves along a line defined by the parametric Eq. ( 13) that is parallel to the line Thanks to a rotation R, the coordinate system (x 1 , . . .x n ) can be changed into a system (ξ 1 , . . .ξ n ) such that the ξ n axis coincides with line (Eq.30).Hence the evolution of the system can be more easily represented in the (n − 1)-dimensional hyperplane ξ n = 0.An example will be shown in Sect.7 for n = 3.

Retrieval of the initial and final states
On the basis of the model, if we observe a seismic sequence, we can retrieve the state vector x(t) at any time during the sequence.In particular, we can calculate the state of the system at the beginning and at the end of the sequence.Suppose that we observe a sequence made up of n events that can be ascribed to the failure of n faults belonging to the same system.Let t 1 , t 2 , . . .t n be the observed occurrence times of the events.From the knowledge of the fault geometry and of the seismic moments, we can calculate the stress transfer matrix x ij .If we know the strain rate ė from geodetic measurements, we can calculate the stress rate ẋ from Eq. ( 10).
Consider the generic fault i k that has produced the kth event in the sequence.It is easy to see that the Coulomb stress of fault i k at the beginning of the sequence is Apart from the signs, the first term in the rhs is the stress accumulated on the fault from the beginning of the sequence up to instant t k and the second term is the sum of stress transfers that fault i k has received from faults i 1 , i 2 , . . .i k−1 that slipped before it.Hence, apart from the sign, the rhs is the total stress accumulated on fault i k since the beginning of the sequence.This stress must cancel the initial Coulomb stress x i k (t 1 −): hence the initial Coulomb stress must be the opposite of the accumulated stress.
As to the final state of fault i k , it is given by Eq. ( 25).If we replace x i k (t 1 −) in Eq. ( 25) with its expression (Eq.31), we obtain Since the Coulomb stress of fault i k was equal to zero at t = t k −, the final stress is equal to the tectonic stress accumulated in the time interval from t k to t n plus the stress drop associated with the failure of fault i k and the stress transfers of the faults i k+1 , i k+2 , . . .i n that have slipped after fault i k .
Hence the initial Coulomb stress of fault i k depends only on what happened before the instant t k , while the final Coulomb stress depends only on what happened after t k .However, the retrieval of the complete state vector requires the knowledge of the entire sequence.In Sect.8, we shall retrieve the initial and final states of a fault system in a real case.
The degree of heterogeneity of the x i can be expressed by their standard deviation where A relevant point for the subsequent evolution is whether the differences between the x i change during a sequence.We define The d ij form an antisymmetric matrix having n(n − 1) nonvanishing components that are related by (n − 1) 2 equations.Therefore d ij is known if we know only n − 1 components, for example the d 1j with j = 2, 3, . . .n. Thanks to Eq. ( 32), we obtain where the rhs is different from zero because the sum of stress transfers received by a fault during the sequence is in general different from that received by the other faults.In particular, faults located at the centre of the system receive a greater total stress than faults located at the ends, if the events have similar magnitudes.For instance, if n = 3, fault 2 receives a greater stress transfer than faults 1 and 3.
5 The order of events Since the components x i of the state vector are always different from each other, they can be ordered according to their magnitudes.Then, at any instant t of time, the set X of the x i (t) is a well-ordered set.This order controls the order of events in the seismic sequence.Let N n be the set of the first n natural numbers.With each state x of the system we can associate a permutation α of N n , expressing the order of faults in relation to the value of their Coulomb stress: so that with k = 2, 3, . . .n. Hence the parallelepiped S can be divided into a number n! of subsets S j corresponding to the n! permutations of N n .During the interseismic intervals, the permutation α j associated with the system does not change, because all the x i increase with the same rate, according to Eq. ( 13).Therefore x remains in the same subset S j .However, when an event occurs, x switches to a different subset S k , characterized by a permutation α k .Suppose that, before the sequence, the permutation associated with the system is implying that the first event is the failure of the i 1 th fault.This event changes the magnitudes of all Coulomb stresses, so that the new state of the system is associated with a different permutation www.nonlin-processes-geophys.net/23/419/2016/ Nonlin.Processes Geophys., 23, 419-433, 2016 implying that the second event is the failure of the j 1 th fault, and so on.After the (n − 1)th event, the permutation is implying that the last event is the failure of the k 1 th fault.Therefore the order of events in the sequence can be expressed as a permutation As to the order of events, the number of possible sequences in a system made up of n faults is equal to n!.Since every fault may slip only once in a sequence (assumption 4), there are n! alternatives for the initial permutation α 0 , but only (n − 1)! for α 1 and (n − k)! for the generic permutation α k .
If the permutation after the nth event is the duration of the interseismic interval preceding the next sequence is and the sequence will start with the failure of the ith fault.In order to find out the relationship between the initial permutation α 0 and the order of events given by α * , it is necessary to examine which are the orders of magnitude of the quantities controlling the evolution of the system.

Discussion
The evolution of the system is controlled by the stress rate σ and by the stress transfer matrix x ij .We calculate the typical values of these quantities for seismic sequences.For many seismogenic regions, typical strain rates are of the order of 10 −15 to 10 −14 s −1 .With µ = 30 GPa and ν = 0.25, Eq. (3) and Eq. ( 6) yield stress rates | σ | of the order of 2 kPa a −1 for the lower value and of 20 kPa a −1 for the higher value of the strain rate.Calculation of x ij requires the knowledge of the seismic moments associated with each fault, of the fault areas and of the distances between them.A typical seismic moment of an event in the sequence can be calculated by assuming an area A i = 100 km 2 and an average slip u i = 0.5 m; hence, m i 10 18 N m.With these values, Eq. ( 19) yields a stress drop σ i 1 MPa.This value corresponds to the stress that is accumulated in a time interval t 0 = 500 a at a rate of 2 kPa a −1 or t 0 = 50 a at a rate of 20 kPa a −1 .From Eq. ( 20) with a typical value = 0.7 (Scholz, 1990), we obtain σ i = 0.6 τ i .
As to the distance between the centres of two neighbouring faults, we may roughly assume that r = 2 √ A i = 20 km.
Then, according to Eqs. ( 16) and ( 17), the stress σ ij transferred from a fault to its first neighbours is of the order of 10 kPa, so that the maximum value of σ ij (i = j ) is of the order of one-hundredth of a stress drop.It must be noted that a greater value of m i does not entail a proportionally greater value of σ ij , because it implies a greater fault area and a greater distance between faults.
An obvious effect of fault interaction is the shortening of time intervals between seismic events: for neighbouring faults, the gained time σ ij / | σ | ranges from 5 to 0.5 a according to the value of σ .Hence the maximum stress transferred by one event is equivalent to several months or several years of tectonic loading.
The magnitude of the differences d ij between Coulomb stresses is critical for the occurrence of a sequence.A lower limit for d ij is set by the magnitude of transferred stress σ ij (i = j ).If d ij is smaller than σ ij , the failure of the ith fault would immediately produce the failure of the j th fault, in contrast to assumption 5. Hence a condition for having a sequence of n distinct events is d ij > σ ij at any time.
An upper limit for d ij is set by the observed durations of seismic sequences.The d ij must be small enough that a sequence is completed within a few months, if we take into account the effect of stress transfer between faults.Hence we may assume as an upper limit for d ij the stress change ẋ δt that tectonic loading produces in a time δt T (assumption 6) plus the sum of transferred stresses x ij (i = j ).A greater value for δt (several decades) can be assumed for lower stress rates, a smaller value (several years) for higher stress rates.
With this premise, we consider how the order α * of events in a sequence is determined.The quantities determining α * are the initial stress state of the fault system, the stress drops and the stress transfers associated with each event.
The simplest case is when the stress drops are about equal to each other and the d ij are always greater than the transferred stresses x ij (i = j ).If these conditions are fulfilled, the only effect of the kth event is to shift the label i k to the last position in the permutation α k , while the stress transfers x i k j do not change the relative positions of the other labels.So we can associate with each event a permutation It follows that the order of events is given by the initial permutation, i.e. α * = α 0 .The final permutation α n is also equal to α 0 , but this does not imply the repetition of the order α * in the following sequence.According to Eq. ( 36), the new sequence will start with different values of d ij that may produce a different order of events.
Apart from the case just described, the order implied by α 0 is generally changed during a sequence, because the d ij have the same order of magnitude as the x ij (i = j ).In addition, if an event j has a stress drop that is considerably greater than the others, the label i j will permanently occupy the last position in the permutation, thus altering the initial order.It follows that the order of events is different from the initial order of stresses, i.e. α * = α 0 .The final permutation α n is also different from α 0 .
7 An example: the case n = 3 As an example, we consider a system made up of three faults with a strike-slip mechanism.This case is considered because it can be illustrated graphically, owing to the small number of variables involved.Cases with n > 3 would require higher-dimensional spaces.The graphical representation allows a better understanding of the evolution of the state of the fault system during a seismic sequence.
For the sake of simplicity, we suppose that the faults are equal to each other, with distances r 12 = r 23 = 20 km between their centres.With a strain rate ė = 10 −14 s −1 , the stress rate calculated from Eq. ( 6) is σ 19 kPa a −1 .For the sake of simplicity, we suppose that the faults have the same static friction τ = 1 MPa and produce events with the same seismic moment m 0 = 10 18 Nm.The stress transfer matrix Eq. ( 15) is symmetric, with nondiagonal components x 12 17 kPa and x 13 2 kPa.With a typical value = 0.7, Eq. (20) yields stress drops x i = 600 kPa.
We consider coordinates x i /τ .The parallelepiped S is a cube with unit edge, defined by the disequalities and can be divided into six subsets S j corresponding to the six permutations of N 3 .During the interseismic intervals, the representative point x(t) of the system draws an orbit parallel to the line An event occurs whenever the point reaches one of the coordinate planes.
As anticipated in Sect.3, we can rotate the coordinate system so that the ξ 3 axis coincides with line (Eq.49).It is easy to see that where The projection of S on the plane ξ 3 = 0 is a regular hexagon H with side a (Fig. 2a).It is divided into six equilateral triangles H k that are the projections of the subsets S k .We can follow the evolution of the system by looking at the projection P of the representative point on H .During the interseismic intervals, P does not change, because the representative point moves on a line perpendicular to H , and is close to the origin, because the d ij are much smaller than τ .Whenever an event takes place, P moves to a different subset of H . Suppose that, at a certain instant t 0 of the interseismic interval preceding a sequence, the state of the system is with x 01 > x 02 > x 03 .The associated permutation is then and the representative point P 0 belongs to the subset of H labelled with 123 in Fig. 2b.We choose the vector x 0 in order that the sequence is made up of three distinct events, occurring in the order given by α 0 .Then d 12 and d 13 must be positive, with d 12 < d 13 .According to the discussion in the previous section, where we assume δt = 1 a, so that σ δt 19 kPa.Then 17 kPa < d 12 < d 13 < 38 kPa and we choose d 12 = 20 kPa and d 13 = 25 kPa.At the beginning of the sequence, the state is The mean is x = −15 kPa, with a standard deviation s = 10 kPa.The state immediately after the first event is According to Eq. ( 47), the associated permutation is and the representative point P 1 belongs to the subset of H labelled with 231 in Fig. 2b.From Eq. ( 29), the time interval between the first and second events is t 1 = 66 d.At the end of this interval, the state is and the second event takes place.The state becomes The associated permutation is α 2 = ηα 1 and the representative point P 2 belongs to the subset of H labelled with 312   in Fig. 2b.The time interval between the second and third events is t 2 = 56 d.At t = t 3 − the state is The third event takes place and the state becomes with α 3 = ηα 2 , which coincides with α 0 .The representative point P 3 belongs to the subset of H labelled with 123 in Fig. 2c.According to Eq. ( 26), the duration of the sequence is t = 122 d.In the absence of interaction, the duration would have been t = 482 d from Eq. ( 28).The evolution of the three components in time is shown in Fig. 3.
The final state is very different from the initial one, with a value of x that is about 40 times greater, while s has not changed.The differences between components are d 12 = 5 kPa and d 13 = 25 kPa, so that x 2 has become closer to x 1 and farther from x 3 .In fact, the state no longer fulfills the condition d 12 > x 12 , entailing that, when the next sequence occurs, the stress x 12 transferred to fault 2 by the slip of fault 1 will induce the immediate failure of fault 2.
According to Eq. ( 45), the duration of the time interval before the next sequence is T 30 a.Even though the order of events in this sequence is the same as in the previous one, the durations of the intervals t k are different.Moreover, the stress distribution is altered in such a way that α 3 = α 0 , entailing that even the order of events will be different in the following sequence.
The example shows that, even in a simple case such as the one considered, the changes in d ij occurring in a sequence inevitably change the characteristics of the following sequences.A further source of change intervenes when the events have different sizes, as shown in the following section.

The 2012 Emilia sequence
We consider the 2012 Emilia (Italy) seismic sequence that was made up of seven events with magnitudes between 5 and 6 et al., 2013).They occurred in the period between 20 May and 3 June 2012, and can be ascribed to a fault system approximately lined up in the west-east direction, with a total length of about 50 km.The faults are all of thrust type and with shallow hypocentres between 5 and 10 km in depth (Fig. 4).
The Emilia sequence has been studied in detail by several authors, who determined the locations and seismic moments of the events (Scognamiglio et al., 2012), the source func- tions and seismic spectra (Castro et al., 2013) and the coseismic deformation (Pezzo et al., 2013).Convertito et al. (2013) suggested that dynamic triggering caused by seismic waves might be the primary factor to explain the evolution of the Emilia sequence, in addition to the variation in permeability and pore-pressure effects due to a massive presence of fluids in the Po Plain basin.
As stated in Sect.2, we neglect the effect of pore fluid diffusion, on the basis of considerations in Appendix B, as well as the possible effect of seismic waves.Our aim is not to simulate the Emilia sequence in detail, but to use it as an example of a complex sequence for which the present model can afford the retrieval of the initial and final stress states.If further effects are relevant and are introduced in the calculations, they may alter the sequence of permutations and yield a different final permutation.However, they would not change the general conclusions of the paper.
According to the model, we approximate the real fault system with a set of n = 7 faults having the same strike and dip angles and the same average depth.If we number the faults from west to east, the order of events is given by α * = 1 2 3 4 5 6 7 5 6 7 4 3 2 1 . (62) Therefore, fault slip started about the middle of the system and propagated eastward up to the end of the system (5, 6, 7); then it propagated from the middle to the west end (4, 3, 2, 1) (Fig. 5).
The data required for retrieving the initial and final states according to Eqs. ( 31) and ( 32) are the elastic and frictional properties of the medium (the rigidity µ, the Poisson ratio ν, the coefficient κ of static friction), the geometry of the faults (the areas A i , the distances r ij , the dip angle δ), the strain rate ė, the occurrence times t i and the seismic moments m i of the events.As to the elastic and frictional properties, we take µ = 30 GPa, ν = 0.25 and an effective coefficient of friction κ = 0.6.The areas and the locations of the faults have been inferred by employing the distances between hypocentres along the strike direction as constraints (Caporali and Ostini, 2012;Serpelloni et al., 2012).The distances between the centres of the faults are r 12 = r 23 = r 67 = 5 km, r 34 = r 56 = 8 km, and r 45 = 12 km.The projection of the faults on a vertical plane is shown in Fig. 6.
We treat all sources as pure reverse, dip-slip faults with δ = 40 • an average of the values given by Convertito et al. (2013).The strain rate is ė = −3 × 10 −15 s −1 (Caporali and Ostini, 2012).The moments m i are derived from the moment magnitudes reported in Tramelli et al. (2014), whereas fault slips u i are calculated from m i and A i .The data are shown in Table 1.From Eq. ( 19) with C = 1, the values of stress drops σ i are in the range between 0.9 and 1.9 MPa, consistent with the range evaluated by Castro et al. (2013) from seismic spectra.
From these data we calculate the rate ẋ of Coulomb stress and the stress transfer matrix x ij .From Eq. (10), ẋ 2 kPa a −1 .The initial state x(t 1 −) and the final state x(t 7 +) are shown in Fig. 7.The initial and final permutations are, respectively, α 0 = 1 2 3 4 5 6 7 5 4 7 1 2 3 6 , (63) α 7 = 1 2 3 4 5 6 7 6 7 1 2 4 3 5 . (64) The evolution of the system during the sequence shows that Eq. ( 47) does not hold for any value of k.Therefore α 7 is different from α 0 and they are both different from α * .This is a consequence of the heterogeneous distribution of seismic moment in the fault system.In particular, the evolution of stress was conditioned by the first and fourth events, due to the failures of faults 5 and 4, respectively, having greater seismic moments than the average.As a consequence of the greater stress drop, fault 5 permanently occupies the last position in all permutations from α 1 to α 7 .The stress transfers associated with each event also play a role in determining the evolution of the sequence, contributing to the rearrangement in the permutations.Due to the many stress transfers to fault 6, Eq. ( 28) yields t 46 a for the duration of the sequence in the absence of fault interaction, a much longer time than the observed duration t 15 d. Figure 6 shows that, at the beginning of the sequence, the mean Coulomb stress was x −0.05 MPa, with a standard deviation s 0.03 MPa, whereas at the end of the sequence x −1.2 MPa and s 0.4 MPa.Therefore Coulomb stresses are much more spread out after the sequence than before, since the standard deviation is 1 order of magnitude larger.
According to Eq. ( 64), the faults with the highest values of x i after the sequence are the sixth, seventh and first, showing that, in the absence of external perturbations, the next sequence would start in the proximity of one end of the system, rather than in the middle as the 2012 sequence.According to Eq. ( 45), the next sequence will take place after an interseismic interval T 440 a.This figure appears to be representative of typical recurrence times of moderate-size earthquakes in this area: the largest event before the 2012 sequence was the M w 5.5 17 November 1570, Ferrara earthquake (Rovida et al., 2011).

Conclusions
The aim of this study was to enlighten the conditions allowing the occurrence of seismic sequences and the processes controlling the order of events in a sequence.When we observe a sequence, we acknowledge that it is due to a system of n faults that fail one after the other.However, we do not know why the faults fail in that particular order.The order must be a consequence of the initial stress state of the fault system and of the mutual interaction between faults.
In order to unravel this problem, we introduced the concept of permutation of the n faults, ranking the faults according to the magnitudes of their Coulomb stresses.Such a permutation describes the state of the system at a given time and changes whenever a fault is activated.The order of activation itself can be described by a particular permutation of the faults.We have shown that the knowledge of the order of ac-Table 1.Data for the seismic events of the 2012 Emilia sequence.The origin times and the seismic moments m i are taken from Pezzo et al. (2013) andTramelli et al. (2014), respectively.The areas A i take into account the analysis of Caporali and Ostini (2012) and Serpelloni et al. (2012).Fault slips u i are calculated from m i and A i .See Fig. 6  tivation of faults in the sequence yields information on the state of the fault system before and after the sequence.The evolution of the fault system during a seismic sequence can be better understood if we consider the state space.Since a permutation can be associated with each state of the system, the state space can be divided into n! subsets corresponding to the n! permutations.The permutation does not change as long as the system is at rest, so that the state remains in the same subset of the state space.Whenever a seismic event takes place, the order of Coulomb stresses is changed and is expressed by a different permutation: an event corresponds to a switch of the state vector to a different subset of the state space.
We are now in a position to answer the questions we asked in the Introduction.
1.The occurrence of seismic sequences requires particular stress conditions.A crucial role is played by the differences between the Coulomb stresses of faults.A lower limit for these differences is set by the magnitude of transferred stresses and an upper limit by the observed durations of seismic sequences.This constrains their values in a narrow interval of the order of tens of kPa.
2. The order of events in a sequence is determined by the initial distribution of Coulomb stresses, by the stress drops and by the stress transfers associated with each event.The order of events that is implicit in the initial state is generally modified by the changes in the state vector intervening during the sequence.The dominant contribution to stress changes is given by stress drops that are typically 100 times greater than stress transfers.However, stress transfers have a major role both in anticipating the occurrence times of the events and in altering the order of events implicit in the initial state.
3. The characteristics of consecutive sequences originated by a fault system are bound to change.The order implicit in the initial stress distribution is generally changed during the sequence, because the stress transfers between faults have the same order of magnitude as the differences in Coulomb stresses and because one or more events in the sequence may have greater stress drops than the others.In all cases, the state of the system at the end of a sequence is different from the initial one, entailing that the durations of the interseismic intervals between consecutive sequences and between events in a sequence are different.
4. The observation of a seismic sequence allows the retrieval of the state of stress at any time during the sequence.In particular, we can calculate the state of stress at the beginning and at the end of the sequence.This has been done as an application for the 2012 Emilia (Italy) seismic sequence that was made up of seven events with magnitudes between 5 and 6.In this case, the evolution of stress was conditioned by the first and fourth events, having greater seismic moments than the average.The model shows the complex interplay between fault dislocations that produced the observed order of events, resulting in a greater stress heterogeneity at the end of the sequence.It predicts that, in the absence of external perturbations, the next sequence will occur after an interseismic interval of a few centuries and will be completely different from the 2012 sequence.We neglected the effects of pore fluid diffusion and of dynamic triggering that may have had a role in this case (Convertito et al., 2013).These effects may change the details of the results, but not the general conclusions illustrated here.

Figure 1 .
Figure1.Sketch of the model with n coplanar faults.The x axis is the strike direction.Distances r ij between the ith and j th faults are computed from the fault centres.

Figure 2 .
Figure 2. Case n = 3: (a) The hexagon H with its subsets labeled by the associated permutation; (b) States of the system during a seismic sequence: each point belongs to the subset labeled by the associated permutation; (c) magnification of H showing the initial and final states.

Figure 2 .
Figure2.Evolution of a system made up of three faults (n = 3), represented in the state space: (a) the hexagon H , defined in Sect.7, with its subsets labelled by the associated permutations; (b) states of the system during a seismic sequence: P 0 is the initial state; P i (i = 1, 2, 3) is the state after the ith event of the sequence; (c) magnification of H showing the initial and final states P 0 and P 3 .

Figure 3 .
Figure 3. Evolution of a system made up of three faults (n = 3): components of the state vector x as functions of time during the seismic sequence shown in Fig. 2 (τ = 1 MPa, δt = 1 a).The steps labelled by i = 1, 2 and 3 correspond to the occurrence of the events in the sequence.

Figure 4 .
Figure 4. Geographic location of the 2012 Emilia seismic sequence (Italy).Stars indicate the epicentres; numbers indicate the order of fault activation.

Figure 5 .
Figure 5. Seismic moments m i of the events in the 2012 Emilia sequence.The upper scale indicates the fault number, the lower scale the order of activation.The two strings of numbers yield the permutation α * in Eq. (62).

Figure 6 .Figure 7 .
Figure 6.Geometry of the model for the 2012 Emilia seismic sequence.The rectangles are the projections of faults on a vertical plane.Faults are numbered from west to east.Stars indicate the hypocentres.
for fault numbers.