Geological earthquake simulations generated by kinematic heterogeneous energy-based method: Self-arrested ruptures and asperity criterion

1 Introduction

The study of spatial and temporal properties of seismic source parameters allows us to assess the radiated energy, amplitude, and frequency of ground motion [1–4]. Two of the most relevant source parameters are the slip rate u ̇ and final slip u f . This is due to the former constraining the radiated energy while the latter controls the spatial variability of self-arrested ruptures.

The slip rate u ̇ is proportional to the seismic moment rate M ̇ 0 which implies that sudden increase in u ̇ generates a burst of released seismic energy [5]. Thus, the consideration of realistic u ̇ is of uttermost importance because its proper consideration in heterogeneous faults implies better predictions of strong ground motion [6–8]. One way to describe u ̇ properly is by means of the Yoffe function [9]. Despite a few problems, such as indetermination, it has been proposed that the truncated version of the Yoffe function is consistent with earthquake dynamics for finite slip-weakening distances [10].

Second, the spatial domain that constrains where u ̇ evolves could be determined by the distribution of u f . Nevertheless, the proper distribution of u f is still a matter of debate. For example, several authors assume that u f is well described by stochastic von Karman correlation function in dynamic simulations [11,12]. Contrarily, other studies show that the complexity of u f can be well described by power law (fractal) distributions [13]. Due to the fact that one of the main features of fractal model is the lack of any characteristic length, and as von Karman model is characterized by it, the latter is commonly used in simulations [14]. Nonetheless, the usage of von Karman in detriment of fractal model is based on oversimplified distribution of u f or non-physical ruptures that cover the domain totally [14,15,16]. In contrast, the spatial complexity of u f from real earthquakes is associated more with the self-similarity of fractal distributions [17,18,19]. Besides, the self-similarity has also been found in the roughness of exhumed faults [20], in experiments describing the relative displacement of rock samples [21], stress states in systems of faults [22], or in the scale invariant path of cracks [23]. Therefore, self-similarity not only describes u f but is also a characteristic of various rupture parameters, such as mechanical properties, friction, shear resistance, fault proportions, source frequency content, initial stress concentration, changes in the b-value, transition to crack-like rupture, or rupture velocity scaling [24–34].

Special attention should be given to the zones that represent large slips, referred to as asperities, as they may also exhibit a distribution following power laws [35]. As the asperity zone accounts for about a quarter to a third of the total rupture area (asperity criterion) [36,37], it is expected that the total u f could be described as a power law as well. Thus, the abovementioned constraints support the idea that u f could be better described as power laws for realistic self-arrested ruptures that match the asperity criteria. Nevertheless, the non-consideration of these spatial heterogeneities leads to simulations of ruptures on flat surfaces which are considerably sensitive to the initial conditions [38]. This increases the difficulty to obtain simulations of self-arrested ruptures characterized by self-similar properties. Some computational simulations of dynamic ruptures have solved the problem of obtaining self-arrested ruptures. Nonetheless, the sensitivity of the initial conditions makes it difficult to prescribe any rupture magnitude. This implies that large self-arrested ruptures tend to be rare. For example, Ide and Aochi [39] have shown that the generation of a single self-arrested rupture of large magnitude requires thousands of small magnitude ruptures. Then, the generation of realistic self-arrested ruptures of specific prescribed magnitudes is also required.

In order to incorporate the abovementioned phenomena, Section 2 analyses the basis statistics of fractal distributions that match the asperity criteria. Section 3 develops an efficient model to generate infinite synthetic ruptures that are statistically equivalent to the dynamic self-arrested scenarios. The incorporation of truncated Yoffe’s function is also given in Section 3. Statistical analysis of the energy-based method constrained with fractal final slip for several examples is described in Section 4. Discussion and conclusion are provided in Sections 5 and 6, respectively.

2 Properties of random power law surfaces

As the u _f distributions could be regarded as fractal, it is worthwhile to describe the statistics of two random normalized fractal surfaces that can be created by using the methodology developed by Chen and Yang [40]. Note that this method allows for the generation of fractal surfaces where patches of maximum values are not necessarily located in the same spatial regions. This indicates that the final distribution of slip can have more than one asperity. The first one (Figure 1a) represents the case when the rupture reaches the domain boundaries (runaway rupture) and where its histogram is shown in Figure 1b. The second one (Figure 1c) represents the case of ruptures that are arrested before any section of the rupture front reaches the boundaries (self-arrested ruptures) and its histogram is shown in Figure 1d. The truncated version of the runaway case involves considering values between 0.5 and 1 (a.u.) from the original runaway case. Any values below 0.5 (a.u.) are set as zero, which is visualized as dark blue areas in Figure 1c. This selection leads to a concentration of positive final slip values in the central region of the domain, without extending to the domain’s boundaries. Note that a.u. stands for arbitrary units as they simplify calculations and allow for relative comparisons between variables. For both cases, the fractal dimension D is 2.2. The minimum value for the runaway-like and self-arrested-like cases are 0.0044 and 0.5 (a.u.), respectively. The histogram of the truncated fractal surface (red solid line) is compared to the runaway-like (black solid line) in Figure 1b. The histogram of the runaway-like surface is obtained by considering 20 bins from 0.0044 up to 1 (a.u.) while the truncated surfaces also consider 20 bins but from 0.5 to 1 (a.u.). The difference in the truncation is manifested as a short horizontal axis (red solid line between dashed magenta vertical lines in Figure 1b). If the horizontal values of the truncated histogram are scaled, that is, expand the range of values from (0.5, 1) to (0.0044, 1) the histogram in Figure 1d is obtained (red solid line). Note that these types of curves resemble second-order power law functions, which are defined as N ∼ a x b + c , where N represents the histogram and a , b , and c are constants that correspond to the best fit of the histogram. This can be seen as a dashed blue line in Figure 1d. The area of large values of Figure 1c can be compared to that used for real earthquakes. For instance, according to Somerville et al. [36], the zone of large slip known as asperity corresponds to an area between fourth and third of the total ruptured area. The definition of large slip corresponds to those places where the slip is larger than 1.5 times the average slip. By considering that the truncated and scaled fractal distribution of Figure 1c resembles the final slip distribution, the area that resembles the asperity corresponds to 23.9% of the total area while that resembling the large slip is larger than 0.5343 (a.u.). By considering the power law, the total ruptured area is proportional to the number of points N of the histogram. That is, A = N δ A 0 , where δ A 0 is a constant. Thus, the ratio between the total ( A T ) and asperity ( A asp ) area is as follows:

where u asp = 1.5 u ̅ and B 0 = a / ( b + 1 ) (Figure 1e). By using the values of a, b, and c (Figure 1d), the ratio represents ∼ 22 % of the total ruptured area which is similar to that described by Somerville et al. [36]. Thus, the second order power law describes a basic feature of ruptures. Then, a physical method for obtaining ruptures statistically equivalent to that shown in Figure 1 is required.

Figure 1

Difference of power law statistics as to whether it is a runaway or self-arrested case by using a fractal surface generated by the code developed by Chen and Yang [40]. (a) 2D Power law normalized distribution. This scenario resembles the runaway ruptures type because it covers the entire domain. (b) Histogram of the distribution shown in (a) (black solid line) and Figure c (red solid line). (c) Self-arrested type is characterized by a distribution that does not reach boundaries. This distribution corresponds to the truncated version of the original distribution shown in (a). That is, only considering those values larger than 0.5 a.u., the second color bar shows the scaled values. Then, scaling the horizontal axis of the truncated histogram (red solid line shown in (b)). This scale reveals that ∼25% of the area is within the zone of large values (>1.5 times the average). (d) Histogram of the distribution for self-arrested rupture shown in (c). The best fit (blue solid line) for the distribution that resembles the self-arrested case which is a second-order power law. (e) schematic representation of the second-order power law and the area that covers the asperity (A _asp).

3 Heterogeneous energy-based (HE-B) method

Dynamic ruptures evolve from no slip up to fractally distributed final slip. Thus, it is relevant to find out the manner to describe the transition from the initial to final state. One way to do this is by considering the kinematic source time functions. For instance, the slip rate function has been described by using the Yoffe function [9]. Nonetheless, the kinetic parameters should be constrained by considering realistic scenarios. One of them could be regarded as the energy-based method [41]. This method considers that the residual energy E res is the result of the available energy Δ W 0 minus the fracture energy G C for linear slip weakening law [41].

where x represents any points within the fault. In addition, those places within the fault where positive values of E res are found are also more prone to generate ruptures. That is, E res ( x ) > 0 implies slip and E res ( x ) ≤ 0 implies no slip.

The available energy Δ W 0 depends on a function that describes the heterogeneous initial stress states S 0 ( x ) [41,42]. Then, the available energy can be written as Δ W 0 ∼ S 0 . The available energy could also be affected by the existence of asperities which generates large stress coupling C [43]. This implies that the available energy can be regarded as:

where C corresponds to the coupling function that incorporates the influence of large slip deficit within faults. Note that the function C describes the mechanical coupling which corresponds to the areas where the stress is larger than the surrounded area [44]. The fracture energy G C means the energy per unit of area required to spread the rupture even further [45,46] and it could be regarded as a material property which is proportional to the geometrical irregularities within faults [47]. In addition, there is ample evidence indicating that the fracture energy is dependent on the material used. This has been demonstrated by various studies on rock samples which have shown that the fracture energy varies based on whether the material is brittle, ductile, or composite in nature [48–50]. Thus, each fault is characterized by their own heterogeneous values of G C . When the other dissipation mechanisms are also included, this energy is known as apparent fracture energy and tends to be related to the final slip [51]. As this geometry of faults has been described as power laws [20] and G C is related to the final slip despite the initial stress states [52], it means that the fracture energy could transfer its spatial properties to u f . That is, the final slip is also a function of the fracture energy and its spatial properties ( u f ∼ F ( G C ) ). Examples of S 0 , C , Δ W 0 , G C , and E res distributions per unit of area can be found in Figure 2a–e. The G C and S 0 distributions are power laws generated by the indications shown in the study by Chen and Yang [40]. Note that the maximum of G C is ∼ 10 9 J / m 2 because this is the expected order of magnitude obtained for largest earthquake with slip larger than 10 m (Figure 3 in Ke et al. [51] and Figure 5 in Nielsen et al. [52]). The domain size of the distribution is 1,000 × 200 and resembles a fault of 100 and 20 km of length and depth, respectively, where d x = d y = 0.1  km. Note that S 0 has been filtered at its shallower and deeper sections (blue colors in Figure 2a). This allows to obtain no rupture at these depths and allows to define a seismogenic zone (orange to red shades in Figure 2a). The coupling functions arise from heterogeneous stresses within the plates, making it impossible to create them analytically. However, since the coupling distributions of real faults tend to be smoothed, a method is proposed to randomly generate them. That is, the coupling function has been obtained by using the sum of seven bimodal Gaussian distributions where the standard deviation depends on the seismic moment M 0 ( σ = σ ( M 0 ) ). The center of each distribution is randomly located within a specific rectangle of length L ∼ M 0 2 / 5 and width W ∼ L β [29], where β is a constant. Thus, the coupling function is defined by the prescribed seismic moment of reference. In Figure 2b, an example of the distribution of C for a reference dimension L and W equivalent to a magnitude of 6.8 can be observed. The white dot in Figure 2b indicates the maximum coupling function. The available energy is shown in Figure 2c, and the residual energy is shown in Figure 2e. Here the negative values of E res has been set as zero. Note that the distribution of E res tends to zero at its boundaries. The final slip also tends to be zero at the rupture area boundaries for self-arrested ruptures [53]. This spatial feature shared between E res and u f implies that u f should also depend on E res ( u f ∼ F ( G C , E res ) ). Additionally, the rupture velocity v r is positively correlated to the peak slip rate u ̇ p while a weak correlation between u f and u ̇ p have been found [54–58]. This also implies that a certain dependence between u f and v r could exist. This leads to the following equation:

Figure 2

Example of the generation of HE-B model. (a) The background stress states S 0 is generated by a random fractal surface [40] and filtered by a function that forces zero values at the upper and lower boundaries. (b) The coupling factor C is considered as a set of normalized bimodal Gaussian distribution where its center is located within a rectangle of width and large governed by the Leonard scaling relations [29]. (c) The available energy Δ W 0 corresponds to the multiplication of S 0 and C . (d) Another normalized random fractal surface is generated. This represents the fracture energy. The difference between Δ W 0 and G C is shown in (e). The negative values are not considered. (f) The image depicts the gradient of residual energy ∣ ∇ E res ∣ . (g) The rupture velocity obtained by using equation (5) for a 0 = b 0 = − c 0 = 1 . The obtained rupture time by using equation (6) is shown in figure (h). Note that the magnitude order of the energy per unit of surface values has been taken from real earthquakes [52] and the coupling function is defined between 0 and 1.

Figure 3

Different snapshots that show how the slip rate evolves within the fault and in time. (a) The image shows that the rupture begins at a specific location (maximum E res ). That is at a length and depth of ∼20 and ∼10 km, respectively. (b), (c), and (d) display the times when the rupture behaves as an elliptically-shaped pulse-like rupture. The rupture tends to accelerate to the right in (e). This is because the rupture has reached the place where the maximum coupling is located (length and depth of ∼40 and ∼10 km, respectively). The burst of slip rate is still present in (f). (g) and (h) show how the rupture tends to arrest.

The v r also should tend to zero when it arrives at the boundaries of E res . This implies that a direct spatial relation between v r and E res could also exist. In order to incorporate the spatial variability of E res , the rupture velocity could be written as follows:

where a 0 , b 0 , and c 0 are constants. Note that the sign of c 0 should be negative as the rupture onset is set at the place where the largest E res is found (blue star in Figure 2). This implies that v r could only decrease towards the rupture boundaries. The gradient of E res is shown in Figure 2f while the rupture velocity is shown in Figure 2g for a 0 = b 0 = 1 and c 0 = − 1 . The initiation of rupture is determined by the rupture time t 0 [59], which is commonly related to the rupture velocity v r through the expression v r ∼ 1 / ∇ t 0 . Therefore, by integrating this expression, the rupture time can be obtained using the following equation:

where d 0 , e 0 , and f ₀ are constants. The function t ( r ) is the initial rupture time at each point within the fault for circular of elliptical ruptures and represent the homogeneous case. This implies that equation (6) could be regarded as a homogeneous rupture time plus the heterogeneities concentrated in the third term that could accelerate or delay the rupture time. For example, the rupture can easily spread if G C is small, while zones with large G C could take more time to rupture. In other words, G C could rush or delay the onset of the rupture. For radial (circular or elliptical) homogeneous rupture, the second term of equation (6) can be written as e 0 t ( r ) ∼ e 0 ′ r . Note that small values of f 0 could generate homogeneous ruptures. In contrast, high values of f 0 can cause certain areas of the fault to initiate rupture before the arrival of the rupture front. Thus, the proper election of this factor could ensure that the heterogeneous rupture front velocity is lower than the super shear regime [60–62]. For example, the factor f 0 can be constrained by the time required for the rupture to spread from the nucleation point to the furthest point of the ruptured area. This could be achieved by f 0 ∼ R L / ( 0.8 v s ) , where R L is the domain length and v s is the shear velocity ( 3 , 464 m / s ) (Table 1). The example of t 0 can be seen in Figure 2h.

Additional considerations should be made in order to obtain realistic kinematic ruptures. For instance, the rupture could be regarded to begin at the maximum value of E res . This consideration set the origin of r . In addition, G C is proportional to the final slip u f which is proportionally related to the rise time τ R [55]. That is: τ R ∝ G C . The final time t f at any point in the fault is defined as t 0 + g 0 × τ R , where g 0 is a dimensionless factor that can dilate or constrict the rupture front. That is, small (large) values of g 0 could describe pulse-like (crack-like) ruptures. Note that g 0 has been constrained by ∼ 1 / τ ¯ R and implies that large τ R can be reduced by its average in order to obtain values of g 0 ∼ 1  s which can be regarded as pulse-like ruptures (Table 1) [63]. The slip rate is described by the truncated Yoffe function which corresponds to the convolution between the normal Yoffe function and a triangular function of width t s [10]. The amplitude of the truncated Yoffe’s function is normalized and then integrated in order to obtain the slip in terms of time at each point within the fault. Then, the slip at each point is set as proportional to the equation as follows:

where F 0 is a constant. Note that equation (7) corresponds to the simplest version of equation (4). Finally, the slip is differentiated in order to obtain the time series of slip rate for each point in the fault. An example of the slip rate distribution can be found in Figure 3a–h. Here the slip rate is shown for different times. The parameters used are shown in Table 1. The rupture begins where E res is maximum. This is shown in Figure 3a where the rupture is concentrated close to the length and depth, 20 and 10 km, respectively. This can be seen as a blue star in Figure 2e. The first stage (Figure 3b–d) shows the ruptured front is mainly elliptically shaped as it was set. Nevertheless, the heterogeneities are relevant for Figure 3e, and f. This is because the rupture is entering the zone where the maximum coupling function is found (length and depth of ∼40 and ∼10 km, respectively, in Figure 3e, and f and in the white point in Figure 2b). The arrest stage is shown in Figure 3g and h where the slip rate tends to zero.

The obtained final slip can be seen in Figure 4a. This final slip matches the asperity criteria because a portion of the rupture area releases most of the total seismic moment. Specifically, 29.1% of the total area (area within the black solid line) releases 70.9% of the total seismic moment. It is possible to observe that the histogram of u f (solid black line in Figure 4b) also tends to be similar to the second-order power laws (red solid line Figure 4b). The seismic moment rate is shown in Figure 4c. It can be seen that the maximum M ̇ 0 is found when t ∼ 12 s. This corresponds to the time when the right side of the ruptured front enters the area characterized by a large coupling (white point in Figure 2b). This can roughly be seen in Figure 3e for the slip rate where a patch of yellow to orange shades cover an area close to a length of 40–47 km at a depth of 9–14 km. Figure 4c shows the impulsive increases in M ̇ 0 for t ∼ 10  s (vertical orange dashed line). The duration of this burst of released seismic moment lasts ∼5 s (Figure 4c). Figure 4d shows the frequency amplitude of M ̇ 0 . This spectrum is characterized by a double corner frequency ( f Cc 1 = 0.019 Hz and f c 2 = 0.21 Hz ). The decay f − 1 is found between f c 1 and f c 2 and f − 2 is found for f > f c 2 . Note that the frequency f ∼ 0.068 Hz represents the projection of the spectrum with a decay of f − 2 towards the plateau level. This feature of double corner frequency has also been observed for real earthquakes [64,65]. These corner frequencies are related to the total rupture time ∼35 s ( Δ t 1 in Figure 4c) and the time required to rupture a coupled zone ( Δ t 2 ∼ 5  s). In terms of frequency, this is 1 Δ t 1 ∼ 0.02 Hz and 1 Δ t 2 ∼ 0.2 Hz which are the corner frequencies shown in Figure 4d. Note that no consideration has directly been taken in order to constrain the decay of this spectrum. That is, the spectrum appears naturally from the considered heterogeneities and rupture size.

Figure 4

Statistical analysis of the example shown in Figures 2 and 3. (a) The final slip shows that the rupture does not reach the domain’s boundaries. The percentage of asperity is ∼29% and releases ∼70% of the total seismic moment. (b) The statistic of the obtained u f reveals the second-order power law tendency. (c) The seismic moment rate M ̇ 0 is also obtained. The maximum value is found when t ∼ 12 s . That is the time when the rupture front enters the main coupling (Figure 3e and f). (d) Spectrum of M ̇ 0 . This is characterized by a plateau (yellow segmented line) and a double corner frequency. The decay is f − 1 (blue segmented line) and f − 2 (red segmented line).

4 Examples

Let us consider two sets of kinematic simulations to account for the effects of changes in fault’s heterogeneities. The first set is that when the distribution of G CC and S 0 are the same and different distributions of C are considered. The second set is when C and S 0 are the same and changes come from the distribution of G C .

4.1 Changes in the coupling distribution

The different magnitude examples can be done by considering an initial magnitude of reference M W ref . This allows to set the boundaries where the center of each bimodal Gaussian function can be located. Specifically, the boundaries are determined by the Leonard relation [29] which established the scaling between the width and length of rupture in terms of a specific magnitude (Table 1). Note that the standard deviation σ also could be regarded as a function of the magnitude ( σ = σ ( M W ) ). Some examples of different distributions of C are shown in Figure 5a–j. Particularly, the increase in the reference magnitude generates large scaling relations between width and large (Figure 5a, c, e, g, and i) which generates large, ruptured areas (Figure 5b, d, f, h, and j). The full results are shown in Tables 2 and 3. On average, the ruptured area that can be regarded as an asperity is 25.6303% and releases 63.0439% of the total seismic moment (Table 2). This confirms that the spatial properties of u f are self-similar regardless of the distribution of C . Regarding the scaling of frequency spectrum, Figure 6a shows how the plateau increases with the magnitude as well as the corner frequencies decrease. Decay found in these spectra is similar to that found in Section 3. That is, there is a decay of f − 1 between the two corner frequencies and a decay of f − 2 for the high frequency content. In reference to the distribution of final slip shown in Figure 6b, it is evident that they display a self-similar characteristic. This implies that each distribution can be described by a second-order power law ( ∼ a u f b + c ), irrespective of the earthquake magnitude. The varying shades from magenta to black in Figure 6b indicate earthquakes of magnitudes close to the range of 5–7. The seismic moment rate is shown in Figure 7a. Both, the peak M ̇ 0 and rupture duration depends on the magnitude. Large magnitude earthquakes (black shades in Figure 7a) tend to have large max M ̇ 0 and last for more than 25 s. This correlation between max M ̇ 0 and M W can be seen in Figure 7b. In addition, this figure also shows that max M ̇ 0 , M W , and the maximum fracture energy within the ruptured area ( max G C ) are well correlated among them. That is, correlations larger than 0.85. An almost perfect correlation exists between the maximum frequency f max and the average fracture energy ( G ¯ C ) as it is shown in Figure 7c. However, no correlation is found between G ¯ C and M W and between f max and M W . This implies that there are differences between the physical interpretation of G ¯ C and max G C . Besides, f max does not depend on the earthquake magnitude. Figure 7d shows the positive correlation between the area of asperity and the amount of seismic moment that is being released for asperity. The larger the area, the larger the released seismic moment. In contrast, no correlation is found between the latter and former to f max . That is, f max does not depend on the magnitude, the asperity area, and the seismic moment released by it.

Figure 5

Example of increasing heterogeneous magnitude. The full details of the results is listed in Tables 2 and 3. It can be seen that the reference magnitude generates a random coupling function that satisfies the scaling relation proposed by Leonard [29]. This implies that each E res is restricted to a specific area and generates a rupture of similar magnitude. (a) Displays the coupling functions for a reference magnitude of Mw5, while (b) shows the resulting final slip. Similarly, (c), (e), (g), and (i) display the couplings for reference earthquakes with magnitudes of 5.5, 6.0, 6.5, and 7.0 respectively. The final slips for each of these coupling functions are shown in (d), (f), (h), and (j) respectively.

Table 2

Values obtained by different simulations of different magnitudes

M Wref	M W	u som   ( m )	Asperity  ( % )	M 0   ( % )	a	b	c
5.0	5.0519	0.1862	24.4611	63.6338	1600.1427	−0.0170	−1614.8813
5.1	5.1103	0.1982	25.5094	63.5403	665.4257	−0.0424	−679.9147
5.2	5.0997	0.1958	26.8170	62.2115	−640.9314	0.0519	623.7332
5.3	5.1963	0.2190	24.5649	59.1244	−683.2097	0.0601	665.5470
5.4	5.3888	0.2716	24.7119	59.3898	4948.3241	−0.0137	−4943.67243
5.5	5.4056	0.2780	25.0309	59.6445	−1937.4294	0.0343	1929.6274
5.6	5.5770	0.3380	25.2357	60.5781	−4680.2396	0.0214	4694.5254
5.7	5.6702	0.3759	25.3745	61.9355	−1205212.6361	0.0001	1205240.3304
5.8	5.8517	0.4618	25.7776	62.6430	−19169.1590	0.0093	19244.7341
5.9	5.9808	0.5357	25.6975	62.6764	165862.0483	−0.0014	−165720.8613
6.0	6.0114	0.5547	24.6447	59.8050	−9725.3506	0.0281	99071.4110
6.1	6.1643	0.6624	27.0784	62.0645	−7185.3387	0.0492	7445.8838
6.2	6.2925	0.7654	25.1499	60.2559	−23764.3912	0.0225	24345.8120
6.3	6.4236	0.8916	24.9902	62.5688	−46234.8098	0.0131	46837.8276
6.4	6.5612	1.0399	24.3248	61.3599	164746.2798	−0.0058	−163456.1167
6.5	6.6942	1.2126	26.0065	65.3109	25741.0613	−0.0420	−24332.0990
6.6	6.8592	1.4651	25.4644	66.4560	23729.0286	−0.0626	−21611.2308
6.7	6.9052	1.5466	24.7817	61.6741	−136758.5092	0.0130	139540.9435
6.8*	7.0814	1.8926	29.0826	71.7570	15015.9276	−0.1293	−11921.6992
6.9	7.1841	2.1309	26.5846	68.8305	51295.4564	−0.0557	−46383.1040
7.0	7.1973	2.1651	26.9484	68.4622	38802.5636	−0.0721	−33986.6888
Average	—	—	25.6303	63.0439	—	—	—

These results are obtained by changing the coupling function and letting the other parameters remain the same. Besides the reference magnitude M Wref is used to constrain Δ W 0 , while the other results listed are the real magnitude M W , the Somerville average slip is defined as the area where u f is larger than 1.5 × u ¯ f [36], percentage of rupture area that corresponds to the asperity, amount of seismic moment released by the asperity and the second-order power law that fit the histogram for each rupture (a, b and c). The *means that this is the simulation used as example in Section 3.

Table 3

Other results obtained by different simulations of different magnitudes (same simulations as those presented in Table 2)

M Wref	M W	u ̇ p   ( m / s )	g 0	f max   ( Hz )	G ¯ C   × 10 8 J m 2	G Cmax × 10 8 J m 2	e 0   ( s )
5.0	5.0519	2.2139	1.0185	8.9552	5.8913	6.9557	1
5.1	5.1103	2.4226	1.1491	7.9372	5.2216	6.7075	1.4500
5.2	5.0997	2.2840	0.8135	11.2113	7.3754	8.3420	1.9000
5.3	5.1963	2.5343	0.7384	12.3519	8.1259	9.1527	2.3500
5.4	5.3888	4.4793	0.9942	9.1739	6.0351	6.9125	2.8000
5.5	5.4056	3.4795	0.7546	12.0857	7.9507	9.1527	3.2500
5.6	5.5770	4.9094	0.8662	10.5293	6.9268	8.3879	3.7000
5.7	5.6702	5.2705	0.9275	9.8335	6.4691	8.6006	4.1500
5.8	5.8517	5.9096	0.9729	9.3749	6.1674	8.0444	4.6000
5.9	5.9808	6.9127	0.9648	9.4533	6.2190	7.7943	5.0500
6.0	6.0114	7.8100	0.9561	9.5391	6.2754	9.2076	5.5000
6.1	6.1643	9.7900	0.9028	10.1028	6.6462	10.3750	5.9500
6.2	6.2925	12.9623	0.9120	10.0000	6.5786	9.1767	6.4000
6.3	6.4236	9.2945	0.8301	10.9877	7.2283	10.1323	6.8500
6.4	6.5612	17.5137	0.8888	10.2613	6.7505	9.5637	7.3000
6.5	6.6942	11.5334	0.8027	11.3621	7.4747	10.9065	7.7500
6.6	6.8592	16.3985	0.8879	10.2718	6.7574	10.0888	8.2000
6.7	6.9052	19.2114	0.9074	10.0508	6.6120	11.4204	8.6500
6.8*	7.0814	15.7329	0.8766	10.4049	6.8450	11.5286	9.1000
6.9	7.1841	16.5019	0.8976	10.1614	6.6848	12.0000	9.5500
7.0	7.1973	18.7485	0.9036	10.0930	6.6398	11.7488	10

This table shows the maximum slip rate of the rupture u ̇ p , g 0 , the maximum frequency f max , the average fracture energy G ¯ C , the maximum fracture energy G Cmax , and e 0 .

Figure 6

(a) The frequency spectra and (b) the power law histogram for ruptures shown in Tables 2 and 3. The only change is in the coupling function. Color gradient means the magnitude. Magenta is for those magnitude close to 5 and the black means magnitude close to 7.2. It can be seen that the basic properties of real earthquakes are obtained. For instance, the amplitude of the plateau is large for large earthquakes and its corner frequency tend to be smaller (a). (b) the histogram for each rupture. The larger the earthquake, the larger the final slip and the number of points (area). However, the same second-order power law is obtained regardless of the magnitude.

Figure 7

Seismic moment rate and correlation among different rupture parameters. (a) Seismic moment rate of different magnitude ruptures where the ruptures are the same that are described in Tables 2 and 3. Ruptures only differ in the coupling function. Large M ̇ 0 (∼7.2) tends to be in black shade and smaller M ̇ 0 (∼5) is in magenta shade. (b) Positive correlation among peak seismic moment rate ( max M ̇ 0 ), maximum fracture energy ( max G C ), and magnitude ( M W ). (c) The maximum frequency ( f max ) is correlated to the average fracture energy ( G ¯ C ) but none of them are correlated to the magnitude. (d) The percentage of ruptured area that can be regarded as an asperity tend to be positively correlated to the amount of the released seismic moment. None of them are correlated to the maximum frequency.

4.2 Changes in the fracture energy distribution

The next set of examples considers that the changes are only performed in the distribution of G C . That is, the same distribution of S 0 and C are used. The results are shown in Tables 4 and 5. Here it is possible to observe that all the cases turned out to be ruptures with similar magnitude (∼7). Figure 8 shows some examples. Note that the number of simulations is located at the upper left corner of each plot in Figure 8 and the full results of the labeled examples can be seen in Tables 4 and 5. It can be seen that all the final slip distributions tend to be quite similar among them. In fact, the asperity area and the seismic moment released there is almost the same for all the cases (average of 28.6848 and 68.9061%, respectively, as shown in Table 4). This indicates that the parameter that drives the size of earthquakes is the coupling function. In addition, note that the maximum coupling (white point in Figure 8), initial rupture point (blue star in Figure 8), and maximum slip (magenta triangle in Figure 8) may not be all located in the same place.

Table 4

The same parameters described in Table 2 are presented. However, these 20 new simulations have the same coupling function and S 0 distribution but different fracture energy distributions G C among them. While the coupling function uses a magnitude 7 M W earthquake as a reference, earthquakes of similar magnitude to the reference are generated regardless of the G C distribution

Number #	M W	u som   ( m )	Asperity  ( % )	M 0   ( % )	a	b	c
1	7.0814	1.8926	29.0826	71.7570	15015.9276	−0.1293	−11921.6992
2	7.1004	1.9392	29.6045	69.4418	35906.5182	−0.0559	−32579.4075
3	7.0583	1.8404	27.9023	69.6037	19477.0796	−0.1054	−16217.7947
4	7.0209	1.7679	28.7229	66.8259	34023.6423	−0.0522	−31205.1155
5	7.0735	1.8747	29.2123	71.9647	10651.2396	−0.1716	−7718.0547
6	7.0745	1.8784	27.5368	68.9274	19379.2027	−0.1044	−16145.9520
7	7.1004	1.9392	29.6045	69.4418	35906.5198	−0.0559	−32579.4090
8	7.0597	1.8488	28.9776	69.3402	21215.8819	−0.0854	−18309.9258
9	7.0977	1.9333	29.0871	68.3850	110937.0692	−0.0208	−107090.9096
10	7.0525	1.8294	27.5905	68.7998	25189.6916	−0.0793	−22026.6679
11	7.0827	1.8897	28.3338	73.4226	11593.7537	−0.1666	−8582.2291
12	7.0593	1.8503	29.8987	68.0159	110799.0800	−0.0172	−107663.9327
13	7.0888	1.9127	27.0004	67.6688	42955.4155	−0.0485	−39575.5718
14	7.0652	1.8588	29.0743	70.5901	15105.8194	−0.1177	−12265.9039
15	7.0847	1.9017	29.8438	68.0794	48611.6106	−0.0429	−45130.5152
16	7.0588	1.8463	28.3138	68.9105	24810.5266	−0.0748	−21842.2727
17	7.0747	1.8857	27.6533	65.7931	−58819.5243	0.0364	62400.5053
18	7.0528	1.8380	28.5759	67.1862	78934.4536	−0.0236	−75905.1812
19	7.0663	1.8680	29.6638	67.8802	1640425.8229	−0.0013	−1637030.0899
20	7.0636	1.8620	28.4546	66.0260	−57646.2064	0.0374	61267.2608
21	7.0884	1.9110	28.2483	68.9678	44698.7289	−0.0519	−40877.8246
Average	7.0716	1.8747	28.6848	68.9061	106151.0600	−0.0634	−102904.7948

Nevertheless, the fracture energy is changed. This implies that the coupling function and S 0 distribution are the same. The changes in G C generate earthquakes of similar magnitude ( M W ∼ 7 ).

Table 5

Same results as that shown in Table 3 but for the examples described in Table 4

Number #	M W	u ̇ p   ( m / s )	g 0	f max (Hz)	G ¯ C   × 10 8 J m 2	G Cmax × 10 8 J m 2
1	7.0814	15.7329	0.8766	10.4049	6.8450	11.5286
2	7.1004	16.7908	1.1711	12.9430	5.1233	8.7511
3	7.0583	13.6012	0.7832	11.4797	7.6608	11.8212
4	7.0209	13.5083	0.7622	8.8000	7.8715	12.0000
5	7.0735	15.4551	0.8202	13.5270	7.3152	12.0000
6	7.0745	15.3344	0.8776	13.5853	6.8372	11.6546
7	7.1004	17.2497	1.1711	12.9430	5.1233	8.7511
8	7.0597	15.6780	0.8515	20.6447	7.0461	11.2354
9	7.0977	16.2587	1.1177	17.8148	5.3681	8.7329
10	7.0525	12.9551	0.8238	19.5617	7.2832	12.0000
11	7.0827	16.6316	0.7622	45.7667	7.8723	12.0000
12	7.0593	16.3859	1.0190	27.1925	5.8880	9.8759
13	7.0888	15.2155	1.1151	25.6441	5.3806	8.4600
14	7.0652	15.3827	0.8254	15.4226	7.2690	11.1246
15	7.0847	16.8932	1.1785	13.7856	5.0912	9.3091
16	7.0588	13.9613	0.8431	11.4010	7.1163	11.3021
17	7.0747	13.7973	1.0753	21.2124	5.5796	9.8109
18	7.0528	13.3773	0.8097	10.4060	7.4103	12.0000
19	7.0663	14.4638	0.8599	17.6595	6.9773	11.8272
20	7.0636	12.4038	0.8883	16.7349	6.7544	10.6858
21	7.0884	13.3394	0.8572	19.7810	6.9995	12.0000

Figure 8

(a)-(j) showcase examples of ruptures appearing in Tables 4 and 5 when different fractal distributions of fracture energy are used alongside C and S 0 distributions that are the same for all cases. (a), (c), (e), (g), and (i) depict the distributions of fracture energy for simulations numbered 2, 6, 11, 16, and 21 (Tables 4 and 5). (b), (d), (f), (h), and (j) display their respective final slips. Note that the change in G C results in locations of maximum coupling (white point), rupture onset (blue star), and maximum displacement potentially differing, even though the shapes of the ruptures remain similar among them.

The change in G C generates frequency content that is almost identical among all the cases (Figure 9a). That is, the same double corner frequency and the double spectrum decay ( f − 1 and f − 2 ). Similarities are also found in the final slip distribution (Figure 9b). The change in G C does not affect the basic statistics of ruptures. Minor differences can be seen in the seismic moment rate (Figure 10a). It can be observed that there are slight differences in M ̇ 0 due to the change in G C . Nevertheless, the main shape of M ̇ 0 is shared among all the cases. Despite these minor changes, there are some evident changes in the correlation among rupture parameters. For instance, there are weak correlations among max G C , max M ̇ 0 , and M W (Figure 10b) while the correlation between f max and G ¯ C vanished (Figure 10c). Here a weak negative correlation is found between G ¯ C and M W . Figure 10d shows that there is no correlation among the asperity area, the seismic moment released, and f max .

Figure 9

(a) Frequency spectrum of the earthquakes listed in Tables 4 and 5. The spectrum behaves similarly across all the cases. That is, the double corner frequency determined by the double decay ( f − 1 and f − 2 ). (b) The histogram reveals that the statistics of final slip remain identical among the cases regardless of the fracture energy distribution.

Figure 10

(a) Seismic moment rate. Regardless of the fracture energy distribution, the main shape of M ̇ 0 remains almost the same. (b) A weak negative correlation exists among max G C and max M ̇ 0 (∼ −0.38). A weak positive (negative) correlation exists between max M ̇ 0 and M W ( max G C and M W ). (c) No correlation is found between f max and G ¯ C and between f max and M W . A weak negative correlation is found between G ¯ C and M W . (d) There is no correlation among the asperity percentage, seismic moment released for that asperity, and f max .

5 Summary and discussions

6 Conclusion

Based on the discussions presented in the previous sections, the utilization of fractal distributions to describe the distribution of final slip and rupture parameters constrained by energy conditions leads to kinematic simulations of earthquakes that can reproduce properties observed in real earthquakes. Additionally, being heterogeneous kinematic simulations, they can contribute to the study of seismic sources and, potentially, improve predictions of strong ground motion with reduced computational resources and increased efficiency. The main conclusions regarding the development of the HE-B model are listed below:

• The final slip distribution of self-arrested ruptures resembles truncated fractal surfaces.

• Somerville et al. [36] asperity criterion is a property of fractal surfaces.

• The energy-based model presented by Noda et al. [41] can be spatially constrained by heterogeneous fractal fracture energy and coupling distributions.

• The work done by Leonard [29] can constrain the coupling function in terms of a required earthquake magnitude.

• The heterogeneous coupling function determines the residual energy which defines the total ruptured area.

• The fractal final slip can be obtained by considering the residual energy, the fracture energy and rupture velocity in order to generate realistic self-arrested ruptures.

• The HE-B method leads to the same spatial statistics in the final slip as those found in real self-arrested ruptures.

  As the HE-B model enables the generation of realistic simulations, it also facilitates the analysis of the rupture process, such as spatial correlations among rupture parameters and the association of corner frequencies with source properties, such as the duration of asperity rupture. These findings can be summarized as follows:

• Double corner frequency is naturally obtained by using the HE-B model.

• The HE-B model also leads to the decay in the frequency content characterized by f − 1 and f − 2 .

• The double corner frequency could be a feature of the rupture of a large coupled area.

• The maximum fracture energy is related to global parameters such as magnitude and peak seismic moment rate for a given fault.

• The average fracture energy is related to the maximum source frequency for a given fault.

• No significant correlation is found among maximum and average fracture energy, maximum source frequency, peak seismic moment rate, and magnitude for different faults even though they share the same magnitude and available energy distribution.

Future works should explore other realistic approaches such as fracture energy dependence on depth or temperature, presence of fluids, or changes in the fault’s strength over time.