Dynamic nucleation as a cascade-up of earthquakes depending on rupture propagation velocity

Earthquake dynamic rupture requires a nucleation process to provide sufficient energy to overcome the fracture energy. Large earthquakes may occur via a cascading rupture process, which includes many triggering processes that cascade from small to large sections of the fault system. During such a process, the nucleation of a large section of the fault plane may occur dynamically via the propagating rupture from a small section of the fault plane. A quasi-static view of seismic nucleation has been widely discussed in earthquake seismology; however, the dynamic nucleation process remains poorly known. Here, we investigate one aspect of the dynamic nucleation process by focusing on the rupture propagation velocity during the nucleation process. We simplify this process as self-similar crack propagation at a constant rupture velocity in a finite nucleation zone within a target region that possesses a uniform fracture energy. We numerically solve this elastodynamic problem in two dimensions for both the anti-plane and in-plane cases using the Boundary Integral Equation method. As the rupture velocity increases, the critical ratio of the fracture energy step to continue the rupture increases and the critical size of the dynamic nucleation zone decreases. The rapid increase in the ratio of the fracture energy step toward infinity could explain why earthquakes never propagate at slow rupture velocities. However, the effect on the size of the nucleation zone is rather limited, with the size of the dynamic nucleation zone decreasing to ~ 70% of the static nucleation zone size. However, such a small difference would result in a significant overall difference if such a dynamic nucleation process repeatedly occurred in the cascading rupture process of a large earthquake, which would be a difficult situation for earthquake early warning.


Introduction
Is there any difference in rupture growth process between small and large earthquake?Is the final size of an earthquake predictable from the very beginning of seismic signals?These are the questions that have been discussed intensively since around 1990, in part owing to their significance in improving earthquake early warning capabilities (Abercrombie and Mori 1994;Ellsworth and Beroza 1995;Iio 1995;Mori and Kanamori 1996;Olson and Allen 2005;Colombelli et al. 2014).Recent observations have indicated that the initial part of seismic signal is statistically universal among different ranges of event magnitudes (Meier et al. 2016;Noda and Ellsworth 2016).Furthermore, many pairs of different-sized earthquakes have been identified as sharing almost identical initial seismic signals, thereby suggesting that these events have undergone almost the same rupture growth process (Okuda and Ide 2018;Ide 2019).
An earthquake is the propagation of fast slip on fault planes that are loaded with high stress.For this fast slip to propagate spontaneously, fracture mechanics dictate that there must be an finite area on the fault plane, where the stress is released in advance by a particular process.This process is called seismic nucleation, and the area in which this process occurs is called the nucleation zone.Therefore, the above questions can be revised to how the nucleation of a large earthquake may differ from that of a small earthquake.However, the understanding of the physical process of seismic nucleation remains limited mainly in some quasi-static processes, such as stress corrosion (Das and Scholz 1981) and slow preslip (Dieterich 1992;Matsu'ura et al. 1992;Uenishi and Rice 2003;Rubin and Ampuero 2005), have been suggested.Nevertheless, it should be emphasized that slow deformation is not always a requirement of seismic nucleation process.It may occur via faster seismic processes, as demonstrated through the heterogeneous distribution of some physical parameters (Ampuero et al. 2006;Ripperger et al. 2007;Galis et al. 2015Galis et al. , 2017)).
The size (radius) of the nucleation zone is proportional to the fracture energy of the surface to be ruptured (e.g.Andrews 1976;Rubin and Ampuero 2005).Although the fracture energy of large earthquakes is larger than that of small earthquakes (Tinti et al. 2005;Abercrombie and Rice 2005;Mai et al. 2006;Viesca and Garagash 2015), fracture energy is usually only measured after the slip has progressed sufficiently, such that this value is uncertain at the beginning of a given earthquake.Rather, earthquakes can grow by sequential fracturing of areas with small, medium, and large fracture energy, which is known as a cascading process and has attracted much attention in earthquake seismology research (Ellsworth and Beroza 1995;McLaskey and Lockner 2014;Ellsworth and Bulut 2018;McLaskey 2019).In such cascading process, any preceding processes before the rupture of a target region should be considered as its own nucleation process.
Cascading rupture growth processes are expected owing to various fractal heterogeneities in the fault, rock properties, and stress and earthquake statistics like the Gutenberg-Righter law.The distribution of fault traces on map is fractal (e.g.Okubo and Aki 1987;Hirata 1989), and the topography of the fault surfaces is also consistently fractal from the sub-millimeter-to hundred-kilometer scale (e.g.Power and Tullis 1991;Candela et al. 2012;Candela and Brodsky 2016).The stress distribution may be also fractal (Day-Lewis et al. 2010), with such a distribution likely arising from the kinematic interaction of slips on a fault with a fractal geometry and/or material properties (Langenbruch and Shapiro 2014).Seismicity also has a fractal distribution (Okubo and Aki 1987;Hirata 1989;Mandal and Rastogi 2005), and the power-law magnitude-frequency distribution, which is known as the Gutenberg-Richter law, is associated with self-organized criticality.Although fault systems have some characteristic size, such as widths of the damage zone and seismogenic zone, which can be characterized by some exponential distributions of material properties (Mitchell and Faulkner 2009;Savage and Brodsky 2011), the fractal or powerlaw structure should exist at scales smaller than these characteristic sizes.
This idea of a cascading rupture process in fractal media has been incorporated into the numerical model of Ide and Aochi (2005), whereby circular patches with different fracture energies are distributed fractally along the fault plane.This model successfully explains the dynamic rupture processes of earthquakes as a scale-free process.It is also consistent with both initial seismogenic signal observations (Yamada and Ide 2008;Ide 2019) and the nature of the source time function catalog obtained for many earthquakes (Renou et al. 2022).The similarity of the initial stage to that for the entire rupture process (Uchide andIde 2007, 2010) also supports this cascading rupture process.
Each step in the cascading dynamic rupture process, whereby a small rupture induces the failure of a larger area, can be considered a building block in the overall rupture process.This small dynamic rupture process acts as a nucleation process for the larger area, which we term the target region.Therefore, it is a problem of nucleation process for the target region.In fact, such nucleation process has been commonly introduced in numerical simulations of earthquake dynamic rupture since the pioneering study of Day (1982).Since the process may sufficiently fast to radiate seismic waves, we call it dynamic nucleation process and distinguish it from quasi-static nucleation process.Although there are a few theoretical studies on dynamic nucleation (Kame and Yamashita 1997;Madariaga and Olsen 2000;Bizzarri 2010), the quantitative difference between dynamic and static nucleation remains largely unknown.Since static nucleation is based on the equilibrium equation of the continuum medium, such that there is no inertia term (i.e.considered an extremely slow nucleation process), an important characteristic that needs to be quantified is the dependence of dynamic nucleation behavior on the rupture propagation velocity.Here, we investigate dynamic nucleation conditions, with a focus on the rupture velocity inside the nucleation zone, to clarify the difference between dynamic and static nucleation processes.

Theoretical background
We first outline the theoretical framework for a simple dynamic nucleation process, which we treat as a crack problem, whereby a small crack propagates at a constant rupture velocity within a target region that possesses a constant fracture energy.Let us consider either an antiplane or in-plane crack on a planar fault in a two-dimensional infinite homogeneous Poissonian elastic medium with rigidity, µ , and shear-wave velocity, v s .Since the absolute value of the residual stress has no effect on these calculations, the shear stress is defined relative to the residual frictional stress, which is assumed to be constant and sufficiently high to prohibit slip in the opposite direction.We assume that the initial shear stress, σ 0 , is also homogeneously loaded on the fault plane, and introduce a linear slip-weakening friction law (Ida 1972) that is characterized by a peak stress, σ p , and slip-weakening distance, D C (Fig. 1a).The uniform stress condition is simple and sufficiently generic, as suggested by the scale invariance of the apparent stress over a broad range of earthquake sizes (M w − 4 to 8; e.g.Ide and Beroza 2001;Kanamori and Rivera 2004;Baltay et al. 2011).This uniform stress condition has been effectively utilized in previous numerical studies (e.g.Ide and Aochi 2005).
We generate spontaneous rupture propagation at a constant rupture velocity in the numerical simulations by assuming that the fracture energy within the nucleation zone increases linearly with the increasing spatial extent of the crack (e.g.Andrews 1976;Aochi and Ide 2004) until the rupture encounters the target region boundary (Fig. 1).The fracture energy (Fig. 1a), G C , can change by several orders of magnitude, as has been suggested by the large discrepancies in the G C estimates for large earthquakes (Beroza and Spudich 1988;Ide 2003;Tinti et al. 2005;Abercrombie and Rice 2005;Mai et al. 2006) and those from laboratory rock experiments (Ohnaka 2003;Scholz 2019).
Our assumption of a homogeneous stress condition allows us to define D C as a function of distance from the initial rupture point, r = |x| (Fig. 1b): where D ′ C is a constant showing the spatial gradient of D C .There is a sudden jump in D C from D ′ C R dyn to D ∞ C , which is constant in the target region, at the edge of the nucleation zone ( r = R dyn ).The effect of this jump on numerical result will be discussed later.The size of this discontinuity is represented by the fracture energy step ratio, γ , which is defined as We recall that the rupture velocity in the nucleation zone, v i , is constant.The rupture velocity is determined entirely by the energy budget beyond the rupture tip in the stationary condition.This equilibrant condition is given analytically [Eq.(6.9.148) in Broberg (1999)]: where R is the half-length of a self-similar crack; K = 1/π ; E(•) is a complete elliptic integral of the second kind; v s is S-wave velocity; and ν is the Poisson ratio.The left-and right-hand sides of Eq. ( 3) are the energy release rate and fracture energy at both crack tips, respectively.
We can define the non-dimensional spatial gradient of the fracture energy, ξ , in the nucleation zone as (3) Since ξ is a monotonic function of v i /v s , Eq. ( 4) can be rewritten as The energy balance equation [Eq.( 3)] is more complicated for an in-plane crack, as given by Eq. (6.9.90) in Broberg (1999).Nevertheless, we can also derive a similar function between ξ and v i /v s for 0 < v i < v R , where v R is the Rayleigh-wave velocity.ξ is directly related to a non-dimensional parameter, κ , that has been intro- duced by Madariaga and Olsen (2000) as ( where L is the characteristic length scale.With respect to Eq. ( 4), we can regard that L = R dyn and D C = R dyn • D ′ C , and therefore, ξ = K /κ .Rupture can only propagate when κ is larger than a critical value, κ c , which is esti- mated as either κ c = 1/π (anti-plane), κ c = 1/{π(1 − ν)} (in-plane), or κ c = 7π/24 (circular crack in three dimensions).
In the outer target region, which is where D ∞ C is assumed, there is a critical size for the static nucleation zone, R sta C , that can be determined via the Griffith's fracture criterion (e.g.Ida 1972) as where K is 1/π [the same as Eq. ( 3)] and 1/{π(1 − ν)} in the anti-plane and in-plane cases, respectively.While R sta C is the size of the nucleation zone for a static problem, where v i = 0 , we expect the critical size of the nuclea- tion zone to change when v i > 0 .We define this as the critical size of dynamic nucleation zone, R dyn C , which we will quantitatively estimate in the following sections.We note that R dyn C is related to the critical size of the fracture energy step ratio as from Eq. ( 2).Furthermore, we can reorganize Eqs.(2), (5), and (7) to obtain the following relationship:

Numerical method
We conducted dynamic rupture simulations along a planar fault in a two-dimensional medium using the Boundary Integral Equation Method (BIEM) (Tada and Madariaga 2000) for both the anti-plane and in-plane cases in a homogeneous, infinite elastic medium.This is a scale-free problem and the fault size can be that of real earthquakes of any size.The fault was discretized by 1024 or 2048 elements.We started with the coarser grid size (e.g.Fig. 1) and then the simulations from two different resolutions confirm that the results are identical and the simulation resolutions are sufficiently good.The following simulation results are shown from the finer resolution simulations.The D C -value distribution was determined via Eq.( 1).The time step, t , was set to �x/2v s for the anti- plane case and reduced to �x/4v s for in-plane cases to resolve P-wave.Rupture was initiated by a sudden reduction in yield strength from σ p to zero at t = 0 in a small region at the center of the model ( |x| < r ini ) . Figure 1c, d illustrates two slip evolutions when slightly different fracture energy values (0.6%) are employed in the target region for the anti-plane case.There was a self-similar increase in slip within the nucleation zone as the rupture propagated from the center of the model space at a constant velocity.The rupture did not abruptly stop when it reached the edge of the nucleation zone and entered the target region, but it also could not maintain a constant rupture velocity.The rupture terminated when there was a sufficiently large fracture energy to stop the propagation (Fig. 1c).Otherwise, the crack accelerated again and ruptured the entire model space (Fig. 1d); we defined this latter simulation result a successful dynamic nucleation process, or simply a single cascade-up process.
We sought the condition that controlled this bifurcation through a bisection method, whereby we changed the model parameter ( γ ) with respect to several given parameters ( D ′ C , D ∞ C , σ p , and σ 0 ).We note that R dyn was determined once γ was given.Once rupture has initiated and stabilized, it will then propagate at a constant velocity owing to the proportional increase in D C with dis- tance (Aochi and Ide 2004); this has been numerically confirmed in Fig. 2. We measured the rupture velocity based on the expansion of the crack tip, with the crack tip at a given time step corresponding to the leading edge, given by discrete units.The rupture velocity, v i , is normalized by shear-wave velocity, v s , in the upper panel.Small fluctuation in the velocity is due to quantization in space and time.The bold black and thick gray curves correspond to the rupture front (leading edge) and cohesive zone end (trailing edge), respectively.The theoretical R sta C is indicated by vertical dashed line.The artificial effect of sudden rupture onset disappears quickly by r/R sta C ∼ 0.2 , and the rupture velocity then remains quasi-constant until r/R sta C ∼ 0.7 .The rupture velocity drops suddenly at around r/R sta C ∼ 0.8 , because we introduce a fracture energy step, which controls the successful/unsuccessful rupture growth on the target region.This size is not always the same as R sta which is at the point of peak stress.We note that it is also possible to define the rupture velocity based on the trailing edge, which is where the stress has decreased to the residual stress level, namely, the end of the cohesive zone (e.g.Andrews 2004), with a high degree of similarity obtained using the leading and trailing edge approaches (Fig. 8).Furthermore, we note that the effects of artificial rupture initiation are negligible, and that the rupture velocity rapidly approaches a constant value.
In practice, the initial crack size, r ini , in our numerical simulations was assumed to be at least 15 times smaller than R dyn , empirically based on our preliminary study.This artificial influence of the initiation process is limited in the beginning of simulations as demonstrated in Fig. 2. The rupture velocity gets stable at a distance of 0.2 × R sta .Although we do not introduce any renormalization process (Aochi and Ide 2004), which showed a constant rupture velocity after twice the renormalizations, our resolutions are entirely sufficient except for the beginning of rupture process.We set R dyn to be smaller than R sta , as we expect energetic impact of the dynamic rupture.The model parameters are summarized in Table 1.The slip amount and D C values were normalized by 10 −3 x , where x is the element size.The assumed range of D ′ C values (0.001-0.1) is consistent with earthquake observations, where a D C value in the ~ 0.1 − 1.0 m range cor- responds to a crack with a half-width of ~ 10 km (e.g., Ide and Takeo 1997;Olsen et al. 1997;Mikumo et al. 2003;Kaneko et al. 2017;Aochi and Twardzik 2020).
For each simulation, the fracture energy step ratio, γ , was calculated via Eq.(2).We sought the critical value, γ c , which denotes the bisection between unsuccessful (Fig. 1c) and successful (Fig. 1d) triggering of the target region.We note that there were some limitations of the available parameter space for two reasons.One was the limitation of the numerical method, as the simulation had to thoroughly resolve the slip-weakening process occurring in the cohesive zone (Kame et al. 2003) during the entire simulation.We required v l i − v t i /v l i to be larger than 0.03 and 0.94 for the in-plane and anti-plane cases, respectively, where l and t refer to the leading and trailing edges, respectively.The other reason was to avoid any potential effects that depended on the details of the friction law.R sta C is the size of a singular crack at a limit of infinite σ p and zero D ∞ C , and a value between R l and R t , which are the sizes measured using the leading and trailing edges, respectively.The region between R l and R t is the cohesive zone.R l and R t depend on the details of the friction law, and R sta C might not be a good reference when R l and R t are very different and the cohesive zone is large.Therefore, we required R t /R l > 0.70 as a condition for simulation.See Appendix for a detailed discussion of this requirement.

Results
Here we aim to analyze the scaling property of the simulation results.Figure 3a shows the D ′ C , parameter space corresponding to various v i values for the anti-plane case.We can explore a wide range of rupture velocities due to the scaling relation in the nucleation zone.Figure 3b shows the simulated rupture velocities within R dyn as a function of ξ for the varied parameter values that were used in Fig. 3a.We observe that there appears to be a unique v i value for a given ξ value, with the numerically obtained relationships being very close to the theoretical curve given by Eq. ( 4).We observe a slight deviation in the simulations that employed a relatively low σ p and a high ξ .The low σ p simulations yield a rupture velocity that is faster than that determined via the analytic prediction, which assumes an infinitely large σ p .The high ξ simulations impose a large cohesive zone and induce a deviation from the analytic prediction, which has no cohesive zone.These deviations also affect other relationships and are discussed below.The parameter space of the simulations for the in-plane case is shown in Fig. 9.
We searched independently for γ c in two param- eter spaces, (σ 0 , D ′ C ) and (σ p , D ′ C ) , and kept the remain- ing parameters invariant from the reference parameters given in Table 1.The obtained γ c values with respect to different parameter combinations are shown for the antiplane case in Fig. 4.This includes the following (σ p , σ 0 ) combinations: (14 MPa, 4 MPa), (14 MPa, 3 MPa), (20 MPa, 4 MPa), and (20 MPa, 3 MPa).Fixing either σ p or σ 0 , we can induce slight variations in γ c .The compre- hensive result highlights that nucleation becomes easier when σ 0 is higher (Fig. 4a) or σ p is lower (Fig. 4b).This can be interpreted exclusively in terms of ξ .However, a comprehensive investigation of the ξ − γ c space (Fig. 4c) indicates that a global trend exists between ξ and γ c , and γ c is well-approximated as a function of ξ alone via Eqs.( 4) and ( 5).Therefore, we do not lose this generality by focusing our study on a limited parameter space (Fig. 3).
We then investigate the influence of v i on γ c .It is expected that γ c = 1 at the limit of a static nucleation Slip weakening rate Dc' (Fig. 3b) 0.001-0.1 process ( v i → 0 ), such that a faster rupture velocity may trigger a region with a larger fracture energy.Figure 5a highlights the relationship between γ c and v i from the previously simulated anti-plane cases (Fig. 3).We observe that γ c increases monotonically with v i .A crack propagating at a slower rupture velocity is, therefore, likely to be terminated by smaller perturbations in the fracture energy.We also calculate the same relationship for the in-plane case in Fig. 5b, with its respective parameter space summarized in Fig. 9.We observe the same tendency in this case and note that the rupture velocity is saturated until the Rayleigh wave velocity reaches v i ∼ 0.92v s .Regardless of the difference in σ p , γ c can be approximated as a unique function of v i for both the anti- plane and in-plane cases.We will further discuss the approximate theoretical background for this relationship in the Discussion section.We can evaluate the minimum rupture propagation velocity that is required to overcome a certain γ value from these results.For example, we consider the rupture velocity to overcome a gap of γ = 2 , which is the ratio used in the Ide and Aochi (2005) fractal patch model.In this case, v i must be larger than 0.4v s and 0.6v s .respectively.The ability to trigger a gap of γ = 4 , which means skipping one discrete scale in the Ide and Aochi (2005) model, is influenced by a more strict condition, and v i must be larger than ~ 0.8v s for both cases.

Comparison with a singular crack model
We have investigated the dependence of γ c on v i using a numerical approach.Here, we take an alternative approach and derive the relationship using the energy balance of an analytical solution.We begin by explaining the anti-plane problem in this context, and then extend the discussion to the in-plane problem.
We find that the function (1 − v i /v s ) −1 is a good approximation of γ c (v i ) , as shown in Fig. 5a, for the following reason.At the edge of the nucleation zone, where R = R dyn , the rupture velocity suddenly decreases to a small value ( < 0.1v s ), because it takes time to slip a longer D ∞ C (Fig. 2).Given that the rupture velocity suddenly decreases to zero, the stress intensity factor for the anti-plane problem changes from K III (v i ) to K III arrest , obeying the following relationship [Eq.( 6.12.39) in Broberg (1999)]: The corresponding energy release rate, G arrest , is where G ss (v i ) is the energy release rate of a self-sim- ilar anti-plane crack propagating at v i [i.e. the left- hand side of Eq. ( 3 (10) this approximate energy release rate must be close to The problem is qualitatively the same for the in-plane problem, but the ratio of the stress intensity factors (Eq.9) is expressed in a complicated form: where: and v R and v p denote the Rayleigh and compressional- wave velocities, respectively [Eq. (6.11.51) in Broberg (1999)].γ c , therefore, becomes [k II (v i )] −2 .This theoretical value is plotted in Fig. 5b and modestly agrees with our numerical results.
Although these analytical approximations explain the behavior of γ c well, it is difficult to explain the qualitative (11 Fig. 5 Plots of the γ c as a function of v i for the a anti-plane and b in-plane cases, respectively.Different σ p values are used for a fixed initial stress of σ 0 = 3 MPa.The curves in (a) and (b) are given by Eq. ( 11) and k II (v i ) −2 in Eq. ( 13), respectively.See Fig. 3 for the parameter space relation between the ratio between dynamic and static nucleation zone sizes, R dyn C /R sta C , and the nucleation rupture velocity.Since R dyn C /R sta C is represented by γ c and ξ using Eq. ( 8), we can predict this ratio using Eqs.( 5) and ( 11); the result is shown in Fig. 6.Although good fits to the curves in Figs.3b [ ξ (v i ) ] and 5a [ γ c (v i ) ] are obtained, as demonstrated for the anti-plane case, our numerical simulations do not match well with this theoretical prediction, which is formulated in the stationary condition for an infinitely strong medium without a cohesive zone.In particular, the theoretical curve predicts a very strong triggering effect (a small R dyn C /R sta C ) for a high rupture velocity ( v i > 0.8v s ), whereas the dynamic influence is still limited down to around R dyn C /R sta C = 0.6 .This is probably due to relatively large discrepancies in approximating both curves at a high rupture propagation velocity.At the limit of v i → v s , ξ (v i ) → 0 and γ c (v i ) → ∞ , which highlights that we need a better approximation for γ c (v i ) than Eq. ( 11).

Effect of smoothing
Our simplified fault model assumes an abrupt increase in fracture energy at the edge of the nucleation zone, as previously suggested in Ide and Aochi (2005) and other studies.However, the question arises of what happens if the heterogeneity varies smoothly.We, therefore, checked the validity of our discussion.Let us consider that the abrupt increase of D C is smoothed by tak- ing the simple central moving average of D C on a fault with L = 2 n + 1[n = 0, . . ., 6] nodes around the edge at |x| = R dyn .We note that the moving average of D C is identical to the original value of D C in the nucleation zone, since D C either increases linearly along the fault or is constant, such that v i is constant in the nucleation zone.Figure 7 shows the effect of different smoothing lengths on D C .The estimated R dyn C changes slightly according to the degree of smoothing (Fig. 7a), and the variation in γ c due to the difference in L is at most 3% at v i = 0.87v s (Fig. 7b).This test, therefore, confirms that the assumption of an abrupt change is sufficiently simple but does not alter the generality of our discussion.

Conclusions
The dynamic rupture process of an earthquake requires some nucleation process that occurs along part of the fault system.Furthermore, there is increasing evidence that an earthquake rupture can be approximated as a cascading process (e.g.Ellsworth and Beroza 1995;McLaskey and Lockner 2014;McLaskey 2019), whereby there are many stages of triggering from small to large faults.In this case, the source area of each stage requires some nucleation process.Such a nucleation process may not occur in a quasi-static way, but rather in a dynamic way at some rupture propagation velocity.While the quasi-static nucleation problem has been well-studied (e.g.Dieterich 1992;Rubin and Ampuero 2005), the behavior of such a dynamic nucleation process has not been comprehensively studied.
We investigated an aspect of the dynamic nucleation process through the implementation of self-similar crack propagation at a constant rupture propagation within a homogeneous target region.We estimated the critical fracture energy step ratio, γ c , which is uniquely related to the size of the critical dynamic nucleation zone, R dyn C , using a range of parameter sets for the initial stress, σ 0 , yield strength, σ p , and gradient of the slip-weakening dis- tance during self-similar rupture growth, D ′ C .We found that γ c increased as a function of the rupture velocity, C space for the anti-plane case v i , within the nucleation zone, whereas R dyn C decreased to ~ 70% of the static nucleation size.
Intuitively, it is not surprising that a rapidly propagating rupture can extend into a region with a high fracture energy, because inertia helps continuing the crack propagation owing to the higher rupture velocity.Nevertheless, it is notable that γ c increases quite rapidly with v i at a rupture velocity of > 0.6v s and approaches infinity at the limit of v i .This strong velocity dependence may explain why there are few earthquakes with slow rupture propagation velocities (< 0.5v s ), given several cascade-up processes are repeated in the growing process of earthquakes.If the likelihood of successful dynamic nucleation depends exponentially on the nucleation velocity, then a slower rupture is more likely to be arrested by fluctuating fracture energy on a fault system that possesses multiscale heterogeneity.Slow rupture propagation is suggested for some tsunami earthquakes (Kanamori 1972), but these events are quite rare.
It is the first time to have estimated the critical size of a dynamic nucleation zone for a rupture propagating at a constant velocity.We found that the dynamic size was smaller than the static size, but the difference was relatively small, with a maximum difference of up to 30%.This maximum difference is within the range of any uncertainty associated with the rupture process estimates from seismological observations.However, given that such a dynamic nucleation process is simply a building block in a cascading rupture process, with similar processes being repeated many times during the rupture process of a large earthquake, a 30% difference becomes significant in terms of forecasting the final size of the earthquake.
These findings provide important information for numerical simulations of earthquake rupture as a cascading process from one scale to another (e.g.Ide and Aochi 2005;Aochi and Ide 2009;Noda et al. 2013).While many simulations of an earthquake scenario generally employ a quasi-dynamic approximation, their difference from a fully dynamic nucleation approach is likely non-negligible, even in an earthquake cycle.This difference has not been welldocumented or discussed in previous studies but should be given more attention to advance our understanding of earthquake nucleation and rupture.

Definition of the nucleation rupture velocity in the numerical simulation
Averaging in both time and space is required in the numerical simulations to estimate v i , because the result- ant crack propagation is discretized in both spatial ( x ) and time ( t ) steps.We calculated the nucleation rup- ture velocity defined by the leading edge, v l i , as where t l (R) denotes the time when the half-width of the leading edge becomes R .For the trailing edge, v t i is defined using the average time for the rupture to move from t l (R 1 ) to t l R dyn : where R t [t] denotes the half-width of the trailing edges at time t.We set R 1 = 2/3R dyn to ensure that the rupture velocity has converged sufficiently.Figures 2 and 8 show the convergence properties of the rupture velocity.The rupture velocity at each spatial grid in Fig. 2a is obtained by dividing x by the time that it took for the crack tip to propagate to the next grid.The rupture velocity is quantized owing to the discretization in time.A moving average over five grids is used.There is good coherence between the two rupture velocities, v l i and v t i , as shown in Fig. 8. v l i tends to be slightly faster than v t i , which suggests that there may be a slight delay in the breakdown zone progress.However, this difference is negligible, such that the presented generality in our discussions in the main text remains valid.Here, we use v l i for v i in our numerical v l i = R dyn − R 1 t l R dyn − t l (R 1 ) , v t i = R t t l R dyn − R t t l (R 1 ) t l R dyn − t l (R 1 ) , simulations, because the leading edge is intuitively easier to recognize.
Since the BIEM requires a cohesive zone that is no smaller than four grid points (Kame et al. 2003) to appropriately resolve the weakening process, in practice, we used parameter sets that were chosen for a rupture velocity ratio of v l i − v t i /v l i > 0.04 in the anti-plane case.This ensures that the numerical expression of the cohesive zone includes more than four grids for a crack that is larger than a halfwidth of 100 grids.For the in-plane problems, we loosened the restriction to v l i − v t i /v l i > 0.03 , since the smaller time step and faster wave propagation improve the convergence of the rupture velocity.We also imposed the restriction that the cohesive zone must remain small enough to ensure that the shape of the slip distribution does not deviate too much from that of the analytical self-similar crack.While this latter restriction is not a physical requirement, it does avoid potential effects that depend on the details of the friction law.In practice, we require a crack size ratio of R t /R l > 0.70.

Fig. 1
Fig. 1 Schematic illustration of the numerical setting of the dynamic rupture simulations.a Slip-weakening friction law used in this study.σ p , σ 0 , and D C are the peak stress, initial stress, and slip-weakening distance, respectively.This study treats the residual friction as zero owing to the symmetry of the problem.The shaded region represents the fracture energy, G C .b Spatial distribution of D C (Eq. 1).The model space is divided into an inner nucleation zone and outer target region, with the two separated at |x| = R dyn .The step in G C at |x| = R dyn is quantified by the proportionality constant, γ .c, d Time evolution of slip (black line) for an anti-plane case, where the seismic nucleation initiated at x = 0 as a crack with a half-width of 20 x , which is negligible, under the D C distribution indicated by the blue line.The following parameters are used in the simulations: σ 0 = 4MPa , σ p = 12MPa , D ′ C = 6 × 10 −4 , and D ∞ C = 0.06 x .The time interval between each adjacent pair of solid lines is 160 t , and each gray dashed line is spaced 80 t from the next solid line to clarify the self-similarity.c Example of a failed nucleation, whereby the rupture decelerates and terminates.d Example of a successful nucleation, whereby rupture propagation initially decelerates until the crack reaches the static critical crack size R sta C , and rupture propagation then accelerates again.The difference between D ∞ C in c and d is 0.6%

Fig. 2
Fig.2Example of the simulated rupture velocity with rupture growth for an anti-plane case.The distance, r, and time are given by discrete units.The rupture velocity, v i , is normalized by shear-wave velocity, v s , in the upper panel.Small fluctuation in the velocity is due to quantization in space and time.The bold black and thick gray curves correspond to the rupture front (leading edge) and cohesive zone end (trailing edge), respectively.The theoretical R sta C is indicated by vertical dashed line.The artificial effect of sudden rupture onset disappears quickly by r/R sta C ∼ 0.2 , and the rupture velocity then remains quasi-constant until r/R sta C ∼ 0.7 .The rupture velocity drops suddenly at around r/R sta C ∼ 0.8 , because we introduce a fracture energy step, which controls the successful/unsuccessful rupture growth on the target region.This size is not always the same as R sta

Fig. 3 aFig. 4
Fig. 3 a Simulated v i in the nucleation area for given D ′ C values and the anti-plane case.D ′ C is normalized by 0.001 x .b Plot of the same v i values versus the normalized fracture energy rate, ξ = K µσ p D ′ C /σ 2 0 , from Eq. (5).The colors indicate different σ p values for a fixed initial value of σ 0 = 3 MPa

Fig. 6
Fig. 6 Ratio of dynamic to static critical nucleation sizes, R dyn C /R sta C , as a function of v i for the a anti-plane and b in-plane cases, respectively.See Figs. 3 and 9 for the parameter spaces, respectively.The function of 1/(ξ(v i ) • γ c (v i )) is plotted in v i /v s − R

Fig. 7 Fig. 8
Fig. 7 Effect of spatial smoothing on D C .A moving average of L = 2 n + 1(n = 0, . . ., 6) is employed, and the results are presented for the anti-plane case.a Estimation of R dyn C .b Estimaton of γ c .Note that v i is invariant, whereas D C increases linearly in the nucleation zone

Table 1
Model parameters