Dependence of Sandpile Avalanche Frequency–size Distribution on Coverage Extent and Compactness of Embedded Toppling Threshold Heterogeneity: Implications for the Variation of Gutenberg–richter B Value

The effects of the spatiotemporal evolution of failure threshold heterogeneity on the dynamics of fault crit-icality, and thus on regional seismogenesis, have attracted strong interest in the field of regional seismotectonics. The heterogeneity might be a manifestation of the macroscopic distribution and multiscale strength variation of asperities, the distinct regional stress level, and (microscopically) heterogeneous fault surface roughness or friction regimes. In this study, rather than attempting to mimic the complex mi-croscale slipping physics on a fault surface, sandpile cellular automata were implemented with a straightforward toppling rule. The objective is to examine the influence of distinct configurations of the embedded heterogeneous toppling threshold field on the global system avalanche event statistics. The examination results revealed that increasing the coverage extent and decreasing the compactness of the heterogeneous failure threshold, rather than the magnitude, range of contrast , diversity, or the geometric configuration of the threshold heterogeneity, leads to a systematic increase in the scaling exponent of the avalanche event power law statistics, implying the importance of mutual interaction among toppling sites with distinct thresholds. For tectonic provinces with differing stress regimes evolving spatiotemporally, it is postulated that the distinct extent and compactness of the heterogeneous failure threshold are critical factors that manifest in the reported dynamic variations of seismicity scaling.


Introduction
Recent studies of repeating earthquakes (e.g., Uchida et al., 2012;Yu, 2013) and mega-earthquakes such as the 11 March 2011 Tohoku-oki earthquake (Lay and Kanamori, 2011;Hashimoto et al., 2012;Tajima and Kennett, 2012) have indicated that heterogeneous slipping thresholds, which refer to the spatiotemporally distinct limits of failure stress that must be overcome to initiate slipping such as that manifested because of variations in the plate coupling strength (e.g., Yu, 2013), introduce dynamic complications in the plate boundary seismogenesis (Chlieh et al., 2011).Brodsky and Kanamori (2001) showed that hydrodynamic lubrication can help in reducing fault friction and that earthquakes can increase the effective permeability of distant aquifers.Furthermore, rate-and state-dependent friction regimes also lead to the manifestation of dynamic failure heterogeneity along a fault surface (Skarbek et al., 2012), whereas Prieto et al. (2013) discussed thermal runaways that are controlled by temperature, pressure, and the degree of partial melting as another potential heterogeneous mechanism.Wang and Bilek (2014) proposed that relatively small events dominate in areas of rough seafloor subduction, whereas a smooth subduction contact interface tends to be associated with a greater number of large earthquakes.In summary, the failure thresholds of fault surfaces are extremely likely to be heterogeneous and evolve dynamically within seismogenic zones.
The beginning of modern modeling of macroscale earthquake physics and thus the scale invariant power law statistics can be attributed to the classical mechanical sliderblock or Burridge-Knopoff (B-K) train model (Burridge and Knopoff, 1967).Subsequent follow-up modeling studies using cellular automata approaches (Carlson and Langer, 1989;Feder and Feder, 1991;Rundle and Klein, 1993;Maya and Stefan, 1996;Gran et al., 2012) as well as the continuous Olami-Feder-Christensen (OFC) simplifications (Olami et al., 1992) have been successful in reproducing the power law and have highlighted distinct characteristics related to the seismic activity.Bak et al. (1987) implemented the Bak-Tang-Wiesenfeld (BTW) sandpile to model complex dynamic systems categorized as self-organized criticality (SOC).The merits of the BTW sandpile model is that with its simple microscopic redistribution rule, it is capable of reproducing the macroscale criticality of complex systems.These complex systems, including regional seismicity, are typically characterized by the independence of the initial configuration, slow external drive, and capability of achieving criticality as manifested in the distinct event frequencysize power law distribution without the presence of additional fine tuning.Bak and Tang (1989) also used a sliderblock-type model, with a redistribution rule very similar to the BTW sandpile model, to study the scale-invariant statistics of regional seismicity by considering the regional seismicity to be an SOC phenomenon.Briefly, cellular automata approaches are usually set up on a two-dimensional square lattice.In a slider-block or spring-block model, the blocks in each lattice cell are subjected to a force exerted by their neighbors (connected by springs) and a constant "tectonic" driving force from the loader plate overlaid and connected to all the blocks by the driving spring.When the accumulated force on a block exceeds a predefined failure threshold, the block slips to a nearby position to release the accumulated strain while redistributing and increasing the force that it exerts on the neighboring blocks.Similarly, the toppling threshold of a sandpile model is the maximum number of sand grains a cell can sustain; external sand grains fall on the lattice slowly and randomly until the total number of grains in a particular cell exceeds the threshold.Toppling then occurs, and, consequently, the cell is emptied and its grains redistributed among neighboring cells.All these redistributions of force or sand grains are potentially capable of inducing large avalanche events involving a considerable number of cells in the cascaded loading-unloading process (e.g., Bak et al., 1987;Bak and Tang, 1989;Carlson and Langer, 1989;Feder and Feder, 1991;Christensen and Olami, 1992;Olami et al., 1992).Interestingly, although with the aforementioned seismic system heterogeneity, most typical modeling efforts had invoked homogeneous failure threshold fields.Rundle and Klein (1993) pioneered the slider-block cellular automaton with heterogeneous failure threshold fields.The heterogeneities are either arranged randomly or in a fractal distribution.They showed that the frequency-size scaling associated with these heterogeneous failure threshold fields behaves differently from that associated with a uniform threshold field.They postulated that when the dispersion in the failure threshold decreases, the critical system with the existence of a spinodal line crosses over to the SOC regime and that pronounced heterogeneous failure threshold enforces the dependence of the critical behavior on the velocity of the loader plate.This loader plate velocity is a crucial parameter in the slider-block model, and it is sometimes envisioned to be analogous to the plate convergence rate that acts as the slow external drive for plate boundary seismogenesis.These findings are intriguing because a heterogeneous failure threshold, as a plausible manifestation of varying roughness and friction regimes on the fault surface, is consistent with field observations and general geological understanding.However, correlation of the seismicity criticality variation with the tectonic plate kinematics, analogous to the loader plate velocity in the slider-block model, has not been empirically established.
Among sandpile models, extensive studies of the deterministic, Abelian BTW model (Dhar, 1999(Dhar, , 2006)), as well as the comparison of several alternative variants to identify common features, have shown the capability of these models to reproduce the scale-invariant frequency-size event distribution (e.g., Chessa et al., 1999;Castellaro and Mulargia, 2002).In addition, variations of the critical behavior, for example, due to stochastic rather than the nearest-neighbor interaction network topology, have also been a focus (Manna, 1991;Chen et al., 2008).The situation is similar for the slider-block models (Rundle and Klein, 1993;Karmakar et al., 2005;Christensen and Olami, 1992).Huang and Turcotte (1988) attributed the variation in the scaling exponent to the fractal distributions of stress and strength, whereas the OFC model implied that the distinctive amount of energy dissipated during a rupture is responsible for the variations in the scaling exponent.
In addition to the numerical modeling efforts, numerous studies have been conducted on complex slip behavior on frictional surfaces in more realistic physical models, either on rock specimens (Dieterich, 1978(Dieterich, , 1979;;Scholz, 1998) or on other soft materials (Yamaguchi et al., 2009(Yamaguchi et al., , 2011)).Yamaguchi et al. (2011) demonstrated that the variation of the scaling exponent is a function of the sliding plate velocity in their physical model.Hallgass et al. (1997) modeled fault friction by using two fractional Brownian profiles that slide over each other and showed that the scaling exponent is a function of the roughness index.Fractal geometry was explicitly introduced in their study, similar to the treatment of Huang and Turcotte (1988).In short, a crucial tendency of studies of these complex systems appears to be the recognition of the distinct variation as a function of the system heterogeneity.
We have briefly reviewed the observational evidence supporting heterogeneous failure threshold fields on realistic faults and the noticeable variation of the G-R scaling exponent, the b value.We have also reviewed both the cellular automata and the physical modeling approaches that typically invoke homogeneous failure threshold fields with different loader plate velocity or implanted fractal heterogeneity to account for the variation of the scaling exponent of the frequency-size power law distribution.We were motivated to investigate the dependence of the variation of the system scaling exponent on the distinct configuration of the toppling threshold heterogeneity.We employed a sandpile model with a highly simplified toppling rule, instead of the complex microphysics of realistic friction and slipping regimes.In other words, we avoided being restricted by our limited knowledge of the complex fault failing mechanisms and concentrated on the feasible macroscale drives associated with heterogeneous system configurations to make it tractable to obtain analogous insights into the reported spatiotemporal variation of the seismic b value.In Sect.2, we briefly introduce the general method of the loading-unloading procedure of the sandpile cellular automaton, and our implementation and comparison of toppling threshold fields with distinct coverage extents of the embedded heterogeneity.In addition, we cross-examine the effects of varying diversity, compactness, and geometry of the heterogeneous threshold configuration in Sect.3. In Sect. 4 we then summarize our findings and raise discussions on the contributions from our study as compared to previous studies.

Method: toppling threshold fields with distinct extents of the embedded heterogeneity
The sandpile cellular automata models are usually carried out on a two-dimensional square lattice.The finite size effect is caused by using small lattice area that leads to bias of the large size events distribution.To avoid the finite size effect as much as possible, all our sandpile models were implemented on a relatively large 256 × 256 lattice.The background toppling threshold of each cell was set to a uniform value of 4, with some of the randomly picked cells being assigned heterogeneous thresholds ranging from 2 to 24 with uniform probability.The coverage, as portions of the entire lattice, of the heterogeneous thresholds was modeled by varying the extent of coverage (e.g., Fig. 1a).We randomly dropped sand grains on the lattice, one at a time, for 10 6 successive iterations.The invoked toppling and relaxation rule is straightforward that when a cell with the accumulated number of sand grains exceed the local threshold, toppling occurs.For the background cells with a normal threshold of 4, each of the closest neighbors to the north, east, south, and west received one grain of sand when a cell toppled.For heterogeneous sites with a threshold that is not a multiple of 4, The scaling exponents (denoted along with the standard errors) of the log-log frequency-size plots of the four cases of distinct heterogeneity coverages (with corresponding color markers) increase systematically for low coverage extent under 30 %, and approach asymptotically toward 1.41 beyond 30 % and up to 99 % (not shown here) coverage (note that the data are simplified by the binning the same as Fig. 1b).
to avoid any directional preference, the redistributed sand grains were transferred to a cell randomly chosen among the closest neighbors and then the subsequent grains were transferred to cells rotating through these closest neighbors.The redistribution can lead to cascaded topples and results in an avalanche event; the size of the avalanche event is determined by the total number of topples or the number of cells involved in an avalanche.In this study, the number of cells involved in an avalanche, denoted by s, was adopted.Sand grains falling out of the lattice along the lattice boundary are permanently lost.The dissipative system reaches a critical state where a quasi-static number of the total sand grains are retained within the lattice.In other words, the number of total sand grains remaining in the lattice fluctuates slightly around an established near-stationary value even when the sand grains are continuously fed into the critical system.Subsequently, every sand grain input can cause avalanche events of all sizes, with the maximum event size being determined by the size of the lattice.We obtained the average of the slightly fluctuating quasi-static number of total grains retained in the lattice when the dissipative system reached criticality.Evidently, the quasi-static number of total grains increased as the extent of the heterogeneity increased, from 1.06 × 10 5 (for 1 %) through 1.44 × 10 5 (10 %) and 2.16 × 10 5 (30 %) up to 2.90 × 10 5 (60 %), because with increasing coverage, there are more cells capable of sustaining more sand grains.To avoid the controversy concerning the actual establishment of the power law, the cumulative distribution function (CDF) of the event frequency was first determined invoking an implementation (Clauset et al., 2009) of the maximum likelihood method (Aki, 1965) to validate the power law scale-invariant frequency-size distribution that is of the form f = a s s −b s (with the frequency, f ; the size, s; and the scaling exponent, b s , where subscript "s" is used to denote the sandpile model).Then, an estimate of the noncumulative scaling exponent, b s , was obtained in the log-log domain by least-squares regression on the linear relationship, log 10 (f ) = log 10 (a s ) − b s log 10 (s).Standard errors were also estimated for the obtained scaling exponents.To further simplify the scattering of the raw data of frequency-size in the large events, we binned the sizes into 1.1 i -1.1 i+1 intervals (Fig. 1b).Representative size of each interval is determined by s i max s i min , whereas the frequency of the ith interval is normalized by the total number of events within that interval divided by (s i max − s i min + 1), where s i max and s i min are the maximum and minimum event sizes of that interval (Christensen and Moloney, 2005).A systematic, statistically significant increase in the scaling exponent (b s = 1.08 ± 0.01 to 1.30 ± 0.03) was observed to accompany an increase in the extent of relatively low degree of heterogeneity below 30 % coverage (Fig. 2).The scaling exponent approached asymptotically toward the value of 1.41 with heterogeneity coverage up to 99 %.It is worth mentioning that the result for the homogeneous case without any heterogeneity is not differentiable with the 1 % heterogeneous case.We used quenched disorder (i.e., the heterogeneity configuration was fixed during the process of modeling in each case) to avoid potential interference that might have obscured the effects of the quantitatively fixed proportion of the heterogeneity.However, dynamic evolution can be induced through invoking the annealed updating while fixing a particular proportion of the heterogeneity.

Effects of varying diversity and compactness
Although it has been shown that a higher extent of the coverage of the heterogeneous threshold systematically increases the scaling exponent of the avalanche event frequency-size power law, all our previously modeled heterogeneous threshold fields were equally diversified within the range of 2 to 24.
To clarify whether the degree of diversity might modify the systematic increase, we conducted a series of simulations in which the diversity was varied for a given coverage extent of the heterogeneous sites.All modeled cases with the same coverage proportion yielded highly similar scaling exponents, irrespective of the imposed diversity.Two extreme cases for the coverage extent of 10 %, with either full diversity through the range of 2 to 24 (Fig. 1a) or a single heterogeneity of 24 (as opposed to the normal background threshold of 4), yielded distinct quasi-static numbers of total grains but almost the same scaling exponents, 1.30 ± 0.03 and 1.29 ± 0.03 (black and red marks in Fig. 3).The result of Nonlin.Processes Geophys., 21, 1185-1193, 2014 www.nonlin-processes-geophys.net/21/1185/2014/ the nondiversified uniform heterogeneity might be related to the extreme value of 24 and transitional behavior might occur if the heterogeneous value is gradually decreased to be close to the background of 4. We thus performed three additional model simulations for the same 10 % of sparse area coverage (similar to the case shown in Fig. 1) but by considering single heterogeneous thresholds of 15, 7, and 5 for each case.The scaling exponents for the five cases with the heterogeneity either diversified uniformly in the range of 2 to 24 (black case in Fig. 3) or set at the value of 24 (red case in Fig. 3), 15, 7, and 5 are all within the range of 1.28 ± 0.03 to 1.31 ± 0.03, demonstrating the independence of the intensity of the heterogeneous thresholds.This reminds us of an earlier set of testing experiments without heterogeneity in which a sandpile model with a uniform threshold of 24 yielded the same scaling exponent as the one from the standard BTW model with the threshold set to 4, indicating similar independence of the intensity of the (homogeneous) thresholds.A heterogeneity field of the uniformly high threshold of 24 might be capable of retaining a larger quasi-static number of grains compared with the diversified heterogeneity fields with generally lower thresholds; however, it is not trivial concerning the reason why their scaling exponents appear the same.The influence of varying diversities has also been examined for various cases with distinct coverage extents and for various different degrees of diversity for a given coverage extent.The systematic trend, as demonstrated in Fig. 3, was the same.In other words, the degree of diversity and the intensity of the heterogeneity do not affect the scaling exponent of the power law event statistics.
Another potentially crucial consideration is that heterogeneous sites in all our previously modeled cases were randomly selected following a uniform distribution for a given coverage extent.In other words, the heterogeneity density in each case was more or less uniform across the study domain.However, observed realistic failure threshold heterogeneities on the fault plane are not uniformly distributed.In another series of models, we therefore varied the compactness of the heterogeneous threshold fields while preserving their coverage extents and degrees of diversity.Examples of the coverage extent of 10 % for the heterogeneous threshold field with a fully diversified range of 2 to 24 were again compared among various compactness configurations.In other words, sites with heterogeneous thresholds can either be sparsely distributed uniformly throughout the entire lattice or be restricted to gradually smaller areas, such as 70, 50, and 25 %, of the lattice until they are squeezed into the 10 % area limit (Fig. 4a).The scaling exponents showed a systematic decrease with an increase in the compactness.A comparison between two extreme cases, compactly configured (Fig. 4a) versus fully sparsely distributed (Fig. 1), showed that the scaling exponents were significantly different, 1.08 ± 0.01 versus 1.30 ± 0.03 (Fig. 3).A nearly identical trend was observed for a distinct set of diversity degrees and geometric configurations (compare Fig. 4a-d).In The log-log frequency-size display of four cases -the black and red cases are sparsely distributed and the blue and green are compactly grouped (Fig. 4a, and b).There are apparently two different groups for the scaling law: the black and red cases with the same sparse 10 % coverage (although the heterogeneous threshold field is diversified for the black case but uniform for the red case; see the text) have similar values of the scaling exponent (approximately 1.30 ± 0.03), whereas the second group with compactly grouped heterogeneity irrespective of the diversity, location, and geometry (the blue and green cases in Fig. 4) all have similar b s values in the range of 1.08 ± 0.01 to 1.10 ± 0.02.summary, a specific compactness yields similar scaling exponents in the range of 1.08 ± 0.02 to 1.10 ± 0.02, irrespective of the degree of diversity, intensity, and geometric configuration of the heterogeneity.Tejedor et al. (2009) reported an investigation of the influence of the geometric configuration of asperities of a fault on the scaling variation by using a modified two-dimensional forest-fire model (2-D FFM); they defined the asperity as the dual of the specifically assigned trigger sites.It is unclear whether the modified 2-D FFM is directly comparable to the modified sandpile model presented in this study and whether the asperity defined in their model is analogous to the heterogeneous thresholds discussed here.However, we noticed that first, unlike the 2-D FFM, there are no subcritical, critical, and supercritical transitions in our models.Second, their type 5 configuration with dispersed asperity appears to be considerably close to our fully sparsely distributed heterogeneity with low coverage extent; moreover, their configuration yields a relatively low scaling exponent that is consistent with our finding.However, it is also noteworthy that, for their type 5 configuration, the power law exhibits an apparent subcritical behavior on the square 50 × 50 lattice.The scaling exponent for their type 4 configuration, with the asperity covering the complete lattice except a central relatively small hole, conforms to the next to 24, are compactly grouped in a small area (the blue case).(b) Heterogeneous threshold fields with the constant value of 24 are similarly compactly grouped (the green case).(c) Diversified threshold heterogeneity that is compactly grouped in a small area (as in a) but is moved to the center (as opposed to the edge) of the board (the magenta case).(d) Single threshold heterogeneity (as in b) that is also compactly grouped but with elongated geometry, resulting in a change in the aspect ratio (the cyan case).All these cases yield a similar scaling exponent of around 1.09 -the magenta (c) and cyan (d) cases had scaling exponents of 1.08 ± 0.02 and 1.10 ± 0.02, and were not shown in Fig. 3 that is different from the value of around 1.30 for sparsely distributed heterogeneities (Fig. 3).
smaller scaling exponent with apparent supercritical behavior on the square 50 × 50 lattice.Although our models did not show either super-or subcritical behavior, if we treat the central hole as a compactly grouped heterogeneity embedded in the background of what they defined as the asperity, then the scaling exponent would be expected to be slightly larger than the previous dispersed asperity but smaller than other configurations, which is consistent with their reports.Overall, their models clearly demonstrated dependency of scaling on the aspect of the fault, which was not considered in our model, and the geometric configuration of the asperity.It is inconclusive concerning whether the effects of coverage extent and compactness of the heterogeneous failure threshold fields in the current study are inconsistent with their findings.The most critical and difficult task that is yet to be completed would be the reasoning and justification of the type of cellular automaton most appropriate for the study of the statistics of seismogenesis.

Conclusions and discussions
We conducted a series of sandpile model simulations using various configurations of distinct coverage proportions and compactness for the heterogeneous toppling threshold field.It is concluded that the systematic variation of the scaling exponent is controlled by the coverage extent of the embedded heterogeneity rather than the intensity level or the contrast range of the threshold (Fig. 2).Further dependency of the scaling exponent variation on the compactness of the heterogeneity was also found (Fig. 3).Karmakar et al. (2005) postulated that the variation of criticality can be attributed to the absence of the toppling balance between all the relaxed grains of a toppled site and all the received grains when all the neighbors topple once.This might appear to be an unspecific criterion because all models with a stochastic toppling rule (e.g., Manna, 1991;Chen et al., 2008) are likely to show absence of that balance.However, in our study, the hypothesis is confirmed in a more systematic way.In other words, in addition to the observation that all our experiments with noticeable scaling variations are characterized by the absence of toppling balance as opposed to the presence of it in the reference BTW deterministic model without heterogeneity, we also found that the greater the number of sites lacking balance, which might be caused by an increase in the threshold extent, the more pronounced the distinct criticality.However, it cannot be the only reason for the systematic trend of the criticality variation.This can be inferred from the result that, although a given coverage proportion and range of contrast of the heterogeneity correspond to a relatively high degree of toppling balance, varying the compactness configuration of the heterogeneity also leads to a systematic variation trend of the scaling exponent.In contrast, varying the diversity degree that perturbs the extent of the toppling balance exerts negligible effects on the scaling variation (Fig. 3).Christensen and Olami (1992) and Olami et al. (1992) indicated that the redistribution of force as a function of the elastic constants in the slider-block model is intrinsically nonconservative locally.In other words, the unloaded amount of force from a slipping cell is not the same as the sum of all the amounts received by its neighbors.They attributed the variation of the scaling exponent to the level of the conservation by adjusting the ratio of the stiffness of the driving spring of the loader plate to the inter-block spring.However, all our models with distinct heterogeneity configurations preserve the local conservation.In other words, unlike the result of Christensen and Olami (1992), we demonstrate that variations of the scaling exponent revealed from our modeling are caused by the distinct configurations of the heterogeneous failure threshold fields rather than the level of the local conservation.By contrast, Rundle and Klein (1993), who pioneered the study of heterogeneous failure threshold of the slider-block-type cellular automaton, indicated that the criticality variation is strongly influenced by the loading, the form of failure threshold, and the mechanism of energy dissipation.Tejedor et al. (2009) also reported the influence of the aspect ratio of the fault plane and the geometric configuration of asperity.We agree with their profound insight that a spatially varying failure threshold is a crucial factor in the development of a realistic model of earthquake rupture along fault surfaces.However, we did not observe either the phenomena of in and out of criticality or transitions among subcritical, critical, and supercritical states in any of our implementations of a heterogeneous toppling threshold.In the sandpile model, the loading is performed by dropping sand grains one at a time until the cascaded toppling rests completely, and therefore there is no comparable appraisal of the strong dependence of the scaling variation on the loading velocity in the slider-block type models.We believe that the regional kinematic loading (e.g., plate convergence) is incredibly small in the timescale of an individual earthquake event, and therefore the analogous loading velocity in the slider-block model should likely be maintained at a relatively low speed.Furthermore, unlike the study of Rundle and Klein (1993), which examined the dependency on the loading velocity based on a random or fractal pattern of a full heterogeneity field, our models considered only the coverage extent and compactness without including other specific configuration patterns.In other words, the revealed systematic trends of criticality variations are not induced by the implanted pattern of the heterogeneity.We postulate that the extent and compactness of the heterogeneous failure thresholds or asperities on the fault surface, as well as the dynamic connection topology of the interaction network before, during, and after a rupture process (Chen et al., 2008;Chiao, 2012), rather than the magnitude or diversity degree of the failure threshold, are plausible dynamic parameters related to the maturity of realistic faults.
Heterogeneous material and rheological properties, distinct fault surface maturity, roughness and friction regimes, and the three-dimensional porosity-permeability structure that controls the effect of the potential hydraulic behavior on the failure threshold might all be factors promoting failure threshold heterogeneity within a complex fault system.Variation in the healing rate is another crucial plausible factor worth mentioning.Depending on the pressure, temperature, and mineralogy, distinct degrees of healing might induce complex interactions between rapid, unstable failures of fault patches and stable, slow creeps of the surrounding region (e.g., Chen and Lapusta, 2009;McLaskey et al., 2012).It is speculated that a tectonic province under a tensional regime endures accumulated historical ruptures on mature faults that are likely to be phyllosilicate-rich, and thus lacks effective healing mechanisms (e.g., Tesei et al., 2012).A higher extent of threshold heterogeneity might manifest and, consequently, a higher b value might be observed.On the other hand, the healing through time of a region under high compression tends to erase and lower the extent of lateral heterogeneity, thereby leading to a lower b value.Furthermore, the frequently observed increase in the b value after major earthquakes (Chiao, 2012) might at least be partially attributable to the fragmentation of the compactly configured asperities.
Figure 1.(a)Example of a case of diversified heterogeneity with toppling thresholds randomly sampled in the range of the value 2 to 24 (color bar); the heterogeneity covers 10 % of the area of the 2-D model lattice with 256 × 256 cells, and it is embedded in a background field of threshold 4 (medium blue), which occupies the remaining 90 % of the lattice area.(b) The power law frequency-size distribution yields an scaling exponent of b s = 1.30 ± 0.03 (denoted along with the standard error).Note that the inevitable scattering across the large size events (marked by black dots) is considerably simplified by regrouping the sizes into bins of 1.1 i -1.1 i+1 size interval (blue circles; see the text).
Figure2.The scaling exponents (denoted along with the standard errors) of the log-log frequency-size plots of the four cases of distinct heterogeneity coverages (with corresponding color markers) increase systematically for low coverage extent under 30 %, and approach asymptotically toward 1.41 beyond 30 % and up to 99 % (not shown here) coverage (note that the data are simplified by the binning the same as Fig.1b).
Figure3.The log-log frequency-size display of four cases -the black and red cases are sparsely distributed and the blue and green are compactly grouped (Fig.4a, and b).There are apparently two different groups for the scaling law: the black and red cases with the same sparse 10 % coverage (although the heterogeneous threshold field is diversified for the black case but uniform for the red case; see the text) have similar values of the scaling exponent (approximately 1.30 ± 0.03), whereas the second group with compactly grouped heterogeneity irrespective of the diversity, location, and geometry (the blue and green cases in Fig.4) all have similar b s values in the range of 1.08 ± 0.01 to 1.10 ± 0.02.

Figure 4 .
Figure 4. Heterogeneous toppling thresholds (2-24) for the 10 % extension and compactness, but with distinct diversity or geometry, embedded in a background threshold field of 4. (a) Diversified heterogeneous thresholds, with values randomly sampled within the full range of 2 to 24, are compactly grouped in a small area (the blue case).(b) Heterogeneous threshold fields with the constant value of 24 are similarly compactly grouped (the green case).(c) Diversified threshold heterogeneity that is compactly grouped in a small area (as in a) but is moved to the center (as opposed to the edge) of the board (the magenta case).(d) Single threshold heterogeneity (as in b) that is also compactly grouped but with elongated geometry, resulting in a change in the aspect ratio (the cyan case).All these cases yield a similar scaling exponent of around 1.09 -the magenta (c) and cyan (d) cases had scaling exponents of 1.08 ± 0.02 and 1.10 ± 0.02, and were not shown in Fig.3that is different from the value of around 1.30 for sparsely distributed heterogeneities (Fig.3).