Modelling of the Kopili Fault based on slip rate, moment rate and seismic activity in Mikir Hills Plateau of Northeastern India

ABSTRACT The recurrence period of maximum magnitude earthquake in a seismogenically active formation along the Kopili Fault has been estimated having adequate dependence on the slip rate, moment rate and its seismicity pattern. Here, the frequency–magnitude cum fault area–maximum magnitude relations play a key role with input parameters pertinent to the Kopili Fault zone. Subsequently, where the fault slip rate estimates are not available, the seismic activity is studied from the seismic moment release. The results of this study show that the return period has a strong relation with the fault length, slip rate, strain drop and rigidity. This study ascertains the activity rate in terms of the return period as ∼50 ± 5 years with the moment release of 2.12E+23 dyne-cm from the most active 80-km fault length considering Mw 5.5 as reference magnitude under the Kopili Fault zone that may produce a maximum magnitude of Mw ∼ 7.6. Finally, we conclude that these models can be used to study the rate of seismicity of the active faults in Northeast India which will provide us prime inputs for seismic hazard analysis of the region and is especially significant for estimating expected return period for poorly known faults or blind faults that lack surface expression.


Introduction
Estimation of the size of the largest earthquake and subsequent prediction of ground motions a particular fault can generate is a prime factor in seismic hazard studies (Anderson et al. 1996). Since the size of an earthquake depends on the return period and the tectonic environment (Kanamori & Allen 1986;Scholz et al. 1986), it allows evaluation of future earthquake potential of a fault mainly from the fault rupture parameters (e.g. Tocher 1958;Iida 1959;Chinnery 1969;Wells & Coppersmith 1994;Hanks & Bakun 2008). Crustal deformation at any specific site is dependent on the slip rate (mm/year) of an active fault and in turn allows to estimate the return period pertinent to the rupture area. The Kopili Fault in Northeast India, being an intraplate earthquake source zone, nearly 300-km long NW-SE trending fault, is characterized as a potential zone for a large impending earthquake. The activity of the Kopili Fault is estimated through seismicity pattern over last few decades, extrapolating this into the future occurrence of probable maximum magnitude earthquake utilizing the statistical models by Bungum (2007). Although this fault generated two large earthquakes besides several intermediate and small-magnitude earthquakes (Baruah et al. 1997), the probable time period for next large earthquake is unknown.
This study tries to quantify the seismic potential and activity rate in terms of the return period of the Kopili Fault based on the relations among earthquake size (magnitude or seismic moment), seismicity pattern, slip rate, moment rate, the rupture length, area and displacement of the Kopili Fault constrained by the maximum magnitude for seismic hazard assessment of the region.

Study area and database
The study area pertains to the Kopili Fault and its vicinity which lies in between the latitudes 25.00 N to 27.80 N and longitudes 91.00 E to 94.00 E (Figure 1) between the Shillong Plateau and the Mikir Hill massif. The Kopili is the main river of the Kopili Valley. Geologically, Kopili Valley area comprises Neogene-Quaternary sediments which were deposited directly over the Archean basement. The Kopili Fault is a NW-SE trending strike-slip fault where intense seismic activity occurs down to a depth of about 50 km beneath the Kopili Fault, and the activity continues to the Figure 1. Location map of the study area (within latitude 25.00 -27.80 N and longitude 91.00 -94.00 E, shaded in grey) showing major earthquakes along with tectonic features of Northeast India region (Murthy et al. 1969;Nandy 2001;Baruah & Hazarika 2008). Two great earthquakes (M s » 8.7) in the region are shown by stars, and the large earthquakes (M b 5.5) by circles. The digital seismic stations are shown by triangles.
Main Central Thrust (MCT) in the Bhutan Himalaya (Baruah et al. 1997;Bhattacharya et al. 2002Bhattacharya et al. , 2008Bhattacharya et al. , 2010Bora & Sitaram 2007;Kayal 2008;Kayal et al. 2006Kayal et al. , 2012. Two large earthquakes occurred in this fault. One of these events occurred in 1869 (M 7.7) in the southeastern end of the fault transgressing Naga-Disang thrust while the other event (1943; M 7.2) occurred farther north of 1869 event within a span of about 75 years (Kayal 2008). The Kopili Fault zone is, under compressional stress from the Indo-Burma arc to the east and from the Himalayan arc to the north, characterized by transverse tectonics.
A database containing 10,109 events for the period 1964-2015 is taken into consideration for the study. The hypocentral database is prepared from the seismological bulletin published annually by CSIR North East Institute of Science and Technology (NEIST)-Jorhat and CSIR National Geophysical Research Institute (NGRI)-Hyderabad. Phase data bulletins of Indian Meterological Department (IMD), Shillong, Manipur University, Mizoram University and Gauhati University are also taken into consideration wherever needed. The threshold magnitude is, however, estimated to be 2.5 for this region. The epicentral plot along with focal mechanism solutions (FMS) (Figure 2 , 29, 20, 26, 9, 13, 2, 24 and 34) of focal mechanism solutions could be associated with the Kopili Fault. Out of these, 7 numbers (event numbers 1, 29, 20, 26, 13, 2 and 34) of events are characterized by strike-slip faulting, 1 number (event number 9) is characterized by thrust with strike-slip component and 1 event (event number 24) indicates thrust component. Event number 26 characterizes the Kopili Fault with strike-slip mechanism having the fault plane very clearly oriented along NW-SE. Similarly, event numbers 2, 13, 34 and 29 clearly ascertain the trend of the fault plane which is pertinent to the regional trend of the Kopili Fault. Although the entire length of the Kopili Fault is 300 km, but primarily, we consider the 80-km fault length to be most active part for this study since the length of the most active part of the fault is a critical input parameter ( Figure 2). In order to identify the most active part, seismicity pattern over this entire fault length, segregated in three different blocks, is observed since the occurrence of past major earthquakes till date, to ascertain the maximum earthquake potential the fault can generate as well ( Figure 2). The middle part of the Kopili indicated as Block 2 (length 80 km and width 40 km) is considered to be the most active. In addition, seismic moment release is relatively large during last two decades  in comparison to earlier period  in the study area ( Figure 2(B-i)).
Step-wise increment in seismic moment followed by abrupt increase in seismic moment is observed. Three rises are very prominent which correspond to the three large events (1998, M w 5.5; 2009, M w 6.0 and 2013, M w 5.9) in the region. The temporal distribution of events (greater than M w » 3.5) (Figure 2(B-ii)) indicates occurrence of three events greater than M w » 5.5 since 1998 till 2015 and last two events (M w 5.9 and 6.0) merely occurred within a span of 4 years. We have studied the cumulative number of events vs. time for a single magnitude range (M w > 3.5) following Sharma et al. 2013 (Figure 2(B-iii)). Negligible changes in the slope of the curve are perceived for the felt events. However, remarkable changes are found for events occurred in 2009 (M w 6.0) and 2013 (M w 5.9). We have considered here that the higher range of the expected magnitude the Kopili Fault can produce (M max ) is about M w 7.6 and a lower range of reference magnitude of about M w 3.5. This lower bound or the reference magnitude is the starting point for seismic hazard integration. The slip rate for the present study is considered to be 1.3 mm/year as per the Northeast India permanent network baselines (Jade et al. 2007). The slip rate associated to the Kopili Fault is contributing to accumulation of seismic moment to the order of »70.74 £ 10 15 N m/year (Barman et al. 2016). The input parameters (Tables 1 and 2) for the fault characteristic are estimated from existing literatures (Hanks & Boore 1984;Wells & Coppersmith 1994;Jade et al. 2007;Angelier & Baruah 2009;Bhattacharya et al. 2010;Kayal et al. 2006Kayal et al. , 2012.

Methodology
The seismic activity rate of the fault has been studied in two ways. Primarily, related to the slip rate estimation from published models; in this method, the size of the fault constrained by occurrence of the maximum magnitude is an essential component. In addition, the frequency magnitude relation along with the moment-magnitude and fault area-maximum magnitude relations are the typical controlling parameters. Second, where direct fault slip rate estimates are not available, the seismic activity is studied from the seismic moment release. Simultaneously, predefined b-value and the magnitude-frequency distribution are two important factors. Since the seismic moment increases strongly with magnitude, the consideration of probable maximum magnitude estimate is a critical factor in this process. We have used the computer algorithm developed by Bungum (2007) to estimate various related parameters in our present study.
While considering the fault activity from slip rates, the seismicity of an area or a fault usually follows the classical cumulative Gutenberg-Richter relationship: where N is the number of earthquakes equal to or above magnitude M, and a and b are constants.
Here, a determines the absolute activity level while b determines the slope of the curve, or the ratio between the smaller and the larger earthquakes. Fault length range: min, max and step length (km) 5.0, 80.0, 2.5 The minimum fault length is taken to be 5 km as a range of fault lengths are tried since the maximum is taken as 80 km The maximum fault length is taken to be 80 km because it is observed to be the most active part of the fault Fault slip to fault length ratio 1.0E¡04 10 Fault length to fault width factor (L/W) 2.0 (since the length of the fault is 80 km as the active part and the width is considered to be 40 km) 11 Shear modulus (GPa) 30  The cumulative occurrence relationship (Chinnery and North 1975) used in this study can be expressed as where H is the Heaviside step function. Normally, the b-value for the source area is considered to be 1.0 which implies a factor of 10 reduction in the number of events per magnitude unit.
The inferences on seismic activity rates (N-values) involves the seismic moment. The seismic moment, M o , is the most physically meaningful way to describe the size of an earthquake in terms of static fault parameters (Youngs & Coppersmith 1985;Aki & Richards 2002).
where m is the rigidity, or shear modulus (usually taken as ̴ 3 £ 10 11 dyne/cm 2 ); D is the total average displacement (or slip) across the fault and A is the rupture area where A = LW, where L is the length of the fault and W is the width of the fault (Aki 1966). The total seismic moment rate M T o or the rate of seismic energy release along a fault which also determines the time for the slip at each point to be completed is estimated by Brune (1968) as where S is the annual slip rate, ignoring possible aseismic creep (Anderson et al. 1993). M T o should be averaged over several cycles of large earthquakes for this relation to be valid (Anderson 1979;Bungum 2007). Based on the Equations (2)-(4), Anderson and Luco (1983) deduced the following relationship for the determination of the number of earthquakes N above the lower bound magnitude: Seismic moment rates determined from fault slip rates in a region may be directly compared with seismic moment rates based on seismicity data (Doser & Smith 1982). Seismic moment is translated to earthquake magnitude according to the following relation (Kanamori & Anderson 1975): where the parameter d is the magnitude scaling coefficient. In addition to Equation (5), Anderson and Luco (1983) as advocated by Bath (1978), Anderson (1979) and Berril and Davis (1980) developed similar relationships for two other recurrence models, first the incremental relation This relation has a more smooth (rounded) transition towards N(M max ) = 0 in its cumulative version.
This is the occurrence relation (2) in Anderson and Luco (1983) or Model 2 of the present paper.
The occurrence relation (3) in Anderson and Luco (1983) or Model 3 of the present paper was proposed by Main and Burton (1981): This relation has even gentler transition towards N(M max ) = 0 than Equation (7) and the corresponding N-relation is thus For Equations (5), (8) and (10) of Anderson and Luco (1983), the annual number of events N above a given magnitude M is the inverse of the return time T in years, for the same magnitude, i.e. T = 1/N (Bungum 2007).
In addition to these three model of Anderson and Luco (1983), a fourth model, Model 4 developed by Youngs and Coppersmith (1985) is considered to find out the annual number of events. This model is very close to Model 2 which is as follows: where m 0 is some arbitrary reference magnitude, m = b ln(10), A f = LW is the fault area and m u is an upper bound magnitude. The seismic impulse response or seismic moment is the most important parameter for expressing the level of seismicity in a particular fault zone. Subsequently, the b-value which is the principal factor in this method has been constrained from the regional seismic activity as deduced from past and recent seismicity. Essentially, this means that any fault is assumed to exhibit a certain moment rate assessed in one way or other based on a predefined b-value (Bungum 2007). Lower b-value indicates lower seismicity rate, and more stress accumulation which is likely to produce a higher magnitude earthquake with higher seismic moment release. The a-value is calculated iteratively based on integrating the seismic moment over the entire magnitude rate. The input parameters for the calculation of the moment rate are given in Table 2.
On the other hand, seismic moment release determined from fault slip rates in a region may be directly compared with seismic moment rates based on seismicity data (Doser & Smith 1982) which is translated into earthquake magnitude according to the relation by Kanamori and Anderson (1975). These models are used when the dimensions of the fault, maximum magnitude, coseismic slip and slip rates are known or are assumed.

Fault activity, return period and slip rate
The four different models (Equation (5), (8) and (10) of Anderson and Luco (1983), Equation 11 of Youngs and Coppersmith (1985)) and the average of these along with error bars for one standard deviation for a fault-producing magnitudes larger than M w 5.0 are plotted (Figure 3). The log-linear plot shows the relation between the return period for reference magnitude of M w 3.5 associated to a fault length of 80 km. It is seen that the fault length and the maximum magnitude show a linear relation with the return period for a reference magnitude of M w 3.5 concerning all the models indicated including the average model. Simultaneously, Figure 3 indicates that the increase in fault length increases the probability of producing a high-magnitude earthquake. Moreover, we have observed that while fault length increases, the return period decreases which indicates an increase in the seismic activity as reflected in Bungum (2007) (inset in Figure 3) as well. Thus, the earthquake frequency increases strongly (and the return period decreases) with fault length keeping the slip rate constant. On the other hand, keeping the parameters constant as in Figure 3 and taking into consideration the average of the four models, i.e. model 5 as a function of the reference magnitude of M w 3.5, 4.5 and 5.5, we observe that with the increase of the reference magnitude, the return period becomes higher ( Figure 4) even though the graph shows decrease in return period with the increase in the fault length. Thus, for the 80-km-long Kopili Fault, the return period for the reference magnitude M w 4.5 is 5 years while the return period for the reference magnitude M w 3.5 is around half a year. However, it is seen that variation of fault length does not show much influence on the return period as per the model considered. In normalcy, the longer the fault length, the higher the magnitude of earthquake it can produce. Here, it may be mentioned that even though the pattern of figures as in Bungum (2007) (Figure 4: inset) are almost similar, the output parameters are absolutely different. Figure 5 indicates the log-linear sensitivity of the activity rate (in terms of the return period) with respect to slip rate (the measure of slip of a fault from geodetic measurement) and change in fault length. In this graph, the return period (activity rate) covers five orders of magnitude which indicates a strong sensitivity both to the fault length and the slip rate. Here, the return period is inversely proportional to the slip rate, i.e. the return period decreases with the increase in the slip rate of the fault. In case of the Kopili Fault, higher slip rates has increased the rate of activity ( Figure 5). This indicates that a 80-km-long Kopili Fault does have higher rate of activity with low return period in comparison to a fault length of 5 km as considered for the same fault. A smaller fault length of 5 km is likely to slip at a higher rate than a longer fault of about 80 km. Slow slip rates in a relatively longer fault implies that seismic moment budget is biased towards the production of large but infrequent earthquakes. Thus, from the model of slip rate, it can be seen that an 80-km length of seismogenic zone of the Kopili Fault having the potential to produce an earthquake of about M w 7.6 with b value of 1.0, slip rate of 1.3 mm/year and a reference magnitude of 5.5 corresponds to a probable return period equivalent to 50 years (Figure 4). On the contrary, the sensitivity of the activity rate (in terms of return period) with respect to the ratio in-between the slip and the fault length or the strain drop for different fault lengths of 20 and 80 km for the Kopili Fault are analyzed ( Figure 6). However, these two curves as in Figure 6 represent fixed values of fault length and M max . Consequently, it is observed that the relation between the return period and the strain drop is not linear in log-log space domain. Notably, the return period or the earthquake frequency is directly proportional to the strain drop. Eventually, the lesser the  strain drop, the more is the seismic activity. Since the slip rate and the seismic moment rate are also constant, the change in the earthquake frequency (or return period) is therefore related directly to the change in the stress drop, being proportional to the strain drop. The higher value of strain drop, i.e. from 10 ¡1 to 10 ¡5 may be due to intraplate conditions and thereby more moment release. The final sensitivity test as derived from the slip rate is shown in Figure 7 which includes the relation between the return period and the rigidity (resistance to the shear deformation) modulus. It can be inferred from the graph that as the rigidity increases, more moment is released, thereby reducing  the return time period between earthquakes of a given magnitude. Thus, a fault of 80-km length will have less rigidity and lesser reduction in the return time because of more seismic activity than a fault of 20-km length.
On the other hand, the relation between the fault seismic activity and the maximum magnitude the fault can produce indicates that as the return period of a higher magnitude earthquake becomes lower, the probability of maximum magnitude the Kopili Fault can produce becomes low, because of the increase in the seismic activity (Figure 8).
A comparison between fault length (km) and expected magnitude (M w ) among Anderson et al. (1996), Wells and Coppersmith (1994), Bungum (2007) and the present study ( Figure 9) indicates Figure 9. Relationship between magnitude, rupture length and slip derived from various models (indicated along heavy and dotted lines). The relationship as observed in the present study shows relatively underestimated parameters for lower rupture length. that as the fault length increases, the magnitude of the expected earthquake also increases. However, the magnitude of the expected earthquake depends on the slip rate. It is seen from Figure 9 that as the slip rate increases, the magnitude of the expected earthquake decreases. Noticeably, this study shows a trend which is in conformity with the global trend done by various workers. This establishes the robustness of the approach and the sensitivity tests made towards validation and understanding the fault kinematics.

Fault activity from seismic moment
The dependence of seismic moment is certainly varying with reference to estimates of maximum magnitude. As we observe, the fault target moment release (total seismic moment radiated from the target fault) and the area annual moment release both increase as the maximum fault magnitude (M max ) increases at constant b and a values (Figure 10), thereby showing the increase in associated fault length (Table 3). A maximum magnitude of M w 7.6 might have influenced a length close to 80 km of the Kopili Fault, considered to be most active, which in turn will release a moment rate of 2.12E+23, area annual moment of 4.23E+23 and fault target moment of 2.12E+23 (Table 3 and Figure 10). On the other hand, when the M max is kept constant, it is seen that moment rate, area annual moment release and the fault target moment release increase with a decrease of b and a values while the fault length remains constant (Table 4). On the contrary, Figure 11 shows the moment rate, area annual moment release and fault target moment release are inversely proportional to the b-value of the area. Thus, the Kopili Fault having a b-value of 0.6 will produce a moment rate of 3.37E+24, while a b-value of 1.00 produces a moment rate of 2.12E+23, which indicates the seismic moment release is dependent on the b-value and on the maximum magnitude the fault can produce.

Conclusion
The database reveals a specific range of input model parameters, e g. fault rupture length, its activity, seismic moment release, reference magnitude and maximum magnitude pertinent to the study region which incidentally influence the estimation of recurrence period and other associated parameters. Over the number of iterations, uncertainty ranges could be specified for all these parameters. In order to minimize the uncertainty, a long-term dataset is used since 1964 till 2015 which allowed us to assess the best possible input parameters. The slip rate 1.3 mm/year has been recognized based on the adequately developed and published geodetic information along the Kopili Fault area. Intense seismotectonics of the Kopili Fault as a function of both space and depth distribution allow us to demarcate the approximate seismogenic zone which is 80 km in length and 40 km in width, oriented transversely to the E-W Himalayan trend. It is inferred from the study that an increase in fault length increases the probability of producing a high-magnitude earthquake. Simultaneously, it is observed that with increase of fault length, the return period decreases which indicates an increase in the seismic activity decreasing the probability of producing a large-magnitude earthquake. The empirical equation involved shows a strong correlation among the reference magnitude, slip-rate and return period which enabled us to use these reference magnitudes and associated parameters. The estimates are well determined after numbers of iterations. Thus, a reference magnitude of M w » 5.5 corresponds to a probable return period of 50 § 5 years (Figure 4) for the same seismogenic zone of the Kopili Fault at a very high confidence level for any relationship having the potential to produce an earthquake of about M w » 7.6 with b value of 1.0 and slip rate of 1.3 mm/year. Moreover, the maximum seismic moment release from 80-km length of the fault validates the return period to a substantial extent.
Higher value of strain drop observed in this study, i.e. from 10 ¡5 to 10 ¡1 , may be due to intraplate conditions and thereby more moment release. As the rigidity increases, more moment is released thereby reducing the return time between earthquakes of a given magnitude. The geodetic velocity with higher strain map around the Kopili Fault ( Figure 12) revalidates the inferred seismic paramters.
In case of the Kopili Fault, the moment rate, area annual moment release and the fault target moment release increase with constant b and a values, thereby showing the increase in associated fault length. Considering the coefficients, it is observed that a maximum magnitude of M w » 7.6 with annual moment release 4.23E+23 may influence a fault length of about 80 km in the Kopili Fault area.
Thus, this result further helps us to study activity rate in terms of the return period and the moment release of the Kopili Fault zone based on the assumption of the maximum magnitude the fault can produce. These models can be used to study the rate of seismicity of the other active faults in Northeast India which will provide us prime inputs for seismic hazard assessment of the region. Finally, we conclude that this dependence may be used for most sites and is especially significant for estimating expected return period for poorly known faults or blind faults that lack surface expression.