BOX COUNTING FRACTAL DIMENSION AND FREQUENCY SIZE DISTRIBUTION OF EARTHQUAKES IN THE CENTRAL HIMALAYA REGION

To establish the relations between b-value and fractal dimension (D0) for the earthquake distribution, we study the regional variations of those parameters in the central Himalaya region. The earthquake catalog of 989 events (Mc = 4.0) from 1994.01.31 to 2020.10.28 was analyzed in the study. The study region is divided into two sub-regions (I) Region A: 27.3°N -30.3°N and 80°E -84.8°E (western Nepal and vicinity) and (II) Region B: 26.4°N -28.6°N and 84.8°E -88.4°E (eastern Nepal and vicinity). The b-value observed is within the range between 0.92 to 1.02 for region A and 0.64 to 0.74 for region B showing the homogeneous nature of the variation. The seismic a-value for those regions ranges respectively between 5.385 to 6.007 and 4.565 to 5.218. The low b-values and low seismicity noted for region B may be related with less heterogeneity and high strength in the crust. The high seismicity with average b-values obtained for region A may be related with high heterogeneity and low strength in the crust. The fractal dimension ≥1.74 for region A and ≥ 1.82 for region B indicate that the earthquakes were distributed over two-dimensional embedding space. The observed correlation between D0 and b is negative for western Nepal and positive for eastern Nepal while the correlation between D0 and a/b value is just opposite for the respective regions. The findings identify both regions as high-stress regions. The results coming from the study agree with the results of the preceding works and reveal information about the local disparity of stress and change in tectonic complexity in the central Himalaya region.


INTRODUCTION
On 25 April of 2015, the central Nepal was struck by Mw 7.8 (6.9 mb) earthquake. The event took place on the Main Himalayan Thrust (MHT) Yue et al., 2017) and was resulted from unfastening of the lower edge of the sealed portion of the MHT, along which the Himalayan hunk is thrust over India . The MHT is the dynamic structure along the Himalayan arc and reaches the surface at Main Frontal Thrust (MFT) as observed by the Japanese Agency (Hubbard et al., 2015). Within the 17 days (on 12 May 2015), the nation was again struck by another earthquake of Mw 7.3 (6.7 mb) . These two large events with their aftershocks have partly released the strain accumulated between seismic events along the MHT (Hayes et al., 2015;Sreejith et al., 2018). The mainshock and prompt running aftershock sequence caused strong ground shaking not only in Nepal but also in some parts of neighboring nations India and China, causing more than 8900 deaths . Ground shaking in Kathmandu basin was strong enough to collapse many historical structures that had persisted the preceding earthquakes (Galetzka et al., 2015).
From the study of literatures, the spreading of the rupture in the past Himalayan earthquakes were understood to be organized by asperities and supplementary structural complexities (Bilham, 2019;Mugnier et al., 2017). Since the Gorkha earthquake of 2015 happened between the rupture areas of the 1505 Ms 8.2 Lo-Mustang earthquake in western Nepal (Ghazoui et al., 2019) and Nepal-Bihar earthquake (Mw 8.2 ) of 1934, it was expected to have the magnitude greater than 8 (Letort et al., 2016). Sreejith et al. (2018) disclose the existence high strain and stress region (asperity), one towards the west and the other towards the east of the 2015 Gorkha epicenter. Although, the asperity found in the eastern side had broken partially in the Gorkha earthquake, the asperity in the western side is still unbroken and preparing to host the future large earthquake. Furthermore, it is known that the elastic potential energy releasing process in the Himalaya is either by the large rupture reaching the surfaces or by incomplete rupture (Scholz & Campos, 2012). The case of partial rupture for Gorkha earthquake indicates the event did not release all the stored strain in the region and the region is overdue for major earthquake Ghazoui et al., 2019). This investigation applies the theory of nonlinear dynamics (Helmstetter et al., 2005;Yin et al., 2019) to identify the stress regions in terms of b-values and box counting fractal dimension of epicenter distribution. Based on the previous studies on the aftermath of the Gorkha earthquake, the study region is divided into two sub-regions ( Fig. 1)

Data
The earthquake events were collected from International Seismological Centre (ISC) catalog and United States Geological Survey (USGS) catalog. While selecting the events priority was given to ISC catalog and events were scrutinized on the basis of locations (latitude and longitude) and date with time of occurrence The Mw scale of USGS catalog were converted into mb scale from the relation given by Das et al. (2011). We found 2457 events for the period 1964-02-01 to 2020-11-23 with in the latitude 26°N to 31°N and longitude 80°E to 89°E. After declustering the catalog (Gardner & Knopoff, 1974) from the software ZMAP (Wiemer, 2001) , we only retained 1185 events. The software ZMAP is the widely accepted tool for the statistical analysis of the seismicity pattern so equally applicable for the study in the Himalayan region. The cumulative curve (Fig. 2) shows straight line with constant slope after January 1994 giving the completeness (Mc=4.0) of catalog, so the catalog starts from 1994 with 989 events. For the analysis of the epicenters, we have selected rectangular areas of dimension 3°×3° in the region A (subregions A1 to A10) and 3°×2.2° in the region B (subregions B1to B4) as mentioned in the Table 1 and Table 2. The variation between b and D0 over time was evaluated using sliding windows of aforementioned areas with window advance of 0.2 degree along the direction of longitude in order to cover region of the study. The time number histograms (Fig. 3) display the high seismic activity in region A compared to region B whereas the increase of seismic activity around the year 2015 is observed for both regions.

METHODOLOGY Earthquake frequency size distribution b-value
The cumulative frequency size dispersal of earthquakes is given by the relation.
=where is the number of earthquakes with magnitude greater or equal to M, a is the constant signifying the level of the productivity and b is the constant known as b-value that signifies magnitude distribution (Gutenberg & Richter, 1944). The frequency-magnitude analysis of the data is done by the maximum curvature technique and the b-value is estimated by maximum likelihood approach (Aki, 1965). b = log 10 . is the average value of the magnitudes, M is the lower limit of the magnitude in the catalog and ΔM is the binning width of the catalog. The uncertainty occurred in the estimation of the b-value is given by the relation below (Shi & Bolt, 1982).
where is the number of earthquakes in the given sample. The magnitude of completeness (Mc) was determined as 4.0 (for 989 events) using the iterative method (Wiemer & Wyss, 2000). The b-value is normally inferred as the sign of the heterogeneity present in the material of the medium and found to fluctuate from 1.0 to 1.6 for the global seismicity (Mogi, 1967;Sobiesiak et al., 2007). Also, it is interpreted as the indicator of applied tangential stresses (Scholz, 1968;Schorlemmer et al., 2005;Singh et al., 2009). For the earthquakes in the California, the b-values found to vary from 0.45 to 1.5 (Gutenberg & Richter, 1944). The b-value of 0.80 ± 0.05 is computed for aftershocks sequences of Gorkha earthquake between 2015 May 25 and June 8 . For the 820 aftershocks (2.7 < M < 7.3) of the Gorkha earthquake, the b-value of 1.11 ± 0.08 was computed during the period April 25 -November 12, 2015 (Nampally et al., 2018). The b-values for the longitude range 80°E to 89°E and latitude range 26°N to 31°N found to vary from 0.5 to 1.6 for the period 1964-2015 (Ghosh, 2020). Furthermore, the b-value contour map depicts the low b-value patch in the western part of the Gorkha region as the possible zone of future robust seismic activity (Tiwari & Paudyal, 2021). Fractal distribution A distribution with fracture or uneven basic structure reiterate in different scales is a fractal. The index providing a statistical value to compare how a pattern alters with the measuring scale is the fractal dimension (Lopes & Betrouni, 2009). Earthquakes are known as self-organized dynamical system progressing spatially and temporarily towards a critical stationary stage, so the fractal theories are beneficial in learning seismic activities (Rodríguez Pascua et al., 2003). The spectrum of the fractal dimension can be obtained from the generalized fractal dimension and it is calculated from the relation (Grassberger & Procaccia, 1983).
is the Heaviside step function, r is the scaling radius, is the total number of pair of data points within the sample volume, and are the location of the epicenter in latitude and longitude of the i th event and the j th event while ( , ) is the q th integral.
The generalized dimension , in terms of generalized correlation sum can be written as, log ( ) log For q = 0, the above formula gives capacity dimension. It gives information dimension for q=1, and correlation dimension for q=2 (Roy & Padhi, 2007). Therefore, for q = 0, the formula reduces to 0 = − lim →0 log ( ) log The above formula is similar to the box counting technique (Aggarwal et al., 2017;Mandal & Rodkin, 2011;Radziminovich et al., 2019) as log ( ) log where (r) is the number of boxes of side r occupied by point sources. With the advantage of simple mathematical formulation and empirical estimation, box-counting dimension becomes the widely used dimension (Falconer, 1997). The box counting technique estimates the capacity dimension by measuring the space-filling characters of a fracture set with the change in grid of the scale (Gospodinov et al., 2010;Mandal & Rodkin, 2011). In the box-counting method, the epicenter distribution is covered by squares or cubes of decreasing size (r). If the epicenter distribution points have fractal character the slope of the log ( ) − log curve, in particular range of r, provides capacity dimension.

RESULTS AND DISCUSSION
The b-value of the region A ranges from 0.91 ± 0.08 to1.02 ± 0.10 (Table 1; Fig. 4a) while b-value of Region B ranges from 0.64 ± 0.04 to 0.78 ± 0.06 (Table 2; Fig. 5a). The relatively high b-value and high seismicity observed in the region A may indicate that the area is linked with high material heterogeneity and low strength in the crust as explained in the literatures Mogi (1967) and Scholz (1968). Many small earthquakes can be expected in regions. The low b-values and relatively low seismicity in the region B may be linked with low degree of heterogeneity and strong rheological strength in the crust that approves the preceding works (Bayrak et al., 2013;Khan et al., 2011;Khattri & Tyagi, 1983;Wyss, 1973). The creeping process occurring in the faults is generally related with higher bvalue, and asperities present in the faults are often related to lower b-value (Kawamura & Chen, 2017;Oncel & Wilson, 2002). Thus, the earthquakes in the region A may be due to the creeping process and the dominating number of aftershock in the region B are related to the asperities found in the region.
The fractal dimension of the region A ranges from 1.79 ± 0.01 to 1.89 ± 0.03 (Table 1; Fig. 4b). Similarly, for the region B it ranges from 1.82 ± 0.02 to 1.86 ± 0.02 (Table  2; Fig. 5b). These non-integers values of the dimension between 1 and 2 indicates that the earthquakes are spread on the boundary between a line and a plane (Rodríguez Pascua et al., 2003). Tosi (1998) explained that possible values of fractal dimensions do not depend on the dimension of the embedding space and the value ranges between 0 and 2. The limiting value of dimension can be interpreted as D0 close to 0 has all events bunched into single point, D0 near to 1 means the source region has linear nature and D0 near to 2 shows that the epicenters are casually or evenly spread over a two-dimensional (2D) implanting space (Roy & Ram, 2006;Tosi, 1998). Thus, the higher D0 (1.79 -1.89) obtained for the regions may be the indication of scattered epicenter distributions. The lowest b-values and the highest D0 computed for region B indicate the vulnerability of the region for occurrence of the large earthquakes (Naimi-Ghassabian et al., 2018). The works similar to the present investigation have been carried out by the previous researchers. The fractal correlation dimension of 1.66 is noted for 820 aftershocks between April 25 -November 12, 2015 for Gorkha earthquake (Nampally et al., 2018). The fractal dimension of 1.69 ± 0.05 is obtained for the epicenter distribution of 2001 Bhuj earthquake (Mw 7.7) in Gujarat, India (Aggarwal et al., 2017).  The D0 ─ b and D0 ─ (a/b) relations could possibly replicate regional seismicity and earthquake threat and believed to be applicable in earthquake hazard related studies (Bayrak & Bayrak, 2012;Pailoplee & Choowong, 2014). The negative correlations between b-value and D0 (r = -0.253) is obtained for the region A (Fig. 6a) while the positive correlation (r = 0.511) is observed in the region B (Fig. 6b).

Table 1. Subdivided regions of A with period, location, earthquakes numbers, frequency magnitude distribution coefficients a-value and b-value, ratio of a and b, fractal dimension (D0) and coefficient of determination ( )
The positive correlations between a/b and D0 (r = 0.829) is observed for the region A (Fig. 6b) and negative correlation (r = -0.454) is observed for the region B (Fig. 7b). The D0 ─ (a/b) correlation was explained to be quite dependable and effective comparing to the D0 ─ b correlation for indicating seismic hazards (Bayrak & Bayrak, 2012), but in this work, it is justified only for the region A showing the significant of the distribution of D0 ─ a/b with r = 0.829. The negative correlations mean the drop in one parameter for rise in another and the positive correlations means parallel increase or decrease in corresponding parameters. The drop in D0 and rise in b-value or parallel rise or fall in D0 and a/b for the region A indicate the presence of dense and complex network of faults having likelihood of occurring earthquake with larger magnitude (Oncel & Wilson, 2007). The positive correlation between D0 and b or negative correlation between D0 and a/b in region B can be interpreted as reduced probability of larger magnitude earthquake because the stress in the region is reducing through lower magnitude events. The preceding study had also noticed negative correlation between b-value and fractal dimension for this region from the earthquake data of seismological catalogs of the USGS and ISC (Ghosh, 2020).
For the region A, there is a decrease in D0 and increase in b-value as the sliding windows move from west to east while for the region B both parameters increase as the window moves from west to east (Fig. 8). Thus, the spatial variation of D0 and b-value is in the negative correlation for the region A and in the positive correlation for the region B. The negative spatial correlation may be responsible to the rise in stress concentration (lower b) and a decline in the clustering of the epicenter (increased D0).

CONCLUSIONS
The analysis is carried out on earthquake catalog (Mc = 4.0) of 989 events from January 1994 to November 2020. The spatial variation of fractal dimension (D0) with b-value shows negative correlation for western Nepal and vicinity while correlation is positive for eastern Nepal and vicinity. The correlation of D0 with a/b ratio is just opposite. From these results, it can be inferred that regions considered are under high stress and posing risk for generating large earthquake in the future. The western Nepal is safe against the seismic hazards as the b-value of the region is getting average b-value of 1.0 through the small earthquakes event.
The observed low b values for eastern Nepal imply that it is highly stressed and could be a source region for future strong events. The high fractal dimension value ≥1.74 obtained for both regions indicate the earthquake events are distributed in two-dimensional embedding space (source zones). This study indicates the regional variation of stress level and difference in tectonic complexity.