Spatial-temporal Distribution Characteristics of Global Seismic Clusters and Associated Spatial Factors

Earthquakes exhibit clear clustering on the earth. It is important to explore the spatial-temporal characteristics of seismicity clusters and their spatial heterogeneity. We analyze effects of plate space, tectonic style, and their interaction on characteristic of cluster. Based on data of earthquakes not less than moment magnitude (MW) 5.6 from 1960 to 2014, this study used the spatial-temporal scan method to identify earthquake clusters. The results indicate that seismic spatial-temporal clusters can be classified into two types based on duration: persistent clusters and burst clusters. Finally, we analysed the spatial heterogeneity of the two types. The main conclusions are as follows: 1) Ninety percent of the persistent clusters last for 22–38 yr and show a high clustering likelihood; ninety percent of the burst clusters last for 1–1.78 yr and show a high relative risk. 2) The persistent clusters are mainly distributed in interplate zones, especially along the western margin of the Pacific Ocean. The burst clusters are distributed in both intraplate and interplate zones, slightly concentrated in the India-Eurasia interaction zone. 3) For the persistent type, plate interaction plays an important role in the distribution of the clusters’ likelihood and relative risk. In addition, the tectonic style further enhances the spatial heterogeneity. 4) For the burst type, neither plate activity nor tectonic style has an obvious effect on the distribution of the clusters’ likelihood and relative risk. Nevertheless, interaction between these two spatial factors enhances the spatial heterogeneity, especially in terms of relative risk.


Introduction
Numerous observations and research studies have suggested earthquakes show obvious clusters on the earth (Rehman et al., 2014;Liu and Stein, 2016). The clusters are heterogeneous in terms of spatial-temporal distribution. From a spatial perspective, cluster size in the African-Arabian rift decreases northward (Hall et al., 2018), and clustering of large earthquakes is evident along the Sumatra, Kuril, and Tonga subduction zones (Lay, 2015). From a temporal perspective, earthquakes are more likely to cluster around a large one (Kagan and Jackson, 2000;Jiang and Wu, 2005). Lower magnitude events commonly present more complex clustering mechanisms (Parsons and Geist, 2014). In addition, spatial-temporal clustering of earthquakes has also attracted research attention (Mukhopadhyaym et al., 2011).
Clusters could be classified into different types, and features of these types correspond to different physical characteristics. For example, seismic clusters of southern California is classified into single burst-like and swarm-like sequences (Zaliapin and Ben-Zion, 2013a), these types were likely associated with different failure process (Zaliapin and Ben-Zion, 2013b). Clusters in northern California show longer duration but fewer events per cluster relative to those in southern California, which could be due to different faults or the influence of the 2004 Parkfield earthquake (Wu et al., 2015). Analysis on seismic cluster and associate factors could promote accuracy of seismic hazard assessment (SHA). Firstly, seismicity is one of important parts in SHA. Recently, the most popular method is Probabilistic Seismic Hazard Assessments (PSHA). In this method, main indexes to express seismicity are probability distribution of magnitude and earthquake occurrence rate (Main, 1996). However, the indexes still cannot perfectly assess seismicity (Mulargia et al., 2017). Clusters could be distinguished in the temporal evolution of seismicity (Hall et al., 2018;Wang et al., 2018). Then their characteristic would reflect seismic activity. Considering the cluster characteristic of earthquake, we use spatial-temporal scan method to analyse the spatial-temporal characteristics of global seismic energy clusters in this paper. The characteristics could be auxiliary to PSHA hopefully. This paper is based on earthquakes data between 1960 and 2014. Clusters could be divided into two types by duration. To detail their seismicity, we analyse their characteristics respectively.
Secondly, analysis on associate factors is also necessary in SHA. It not only could be conducive to understand seismicity process; but also could play as auxiliary parameters in SHA, and then revises the assessment result. Existing studies still considered the effect of a single spatial factor on heterogeneity, including heat flow, geological activity, and tectonic style (Zaliapin and Ben-Zion, 2016;Jackson, 2010, 2012;Berryman et al, 2012;Faenza et al., 2008). However, multiple factors integration effect should also be considered (Gao, et al., 2011). Earthquake as a complex event (Sharma et al, 2013;Shen et al., 2018), its distribution may be influenced by more than one factor (Stadler et al, 2010;Lay, 2015). Multi-factor interaction might have a more significant influence on the cluster's feature. In this paper, we use geographical detector method to analyse cluster's response on plate space, tectonic style and their interaction respectively.

Data and Method
Earthquake data were from International Seismic Center (http://www.isc.ac.uk/iscgem/; Storchak, 2013). To ensure that the data was completely recorded, this experiment extracted earthquakes not less than moment magnitude (M w ) 5.6 from all data records between 1960 to 2014, totally 16 905 events. Gutenberg-Richter law is used for completeness test (Schorlemmer et al., 2018). It is based on suppose that relationship between M w and seismic count above the M w value would present power-law if data is complete (Main, 1996;Faenza et al., 2008). Observed deviations from the distribution within the lower magnitude range indicate a loss in completeness. According to this law, the frequencymagnitude distribution (log linear plot) of all global earthquakes records from 1960 to 2014 is shown in Fig. 1, in which the x axis is M w and the y axis is the proportion of events above an M w value in all the data, namely, the complementary cumulative distribution function (CCDF). The graph shows the minimum magnitude of completeness M c = M w 5.6, that is, earthquakes not less than M w 5.6 are completely recorded in the ISC catalogue and can be applied to statistical analysis.
The spatial-temporal scan method (Kulldorff et al., 2005; http://www.satscan.org/) is used to detect seismic energy clusters and to measure the cluster likelihood and relative risk. Based on grid point data, it utilizes thousands or millions of overlapping cylinders to define the scanning window, each of which is a possible energy cluster candidate. The steps of this process are as follows: 1) considering the research area as a spatial-emporal cube, the cube base represents the whole geographical area, and height is the research period. 2) A cylinder is used to detect a cluster, and its base and height are the spatial scope and temporal range, respectively. Thus, to detect the entire area, the base's center dynamically changes. 3) Then, the base and height of the cylinder are continuously increased during the detection process until a threshold is reached. 4) All possible energy clusters are detected, and a Monte Carlo simulation is used to evaluate the statistical significance of the cluster likelihood. 5) Finally, the relative risk of each cluster is calculated.
We utilizes Poisson log-likelihood ratio (LLR) to measure the likelihood of clustering in this study: Assuming seismic energy exhibits a Poisson distribution, c is the actual total seismic energy in the cylinder, μ is the expected energy in the cylinder based on the Poisson distribution, and C is the total energy worldwide. The method calculates the LLR of each cylinder. A larger LLR value means cylinder's energy presents higher probability of clustering than other cylinders. Relative risk (RR) evaluates the unexpected degree of actual seismic energy. It is computed as a ratio of actual energy to expected energy inside the cluster versus outside, as shown in Equation (2).
The parameters in the formula are the same as those in Equation (1). A larger RR value means that the cylinder's energy is higher than the expected degree, and its risk is higher than that of other cylinders.
In this article, the geographical detector method (Wang et al., 2010) is used to measure the spatial heterogeneity of the LLR and RR for both a single spatial factor and multiple spatial factors. This method is based on the hypothesis that if the spatial heterogeneity of the variable is dependent on the space, then their distributions would be similar. Wang et al. (2010) defined the operation that measures the spatial heterogeneity in intersection space as the interaction detector. The interaction detector estimates the interaction between two spatial factors by comparing the spatial heterogeneity in a single space with the intersection space. Let the space variable be X and the dependent variable be Y. Overlay-ing X and Y, X subdivides Y according to zones of X. In this paper, X represents the spatial factors, and Y represents the LLR or RR. According to plate space, we divide the clusters into the zones according to their centre locations. Thus, all the clusters are separated into 22 regions. According to tectonic style, all the clusters are separated into four groups according to their centre's distance from a tectonic boundary. Of these, three of the groups correspond to three tectonic styles, and the fourth does not belong to any style. After numerous tests, when a cluster is approximately 10º (Geographic Coordinate System) away from a tectonic boundary, the distribution of seismic events becomes discrete instead of strip-like. Therefore, 10º is set as the threshold for classification. If the distance is greater than 10º, the cluster is classified as the fourth style; otherwise, it is classified as the nearest tectonic style.
To assess the spatial dependence of X and Y, the geographical detector method compares the variance in Y before and after division. If X completely determines Y, then Y's variance decreases markedly after division and approaches zero within all sub-parts of Y. If X does not determine Y, then Y's variance does not differ between before and after division. The power of determinant (q) is shown in Equation (3).
where 2 y  is the variance within sub-part y, N y is the cluster count in zone y,  2 is the global variance of Y, and N is the number of all clusters worldwide.
To test whether the spatial heterogeneity of Y increases at the intersection space, the interaction detector compares values of q x1 ∩ x2 with values of q x1 and q x2 . The term q x1 ∩ x2 indicates the determinant of the intersection space by overlapping plate space (x1) and tectonic space (x2). If q x1 ∩ x2 > max(q x1 , q x2 ), the spatial factors enhance each other; if q x1 ∩ x2 = (q x1 +q x2 ), the two spatial factors are independent of each other; and if q x1 ∩ x2 > (q x1 +q x2 ), the factors nonlinearly enhance each other (Luo et al., 2015).
Technology roadmap in this paper is shown in Fig. 2. 1) Dividing the spatial zones, including plate space and tectonic style. The plate space contains interplate zones and intraplate zones, which are generated based onseismic events. Tectonic style contain trench, ridge and transform margins. 2) Using the spatial-temporal scan method to detect global seismic energy clusters. The characteristics of each cluster are calculated, including duration, cluster likelihood, and relative risk. 3) Dividing clusters into different types by duration. 4) Using the geographical detector method to explore spatial heterogeneity of the two types.

Geographical Spaces Related to the Distribution of Seismic Clusters
Different plates move with different speeds and directions (Copley et al., 2010;Jagoutz et. al., 2015). These differences may cause cluster distribution differences at interplate boundaries or in intraplate zones, and the long-term movement may result in different tectonic styles, which may also influence the spatial heterogeneity of the clusters. Fig. 3 shows spaces constructed by interaction between plates and tectonic styles respectively. First,  based on the 16 905 earthquake events, we divided global space into interplate and intraplate spaces, totalling 22 spaces that are collectively referred to as plate space in the text, shown with Arabic numerals in the figure. Seven plates were considered in the division: Eurasia, Pacific Ocean, India, Africa, America, Antarctica and Nazca.
The division was based on the idea that global earthquakes tend to occur along plate boundaries (Scholz, 2002) and became dispersed with increasing distance from the boundary. Furthermore, the earthquake density differs among each interplate space. Thus, plate space can be divided according to seismic distribution, but a uniform buffer width is not applicable. This article used the following steps to define the space: 1) the seismic density ρ was calculated for a 0.5º × 0.5º global grid. 2) After several trials, spaces with ρ > 0.2 were extracted. We tested several extraction results with density thresholds ρ = [0.1, 0.2, 0.3], as shown in Fig. 4. When ρ > 0.2, the interplate spaces separated well with intraplate spaces. 3) Based on empirical knowledge, zones that belong to the same plate boundary are connected, and regions are separated into independent spaces when divided by a plate interaction zone.
Then, tectonic style was the other spatial aspect considering in this stuty. The data was from University of Texas Institute for Geophysics (UTIG) PLATES program (http://ig.utexas.edu/marine-and-tectonics/plates-project). It contains trench, ridge and transform, shown with coloured lines in Fig. 3. Distribution of tectonic style is similar with the plate space. But levels of detail are different, for example, plate spaces No. 8 and No. 11 are different plate spaces but all belong to ridge in Fig. 3. However, plate space No. 1 contains 3 kinds of tectonic styles.

Spatial-temporal Clusters of Global Seismic Energy and Their Features
We detect seismic energy clusters based on global 0.5º × 0.5º grid data in this study. The seismic energy in each grid is the total energy of the earthquakes located in the grid. Itis calculated according to seismic magnitude (Hanks and Kanamori, 1979;Hall et al, 2018). The relationship is shown in Equation (4): where E is seismic energy and M w is moment magnitude. When implementing the spatial-temporal scan, this article considers the maximum spatial-temporal scope of the scan cylinder, the expected energy in each grid, and the statistical significance of cluster likelihood. Firstly, after numerous tests, this work set the maximum scan radius to 1000 km and the longest duration to 70% of the entire study period. Under these conditions, the detection results gradually stabilize. Secondly, due to stable geological conditions, it is almost impossible to have an earthquake in certain regions. This paper only considers grids in which an earthquake has occurred. Under the Poisson assumption, the expected energy in each grid is the average of the global energy on these grids.
Lastly, after detection, the experiment only considers clusters with cluster likelihoods with high statistical significance, P ≤ 0.01.
Finally, the paper retains a total of 390 clusters. They are shown in Fig. 5. Each of them has LLR and RR attributes, which could reflect the physical processes of the earth.

Different cluster types based on duration and their characteristics
From the perspective of the duration of seismic energy, all clusters are separated into two types: persistent Fig. 5 Spatial-temporal clusters of global seismic energy. Circles and columns represent the spatial scope and temporal information, respectively. Column height is used to express duration; column colour corresponds to the temporal evolution of a cluster. Therefore, a taller column corresponds to a longer duration, and blue and red correspond to earlier and later periods, respectively clusters and burst clusters. The cluster characteristics, including duration, LLR and RR, are shown in Fig. 6. Every point denotes a cluster, and the x and y axes are the start time and end time of the cluster, respectively. Thus, a point far from the diagonal has a longer duration than a point close to the diagonal. The duration of a point on the diagonal is 1 yr. According to the duration data, the clusters show either long durations or short durations. The clusters group together at approximately 30 yr and approximately 1 yr. Taking 10 yr as a critical value here, the clusters are classified as persistent clusters or burst clusters. There are 21 persistent clusters, and 369 burst clusters.
After analysing the LLR and RR values for each type, the results show that seismic energy is more likely to gather in persistent clusters, but burst clusters have greater relative risk. In Fig. 6, the size of a point represents the cluster likelihood. Because the point size of persistent clusters is generally larger than that of burst clusters, the figure indicates that persistent clusters exhibit greater seismic energy. The different point colours represent different RR values. Because the RR values of the burst clusters are generally higher than those of the persistent clusters, burst clusters have a greater relative risk. Further examining the duration distribution of the two types, the experiment results show that usually the persistent clusters last for 22−38 yr, and the burst clusters last for 1.00−1.78 yr, as shown in Fig. 7. Fig. 7a shows the duration distribution of the persistent clusters. The cluster duration ranges from 10−38 yr, and 90% of the clusters located in the range 22−38 yr. Fig. 7b shows the duration distribution of the burst clusters. The cluster duration ranges from 1 to 8 yr, and 90% of the clusterslocated in the range 1−1.78 yr. Fig. 6 Duration, poisson log-likelihood ratio (LLR) and relative risk (RR) of the clusters.In this paper, clusters with duration longer than 10 yr are persistent type; and clusters with duration shorter than 10 yr are burst type. The size of a point represents the cluster likelihood, with larger sizes representing higher likelihoods; the different point colours represent different RR values, with the transition from blue to red representing low to high RR values Fig. 7 Distribution of duration for each type. The blue line corresponds to cluster count for each duration. The red points represent measurements of cumulative distribution function (CDF) and duration, and the red line is a curve fitted to this measurements. The cross point of dashed represents upper (lower) duration 90% of the clusters located. The shaded part indicates the distribution of 90% of the persistent clusters

Spatial heterogeneity of clusters in each type
The spatial distribution of the persistent clusters is shown in Fig. 8. The persistent clusters tend to occur in interplate zones, mainly along the western margin of the Pacific Ocean, 30% of the clusters located in this area (zones No. 1 and No. 2). The area is characterized by a trench. The northern and southern plate boundaries of India (zones No. 6 and No. 9), the Nazca boundary (zones No. 4,No. 13 and No. 15), and southern boundary of the American plate (zone No. 14) each contain 10% of the clusters. In addition, there are no persistent clusters along the boundaries of the southern Pacific Ocean (zone No. 5), northern and southern African plate (zones No. 7 and No. 12), northern Eurasian plate (zone No. 8), and the intraplate zones .
The distribution of the burst type is comparatively dispersed. Burst clusters are distributed in both intraplate and interplate areas, as shown in Fig. 9. In terms of interplate boundaries, the clusters are somewhat concentrated along the India-Eurasia boundary, which is characterized by subduction. Compared with the persistent clusters, colour of space No. 6 (India-Eurasia boundary) is darker, corresponding to an increase in the cluster ratio from 10% to 21%. Similarly, the colours of zones No. 1 and No. 2 (western boundary of the Pacific Ocean) are lighter, corresponding to a decrease in the cluster ratio from 30% to 11%. In the intraplate regions, except zones No. 20 and No. 22 (Antarctic and Nazca), clusters are distributed evenly in the intraplate zones, all with ratios of approximately 4%.

Spatial heterogeneity in the LLR and RR of the persistent type
Based on the geographical detector method, the spatial distribution and heterogeneity of the LLR and RR in the persistent clusters is shown in Fig. 10. For the LLR of the persistent type, the spatial heterogeneity is more obvious in plate space than in tectonic style: seismic energy is more likely to cluster in the western Pacific Ocean. The spatial heterogeneity is enhanced in the intersection space, and seismic energy is more likely to cluster in the intersection space between the western Pacific Ocean and the trench. In plate space, the spatial heterogeneity is 0.8442, as shown in Fig. 10a. The LLR values are significantly separated. They are high overall in plate space No. 1 (the Pacific Ocean-India boundary), Fig. 8 Spatial distribution of persistent clusters. Every cluster center is marked with a red triangle and serial number C i (i = 1, 2,…, 21). Each zone is marked with serial number j (1, 2,…, 22). The colour gradient of the interplate zone corresponds to the percentage of clusters located in that zone relative to the total number of persistent clusters.The columnar graph shows accurate percentage values of cluster count in each interplate and intraplate space, for example, there are 15% clusters located in zone No. 1 Fig. 9 Spatial distribution of burst clusters. Each zone is marked with serial number j (1, 2, …, 22). The colour gradient of the interplate zone corresponds to the percentage of clusters located in that zone relative to the total number of burst clusters. The columnar graph shows accurate percentage values of cluster count in each interplate and intraplate space followed by plate space No. 2 (the Pacific Ocean-Eurasia boundary). In terms of tectonic style, the spatial heterogeneity decreases to 0.3263, as shown in Fig. 10b. This value overlaps the LLR values in other zones. In addition, the values in the same zone show discrete distributions. Therefore, the LLR values do not exhibit sufficient separation in terms of tectonic style. Finally, in intersection space, the spatial heterogeneity is 0.8821, as shown in Fig. 10c. This value is larger than the value of either space but lower than the sum of both spaces. Accordingly, the difference is enhanced compared with that of a single space. The LLR values are high overall in the zone 2a (Pacific Ocean-India boundary and trench) and in the zone 1a (Pacific Ocean-Eurasia boundary and trench).
The spatial heterogeneity in the RR of the persistent type is more obvious in plate space than in tectonic style, and the relative risk is high along the Eurasia-India boundary. The spatial heterogeneity is enhanced in the intersection space, and the relative risk is high in the intersection space between the Eurasia-India boundary and trench style. In plate space, the spatial heterogeneity is 0.5416, as shown in Fig. 10d. The RR values are high overall in plate space No. 6 (Eurasia-India boundary). In tectonic style, the spatial heterogeneity is 0.1033, as shown in Fig. 10e. This value overlaps the RR values in other zones. In addition, the values in the same zone show discrete distribution. Therefore, the RR values do not exhibit sufficient separation in tectonic style. Finally, in intersection space, the spatial heterogeneity is 0.5661, as shown in Fig. 10f. This value is larger than the value of either space but is lower than the sum of both spaces. Accordingly, the difference is greater than that associated with a single space. The RR values are high overall in zone No. 6a (Eurasia-India boundary and trench).

Spatial heterogeneity in the LLR and RR of the burst type
Based on the geographical detector method, the spatial distributions of the LLR and RR values for the burst clusters are shown in Table 1. The LLR values for the burst type exhibit no obvious spatial heterogeneity in plate space or in tectonic style. The spatial heterogeneity remains weak in the intersection space but is higher. Inplate space, the spatial heterogeneity is 0.1027 and is very weak compared to that of the persistent clusters. In tectonic style, the spatial heterogeneity is 0.0605 and is also weak. Therefore, the LLR values do not exhibit sufficient separation either in plate space or tectonic style. Finally, in intersection space, the spatial heterogeneity is0.1334. This value is still low but has increased; it is Spatial distributions and heterogeneity of poisson log-likelihood ratio (LLR) and relative risk (RR) values for the persistent clusters. The x axis is geographical space, and the y axis is value of LLR or RR. The numbers 1-15 correspond to the different plate spaces, and trench, ridge and transform correspond to the tectonic styles. In intersection space, an x-tick label indicates a combination of two kinds of spaces. Arabic numerals correspond to the plate spaces, and a, b, and c correspond to trench, ridge, and transform, respectively. For example, 1a represents the interaction of plate space No. 1 and a trench. The term q is the degree of spatial heterogeneity larger than that of either space but lower than the sum of both spaces. Accordingly, the difference is enhanced compared with that of a single space.
The RR values for the burst type also exhibit no obvious spatial heterogeneity either in plate space or in tectonic style. The spatial heterogeneity remains weak in the intersection space but has been improved In plate space, the spatial heterogeneity is 0.1564. This value is very weak compared to that of the persistent clusters. In tectonic style, the spatial heterogeneity is 0.0517. This value is also weak. Therefore, the RR values do not exhibit sufficient separation either in plate space or tectonic style. Finally, in intersection space, the spatial het-erogeneity is 0.2599. This value is still low but has increased, and it is larger than the sum value of both spaces. Accordingly, the difference is markedly enhanced compared with that of a single space.

Discussion
The clusters are classified as persistent and burst. In the examined dataset, 90% of the persistent clusters last for 22-38 yr, and 90% of the burst clusters last for 1.00-1.78 yr. The persistent clusters are associated with high clustering likelihood than the burst clusters, whereas the burst clusters have a higher relative risk than the persistent clusters. Active geological conditions facilitate the seismic energy clustering. Compared with the burst clusters, persistent clusters are located in geologically active regions. As shown in Fig. 8, persistent clusters are distributed along plate boundaries, including western margin of the Pacific Ocean, northern Nazca, India, and South America. The Pacific Ocean plate is an especially active plate (Stadler et al., 2010;Loveless et al., 2010). It is moving northwestward continuously, converging with the Eurasian and Indian plates. Therefore, earthquakes are frequent in this region and exhibit a high likelihood of energy clustering. For example, the clusters (C1, C17 and C18) associated with Tonga, Kuril Islands andSolomon Islands respectively all last for 30 yr or more, as shown in Fig. 8. The northern Nazca plate interacts with the Cocos and Pacific Ocean plates. Therefore, the geological conditions of this region are highly active and complex (Zhang et al, 2017). For example, clusters that C3, C4, C7, and C8 all last for 22-38 yr. The Antarctic subduction region is also active (Leat et al, 2018). In this region, cluster C9 and C10 located in South Sandwich Trench, last for 27 and 37 yr, respectively. In terms of the RR, relative to the persistent clusters, the burst clusters feature less predictable and more intense energy. Thus, the RR of the burst type is higher.
The spatial distribution of energy clustering shows that the persistent clusters are concentrated in interplate zones, mainly along the western margin of the Pacific Ocean, which is characterized by a trench. The burst clusters are dispersed across the intraplate and the interplate regions. Along interplate boundaries, clusters are slightly concentrated in the India-Eurasia interaction zone, which is characterized by a trench-type margin. In the intraplate regions, except for Antarctica and Nazca, the clusters are distributed uniformly in other intraplate zones. Compared with the interplate earthquakes, the intraplate earthquakes have longer recurrence times (Liu and Stein, 2016). The western Pacific Ocean and the India-Eurasia interaction zone are both characterized by trench-type boundaries but show different proportions of cluster types. This difference might be because the consistent uplift of the Himalaya-Tibet Mountains hinders the northern movement of the Indian plate (Copley et al., 2010;Iaffaldano et al., 2011), further weakening the India-Eurasia interaction. Finally, most clusters in this zone present short durations. Clusters are scarce in the intraplate regions of Antarctica and Nazca. This scarcity could be the result of their stable geological conditions or a biased spatial distribution. In addition, geological activity is almost the same in other intraplate zones.
For persistent clusters, plate interaction plays an important role in their distributions of cluster likelihood and relative risk. The combination of plate interaction and tectonic style further enhances the spatial heteroge-neity. Seismic energy is more likely to cluster in intersection space between the India-Pacific Ocean zone and the trench style, and clusters have higher relative risk in the intersection space between the Eurasian-India zone and the trench style. For cluster likelihood, on the one hand, spatial heterogeneity is more obvious in plate space than for tectonic style. This difference may be because a tectonic style may include multiple plate zones. If the experiment only considers tectonic style, then the spatial heterogeneity of the LLR is poorly reflected. For example, zone No. 6 and No. 13 in Fig. 8 are both mainly associated with the trench style, but their LLR values are different. Therefore, greater difference is present in plate space. On the other hand, spatial heterogeneity is enhanced in the intersection space. This finding illustrates that the spatial heterogeneity could be related to the intersection space between tectonic style and plate space. For relative risk, the spatial heterogeneity mainly depends on plate space. However, with the weak dependence on tectonic style, the spatial heterogeneity is enhanced in the intersection space.
For burst clusters, no spatial factor has an obvious effect on the distribution of clustering likelihood and relative risk. Nevertheless, the interaction between plate activity and tectonic style enhances the spatial heterogeneity, especially in terms of relative risk. The complex distribution shows that, in addition to the two spaces considered in the paper, the spatial heterogeneity may be related to other spatial factors. Thus, it is necessary to explore further. In addition, the heterogeneity is enhanced in intersection space, and RR in particular is enhanced nonlinearly. This finding indicates that the spatial heterogeneity may be result of interaction among multiple spatial factors. Further research should consider the interaction among more spaces.

Conclusions
This article explores spatial-temporal characteristics of seismicity clusters and their spatial heterogeneity. In addition to a single spatial factor, we considered the effects of multiple spatial factors. Main conclusions are as follows: The duration of seismic energy clusters worldwide shows obvious differences. Different durations are also associated with different cluster likelihood and relative risk values. In addition, spatial heterogeneity is more obvious in the persistent type. This pattern is closely related to plate interaction, although tectonic style plays an additional role. Therefore, the plate movement would be conductive to explain the clustering mechanism. The spatial heterogeneity is relatively weak in the burst type but also shows enhancement in the interaction spatial factor, especially in terms of relative risk. Therefore, future research should analyse more spaces to discover the spatial distributions of burst clusters. Moreover, in addition to a single spatial factor, research should pay more attention to effect of interaction among multiple spatial factors.
The conclusions could have some inspirations for earthquake hazard research. Firstly, different evaluation programmes should be developed for different clusters and zones. Cluster characteristics would provide positive information hazard evaluation. From the perspective of duration and LLR, persistent cluster would have higher hazard than the burst; from the perspective of RR, the burst cluster would be more dangerous. According to spatial distribution of cluster characteristics, zones would be developed corresponding evaluation measures. Then difference of influence factors on cluster characteristics also brings insight in hazard evaluation. Interaction between plates and tectonic style plays more important role on persistent cluster than the burst.