The Growth of Earthquake Clusters

Migration of hypocenters is a common attribute of induced injection seismicity and of earthquake swarms, which distinguishes them from aftershock sequences. Spreading of the triggering front is often examined by fitting the time dependence of hypocenter distances from the origin by the pore pressure diffusion model. The earthquake migration patterns however often exhibit not only spreading envelopes but also fast-growing streaks embedded in the overall migration trends. We review the observed migration patterns and show that in the case of earthquake-driven migration, where the new ruptures are triggered at the edge of previous ruptures, it is more suitable to examine the cluster growth as a function of the event index instead of time. We propose a model that relates the speed of seismicity spreading to the average rupture area and the effective magnitude of the hypocenter cluster. Application of the model to selected linearly growing clusters of the 2008 West Bohemia swarm gives an almost linear increase of the measured total rupture area with the event index, which fits the proposed model. This is confirmed by a self-similar scaling of the average rupture area with the effective magnitude for stress drops ranging from 0.1 to 1 MPa. The relatively small stress drop level indicates the presence of voids along the fault plane and a possible role of aseismic deformation.


INTRODUCTION
The spreading of earthquake hypocenters is generally considered to be characteristic of earthquake swarms, which distinguishes them from aftershock sequences. While swarms typically show a clear hypocenter migration starting from a single point of the fault plane, aftershocks usually occur immediately all over the fault plane and along the edges of mainshock rupture. First aftershocks are located close to the rupture plane of the mainshock, while later aftershocks sometimes seem to migrate away from the mainshock (Rydelek and Sacks, 2001) at velocities smaller than 1 km/h. A significant expansion of aftershock clouds is only expected in the presence of aseismic processes such as mainshock-induced afterslip (Peng and Zhao, 2010;Perfettini et al., 2019), while coseismic stress transfer of the aftershocks itself can only explain a minor spreading (Helmstetter et al., 2003). Hypocenter migration of earthquake swarms is generally attributed to fluid flow along the fault plane. Number of examples of swarms from various tectonic environments show this feature, e.g. in Iceland (Woods et al., 2019), Japan (Yoshida and Hasegawa, 2018a), Afar Rift (Wright et al., 2012), Yellowstone (Massin et al., 2013), or West Bohemia/Vogtland (Fischer et al., 2014). Different migration scenarios are observed from monotonous migration in the horizontal or vertical direction or migration with direction changes within the evolution of the seismic cluster. The migration is often attributed to dyke intrusion, hydrofracture propagation, slow slip events, or pore pressure diffusion. Those models will be discussed below. Monotonous hypocenter migration is also observed for microearthquakes accompanying hydraulic fracturing, which is commonly used in the industry to create a permeable channel in tight rocks in order to increase the production of oil and gas or facilitate water circulation in the geothermal heat exchanger. The hydraulic fracture is formed by massive water injection at pressures exceeding the minimum principal stress. The growth of the edge of the hydrofracture (socalled rupture triggering front) can be explained by the conservation law relating the injected fluid volume to the volume of the fracture (Fischer et al., 2008;Shapiro, 2015). Accordingly, depending on the hydrofracture geometry, either linear or parabolic growth is observed. Pronounced migration of epicenters of nonvolcanic tremors, the weak seismic signals associated with slow slip on fault planes, has been observed in subduction zones along distances above 100 km, e.g. (Houston et al., 2011). The observed hypocenter patterns are explained by slow slip initiated by pressurization by fluids interacting with small asperities, which is reproduced by model simulations (Luo and Ampuero, 2017). In many cases, the coordinate-time plots of the migrating hypocenters show two features: slow nonlinear growth of the overall envelope and rapid growth of linear streaks during bursts of activity as observed e.g. within the 2008 West Bohemia swarm ( Figure 1A). Here, a streak is defined as a narrow linear belt of points in the coordinate-time plot, which is usually embedded in the overall large-scale migration of the seismic cloud. We identified similar rapidly propagating streaks also in other swarms from West Bohemia and other areas as Yellowstone (Shelly et al., 2016) or Corinth Rift (Kapetanidis et al., 2015). Interestingly, similar streaks also occur during the nonvolcanic tremors showing both forward and reverse motion of epicenters at velocities from 7 to 100 km/day (Ghosh et al., 2010;Houston et al., 2011) (see Figure 2 in Ghosh et al., 2010). In contrast to the natural earthquakes, the seismicity evolution with time is quite smooth and in most cases no streaks of hypocenters are apparent for the injection-induced seismicity ( Figure 1G). The only observation of a feature similar to the rapid streaks is that of Shapiro and Dinske, 2009, who attributed the linear fronts in the distance-time plots to hydraulic fracture opening. In summary, seismic swarms and nonvolcanic tremors show migration patterns on two scales: 1) an (asymmetric) overall growth over days to weeks 2) rapid migration on small spatial and temporal scales, partially with reversed trends. While the overall growth is to be explained either by pore pressure diffusion or hydraulic fracture growth (possibly also by slow slip, but not analyzed in detail yet), the embedded rapid migration have been studied only rarely in detail. Yoshida and Hasegawa (2018b) analyzed the triggering of the Sendai-Okura earthquake swarm by the 2011 Tohoku-Oki earthquake. And recently, Dublanchet and Barros (2020) proposed a model that reproduces the dual migration speeds observed in real swarms by coupling rateand-state friction, non-linear diffusivity and elasticity along the fault. Regarding the control of seismicity migration, two end-member cases can be distinguished: externally (aseismic) and internally (seismic) controlled propagation. In the first case an external physical process drives and controls the migration, where earthquakes are only passive markers of this usually timedependent process. This could be pore-pressure diffusion (Shapiro et al., 1997), hydraulic fracture growth (Dahm et al., 2010), or slow slip (Schwartz and Rokosky, 2007). In the second case, the migration is self-controlled by the earthquake ruptures themselves in the way that the ruptures allow or facilitate nucleation of adjacent new ruptures. Such an earthquakeinduced effect can be related to coseismically induced elastic stress changes as well as postseismic effects such as creation of pore-space enabling subsequent fluid flow. Both increase the stress state at the rupture borders facilitating subsequent adjacent ruptures. In such a case, the internal processes and conditions control the "time" and no general time dependence of the seismicity migration is necessarily expected. In this paper, we review the existing approaches to understanding the earthquake migration. We show examples of migration of earthquake swarms, where both the slow growth of the triggering front and the fast propagation of streaks occur, in particular when displayed using the coordinate -event order plot. Provided that every rupture opens the way for nucleation of the further rupture, we develop a theoretical model of the seismicity spreading related to the source parameters of the earthquakes suitable for the physical interpretation of the observed migration patterns.

MODELS AND INTERPRETATION OF EARTHQUAKE MIGRATION
In the following, we review the existing approaches for interpreting and employing the spreading of earthquake hypocenters.

Pore Pressure Diffusion
The observed hypocenter migration shows, in many cases, a nonlinear growth of envelope in the distance-time plot, which is commonly interpreted by pore pressure diffusion in a homogeneous isotropic medium. The model was originally proposed to interpret the spreading of hypocenters of injection-induced seismicity (Shapiro et al., 1997). It is based on the assumption that earthquakes map the advance of the pore pressure front. It is also assumed that the diffusivity of the rock is constant and that the pore pressure field is decoupled from seismic rupturing and related stress changes. This is in contrast to, for example, the model of Yamashita (1999) where the effect of dynamic pore creation is considered. Ignoring such coupling effects, the evolution of the pore-pressure field can be described by the diffusion equation whose solution for a point pressure source gives the distance r of the propagating pore pressure front as r 4πDt √ (Shapiro et al., 1997), where t is the time since the first contact of the pore pressure source with the host rock and D is the hydraulic diffusivity. Observations of injection-induced seismicity also show that after stopping the injection the induced earthquakes close to the injection point gradually cease to occur. A triggering backfront is being formed that becomes more distant from the injection point; its distance is also proportional to Dt √ (Parotidis et al., 2004). Rock properties are, however, anisotropic because of strong heterogeneities like faults acting as permeable channels on the one hand, or as impermeable barriers on the other hand. Furthermore, the stress field is heterogeneous due to the stress gradients originating from density contrasts between fluid and rock or Soultz-sous-Forets injection in 2000. Note the episodic occurrence of activity in the coordinate-time plots compared to the overall continuous spreading of activity in the coordinate-event-index plots. Note also that fast propagating streaks of hypocenters are observed in swarms but are missing in the injection induced seismicity. The color of symbols from dark blue to dark red is proportional to the event magnitude. In (A) two spreading episodes are highlighted in black in both coordinate-time and coordinate-event-index plots.
heterogeneous tectonic stress. The combination of both types of heterogeneities necessarily results in the preferred orientation of the fluid flow. Consequently, the seismicity front is expected to depend on the direction. The data of West-Bohemia/Vogtland swarms show a clear directional dependence of earthquake migration. In the case of the 2008 swarm the migration in dip-direction could be very well fitted by the pore pressure front curves. However, it requires different hydraulic diffusivities for up-dip and down-dip directions. The estimated down-dip diffusivity (D 0.1 m 2 s −1 ) is only one-third of the updip value (D 0.3 m 2 s −1 ) (Hainzl et al., 2012). Furthermore, in the down-dip direction, the migration stops completely after approximately nine days, whereas the migration continues in the up-dip direction. These discrepancies indicate that the diffusion model is often not able to explain earthquake migration in more detail and reveals the need for an alternative model in the case of tectonic earthquakes.

Hydraulic Fracture Growth
The pore pressure diffusion model is only appropriate in the limit of poroelasticity, which means, if non-linear effects such as the creation of pore-volume can be neglected and hydraulic properties remain unchanged. If the involved fluid pressure exceeds the minimum confining pressure, hydrofracturing of intact rock takes place. This leads locally to open channels for fluid transport and thus strongly changes the local permeability and hydraulic properties. As a result, instead of pore pressure diffusion in isotropic intact medium, the high-pressure fluid injection results in a subsequent opening of the fault zone and related fluid flows that are driven, in addition to the injection itself, by stress gradients. In a homogeneous rock a hydraulic fracture should open against the direction of the minimum stress component σ 3 . The model of hydrofracture growth of Fischer et al., (2009) and Dahm et al., (2010) can be then used to fit the migration pattern of the swarm activity. The resulting relation between the distance r of the triggering front and the time t since the start of the injection depends on the fracture geometry. Using the mass conservation of the injected fluid and assuming a constant injection rate and thickness of the hydraulic fracture without fluid penetration in the rock matrix, a linear relation of r and t is predicted for fracture propagation in a single dimension (1D). Similarly, for unlimited two-dimensional fracture propagation, a square-root growth of r with t is predicted. If asymmetric growth similar to that shown in Figures 1A, 2 is observed, the relation between the length of the shorter and longer wing can be used to estimate the driving stress gradient responsible for the asymmetry as presented in our previous studies (Fischer et al., 2009;Hainzl et al., 2012;Mls and Fischer, 2017).

Slow Slip Driven Seismicity
In contrast to the two other mechanisms, slow slip describes transient aseismic shear dislocation (slip) on faults. By means of new geodetic instrumentation, slow slip events have been identified to occur frequently on faults (Schwartz and Rokosky, 2007). Particularly, it is expected to occur in fault zones, where frictional behavior is partly characterized by velocity strengthening. However, it can also occur in velocity weakening zones, if the nucleation dimension exceeds the asperity size. The nucleation size is itself inversely proportional to the effective normal stress, thus the evaluated pore pressure favors aseismic vs. coseismic slip (Dieterich, 2007). While aseismic slip events can nucleate spontaneously, they can also be triggered by fluid injection (Guglielmi et al., 2015) or coseismic slip. The latter is known as afterslip which is believed to be one of the main drivers for aftershocks (Perfettini et al., 2019). In particular, a log(t) expansion of the aftershock zone has been forecasted for faults governed by rate-state-dependent friction, where the rupture of the mainshock asperity triggers postseismic creep in the surrounding velocity-strengthening region (creep zone). Simulations show that this creep leads to a log(t) migration front of aftershocks which are triggered in smaller embedded asperities (Kato, 2007). Such a r ∝ log(t) expansion of the aftershock zone has been observed for the early aftershocks of the 2014 M6 Parkfield event (Peng and Zhao, 2010) and for the aftershocks of the 2014 mainshock in Western Bohemia (Hainzl et al., 2016). Fischer and Hainzl (2017) proposed to employ the spreading of hypocenters of earthquake clusters for interpreting the stress release of the involved seismicity. The basic idea is that if a large earthquake rupture is replaced by numerous small ruptures, which is the mainshock on the one hand and the earthquake swarm on the other hand, similar strain should be released in both cases. Then the total seismic moment of the cluster should equal the seismic moment of a single large event and the familiar formula for seismic moment of a circular rupture M 0 7/16ΔσR 3 could be used to relate the total seismic moment ΣM 0 , radius of the cluster R and stress drop Δσ. To distinguish it from the static stress drop, we name it the effective stress drop. The validity of this approach however depends on the fact that all stress has been released during the sequence, which is pronounced in the power of M 0 (R) scaling. For this purpose, we compare the spatiotemporal evolution of the seismic moment release by projecting the hypocenters to a best-fitted plane and measuring the radius of the convex envelope of the cluster. By analyzing the uncertainties of the estimated effective stress drops it was found that the effective stress drop is only comparable to earthquake stress drops in specific cases. In particular, the effective stress-drop values underestimate the earthquake stress drops in the presence of aseismic deformation significantly. The values are only scale-independent if pre-stress and post-stress conditions are uniform in space. The analysis of data from injection-induced seismicity, natural earthquake swarms, and aftershock sequences showed that in most cases the effective stress-drop estimate is rather stable during the cluster evolution. The effective stress drop concept, however, explores only the overall growth of the earthquake cluster. It does not explore the preferred direction of the triggering front advance. We hypothesize that the advance of the overall cluster envelope and, in particular, of the embedded rapid migration patterns are closely related to the triggering forces and can be possibly used for discrimination of the underlying mechanism. In particular, we aim at modifying the effective stress drop concept in order to reflect the specific character of rapidly propagating streaks to learn more about their origin.

IDENTIFICATION OF MIGRATION PATTERNS
Typically, migration patterns are analyzed in the distance-time domain, where the form of observed migration patterns also depends on definition of distance. Traditionally, the dependence of unsigned distance on time, r(t), is explored, which masks however the possible asymmetry of the cluster growth. Such type of plot is suitable for determining the hydraulic diffusivity, which is independent of direction. To explore the asymmetry of the fracture propagation a coordinate-time plot x(t), where the coordinates x are taken including their sign, is more suitable (left panels of Figure 1). Concerning the independent variable, it is usually the time, which is used. This however postulates that it is the time, which controls the seismogenic process. Such plots are suitable to identify and characterize the external driving force, such as e.g. to estimate the hydraulic diffusivity D in the case of seismicity driven by pore-pressure diffusion.
An additional analysis tool is to use the event order as the independent variable, which is also termed natural time (Rundle et al., 2018). In such a case, a coordinate-event-index plot x(N) is produced. This postulates that the seismogenic process itself controls the seismicity. In other words, every rupture opens the way for nucleation of the further rupture, which is also the basis of the Yamashita (1999) model where the effect of dynamic pore creation is considered. Exploring the seismicity spreading dependent on the event order ignores the influence of time on the seismicity and allows for investigating the dynamics of the cluster formation. As a result, it is the seismogenic process that measures the new natural time; the pendulum stops in the period of quiescence and restarts with the new onset of seismic activity. Thus, the combined analysis of both plots, the coordinate-time and the coordinate-event-index plot, offers the possibility to understand the processes in more detail.
The form of migration diagrams using the coordinate-time plot and coordinate-event-index plot is compared in the left and right panels of Figure 1 for selected swarm activities. In the coordinate-event-index plots the periods of quiescence disappear and the original temporal clusters apparent in the coordinate-time plot merge into a continuous seismicity sequence migrating in a single dominating direction. The triggering front shows a linear envelope whose slope (migration velocity measured in meters per event) stays almost constant for longer periods. Note the continuous seismicity growth of the 2008 West-Bohemia swarm in the event-index plot ( Figure 1B) compared to the heavily interrupted spreading in the time plot ( Figure 1A) of the same activity. The change of growth character is highlighted in two rupture front clusters, which are continuous in the coordinate-event-index plot, but strongly discontinuous in the coordinate-time plot. Similar characteristics are illustrated in Figures 1E,F for the 2014 swarm in the Long Valley Caldera (Shelly et al., 2016), which showed, compared to the 2008 West-Bohemia swarm, almost simple unilateral spreading. The plots in Figure 1 also illustrate the two aforementioned phenomena: the overall growth of the seismic clouds manifested in the triggering front envelope and embedded rapid streaks that propagate both in the same and opposite directions as the triggering envelope.

THEORETICAL MODEL
One end-member model considers earthquakes only as a passive result of an underlying driving force, e.g. pore-pressure diffusion or slow slip which is not influenced by the triggered earthquakes. In this case, the analysis of the distance-time or coordinate-time plots is sufficient to characterize the migration pattern and the underlying process. However, in the other end-member model, earthquake migration is only possible due to the occurrence of earthquakes, e.g. by the creation of pore-space during ruptures. The presence of clear triggering fronts and embedded linear streaks in the coordinate-event-index plots, while no clear migration pattern can be observed in the coordinate-time plot, points to the latter case. In particular, in such a case, the speed of growth should increase with the size of the event ruptures. Developing a theoretical model for the earthquake generation in the case of a self-driven seismicity migration would help understand the generation mechanism of earthquake swarms. In accordance with the Yamashita (1999) model we suppose that every rupture facilitates the nucleation of adjacent ruptures. This can be simplified by assuming that subsequent events are without overlap initiated at the outer border of the previous rupture. In two dimensions, two geometrical concepts are considered, the channel (C) model describing a unilateral growth along a channel of width W and the sector (S) model describing a two-dimensional sectorial growth with angle θ (Figure 2).
The speed of the cluster growth is then related to the average position x of the rupture front, which is after the occurrence of N events (as observed in Figure 1) simply given by where both constants ν and D are proportional to the average rupture area 〈A〉 of the earthquakes in the cluster. In particular, model C predicts a linear growth of the rupture front x with increasing event index N according to a velocity ] 〈A〉/W, which is consistent to a constant average advance per observed event. On the other hand, model S predicts a square-root growth of the rupture front x with event index N, with a local velocity 0.5 D/N √ 0.5 360/πθ〈A〉/N per event. In both cases, the total rupture area grows linearly with event index as Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 638336 5

Fischer and Hainzl
The Growth of Earthquake Clusters A 〈A〉 · N. (2) In order to express 〈A〉 we suppose that the frequencymagnitude distribution of the events follows a doubly truncated Gutenberg-Richter distribution with exponent b and truncations at a minimum magnitude M 1 and maximum magnitude M 2 . The corresponding probability density function for the event magnitudes is in this case given by Marzocchi and Sandri (2003).
p(m) ln(10)b 10 − bm 10 −bM1 − 10 −bM2 . ( We also use the relation between the area A of a circular rupture and seismic moment M 0 f ΔσA 3/2 (see, e.g., review of Madariaga and Ruiz (2016) where f is a geometric factor of the source model, Δσ is the static stress drop, c and d are 9.1 and 1.5, respectively, for the moment magnitude m W (Hanks and Kanamori, 1979).
Then, the average rupture area 〈A〉 can be determined (see In the above equations, M 1 is the true physical minimum magnitude of the earthquakes, while we usually have a larger observational cutoff magnitude M c (completeness magnitude). Assuming that magnitude correlations do not exist (magnitudes are random within the sequence), the observed m ≥ M c events are mapping the true advance, but missing on average 10 b(Mc− M1) events between two successive m ≥ M c events. Consequently, the event index is Then Eq. 1 describes the average position of the seismicity as a function of observed events N obs , if the parameters are used instead of ] and D, respectively.

APPLICATION TO THE SEISMICITY DATA
Equation (6) predicts the average rupture area 〈A〉 from the seismicity parameters (b, M 1 , M 2 ). These parameters can be related to the average seismic moment per event 〈M 0 〉 and its corresponding effective magnitude M eff using Eqs. 4, 5. The empirical value of 〈M 0 〉 can easily be calculated for an observed sequence, while the theoretical value is given by Zakharova et al., (2013).
Based on 〈M 0 〉, the effective magnitude of the average seismic moment can be determined by As follows from Eq. 6, the average rupture 〈A〉 area depends on the stress drop and the magnitude distribution of the analyzed sequence reflected in its effective magnitude. This dependence is illustrated in Figure 3 showing the average rupture area as a function of the maximum and effective magnitude for a set of b-values ranging from 0.8 to 1.2 and stress drops between 0.1 and 100 MPa.
As predicted by Eq. 6, except the case b 2/3d, the average rupture area 〈A〉 increases non-linearly with the maximum magnitude M 2 . The growth rate decreases with increasing b-value ( Figure 3A), which is related to the decreasing number of large events. Similarly, larger 〈A〉-values are observed for smaller stress drops ( Figure 3B). Equation 6 further predicts that the slope increases with the moment-magnitude scaling factor d and that the curves shift up with the constant c. Figure 3 predict the average rupture area 〈A〉 growth in the case that the seismic activity is driven internally by the earthquake ruptures themselves. Then, 〈A〉 equals the average advance of the cluster growth per event (Eq. 1) directly. For real data, comparing the observed average rupture area 〈A〉 with the theoretical rupture area for the same effective magnitude would allow us to compare the scaling of the growth velocity with magnitude and verify the validity of our model for real data. Additionally, an estimate of the average stress drop of the cluster could be obtained.

The curves in
To test the validity of the presented models on data we measure the growth of selected clusters of the 2008 swarm in West Bohemia, whose migration patterns are presented in Figure 1 and compare it with theoretical 〈A〉 obtained using Eq. 6. The strike and dip coordinates are obtained by projecting the hypocenters to the best-fitting plane, which is justified by the planar character of the activated area ( Figure 4C). As follows from the diagrams, the 2008 swarm experienced prevailing upward and strike-ward (south-north) migration interrupted by several short reverse trends. The same clusters as indicated in Figures  1A,B were chosen for the analysis; both occurred at the migration front in the middle of the swarm, see also two projections of the x(N) diagram in Figure 4 and panels a and e of Figure 5. The earlier cluster resembles the sector type and the later one the channel type. The events close to the rupture front (the envelope of hypocenters) were selected in the strike projection; outliers in the other projection were removed using the criterion that at least 10% neighboring hypocenters have to occur within the distance equal to the length of the rupture front.
The total rupture area was measured by the convex envelope algorithm (convhull) implemented in Matlab, which determines the smallest total area A that contains the hypocenter projections within a convex hull. To track the cluster growth, we determine the ruptured area for an increasing number of events; 10 equidistant steps were used. The panels b and f of Figure 5 show the growing convex hull of the clusters indicating a mixed channel and sector growth. The dependence of the measured total rupture area on the number of included events is shown in panels c and g of Figure 5 indicating an approximately linear increase of the area with event index, which is in agreement with the model prediction in Eq. 2. The average rupture area per event 〈A〉 is 1750 m 2 for the earlier cluster and 2,870 m 2 for the later cluster. Considering the approximate width W 740 ± 100 m for the earlier and 630 ± 100 m for the later cluster, this corresponds to the speed of 2.4 and 4.6 m per event for the earlier and later clusters, respectively. The difference of 〈A〉 and of the growth speeds is clearly related to the different 〈M obs eff 〉, which amount to 1.1 for the earlier and 1.4 for the later cluster. Using Eq. 6 to estimate the theoretical ratio of 〈A〉 of events with 0.3 magnitude difference gives 1.64, which equals the ratio of measured 〈A〉 amounting to 1.66.
Scaling between the average rupture area and the measured average effective magnitude of the two clusters is analyzed in panels d and h of Figure 5. The measured 〈A〉 is shown (crosses) in dependence on the 〈M obs eff 〉 of each subset of events defined by the expanding window mentioned above. Curves of theoretical 〈A〉 are shown for comparison with the measured 〈A〉 (Eq. 6). We used b 0.9 as determined by Hiemer et al., (2012) for the 2008 swarm, M 1 0 according to the empirical cutoff magnitude minimum, and varied M 2 so that the 〈M eff 〉 (Eq. 10) covers the measured 〈M obs eff 〉 range. It is worth to note that the choice of both M 1 and M 2 influences the theoretical curves. While the range of M 2 -values (for given M 1 ) only determine the 〈M eff 〉-interval of each theoretical curve but not its shape, M 1 slightly affects the vertical shift and curvature of the lines (see Supplementary Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 638336 Figure S1). Accordingly, M 1 should be chosen with care to correspond to the real minimum magnitude. For the moment-magnitude relation (Eq. 4) we used constants c 10.09, d 1.1 obtained for the West Bohemian swarms by Jakoubkova et al., (2017) and f (7/16)π 3/2 for the circular rupture. It turns out that for the earlier cluster ( Figure 5D), 〈A〉 scales with 〈M obs eff 〉 with rather constant stress drop Δσ in the range 0.2-0.8 MPa. Similar behavior is found for the later cluster with a slightly smaller stress drop from 0.1 to 0.5 MPa ( Figure 5H).

DISCUSSION AND CONCLUSIONS
Our analysis shows that the usage of coordinate-event-index plots can yield valuable insights into the earthquake triggering process. However, it is important to note that this analysis should always be performed in conjunction with an analysis of traditional coordinate-time plots. A linear growth with the event-index does not necessarily mean that the migration process is internally driven. In the case that earthquakes are triggered in a homogeneously pre-stressed brittle fault segment by external forcing such as pore-pressure diffusion, the earthquake number would be expected to increase with the pressurized area. Thus, even if the events do not influence each other, this would lead to a linear or square-root growth with earthquake index, similarly to the case of self-driven migration. However, in the case of diffusion-driven seismicity, the space-time plots would also show a clear and continuous migration pattern. In contrast, the spacetime plots are expected to show a less clear picture with discontinuous spreading in the case of self-driven migration, as observed for our analyzed cases (see Figure 1). Only the combined analysis of space-time and space-index plots together can resolve this issue.
The graphs in Figure 5 show that the total ruptured area grows almost linearly with the event index, as expected for the suggested model supposing that every rupture facilitates FIGURE 4 | The 2008 earthquake swarm in West-Bohemia showed prevailing unidirectional growth of the seismicity front in the strike-ward (A) and upward (B) directions. Two linearly growing clusters were selected in the strike direction for the analysis. The triggering front of the earlier cluster (black) growed slower than that of the later cluster (blue). To measure the distance between the points needed for removing outliers, the event index was transformed into a distance in km so that the aspect ratio of the diagram was 1:1. In (C), two vertical sections show the planar character of the swarm activity; the symbol color is proportional to the event origin time.
Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 638336 the nucleation of a new adjacent rupture making the seismogenic process self-organized. The local shape of the x(N) curve is however influenced by the actual magnitude distribution of the running activity; increased magnitudes will also increase the average rupture area and raise the slope of the curve. In general, the linearity of the migration in the  Figure 4: The activated area is marked in colors, which are proportional to the event index (A,E); its convex hull was determined (B,F) for extending equidistant windows that included from 10 to 100% of events; the scaling of the total rupture area with the number of events is almost linear (C,D). The dependence of the average rupture area per event 〈A〉 on the average observed effective magnitude 〈M eff 〉 is examined (D,H).

Fischer and Hainzl
The Growth of Earthquake Clusters coordinate-event-index plots might be distorted by the occurrence of single earthquakes with large magnitudes departing from the trend of the Gutenberg-Richter distribution. However, this is usually only observed in mainshock-aftershock sequences, but not in swarm activity as considered here. Another worth-noting point is the estimated stress drop of less than 1 MPa, which appears rather small compared to the estimated static stress drops of individual events of the West-Bohemia swarms in the order of 10 MPa (Michalek and Fischer, 2013). It is of interest to note that using the concept of effective stress drop (Fischer and Hainzl (2017)), we arrived at a similar value of Δσ eff 0.23 MPa for the whole 2008 swarm. Both approaches are of similar nature and exploit the relation between the rupture area and the seismic moment. The difference to the static stress drop estimates of more than one order of magnitude could be partially related to the usage of the Madariaga's model (Madariaga (1976)) for the static stress drop values, which gives five times higher estimates than the alternative Brune's model (Brune (1970)). Another reason might be related to our assumption that subsequent events are without overlap initiated at the outer border of the previous rupture. If the deformation process involves triggered aseismic afterslip, void areas will remain among adjacent ruptures and the resulting stress drop estimate would decrease accordingly. A direct proof is difficult because of the large uncertainties involved in the estimation of hypocenters, stress drops, and rupture geometries. For illustration, Supplementary Figure  S2 shows the estimated rupture extensions for the two analyzed clusters assuming circular ruptures centered at the estimated and projected hypocenters, where the rupture radius depends on the stress drop. Using a constant stress drop of 1 MPa for all events leads to a significant overlap of the activity, while a stress drop of 10 MPa leads to adjacent ruptures with some void areas.
Migration characteristics are used to typify seismic activity. In contrast to the random occurrence of aftershock hypocenters along the mainshock fault plane, earthquake swarms and injection-induced seismicity outstand by monotonous migration resulting in the spreading of the hypocenter cloud. This is usually attributed to fluid flow along the fault plane, which could be modeled by porepressure diffusion or by a propagating hydraulic fracture. In this endmember case, the aseismic load leads to an increase of the effective stress in an expanding volume which triggers the earthquake migration as function of time, independently of precursory seismicity. Earthquakes are assumed to be triggered in this case in critically loaded areas, where stresses or pore pressure is elevated due to the aseismic driving process. As a result, void areas may occur in the diffusion model, which are bridged by the aseismic pore pressure diffusion or hydrofracture growth. Accordingly, the event coordinates would grow non-linear, respectively jump with increasing event index, which would also apply to the growth of average rupture area with the event index. This would result in low values of the estimated stress drop due to the large activated area and small strain released in a brittle way, as discussed in detail in our previous study (Fischer and Hainzl, 2017). However, the suggested analysis of the coordinate-event-index plots would indicate that aseismic deformation is present in this case due to the non-linear, jumptype shape of the rupture front.
In contrast to the models where the seismicity is externally controlled by the underlying time-dependent aseismic process, the seismicity is self-controlled by the rupture occurrences in the other end-member case of self-driven migration. In that case, fluid flow is only a secondary process where seismic slip allows the fluid to migrate into the reactivated fracture through dilatancy, and once the pressure gets high enough on the rupture edge, the subsequent rupture occurs. Consequently, it is the rupture occurrence that measures the time and, hence, it is more suitable to use the event index as an independent variable when exploring the seismicity migration. A continuous, and in most cases linear advance of rupture front is predicted in the coordinate-event-index plot in the case of self-driven earthquake migration.
For several examples of seismic swarms, we show how the form of migration plots changes when the event index is used instead of time. Discontinuous activities convert to continuous ones and nonlinear growth with time becomes linear when the event index is displayed on the horizontal axis. We interpret this behavior by a model supposing that every rupture opens the way for nucleation of the further one. Then the speed of growth is proportional to the average rupture area of the earthquakes that is related to the magnitude-frequency distribution of the events and the stress drop.
Our test using relocations of the 2008 earthquake swarm in West-Bohemia confirms the validity of the presented model. The two selected clusters, which represent the linearly spreading triggering front in the coordinate-event-index plot, show an almost linear growth of the total rupture area with the event index. Based on the envelopes of the observed hypocenters, the estimated speed of growth is 2.4 and 4.6 m per event and the average rupture area of the cluster events is 1750 m 2 and 2,870 m 2 for the two clusters, respectively. The observed difference in the spatial spreading speed is explained by differences in the frequency-magnitude distribution, which is expressed by the different effective magnitudes of 1.1 and 1.4 for the two clusters, respectively. The average speed measured in time domain is 340 m/day for the first cluster and 720 m/day for the second cluster. The dependence of average rupture areas on the effective magnitude points to self-similar scaling with stress drops ranging between 0.1 and 1 MPa. Although our analysis indicates a self-driven cluster growth, the low stress drops might indicate that parts of the rupture area deformed aseismically either by porepressure induced or seismically triggered creep and did not contribute to the seismic moment release. In future studies, we will address the role of fluids and earthquake ruptures for specific migration patterns such as triggering fronts and embedded linear streaks in the case of various seismic activities in more detail.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://zenodo.org/ record/4294887.

AUTHOR CONTRIBUTIONS
TF and SH developed the idea of linear streaks and wrote the manuscript. TF contributed to the theory developed by SH and carried out the data analysis.

FUNDING
The work of TF was supported by the Czech Funding Agency under the grant 20-26018S.