Self-Organized Fractal Seismicity and b-Value of Aftershocks of the 2015 Gorkha Earthquake, Nepal

Sanchit Minocha¹, Imtiyaz A. Parvez²

¹Indian Institute of Technology Roorkee (IITR), Roorkee, India.
²CSIR Fourth Paradigm Institute (Formerly CSIR C-MMACS), Bangalore, India.

DOI: 10.4236/ijg.2020.118030                

Abstract

The Gorkha Earthquake that occurred on 25^th April 2015 was a long anticipated, low angle thrust-faulting shallow event in Central Nepal that devastated the mountainous southern rim of the High Himalayan range. The earthquake was felt throughout central and eastern Nepal, much of the Ganges River plain in northern India, and northwestern Bangladesh, as well as in the southern parts of the Plateau of Tibet and western Bhutan. Two large aftershocks, with magnitudes 6.6 and 6.7, occurred in the region within one day of the main event, and several dozen smaller aftershocks occurred in the region during the succeeding days. In this study, we have analyzed the 350 aftershocks of the 2015 Gorkha Earthquake of M_w 7.8 to understand the spatial and temporal distribution of b-value and the fractal correlation dimension. The b-value is found to be 0.833 ± 0.035 from the Gutenberg-Richter relation by the least squares method and 0.95 ± 0.05 by the maximum likelihood method, indicating high stress bearing source zone. The spatial and temporal correlation dimension is estimated to be 1.07 ± 0.028 and 0.395 ± 0.0027 respectively. Spatial correlation dimension suggests a heterogeneous distribution of earthquake epicenters over a linear structure in space, while the temporal correlation dimension suggests clustering of aftershock activity in the time domain. The spatial variation of the b-value reveals that the b-value is high in the vicinity of the mainshock which is due to the sudden release of stress energy in the form of seismic waves. The spatial distribution of correlation dimension further confirms a linear source in the source zone as it varies from 0.8-1.0 in most of the region. We have also studied the temporal variation of b-value and correlation dimension that shows positive correlation for about first 15 days, then a negative correlation for next 45 days and after that, a positive correlation. The positive correlation suggests that the probability of large magnitude earthquakes decreases in response to increased fragmentation of the fault zone. The negative correlation means that there is a considerable probability of occurrences of large magnitude earthquakes, indicating stress release along the faults of a larger surface area [1]. The correlation coefficient between b-value and the correlation dimension is estimated to be 0.26, which shows that there is no significant relation between them.

Keywords

Nepal Earthquake, Gorkha Earthquake, b-Value, Fractal Dimension

1. Introduction

An earthquake is a sudden violent shaking of the surface of the earth, resulting from the sudden release of energy within Earth’s crust in the form of seismic waves. An aftershock is an earthquake that occurs after the mainshock in the same region but generally of a smaller magnitude. Aftershocks are the results of the adjustment of crust around the displaced fault plane which is caused by the main shock. The earthquakes occur along the fault zones or material heterogeneities present in the crust. The distribution of faults or fractures in the crust possesses a fractal structure [2]. These fractal structures in time, space and magnitude dimensions are observed not only in the distribution of faults but also in the seismicity of earthquakes [3]. In 1993, Kagan suggested that earthquakes do not occur on a single surface, but rather on a fractal structure of many closely related faults [4]. Also, the number of aftershocks decay follows a power law [5], suggesting a scaling property between the main shock and its aftershocks using the relation $f\left(m\right)=m^{-D}$ (where m is the magnitude and D is a kind of fractal dimension) which further indicates scale invariance property. A large number of studies have shown that scale invariance is the only property to characterize a fractal structure [2] [6]. The relationship between the magnitude and the total released seismic energy ( $N\left(E\right)=E^{-2b/3}$ , where E is the total released seismic energy) also confirms the fractal nature of earthquakes.

The Gutenberg-Richter law ( ${\log }_{10}N=a-bM$ , where M is the magnitude and N is the number of earthquakes of magnitude greater than or equal to M) relates the magnitude of the earthquake with the total number of earthquakes for a given region and time period of at least that magnitude [7]. The constant b is known as the b-value, which represents stress regime and material heterogeneity of the earth’s crust. The b-value is usually close to 1 over long time periods and large regions but statistically significant variations have been found for earthquake data for a short period and small geographical regions. [8] has derived a simple relation between the fractal dimension of the earthquake seismicity (D) and the b-value (D = 2b) which shows a positive correlation between the two. However, in contrary to this, a negative correlation had been noticed between the b-value and the fractal dimension by [9] for the earthquakes of the Tohoku region, Japan. The fact, which the fractal dimensions are usually fractional, is a bit interesting. The fractal dimension represents the measure of crustal deformation and a degree of clustering in space and time [6] [8] [10] [11].

There are different types of fractal dimension depending on the various methods of measurement, namely, Euclidian dimension (d), similarity dimension (D_s), Hausdorff dimension (D_H), box-counting dimension (D₀), information dimension (D₁) and correlation dimension (D₂). The box-counting dimension (D₀), also known as the capacity dimension, measures the space-filling properties of a fracture with respect to changes in grid scale [9]. On the other hand, the correlation dimension (D₂) measures the degree of clustering of a set of points [12]. The correlation dimension has been found for both earthquake epicenters and hypocenters in numerous studies [13]. The correlation dimension is usually found by using the correlation integral method, as it is the simplest and most naïve way to estimate D [14]. In general, D₂ varies from 0 - 2 [15]. When it is 0, all events cluster into a single point; D₂ = 1 indicates a line source or linear structure and when all events are homogeneously distributed over a two-dimensional embedding space, D₂ = 2 [16].

On 25^th April 2015, a large shallow earthquake of magnitude 7.8 occurred in the Gorkha region of central Nepal. The aim of this research is to find the b-value and the fractal correlation dimension in space and time for the aftershocks of the 2015 Nepal Earthquake and to understand the correlation between both of them. This study helps to understand the variation in stress accumulation in different regions of Nepal. For this, we have analyzed the 350 aftershocks of magnitude greater than or equal to 4.0, compiled by the United States Geological Survey (USGS), occurred in the period 25^th April 2015 to 28^th February 2017. The spatial and temporal fractal correlation dimension is estimated using the [14] algorithm while the b-value is estimated using the least squares method. The temporal variation of b and D₂ after the occurrence of the 2015 Gorkha earthquake of M_w 7.8 has been examined, and thus the correlation between the two of them is analyzed. We have also studied spatial variation in b and D₂.

2. Tectonic Setting of the Region

Figure 1. Geological map of Nepal showing the five latitudinal morpho-tectonic zones [22].

Heim and Gansser [23] divided the Himalayan belt into four tectonostratigraphic zones, which can further be divided into five latitudinal morpho-tectonic zones namely (Figure 1):

1) Terai Zone (Gangetic Plain)

Figure 1. Geological map of Nepal showing the five latitudinal morpho-tectonic zones [22].

2) Sub-Himalayan Zone (Siwaliks)

3) Lesser Himalayan Zone

4) Higher Himalayan Zone (Greater Himalaya)

5) Tibetan-Tethys Zone

In between the Lesser Himalaya and the Greater Himalaya, there is a Main Central Thrust Zone, which is bounded by Main Central Thrust (MCT) I and Main Central Thrust (MCT) II [24]. The metamorphic grade, in central Nepal, increases from low to medium towards the MCT as one move from north to south.

Due to the presence of the Main Frontal Thrust (MFT), Main Boundary Thrust (MBT) and Main Central Thrust (MCT), Nepal had experienced several damaging earthquakes. The earthquakes of 1833, 1869 and 1988 occurred along the same segment of the Main Himalayan Thrust (MHT) as that of the 2015 earthquake (Figure 5). This segment is MCT and is adjacent to segments northwest that ruptured in 1505 [25] and to the southeast that ruptured in the 1934 Nepal-Bihar earthquake [26]. The 1833 Earthquake that took place in Kathmandu was of magnitude (M_s) 8.0. The 1988 Earthquake that occurred in Nepal was of magnitude (M_w) 6.9 and caused the death of 709 persons while thousands got injured.

3. Source Parameters of the Main Shock

On 25th April 2015, a severe earthquake struck near the city of Kathmandu in central Nepal. It killed about 9000 people, and thousands were injured. At 6:11 am (UTC time), the initial shock of magnitude (Mw) 7.8 was registered. Its epicenter was located about 34 Km east-southeast of Lamjung and 77 km northwest of Kathmandu. The depth of the focus was 15 km and the cause of the earthquake was the thrust faulting in the Indus-Yarlung suture zone. This earthquake released the seismic energy accumulated because of the compressional pressure between the Eurasian plate and the Indo-Australian plate. As the Indian plate subducts the Eurasian plate at an average rate of 4 - 5 cm annually, the height of the Himalayan mountain is increased by more than 1 cm every year.

A total of 361 aftershocks of the 2015 Nepal earthquake were recorded on the United States Geological Survey (USGS) Global Seismic Network (GSN). The geographical area considered for the aftershocks was from 26.096˚ to 30.581˚ in latitude and from 79.849˚ to 88.330˚ in longitude and the period considered was from 25th April 2015 to 28th February 2017. There are however, we analyzed only 350 aftershocks which were having magnitudes greater than or equal to 4.0.

4. Data Analysis

First, we have made a magnitude histogram to show the decay of the frequency of occurrences of magnitude ≥ 4 aftershocks (Figure 2). It shows that the aftershocks data is by the Gutenberg-Richter law. A Temporal distribution of aftershocks (Figure 3) classifies them on the basis of month and year. It shows that almost all the aftershocks of magnitude greater than 5.5 occurred in the April and May of 2015 and more than 75% of the total aftershocks had occurred in 2015. Thus, the aftershocks have decayed significantly with time. Secondly, to understand the spatial distribution of the aftershocks, 3-D plots of their hypocenters, from three different viewpoints, are made (Figure 4). These figures show that the main shock is surrounded by a large number of aftershocks from all the directions and thus, there is a minimal variation in the depth of these hypocenters. We have prepared the seismicity map of Nepal to examine the spatial distribution of the epicenters of these aftershocks (Figure 5). As all the major faults are drawn on this map, one can easily see that the mainshock occurred near the Main Central Thrust Zone.

4.1. Estimation of b-Value

The b-value of these aftershocks is estimated using the Guttenberg and Richter (1944) relation:

${\log }_{10}N=a-bM$ (1)

where N is the number of earthquakes of magnitude greater than or equal to M, “a” and “b” are constants. Seismicity level is indicated by the constant “a”, and the amount of stress accumulated in a region is indicated by the constant “b” known as b-value. These constants are determined by applying linear regression over the log₁₀N Vs M, using the method of least squares (Figure 6).

Figure 2. A histogram showing the decay of the frequency of the occurrences of aftershocks, having magnitude ≥ 4.

Figure 3. Temporal distribution of the aftershocks (Month-wise and Year-wise). Most of the aftershocks occurred in April and May 2015. It also shows that the aftershock activity decayed significantly with time.

A higher b-value means that there is less number of larger earthquakes as compared to smaller earthquakes, whereas a lower b-value means that there are more numbers of larger earthquakes as compared to smaller ones. Thus, a low b-value indicates higher stress in a region [27].

[28] proposed the following relation to find b-value by the maximum likelihood method:

$b=\left({\log }_{10}e\right)/\left(M_a-M_c\right)$ (2)

where M_a is the average magnitude of all events, and M_c is the minimum magnitude of all the events. But since, the number of aftershocks is quite small, the b-value calculated by the maximum likelihood method is not a good estimate.

The b-value is estimated to be 0.833 ± 0.035 by the method of least squares whereas 0.95 ± 0.05 by the maximum likelihood method.

(a)(b)(c)

Figure 4. Spatial distribution of the hypocenters from three different angles and elevation. The azimuthal angles of (a)-(c) are −60˚, 30˚ and −150˚ while the elevations are 20, 60 and 60 respectively. Aftershock activity is highest in the vicinity of the mainshock.

Figure 5. Map of Nepal showing the distribution of epicentres of all aftershocks along with the major faults. MCT stands for Main Central Thrust, MFT stands for Main Frontal Thrust, and MBT stands for Main Boundary Thrust. It also shows the major and great earthquakes occurred in Nepal in the past.

Figure 6. The b-value is estimated by the slope of a best fit line using least squares method. The b-value is determined to be 0.833 ± 0.035.

4.2. Depth and Spatial Variation of the b-Value

To study the spatial variation in the b-value, we have chosen a small geographical area (84.3˚ to 87˚ longitude and 27.3˚ to 28.5˚ latitude). This area is divided into 36 grids of 0.3˚ × 0.3˚ and b-value is found for each of the grids having a minimum occurrence of 10 events. Mapping of b-value is done by interpolating the b-values (Figure 7). This map shows a spatial variation of the b-value from 0.49 to 1.09, higher values to the west indicates lower stress region whereas lower values to the east indicate higher stress region. This suggests a high probability of the occurrence of a large earthquake in Bihar, Sikkim and Assam Himalayan regions which are to the east of the mainshock and a low possibility of the occurrence of large earthquake in Himalayas west to Pokhara region. Also, the b-values are relatively higher in the vicinity of the mainshock which is due to the release of the accumulated stress in this region. The variation in the b-values in different zones may be related to structural heterogeneity and stress distribution in space [3] [29]. Since out of 350 aftershocks, 282 aftershocks are having a depth of 10 Km, the b-value variation with depth could not be studied.

4.3. Estimation of Spatial and Temporal Fractal Correlation Dimension

Spatial distribution of earthquake epicenters can be represented by a self-similar fractal structure, and the scaling parameter is known as the fractal dimension D [2]. Though many fractal models have been proposed for the distribution of epicenters, the correlation dimension is the easiest and the efficient method to estimate the fractal dimension. The correlation dimension is calculated using the correlation integral algorithm proposed by [14]. The fractal correlation is written as [30]:

$Cp\left(r\right)={\displaystyle {\sum }_{i=1}^{k}H\left(r-\Vert x_p-x_i\Vert \right)}$ (3)

where k is the total number of events in the catalogue, r is the radial distance to

Figure 7. The b-value mapping is done by dividing the area into 36 grids of 0.3˚ × 0.3˚. b-values are higher in the vicinity of the mainshock and are lower to the east of the mainshock.

be considered, x_i is the location of i^th earthquake event, and thus x_p is the location of the point at which correlation sum is evaluated. The double bars (||) represent the scalar distance which is to be considered between two events. H is the Heaviside function that is 0 when the argument is negative and is 1 otherwise. Therefore, correlation sum simply tells the number of points within a distance “r” from the point x_p. The correlation integral is related to the standard correlation function as [12]:

$C\left(r\right)~r^D$ (4)

Thus, the following fractal model can be used to characterize the spatial distribution of earthquakes (using Equation (4)):

${\log }_{10}C\left(r\right)=C+D{\log }_{10}r$ (5)

where C and D are constants and D is the fractal dimension, more strictly, correlation dimension. For the estimation of the correlation fractal dimension (D₂), the generalized correlation integral C(r) is used [14]:

$D_2={\lim }_{r\to 0}\frac{{\log }_{10}C\left(r\right)}{{\log }_{10}r}$ (6)

where C(r) is the correlation integral and r is the radius of sphere of investigation. The correlation integral is given by:

$C\left(r\right)={\lim }_{N\to \infty }\frac{1}{N^2}{\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^{N}H\left(r-\Vert x_i-x_j\Vert \right)}}$ (7)

The estimator of correlation integral is correlation sum and is calculated as shown:

$C\left(r\right)=\frac{1}{N^2}{\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^{N}H\left(r-\Vert x_i-x_j\Vert \right)}}$ (8)

where N is the total number of earthquakes considered, x are the coordinates of the epicenters and H is the Heaviside step function.

$H\left(x\right)=\{\begin{array}{l}0\text{for}x<0\\ 1\text{otherwise}\end{array}$ (9)

Thus, the function H evaluates the number of points which have an inter-event distance less than or equal to r and C(r) evaluates the probability that a distance will separate two earthquake events less than r.

In case of an infinite fractal distribution, the ${\log }_{10}C\left(r\right)$ Vs ${\log }_{10}r$ plot is a straight line whose gradient is the fractal correlation dimension. It has been found that for large values of “r”, the gradient is quite low. On the other hand, for small values of r, the gradient is artificially high. These two conditions are known as “saturation” and “depopulation” respectively. So the fractal dimension is estimated by fitting a straight line to a subjectively chosen straight part of the curve. If r_n and r_s are the distances of depopulation and saturation respectively, then they can be determined by the following formulae:

$r_n=R{\left[\frac{1}{N}\right]}^{1/d},r_s=\frac{R}{2\left(d+1\right)}$ (10)

where R is the maximum inter-event distance, N is the total number of events and d is the embedding dimension. We can also safely start the scaling range from even r_n/3 [31]. For estimating the spatial correlation dimension ( $D_2^s$ ), the angular distance between two events is calculated by using a spherical triangle [9].

$\left|X_i-X_j\right|={\cos }^{-1}\left(\cos {\theta }_i\cos {\theta }_j+\sin {\theta }_i\sin {\theta }_j\cos \left({\varphi }_i-{\varphi }_j\right)\right)$ (11)

where ${\theta }_i$ and ${\theta }_j$ are the latitude of the i^th and j^th event and ${\varphi }_i$ and ${\varphi }_j$ are the longitude of the i^th and j^th event respectively whereas for the estimation of the temporal correlation dimension ( $D_2^t$ ), inter-event time between two aftershocks is computed as:

$\tau =t_i-t_j$ (12)

where t_i and t_j are the origin time of the i^th and j^th event respectively.

In this case, R is estimated to be 5.085˚ while r_n and r_s are found to be 0.272˚ and 0.84˚. Similarly the maximum inter-event time (T) is 730.37 days and t_n and t_s values considered are 2.087 days and 182.6 days respectively. For calculating R and T, the embedding dimensions are taken as 2 and 1 respectively.

By estimating the slopes of the ( ${\log }_{10}C\left(r\right)$ Vs ${\log }_{10}r$ ) and ( ${\log }_{10}C\left(\tau \right)$ Vs ${\log }_{10}\tau$ ) in their scaling range (Figure 8 and Figure 9), the spatial ( $D_2^s$ ) and temporal ( $D_2^t$ ) correlation dimensions are estimated to be 1.07 ± 0.028 and 0.395 ± 0.0027 respectively. The inter-event time has been taken in seconds for plotting ( ${\log }_{10}C\left(\tau \right)$ Vs ${\log }_{10}\tau$ ). We have also estimated the gradient for smaller values of r and $\tau$ to show the depopulation effect. To fit the best line over these curves, we have used the method of least squares.

4.4. Spatial Variation of the Fractal Correlation Dimension

In order to study the spatial variation of the fractal correlation dimension, we have chosen the same geographical area (84.3˚ to 87˚ longitude and 27.3˚ to 28.5˚ latitude) in which the spatial variation of the b-value was analyzed. Here also, the area is divided into 36 grids of 0.3˚ × 0.3˚ and $D_2^s$ is found for each of

(a)(b)

Figure 8. (a) “r” is the angular distance between two earthquake events. The gradient for lower values of r is high while that for the higher values of r is low. These conditions are known as depopulation and saturation; (b) The slope of the best fit line over the curve ( ${\log }_{10}C\left(r\right)$ Vs ${\log }_{10}r$ ) is $D_2^s=1.07\pm 0.028$ in the scaling range of 0.272˚ ≤ r ≤0.84˚.

(a)(b)

Figure 9. (a) “t” or “τ” is the time interval between two earthquake events. The gradient for lower values of t is high while that for the higher values of t is low. These conditions are known as depopulation and saturation. (b) The slope of the best fit line over the curve ( ${\log }_{10}C\left(\tau \right)$ Vs ${\log }_{10}\tau$ ) is $D_2^t=0.395\pm 0.0027$ in the scaling range of 2.087 days ≤ τ ≤ 182.6 days.

the grids having a minimum occurrence of 10 events. Mapping of spatial correlation dimension is done, in the similar manner as it was done for b-values, by interpolating the $D_2^s$ -values (Figure 10). The $D_2^s$ varies from 0.21 to 1.06. Though the lower values of $D_2^s$ are away from the mainshock, in most of the region it varies from 0.8 to 1.0 revealing that the epicenters are distributed in a linear trend. Therefore, the aftershocks are neither clustered at a single point nor homogenously distributed in the region, indicating a non-continuous aftershock activity that can be attributed to the heterogeneous structures, complex geology and tectonics of the area.

4.5. Temporal Variation in the b-Value and the Correlation Dimension

We have divided the whole earthquake catalogue, sorted by date and time, into 25 subsets between April 2015 and February 2017 and each subset consists of 14 aftershocks.

The b-value and spatial fractal correlation dimension is estimated for each subset. The time change of b-value and fractal dimension is plotted (Figure 11(a)). The mean of the first event’s year and the last event’s year in each subset is taken for the value on the time axis.

The b-value is proportional to the fractal dimension and is related to each other by the following equation [8]:

$D=2b$ (13)

Though Aki’s relation suggests a positive correlation between the fractal dimension and the b-value, some studies have shown a negative correlation between the two [9] [32] [33] [34]. Also, [9] showed that the fractal dimension involved in the Aki’s relation is capacity dimension.

To see the correlation between the b-value and the fractal dimension ( $D_2^s$ ), the fractal dimension is plotted against the b-value for 25 subsets (Figure 12) and with the help of the regression analysis by least squares method, the following first order relation between the b-value and fractal dimension is obtained:

Figure 10. The $D_2^s$ -value mapping is done by dividing the area into 36 grids of 0.3˚ × 0.3˚. In most of the region, $D_2^s$ varies from around 0.8 to 1.0 indicating a linear source in the source. zone.

(a)(b)

Figure 11. (a) The temporal variation of b and $D_2^s$ value is plotted by dividing the data into 25 subsets with 14 aftershocks in each subset. The value on the time axis is taken by taking a mean of the year of the first event in the subset and year of the last event in the subset; (b) The temporal variation of b and value for about first 60 days. The correlation between b and for the first 10 days has been positive and after that for about 50 days, it has been negative.

Figure 12. The values are plotted against b-values and the equation of the best fit line is determined:. The correlation coefficient is 0.26.

(14)

The correlation coefficient is 0.26 which indicates that the correlation is not very significant. The mean value of the is estimated to be 1.9 which is in support with the Aki’s relation. But it is statistically not a correct method to find the relation between any two quantities whose uncertainties vary and are not constant. Also the correlation coefficient is found to be decreasing as the number of subsets considered increases. It means that the correlation is significant for a shorter duration of aftershock activities just after the mainshock.

For the first 10 days the correlation is positive (Figure 11(b)) but after that for about 50 days the correlation is negative. The negative correlation between and b-value suggests increased complexity in the active fault network and thus, the rupture is accommodated along the fault planes of relatively larger surface area [1]. That is why, all the aftershocks of magnitude greater than 5.5 occurred within the first two months of the mainshock (Figure 3). After about 60 days the correlation has become positive (Figure 11(a)). Positive correlation suggests that the probability of large magnitude earthquakes decreases due to the increased fragmentation of the fault zone and thus, stress is released along the faults of smaller surface area [5].

5. Discussion and Conclusions

In this research paper, we have analyzed the 2015 Nepal Earthquake Aftershock activity. A total of 350 aftershocks with magnitude greater than or equal to 4.0 have been considered for the analysis. The geographical region for the analysis is from 26.096˚ to 30.581˚ in latitude and from 79.849˚ to 88.330˚ in longitude. The aftershock data encompasses a period from 25^th April, 2015 to 28^th February, 2017. The b-value is computed by the maximum likelihood method as well as the least squares method. But since the number of aftershocks is quite less, the maximum likelihood method is not very reliable.

The b-value is estimated to be 0.833 ± 0.035 by least squares method and 0.95 ± 0.05 by maximum likelihood method. The b-value is identical within uncertainties to the b-value 0.83 ± 0.05 estimated from 2102 events of magnitude greater than or equal to 2.4, between 1995 and 2015, in the same area [35]. A low b-value suggests that large amount of stress accumulates in this region. An alternative view concerning the b variation before a major earthquake results from natural time analysis, which is a recent procedure to analyze complex time series such as seismicity and other complex systems like heart rate fluctuations displaying 1/f noise and fractal dynamics [36] [37]. Such an analysis of seismicity reveals [38] that the decrease of the b-value reflects an increase of the order parameter fluctuations of seismicity.

The spatial and temporal correlation dimensions have been estimated using the correlation integral method. For the scaling range, 0.272˚ ≤ r ≤ 0.84˚, the slope of the regression line is. Thus, the fractal structure of the earthquake epicenters is characterized by a linear source consisting of heterogeneous distribution of points in the source zone. The temporal correlation dimension () is 0.395 ± 0.0027, produced by the linearity in the time scale of 2.087 days ≤ τ ≤ 182.6 days. It suggests a strong clustering in the time domain which further indicates a continuous aftershock activity that may be due to the large stress accumulation in this region.

The b-value and -value mapping is done to understand the spatial variation of these parameters. The b-value varies from 0.49 to 1.09 with higher values in the vicinity of the mainshock indicating low stress as much of the stress is relieved in the form of the seismic energy released by the mainshock. The value varies from 0.21 to 1.06, but in most of the region, it is in the range 0.8 - 1.0 which gives an idea that the distribution of the aftershocks is in the form of a linear structure. The correlation coefficient for b and values is 0.26 which indicates that there is no relation between them. But, the temporal variation of the b and -values suggests a positive correlation between the two, for the first 10 - 15 days and after that there is a negative correlation for about 45 - 50 days. The negative correlation indicates a larger possibility of the occurrence of large earthquakes [1]. That is why, almost all aftershocks of magnitude greater than or equal to 5.5 have occurred in the first 2 months of the mainshock. After about 60 days, correlation is almost positive and thus, the aftershock activity includes occasional occurrences of magnitude 4.0 - 5.0 which shows the decreased probability of the aftershocks.

The variation in the b-value depends on the levels of stress accumulated in the source zone or the degree of heterogeneity of the source medium [39], whereas the change in value depends on the degree of clustering of the epicenters or the complexity of the fault systems. The degree of heterogeneity of the fault structure decides the fractal dimension. In general, the b-value changes from b < 1 to b > 1 immediately before and after the occurrence of an earthquake of large magnitude (Von Seggern, 1980). It is mainly because of the increasing stress release along the main branch of the self-similar fault system, thus, leading to (,) before rupture [40]. Then, stress energy is subsequently released along the sub-branches of the main fault leading to (,). Figure 11 reveals that and immediately after the mainshock.

Acknowledgements

The author S.M. would like to thank the Science Academies’ Summer Research Fellowship Programme, which is jointly conducted by IASc (Bengaluru), INSA (New Delhi) and NASI (Allahabad), for giving an opportunity to work as a Summer Research Fellow. This work is a part of the Summer Research performed by S.M. at the CSIR Fourth Paradigm Institute, Bengaluru.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Mandal, P. and Rastogi, B.K. (2005) Self-Organized Fractal Seismicity and B Value of Aftershocks of the 2001 Bhuj Earthquake in Kutch (India). Pure and Applied Geophysics, 162, 53-72. https://doi.org/10.1007/s00024-004-2579-1
[2]	Mandelbrot, B.B. (1982) The Fractal Geometry of Nature. W. H. Freeman, New York, 468 p.
[3]	Mogi, K. (1962) Magnitude-Frequency Relation for Elastic Shocks Accompanying Fractures of Various Materials and Some Related Problems in Earthquakes. Bulletin of the Earthquake Research Institute, 40, 831-853.
[4]	Kagan, Y. (1993) Statistics of Characteristic Earthquakes. Bulletin of the Seismological Society of America, 83, 7-24.
[5]	Omori, F. (1894) On Aftershocks of Earthquakes. Journal of the College of Science, Imperial University, 7, 111-200.
[6]	Turcotte, D.L. (1986) A Fractal Model for Crustal Deformation. Tectonophysics, 132, 261-269. https://doi.org/10.1016/0040-1951(86)90036-3
[7]	Gutenberg, B. and Richter, C.F. (1944) Frequency of Earthquakes in California. Bulletin of the Seismological Society of America, 34, 185-188.
[8]	Aki, K. (1981) A Probabilistic Synthesis of Precursory Phenomena. In: Simpson, D.W. and Richards, P.G., Eds., Earthquake Prediction: An International Review, Maurice Ewing Ser. 4, AGU, Washington DC, 566-574.
[9]	Hirata, T. (1989) A Correlation between the B Value and the Fractal Dimension of Earthquake. Journal of Geophysical Research, 94, 7507-7514. https://doi.org/10.1029/JB094iB06p07507
[10]	King, G. (1983) The Accommodation of Large Strains in the Upper Lithosphere of the Earth and Other Solids by Self-Similar Fault Systems: The Geometrical Origin of b Value. Pure and Applied Geophysics, 121, 761-815. https://doi.org/10.1007/BF02590182
[11]	Main, I. (1996) Statistical Physics, Seismogenesis and Seismic Hazard. Reviews of Geophysics, 34, 433-462. https://doi.org/10.1029/96RG02808
[12]	Kagan, Y.Y. and Knopoff, L. (1980) Spatial Distribution of Earthquakes: The Two Point Correlation Function. Geophysical Journal of the Royal Astronomical Society, 62, 303-320. https://doi.org/10.1111/j.1365-246X.1980.tb04857.x
[13]	Hirata, T., Satoh, T. and Ito, K. (1987) Fractal Structure of Spatial Distribution of Microfracturing in Rock. Geophysical Journal International, 90, 369-374. https://doi.org/10.1111/j.1365-246X.1987.tb00732.x
[14]	Grassberger, P. and Procaccia, I. (1983) Measuring the Strangeness of Strange Attractors. Physica D, 9, 189-208. https://doi.org/10.1016/0167-2789(83)90298-1
[15]	Tosi, P. (1998) Seismogenic Structure Behaviour Revealed by Spatial Clustering of Seismicity in the Umbria-Marche Region (Central Italy). Annals of Geophysics, 41, 215-224.
[16]	Ghosal, A., Ghosh, U. and Kayal, J.R. (2012) A Detailed B-Value and Fractal Dimension Study of the March 1999 Chamoli Earthquake (Ms 6.6) Aftershock Sequence in Western Himalaya. Geomatics, Natural Hazards and Risk, 3, 271-278. https://doi.org/10.1080/19475705.2011.627380
[17]	Dewey, J.F., Shackleton, R.M., Chengfa, C. and Yiyin, S. (1988) The Tectonic Evolution of the Tibetan Plateau. Philosophical Transactions of the Royal Society A, 327, 379-413. https://doi.org/10.1098/rsta.1988.0135
[18]	Dewey, J.F. and Burke, K.C. (1973) Tibetan, Variscan, and Precambrian Basement Reactivation: Products of Continental Collision. The Journal of Geology, 81, 683-692. https://doi.org/10.1086/627920
[19]	Bilham, R., Larson, K. and Freymueller, J. (1997) GPS Measurements of Present-Day Convergence across the Nepal Himalaya. Nature, 386, 61. https://doi.org/10.1038/386061a0
[20]	McNamara, D.E., Yeck, W.L., Barnhart, W.D., Schulte-Pelkum, V., Bergman, E., Adhikari, L.B., Earle, P.S., et al. (2016) Source Modeling of the 2015 Mw 7.8 Nepal (Gorkha) Earthquake Sequence: Implications for Geodynamics and Earthquake Hazards. Tectonophysics, 714-715, 21-30. https://doi.org/10.1016/j.tecto.2016.08.004
[21]	Molnar, P. (1988) A Review of Geophysical Constraints on the Deep Structure of the Tibetan Plateau, the Himalaya and the Karakoram, and Their Tectonic Implications. Philosophical Transactions of the Royal Society A, 326, 33-88. https://doi.org/10.1098/rsta.1988.0080
[22]	Dahal, R.K. (2006) Geology for Technical Students.
[23]	Heim, A. and Gansser, A. (1939) Central Himalaya. Hindustan Publishing, Delhi.
[24]	Arita, K. (1983) Origin of the Inverted Metamorphism of the Lower Himalayas, Central Nepal. Tectonophysics, 95, 43-60. https://doi.org/10.1016/0040-1951(83)90258-5
[25]	Bilham, R. (1995) Location and Magnitude of the 1833 Nepal Earthquake and Its Relation to the Rupture Zones of Contiguous Great Himalayan Earthquakes. Current Science, 69, 101-128.
[26]	Sapkota, S.N., Bollinger, L., Klinger, Y., Tapponnier, P., Gaudemer, Y. and Tiwari, D. (2013) Primary Surface Ruptures of the Great Himalayan Earthquakes in 1934 and 1255. Nature Geoscience, 6, 71-76. https://doi.org/10.1038/ngeo1669
[27]	Wyss, M. (1973) Towards a Physical Understanding of the Earthquake Frequency Distribution. Geophysical Journal International, 31, 341-359. https://doi.org/10.1111/j.1365-246X.1973.tb06506.x
[28]	Aki, K. (1965) Maximum Likelihood Estimate of B in the Formula Log N = a_bM and Its Confidence Limits. Earthquake Research Institute, the University of Tokyo, 43, 237-239.
[29]	Scholz, C.H. (1968) The Frequency-Magnitude Relation of Microfracturing in Rock and Its Relation to Earthquakes. Bulletin of the Seismological Society of America, 58, 399-415.
[30]	Henderson, J.R., Barton, D.J. and Foulger, G.R. (1999) Fractal Clustering of Induced Seismicity in the Geysers Geothermal Area, California. Geophysical Journal International, 139, 317-324. https://doi.org/10.1046/j.1365-246x.1999.00939.x
[31]	Narenberg, M.A.A. and Essex, C. (1990) Correlation Dimension and Systematic Geometric Effects. Physical Review A, 42, 7065-7074. https://doi.org/10.1103/PhysRevA.42.7065
[32]	Henderson, J., Main, I.G., Pearce, R.G. and Takeya, M. (1994) Seismicity in North-Eastern Brazil: Fractal Clustering and the Evolution of the B Value. Geophysical Journal International, 116, 217-226. https://doi.org/10.1111/j.1365-246X.1994.tb02138.x
[33]	Öncel, A.O., Main, I., Alptekin, Ö. and Cowie, P. (1996) Temporal Variations in the Fractal Properties of Seismicity in the North Anatolian Fault Zone between 31 E and 41 E. Pure and Applied Geophysics, 147, 147-159. https://doi.org/10.1007/BF00876441
[34]	Wang, J.H. and Lee, C.W. (1996) Multifractal Measures of Earthquakes in West Taiwan. Pure and Applied Geophysics, 146, 131-145. https://doi.org/10.1007/BF00876673
[35]	Adhikari, L.B., Gautam, U.P., Koirala, B.P., Bhattarai, M., Kandel, T., Gupta, R.M., Hoste-Colomer, R., et al. (2015) The Aftershock Sequence of the 2015 April 25 Gorkha-Nepal Earthquake. Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society, 203, 2119-2124. https://doi.org/10.1093/gji/ggv412
[36]	Varotsos, P.A., Sarlis, N.V., Skordas, E.S. and Lazaridou, M.S. (2007) Identifying Sudden Cardiac Death Risk and Specifying Its Occurrence Time by Analyzing Electrocardiograms in Natural Time. Applied Physics Letters, 91, Article ID: 064106. https://doi.org/10.1063/1.2768928
[37]	Varotsos, P.A., Sarlis, N.V. and Skordas, E.S. (2011) Scale-Specific Order Parameter Fluctuations of Seismicity in Natural Time before Mainshocks. EPL (Europhysics Letters), 96, Article No. 59002. https://doi.org/10.1209/0295-5075/96/59002
[38]	Varotsos, P.A., Sarlis, N.V. and Skordas, E.S. (2012) Order Parameter Fluctuations in Natural Time and B-Value Variation before Large Earthquakes. Natural Hazards and Earth System Sciences, 12, 3473-3481. https://doi.org/10.5194/nhess-12-3473-2012
[39]	Mogi, K. (1967) Earthquakes and Fractures. Tectonophysics, 5, 35-55. https://doi.org/10.1016/0040-1951(67)90043-1
[40]	Main, I.G., Meredith, P.G. and Jones, C. (1989) A Reinterpretation of the Precursory Seismic b Value Anomaly from Fracture Mechanics. Geophysical Journal International, 96, 131-138. https://doi.org/10.1111/j.1365-246X.1989.tb05255.x