Multifractal Characteristics of Seismogenic Systems and b Values in the Taiwan Seismic Region

Seismically active fault zones are complex natural systems and they exhibit multifractal correlation between earthquakes in space and time. In this paper, the seismicity of the Taiwan seismic region was studied through the multifractal characteristics of the spatial-temporal distribution of earthquakes from 1st January 1995 to 1st January 2019. We quantified the multifractal characteristics of Taiwan at different scales and defined them as ΔD values. Furthermore, we studied the relationship between the ΔD and b values, which signifies the average size distribution of those earthquakes. The results are as follows. (1) The temporal multifractal curve changes substantially before and after the strong earthquakes. (2) The maximum ΔD value of the seismic region in Taiwan occurs at depths of 0~9 km, indicating that geological structures and focal mechanisms is the most complex at these depths compared with other depths. (3) ΔD values for different regions range from 0.2~1.5, and b values range from 0.65~1.3, with a significant positive correlation between them (ΔD = 1.5 × b − 0.68). For this purpose, a statistical relationship is developed between b and ΔD values, and regional and temporal changes of these parameters are analyzed in order to reveal the potential of future earthquakes in the study region.


Introduction
Earthquakes are among the most serious natural disasters and cause great losses to humanity. Therefore, scientifically understanding seismic hazards and designing different seismic fortification intensity for different regions according to the degree of seismic hazard are of great practical significance. To evaluate the seismic hazard of an area, the seismicity of the area must first be analyzed. Numerous statistical models have been proposed to describe the regional and temporal variations of seismicity by using scaling laws in seismology [1][2][3][4]. Different seismic parameters can be used to analyze the characteristic of seismicity in these models, including the 1) magnitude of completeness (MC), which is the lowest magnitude at which all earthquakes in a space-time volume are credibly detected [5]; 2) b value, which defines the frequency-magnitude distribution of earthquakes, and 3) D value, which indicates the number of objects greater than a specified size with a power law dependent on size [6].
Fractal characteristics describe the processes of seismic systems from gestation to the selforganizing critical state. Changes in the heterogeneity degree of seismicity, seismotectonic structures, mechanical rock properties can be described by using the fractal dimension D value [3]. Wang et al. [7] applied fractal theory to the study of the spatial distributions of seismic sources, and the results show significant multifractality and a spatial variation in multifractality for epicentral distributions of earthquakes in west Taiwan. Chen et al. [8] proposed the Hill multifractal estimator to compare the generalized strain energy of seismic activities in New Zealand, Japan and the eastern and western regions of China and analyzed the multifractal characteristics of the seismic activities in each region. The results showed that the multifractal characteristics of seismic activities are closely related to the complexity of the geological structures. Pastén et al. [9] also used multifractal theory to analyze the change sequence of seismic electromagnetic currents and noted that the dimension of seismic electromagnetic currents has an obvious abnormal value before a large earthquake. Balcerak [10] analyzed the seismic time series of the 2003 Tokachi-Oki earthquake in Japan and noted that the multifractal dimension of the seismic activity would also have outliers before the next large earthquake.
Another fundamental parameter in seismology is the b value from the Gutenberg-Richter frequency-magnitude relation. According to global statistics, the b value of a large seismic area is generally close to 1.0 [11][12][13]. However, variations in b values do occur regionally, these variations are highly correlated with the regional characteristics of seismic activity and may be significant, for example, ranging from 0.5 to 1.3 near Parkfield on the San Andreas fault or from 0.5 to 1.5 in California and Japan [14,15]. The variation in the b value is sensitive to the differential stress [16]. Previous studies have confirmed that the b value is inversely dependent on the differential stress in both the laboratory [17,18] and the field [19]. A region with a low b value is implied to exhibit a large differential stress, suggesting its being toward the end of the seismic cycle [14]. Such a relationship can be used as a kind of precursory information for earthquake forecasting [16,20]. In probabilistic seismic hazard analysis (PSHA) studies, the b value is used to reflect the magnitude distribution of earthquakes in a certain area, which is crucial for the accuracy of the final evaluation results.
The relationship between the b value and the fractal dimension D value has been widely discussed in many studies [21][22][23][24][25]. Aki [21] proposed a simple relationship between the b value and the fractal dimension D and identified a positive correlation D=3b/c (where c is the slope of the log moment versus magnitude relation and equal to approximately 1.5). Wang et al. [22] showed a negative correlation between the b value and the D value in western Taiwan and indicated that a negative correlation for the observed seismicity might be reasonable due to the isolation of asperities or barriers on the fault planes. Wang et al. [25] obtained a similar result for earthquakes in Taiwan and showed a negative correlation between the b value and the D value. Mandal et al. [26] showed the correlation between those two scaling exponents could change from a negative one to positive in India.
Most previous studies use fractal theory to describe the fractal features of the study region on a large scale and macroscopically. We focus on the fine multifractal features in the study region. In the present work, the study area was divided into different scales (i.e., regions) and the depths of the regions are divided according to the results of the latest geological structure research, thereby providing the three-dimensional characteristics of the b value and D value of the Taiwan earthquake area. The results of this study are expected to be helpful in quantitatively describing the material composition of the crust and the degree of heterogeneity of geological structures, that is, the complexity of the geological structure environment and the determination of the PSHA.

Study Area
Taiwan island is located in the seismically active region at the junction of the Eurasian plate in the west and the Philippine sea plate in the east, as shown in Figure 1. The Philippine sea plate in northeastern Taiwan is being pushed under the Eurasian plate, while the Eurasia plate is going beneath the Philippine Sea Plate in south Taiwan. As a result of regional tectonic motion, the subduction zone, compressional folds and thrust belts lift the surface of Taiwan Island and most of Taiwan is under a northwest-southeast contraction with the convergence rate of about 80mm/year [27][28][29][30][31]. The unique location of Taiwan at the conjunction point between the Ryukyu arc and the Luzon arc leads to some geological structures at Taiwan such as folds, faults, and volcanoes. Thus, numerous earthquakes have been taking place inland. The Central Weather Bureau Seismic Network (CWBSN), which is responsible for monitoring regional earthquakes in Taiwan. This network currently consists of a central recording system with 71 telemetered stations that are equipped with three-component Teledyne/Geotech S13 seismometers [31]. According to previous research on the seismicity of the Taiwan [32], the CWBSN instruments were operating in a triggered recording mode prior to 1994 when continuous recoding began.
Many devastating earthquakes which have occurred in the past decade ( Figure 1). such as the 1999 Chi-Chi MW7.6 earthquake, which is the largest inland earthquake ever recorded in Taiwan [33][34][35][36][37], the 2002 Hualien MW 7.1 earthquake [38], the 2003 Chengkung MW 6.8 earthquake [39], and the Pingtung MW 7.0 earthquake in December 2006 [40]. The significant events were well recorded in Taiwan provide a valuable database to study both the lithosphere structure, earthquake sources, background stress state and stress perturbation.
Previous studies have inferred crustal strain or stress states in Taiwan from different data sets. Angelier et al. [33], Suppe et al. [34], and Chang et al. [35] investigated the stress state and stress perturbation according to fault-slip data sets, borehole breakout and elongation data GPS observations and earthquake focal mechanisms. They found that the orientations of the maximum horizontal compressive stress are generally NW-SE directed, consistent with the orientation of plate motion. More detailed analysis of spatial variations of stress or strain states after the Chi-Chi earthquake. Hsu et al. [29] discussed the deviatoric stress in the crust from GPS observations. They observed a 10°-20° rotation of fault slip vectors after the Chi-Chi mainshock and suggested a low deviatoric stress of about 1-3 MPa on the décollement beneath the Central Range. Wu et al. [37] evaluated spatial and temporal variations of stress orientations by using the stress tensor inversion methodology. They found significant rotations ( > 20°) of the maximum horizontal compressive axis immediately after the Chi-Chi earthquake. Chan et al. [30] investigated the variations of stress states associated with the 1999 Chi-Chi earthquake in different spatial scales and given the background deviatoric stress in the range of 10-50 MPa and the horizontal NW-SE directed maximum principal stress axis, the predicted fault types show strike-slip and normal faulting near the coseismic surface rupture and thrust and strike-slip faulting in central Taiwan. In general, early works attributed such stress heterogeneity to the existence of fault geometry complexity or pre-existing spatial stress heterogeneity on a local scale by various factors such as existing fault networks, lithology, geologic history, and so on [30].

Earthquake Catalog
In seismicity studies, the integrity of the seismic catalog in the region must be evaluated; otherwise, the seismicity parameter estimates would be biased. In this work, we used catalog from the Central Weather Bureau Seismic Network (CWBSN) for the period 1 st January 1995-1 st January 2019. To evaluate the reliability of the catalog, the completeness of the earthquake catalog was tested by calculating the magnitude of completeness, or MC, is defined as the minimum magnitude at which the cumulative Frequency Magnitude Distribution (FMD) departs from the decay [32]. Most techniques that focus on estimating MC follow the validity of the Gutenberg-Richter law, such as the Goodness-of-Fit Test [41], the Entire Magnitude Range [15], and the maximum curvature (MAXC) method [42].
MC usually shows a non-stable value in time. In the present work, the change of MC values as a function of time is estimated using a mowing window approach with the MAXC technique. Because it is a robust and straightforward method of estimating MC based on finding the magnitude bin with the highest frequency of events in the frequency magnitude distribution plot. The whole earthquake catalog with M ≥ 1.0 is included in the estimation of MC value and it is plotted with its standard deviation for samples of 250 events/windows. Figure Figure 3 that declustering algorithm has removed dependent events from the original catalog and after processing the declustering algorithm, a more homogeneous, reliable and robust earthquake catalog has been obtained. We select events from the catalog with focal depths less than 40 km, which is comparable to the thickness of the seismogenic zone in Taiwan [43,44].

Multifractal Methods
Many irregular forms and phenomena exist in nature, and self-similarity is observable at different scales. Patterns with scale invariance and self-similar irregularity are called fractals. In general, fractals can be divided into two categories: uniform fractals, which are based on strict geometric self-similarity on infinitely many scales, such as the famous Koch curve; and multifractals, which are statistically self-similar and non-uniform. Most fractals in nature are multifractals with statistical self-similarity at a certain scale. Fractal characteristics can be quantitatively characterized by fractal dimension values.
The Renyi generalized fractal dimension is the most common method for describing multifractals. Based on information theory, it defines the q-order probability moment: where Dq is the generalized dimension.
where Pi is the probability that events fall into a box with length ε; q is the estimated order of the moment and can take any real number ranging from -∞ to +∞. Dq for larger positive q shows the fractal property of dense regions, where Pi is large, and Dq for negative q displays that of thin regions, where Pi is small. The probability Pi can be estimated by the box-counting method for the observed data.
Using the box-counting method, the study researches fractal characteristic of the seismic spatiotemporal space by calculating generalized fractal dimension of the time series and spatial distribution in the seismic activity. Each seismic event is set as a data point in the seismic dataset, ε as minimum spatio-temporal interval. With the minimum spatio-temporal scale ε, N(ε) is number of non-empty subsets under the whole seismic dataset. As the ε value head towards zero, the fractal dimension is calculated. The capacity dimension Dc set as [45]: When the i th component contains the number of data points as Ni(ε), the probability for a point data belongs to the i th subset with the spatiotemporal scale ε, is defined as: The information dimension Di is 0 ( ) log ( ) lim log A single value of fractal dimension (Dc and Di) is not enough to characterize the fractal properties of multifractal objects. Therefore, the idea of fractal dimension has been extended to the generalized fractal dimension Dq. Generally speaking, D0 (q = 0) relates to the capacity dimension DC, D1 (q = 1) to the information dimension Di. The value of Dq can be estimated from the linear portion of the plot of log X(ε) versus log ε.

Estimation of b Value
In their 1942 article "Global Seismicity", which was published in the Journal of the American Geological Society, Gutenberg and Richter noted that seismic activity can be characterized by the following relationship [46]: where N(m) is the frequency of earthquakes with a magnitude greater than or equal to m. Parameters a and b can be obtained by regression on the actual earthquake catalog. In the PSHA, b is an important parameter that describes the distributions of earthquake magnitudes and frequencies, and a is an important parameter that actually measure of the level of seismic activity in a certain region. The most common methods for estimating the b value are maximum likelihood and least squares estimation. However, researchers widely recognize that the least squares method should not be used because the cumulative frequencies of earthquakes violate the basic assumption of independence of the data [47].
In the present study, we use a maximum likelihood method to estimate the b value and the standard deviation [48,49]. The b value and the standard deviation (σ) of b value is as follows: where is the average of magnitudes in the earthquake catalog; M0 is the initial magnitude; N is the number of events used for the b value estimation.

Multifractal Characteristics of theTime Distribution
A time series for earthquakes is a finite and discrete point set. To ensure reliable results, this paper selects historical seismic data for the Taiwan seismic province through integrity testing, and the time distribution range is 1 st January 1995 to 1 st January 2019. In a multifractal analysis, the number of point sets to be studied should be as large as possible. When the sample size is too small, the obtained fractal dimension may have pseudo-multifractal characteristics with parameter changes. Generally, the sample size should be no less than 200 [50].
To ensure the reliability of the results, the calculation windows were 1 year long, except for the 1995-1998 calculation window. The results are shown in Figure 4, where the value of q is unified from [-5,5]. In the figure, the curves from top to bottom are the multifractal evolution states over time of q from 5 to 5. In order to make an assessment on the seismic activity with time, the temporal distribution of all M ≥ 1.0 earthquakes from 1995 to 2019 is plotted as seen in Figure 5. There are a few strong earthquakes larger than 6.5 from 1995 to 2019. Many strong earthquakes with magnitudes greater than 6.5 are accompanied by dramatic changes in the corresponding curves in Figure 4. Specifically, most of the curves are in a state of loose tension before strong earthquakes occur; that is, the Dq value in the negative q region increases sharply or the Dq value in the positive q region decreases sharply or both.

Multifractal Characteristics of the Spatial Distribution
Similar to the characteristics of seismic activities in a time series, the spatial distribution can also be regarded as a series of discontinuous point sets in a multifractal analysis. In addition, to further describe the seismic activity characteristics reflected by the multifractals at different depths in the Taiwan seismic province, we first use k-means clustering to stratify all the earthquakes in the study area (focal depths ≤ 40 km). According to the sample size requirements for a multifractal analysis, we divide the study area into 5 layers: H1 (0 ~ 9 km), H2 (9 ~ 14.5 km), H3 (14.5 ~ 21.5 km), H4 (21.5 ~ 29.5 km), and H5 (29.5 ~ 40 km) ( Figure 6). Each layer is latticed, a spatial sequence reflecting the integrity of its spatial distribution is established, and the generalized fractal dimension Dq of each layer is calculated. For the value of q, the entire set of real numbers can theoretically be obtained; the larger the range is, the better the characteristics of the multifractal structure can be displayed [51]. Considering that an increase in q leads to an exponential increase in the calculation results and an overflow, after many experiments, q is set to [-10,10] in this paper. The calculation results are shown in Figure 7. The Dq values at different depths all monotonically decrease with increasing q and are nonlinear, subtractive functions. The value of q gradually increases from negative to zero and then to a positive value, and this change is analogous to that from describing the complexity of seismic system structures in sparse areas to describing the complexity of structures in clustered areas. The curve first decreases sharply and then slowly decreases, which indicates that the spatial multifractal complexity increases and is consistent with the complexity of the seismotectonic structure. To quantitatively describe this complexity, we adopt an area the size of the multifractal complexity of the parameter ΔD (the difference between D10 and D-10), reflecting the complexity of the regional seismotectonic structure [52,53]. Based on the characteristics of the seismicity and seismotectonic structure and zonation in the Taiwan region, the study area is first divided into east and west seismic tectonic zones. Furthermore, the eastern seismic tectonic zone is further divided into the Ryukyu arc (RK), TaiDong valley (TD), and Luzon arc (LZ). The western seismic tectonic zone is divided into the XinZhu fault (XZ), BeiGang uplift (BG), and Southwest basin (SB). The multifractals of each depth layer in different regions are calculated, and the results are shown in Figure 8.  Figure 8 attempts to show the distribution of ΔD values at different spatial scales, and the results show the distribution of ΔD distribution at different depth layers and at different regions of the same depth layer is highly uneven. This phenomenon implies the complexity of regional geological structures.
The ΔD for the study area (b) in H1 (0 ~ 9 km) is the largest, and the value tends to decrease with depth until H4 (21.5 ~ 29.5 km). The ΔD of the last depth layer H5 (29.5 ~ 40 km) increases instead, which indicates that the multifractal structure of this layer becomes relatively complex again. The ΔD results for the east and west seismic tectonic zones (c) show that eastern Taiwan has higher values than western Taiwan at all layers except H3 (14.5 ~ 21.5 km). In the zones of tectonic features (d), different regions have different ΔD values. Compared with those in other zones, the ΔD value in the TaiDong valley (TD) is the largest at all layers. In western Taiwan, the ΔD values of BeiGang uplift (BG) and Southwest basin (SB) are similar and slightly less than that of XinZhu fault (XZ) in the first layer (H1) but greater in the second layer (H2).

Relationship between ΔD and b Values
We consider that seismic activity has the multifractal characteristics, and a single fractal dimension D value may not be able to accurately describe the relationship between spatial fractal characteristics and b values. In this section, we calculate the b value and its standard deviation for each region by using the maximum likelihood estimation. Then, a correlation analysis is carried out based on the above results for the multifractals. Figure 9a shows the relationship of regression fit between b and ΔD values. Negative correlation coefficient (r) is calculated as 0.84 and it can be accepted as a strong enough. Linear regression fit is used and following relationship is obtained: The results also show that the ΔD value ranges between 0.2 and 1.

Discussion
Seismically active fault regions are complex natural systems and they exhibit scale-invariant or fractal correlation between earthquakes in space and time [53]. By understanding the multifractal features of the time distribution of seismic events, the temporal variations of some characteristics of seismotectonic structure can be reflected in detail. Multifractal characteristics of the time distribution results show that the curve oscillates between compression and tension. Most of the curves are in a state of loose tension before strong earthquakes occur. This finding is also consistent with previous research results [54,55]. We also find that the temporal multifractal generalized fractal dimensions of the Dq vs. q curve usually changes abruptly before and after strong earthquakes.
The spatial distribution of seismic events is controlled by the geological and tectonic environments, and the multifractal features of a seismic event's spatial distribution are closely related to the complexity of the seismotectonic structure. To describe this complexity accurately and quantitatively, we first divide the research area into different scales and depths and then use the size of the ΔD values to represent the complexity of the seismotectonic structure in different regions. In the Gutenberg-Richter relation, which is widely used in seismological studies, the b value represents the proportion of all earthquakes of different magnitudes and ranges in a certain area and is sensitive to the differential stress. To date, most studies on b values focus on their temporal and/or spatial variations. Varotsos et al. [56] based on the natural time analysis of earthquake data revealed that a fall of the b value observed before large earthquakes reflects an increase in the order parameter fluctuations upon approaching the critical point (mainshock) and stems from both origins of selfsimilarity. Recently, several studies have proved that b value may decrease before a large earthquake and explained that the variation of the b value could reflect the physical process of stress evolution and crack growth [16,20,57]. laboratory measurements have established that the higher the differential stress applied to the rock sample, the lower the b value [17,18]. Some authors [58,59] pointed out the correlation between b value and tectonic regimes, and the final results are widely consistent: as power laws indicate scale invariance, inverse dependence of the b value on differential stress. Spada et al. [60] and Wiemer and Wyss [61] proved that the b value would decrease with increasing depth and pressure. Grassberger et al. [62] showed that the physical reasons for the decrease in b value with increasing depth are mainly stress changes and anisotropy of the medium. Chan et al. [43] showed that large earthquakes in Taiwan are located predominately in regions with low b values. Table 1 shows the ΔD and b values of different depth layers and regions. The decreasing trend of ΔD from H1 to H5 indicates that the heterogeneity of the crustal in Taiwan crust is decreasing. According to column 3, the maximum b value of the crust in Taiwan is relatively large, which means that the stress level in the upper crust is relatively low and may also indicate that the strength of the rocks constituting the upper crust is relatively low [63]. The decreasing trend of b values from H2 to H5 clearly shows that the pressure increases with increasing depth, leading to a change in the stress state. The b value is a minimum at H5 (29.5 ~ 40 km), indicating that the crustal intermediate stress level is the highest at this depth compared to the other depth layers. Many studies on the relationship between multifractals of seismic spatial sequences and b values have been reported [21][22][23][24][25]50]. Aki [21] proposed a simple relation between the b value and the D value: D = 3b/c (c = 1.5). However, in follow-up studies, the observed results did not satisfy this relationship [52,64]. We argue that the multifractal characteristics of seismic activity leading to a single fractal D value and the b value relationship may not be able to accurately determine the spatial multifractal characteristics of the ΔD and b values. We obtain a significant positive correlation between them, and the range of ΔD values is larger than that of b values.
Based on the complexity characteristics of seismically active fault and fractal theory can be used to reveal the characteristics of non-stability and discontinuity in the complex system. Spatial and temporal analyses of the current seismic behaviors of earthquakes provide valuable evidence of future earthquake hazards in the Taiwan region. For this reason, seismotectonic parameters such as ΔD value, b value, and their relations must be more carefully estimated in order to reveal the significant anomaly regions before strong earthquakes in the future.

Conclusions
Based on data from 1 st January 1995 to 1 st January 2019 in the Taiwan seismic region, we analyzed the spatial and temporal multifractal characteristics of seismicity in this region. The results show that the temporal multifractal generalized fractal dimension of the Dq vs. q curve for the Taiwan earthquake zone oscillates between compression and tension and reflects the repeated changes in the seismogenic system from steady to unsteady, with the curve usually changing abruptly before and after strong earthquakes. Specifically, the Dq value in the negative q region increases sharply, or the Dq value in the positive q region decreases sharply, or both.
Regarding the spatial multifractal characteristics, we calculated the generalized fractal dimension Dq for five depth layers of the shallow source tectonic system. The results show that the ΔD of the H2 layer is the largest, suggesting that the heterogeneity of the stratum in this layer is the most prominent. Thus, the geological and tectonic environment in this layer is the most complex.
Meanwhile, the seismic catalog shows that this layer has the most earthquake events. In addition to the first layer, the ΔD value in H3 (14.5 ~ 21.5 km) in the western area of Taiwan is larger than that in the eastern region, and the ΔD values in other deep layers in the eastern region are larger than those in the west. This finding indicates that the geological and tectonic environment in eastern Taiwan is more complex than that in the west and that the level of seismic activity is also higher in eastern Taiwan.
To estimate a more up-to-date and reliable statistical relation between two seismotectonic parameters b and ΔD values, linear regression is preferred. The relationship of ΔD = 1.5×b-0.68 is suggested with a strong enough positive correlation (r = 0.84) for Taiwanese earthquake distributions. For the Taiwan earthquake area at different depths, the b values range from 0.65 ~ 1.3, while the ΔD values range from 0.2 ~ 1.5.
Author Contributions: The first author C.H. is a major contributor to the paper and responsible for the core writing of the paper. The corresponding author C.C. is another major contributor to the paper and mainly responsible for the direction of research and modifying the paper. The other two authors, L.N. and J.Y., provided opinions on the article's revision and collected some of the data. All authors have read and agreed to the published version of the manuscript.