Artificial neural networks ( ANN ) and stochastic techniques to estimate earthquake occurrences in Northeast region of India

The paper presents the probability of earthquake occurrences and forecasting of earthquake magnitudes size in northeast India, using four stochastic models (Gamma, Lognormal, Weibull and Log-logistic) and artificial neural networks, respectively considering updated earthquake catalogue of magnitude Mw ≥ 6.0 that occurred from year 1737 to 2015 in the study area. On the basis of past seismicity of the region, the conditional probabilities for the identified seismic source zones (12 sources) have been estimated using their best fit model and respective model parameters for various combinations of elapsed time (T) and time interval (t). It is observed that for elapsed time T=0 years, EBT and Kabaw zone shows highest conditional probability and it reaches 0.7 to 0.91 after about small time interval of 3-6 years (2014-2017; since last earthquake of Mw ≥ 6.0 occurred in the year 2011) for an earthquake magnitude Mw ≥ 6.0.Whereas, Sylhet zone shows lowest value of conditional probability among all twelve seismic source zones and it reaches 0.7 after about large time interval of 48 years (year 2045, since last event of Mw ≥ 6.0 occurred in the year 1999). While for elapsed time up to 2016 from the occurrence of the last earthquake of magnitude Mw ≥ 6.0, the MBT and MCT region shows highest conditional probability among all twelve seismic source zones and it reaches 0.88 to 0.91 after about 6-7 (2022-2023) years and in the same year (2022-2023) Sylhet zone shows lowest conditional probability and it reaches 0.14-0.17. However, we have proposed Artificial Neural Network (ANN) technique, which is based on feedforward backpropagation neural network model with single hidden layer to estimate the possible magnitude of future earthquake in the identified seismic source zones. For conditional probability of earthquake occurrence above 0.8, the neural network gives the probable magnitude of future earthquake as Mw 6.6 in Churachandpur-Mao fault (CMF) region in the years 2014 to 2017, Mw 6.8 in Myanmar Central Basin (MCB) region in the years 2013 to 2016, Mw 6.5 in Eastern Boundary Thrust (EBT) and Kabaw region in the years 2015-2018, respectively. The epicentre of recently occurred 4 January 2016 Manipur earthquake (M 6.7), 13 April 2016 Myanmar earthquake (M 6.9) and the 24 August 2016 Myanmar earthquake (M 6.8) are located in Churachandpur-Mao fault (CMF) region Myanmar Central Basin (MCB) region and EBT and Kabaw region, respectively and that are the identified seismic source zones in the study area which show that the ANN model yields good forecasting accuracy.


Introduction
The earthquake forecasting gives the probability of time, location and magnitude of occurrence of next earthquake, which is necessary to understand the seismic hazard of any region [Parvez and Ram, 1997].An earthquake is one of the most precarious natural hazards, which causes the sudden violent movement of the earth's crust, resulting in huge damage to structures and loss of human lives [Milne 1896].The statistical approach based on trends of earthquake such as seismicity patterns, seismic gaps, and characteristics of earthquake is the most appropriate and widely used method for estimation of the seismic hazard in any region.The reliable estimation of seismic hazard requires the prediction of time, location and magnitude of future earthquake events [Anagnos and Kiremidjian, 1988].
The probabilities of occurrence of future earthquake can be estimated by using the different stochastic models.Utsu [1972], Hagiwara [1974], and Rikitake [1974] have proposed a statistical probabilistic approach for forecasting of future earthquake for a particular region.Their models are based on the concept that the earthquake is a renewal process, in which just after an earthquake event the probability of occurrences of next earthquake is initially low.And it goes on increasing over a long period of time until the occurrence of any next earthquake event.Over this long period of time the accumulation of strain energy prepares the fault to release energy in the next earthquake.However, the Poisson distribution model is more appropriate and widely used for seismicity studies [Cornell 1968, Caputo 1974, Gardner and Knopoff 1974, Shah and Movassate 1975, Cluff et. al. 1980] with assumption that the earthquake occurrences are independent in space and time.But the probabilities of occurrence are dependent on the size and time elapsed since the last major earthquake.In the past, several statistical distri-bution models have been proposed for forecasting of future earthquakes including double exponential [Utsu 1972b], Gaussian [Rikitake 1974], Weibull [Hagiwara 1974, Rikitake 1974], Gamma [Utsu 1984] and Lognormal [Nishenko and Buland 1987].This type of study has been carried by many researchers for different areas [Parvez andRam 1997, 1999;Yilmaz et al. 2004, Yilmaz and Celik 2008, Yadav et al. 2008, 2010, Sil et al. 2015].
Recently, Sil et al. [2015] have estimated conditional probabilities for forecasting of future earthquakes (Mw > 6.0) in the northeast region of India using earthquake catalogue from year 1737 to 2011.They have divided their study area (longitudes 86.20°E-97.30°Eand latitudes 18.40°N-29.00°N)into six regions based on seismic event distribution pattern and orientation of all seismic sources or faults.They have suggested lognormal model is the most suitable model for analyzing earthquake occurrence in northeast India.In their study Indo-Burma Range (IBR) and Eastern Himalaya (EH) show > 90 % chances of occurrence for an earthquake Mw > 6.0 in the 5 years period (2012)(2013)(2014)(2015)(2016)(2017).Hence, it needs to be reworked precisely with an updated earthquake catalogue occurred after year 2011.
However, the earthquake magnitude is the most important parameter which describes the severity of an earthquake.It has been observed that higher magnitude earthquakes cause greater amount of fatalities, injuries and devastations.Therefore, estimation of magnitude size of future earthquake is important.Many researchers [Wong et al. 1992, Alarifi et al. 2012, Niksarlioglu and Kulahci 2013, Mahmoudi et al. 2016, Narayanakumar and Raja 2016] have used Artificial Neural Network (ANN) technique to estimate the intensity and magnitude of earthquakes because of its capability of providing higher forecasting accuracy and capturing non-linear relationship than other proposed models.Wong et al. [1992] have presented Modified Mercalli Intensity (MMI) forecast model for regional seismic hazard assessment in California using Artificial Neural Network (ANN) method.The authors have used 'three-layer' back propagation neural network with Hyperbolic Tangent Sigmoid transfer function and Normalized Cumulative-Delta learning rule.On the other hand, Alarifi et al. [2012] have used Artificial Neural Network (ANN) method to estimate the magnitude of future earthquakes in northern Red Sea area.The authors have used different feed forward neural network configurations with multi-hidden layers along with two transfer functions namely Log sigmoid and Hyperbolic Tangent Sigmoid.To measure the performance of the neural networks, authors have estimated Mean Absolute Error (MAE) and Mean Squared Error (MSE).Niksarlioglu and Kulahci [2013] have proposed a 'three-layer' neural network model with Levenberg-Marquardt learning to estimate the earthquake magnitude for East Anatolian Fault Zone (EAFZ), Turkiye using location, radon emission and environmental parameters as input neurons for the neural network.Similarly, Narayanakumar and Raja [2016] have proposed a 'three-layer' feedforward backpropagation neural network model to predict magnitude of earthquakes in the region of Himalayan belt.Location, energy released and seismicity parameters have been taken as input elements and magnitude as an output element to prepare the architecture of neural network.
In the present study, the North East India is being considered to forecast the future earthquake occurrences considering the past historical event data collected since 1737-2015 from available national and international various seismological agencies such as USGS, NGRI, and IMD.In the present work, the study area has been divided into 29 active seismic source zones based on seismicity considering the latest catalogue, distribution pattern of seismic events and orientation of seismic sources or faults.In this work probabilistic recurrence of earthquake (Mw > 6.0) has been estimated using four probability distribution models, namely, Gamma, Lognormal, Weibull and Log-logistic in northeast India.However, we have used the Artificial Neural Network (ANN) technique to estimate the magnitude of future earthquakes in the study region.Future earthquakes have been estimated using different stochastic models for different seismic source zones in this research work.The proposed Artificial Neural Network (ANN) technique is based on feedforward backpropagation neural network model with single hidden layer.The results presented in this study would be helpful for long-term disaster mitigation planning of the region in future.

Study area
For forecasting of future earthquakes, the study area has been chosen between latitudes 19.345° to 29.431°a nd longitudes 87.590° to 98.461°.However, the northeast region of the India is the most seismically active region in the world.This region is situated in the zone-V with a zone factor 0.36g on the latest version of seismic zoning map of India, given in the earthquake resistant design code of India [IS 1893:2002 (Part 1)], which expects the highest level of seismicity in the country.There are many identified geologically active faults in this region, whose activities make the region seismically vulnerable and cause huge destruction of buildings and other structures referred as seismic risk.To overcome the seismic risk first the systematic evaluation of the seismic hazard ZAROLA AND SIL 2 tudes 19.345° to 29.431°.Sil et al. [2015] have collected preinstrumental (historic) seismic events from different literatures [Oldham 1883, Basu 1964, Rastogi 1974, Chandra 1977, Dunbar et al. 1992, Bilham 2004] and instrumental data from various national and international seismological agencies, such as IMD, BARC, NGRI (National agencies) and USGS, ISC, COSMO Virtual Data Center (International agencies).They have removed all the repeated events (from different seismological agencies) and declustered the foreshocks and aftershocks using methodology given by Gardner and Knopoff [1974] and modified by Uhrhammer [1986].This methodology assumes that the temporal and spatial distribution of foreshocks and aftershocks are dependent on the size of the main event.The size of space and time window to decluster foreshocks and aftershocks are given as e-1.024+0.804*Mw(km) and e -2.87+1.235*Mw(days) respectively.Sil et al. [2015] have converted different types of magnitude scale into a common magnitude scale i.e. moment magnitude (Mw) using relationships proposed   by Sitharam and Sil [2014].With addition to this earthquake data base, we have collected earthquake events (magnitude ≥ 4.0) from year 2012 to 2015 (4 years) in the study area from U.S. Geological Survey.All collected events (from year 2012-2015) have been homogenized using the relationships proposed by Sitharam and Sil [2014].
Where Mw is the moment magnitude, Mb is the body wave magnitude, Ml is the local magnitude and Ms is the surface wave magnitude.
This newly collected data set (from year 2012-2015) has been declustered using the same method as used by Sil et al. [2015].A complete earthquake catalogue from year 1737 to 2015 has been prepared after adding two earthquake catalogues, first from year 1737 to 2011 [Sil et al. 2015] and second from year 2012 to 2015 (prepared in this study).Some scattered events are outside the proposed source zones and study area and hence are not taken into account for further study.Total 2508 earthquake events (Mw ≥ 4.0) with range of longitudes 87.599° to 97.918° and latitudes 19.392° to 29.026° have been compiled in data base.
A histogram of earthquake data (Mw ≥ 4.0) that occurred from year 1737 to 2015 in the study area has been prepared (see Figure 2).The study area has been divided into 29 seismic source zones (see Figure 1.a and 1.b ) on the basis of distribution patterns of seismic events (Mw ≥ 4.0) and orientation of seismic sources or faults after superimposing of total 2508 earthquake events (Mw ≥ 4.0) on the digitized tectonic map (Seismotectonic Atlas of India [SEISAT], 2000).A total of 52 scattered events are outside the proposed source zones but within the study area, hence are not taken into account for further study.
The general property of the size distribution of earthquakes that smaller earthquakes occur more frequently than the bigger earthquakes is well known.It is observed that model parameters give a higher probability of occurrence caused by an extremely short recurrence interval (for smaller earthquakes).Therefore, for long term forecasting only main events of magnitude Mw ≥ 6.0 have been taken.
After declustering the catalog, Sil et al. [2015] have observed that some of the main shock events (Mw ≥ 6.0) are available in the same year (only there is a month's difference) in the same identified seismic source zone.Thus, an extremely short recurrence interval makes the same problem as a problem of higher probability of occurrence.Hence, in the same year only the earthquake event having the highest magnitude amongst all has been considered for further study and all other events have been removed from the data base.
Therefore, for further study, they used total 148 events (Mw ≥ 6.0) that occurred from 1737 to 2011.But in the present study, after declustering the catalog, without deleting any real data (even if, there is a month's difference), all the mainshocks (Mw ≥ 6.0) have been considered for further study.Therefore, total 191 earthquake events (Mw ≥ 6.0) that occurred from 1737 to 2015 have been used in this study and superimposed on the digitized tectonic map (Seismotectonic Atlas of India [SEISAT], 2000; shown in Figure 3a).

Identified seismic sources in the study area
In order to carry out forecasting of future earthquake occurrence of magnitude size (Mw ≥ 6.0) in the region that may generate a most catastrophic scenario in the region, the fault data have been collected from SESAT-2000, remote sensing image interpretation through image processing techniques, and published available literature.However, considering the seismicity pattern, epicentre distribution of past events, seismicity parameters and orientation of faults, the following lineaments are identified as the most active faults in the region for further study.

Main Boundary Thrust (MBT) and Main Central Thrust (MCT)
The Eastern Himalayan range is characterized by number of north heading thrusts.Main boundary thrust (MBT), main central thrust (MCT) and main frontal thrust (MFT) are very important amongst them.The entire NE Himalayan belt shows transverse tectonics [Kayal 2001].This region is seismically very active and responsible for 1947 earthquake, Mw 7.8 and recent September 17, 2011 Sikkim earthquake, Mw 6.9 [Kumar et al. 2012].

Mishmi and Lohiti thrust
Mishmi thrust, Lohiti thrust and end part of the Disang thrust is present in the syntaxis zone which is the meeting place of the Himalayan and the Indo-Burmese arc.Molnar [1990], Seeber and Armbruster [1981] suggested that ongoing subduction of India plate under the Eurasian plate causes great and major earthquakes in this zone.This Assam syntaxis is responsible for 1950, great Assam earthquake.

Oldham fault
As per Bilham and England [2001] the 110 km long Oldham fault is a steep, ESE striking fault dipping SSW at 57° present in the northern part of the Shillong plateau.Which demarcates the northern boundary of the Shillong plateau and it is responsible for 1897 earthquake as suggested by Bilham and England [2001], Nayak et al. [2008] and Saha et al. [2007].But Rajendran et al. [2004] suggested that the great earthquake of 1897 had occurred on the Brahmaputra fault and it is the actual northern boundary of the Shillong plateau.Olympa and Kumar [2015) suggested that the Dauki fault is an E-W trending fault to the south of the Shillong plateau having length approximately 320 km.While Evans and Nandy suggested that Dauki fault is a near vertical or a south dipping strike slip/normal fault.1923 Meghalaya earthquake occurred on Dauki fault.But Morino et al. [2011] suggested that among the historic earthquakes the 1548 earthquake, the 1664 earthquake and 1897 earthquake occurred due to the activity of the Dauki fault.

Dapsi thrust
Olympa and Kumar, 2015 suggested that the Dapsi thrust is the extension of the Dauki fault along NW-SE direction.It is 90~100 km long reverse fault with strike slip component.Kayal and De [1991], Kayal [2001] Suggested that this thrust truncates the maximum isoseismal of the 1897 Shillong earthquake in west of the Shillong plateau.4.6 Chedrang fault, Samin fault and Dudhnoi fault Rajendran et al. [2004] suggested that the 20 km long NW-SE trending Chedrang fault and 4 km long E-W trending Samin fault were developed after the 1897 Shillong earthquake near the Shillong plateau.Baruah et al. [2010] studied the tectonics of the Chedrang valley and its vicinity area which is the western part of the Shillong plateau.They suggested the tectonics of this region is influenced by NW-SE oriented Chedrang fault, N-S oriented Dudhnoi fault, NW-SE oriented Samin fault and E-W oriented Dapsi thrust.

Dhubri fault, Tista fault, Padma fault and Madhupur Blind Fault (MBF)
Dhubri fault is a N-S trending fault which lies to the west of the Shillong plateau and separating the plateau from the Indian subcontinent.It is responsible for 1930 Dhubri earthquake [Olympa and Kumar, 2015].Rao and Kumar [1997] first proposed the pop-up tectonics for Shillong plateau and they suggested that the pop-up mechanism was facilitated by Disang thrust in the east, Dhubri fault in the west, Brahmaputra fault in the north and Dauki fault in the south of the Shillong plateau.The Tista fault is a NW-SE trending fault.Padma lineament and MBF are also traversing in the same direction.

Kopili fault
Kopili fault zone is also a main intra plate fault zone in the region as Shillong plateau which transgresses into the Himalayan upto the main central thrust (MCT).Kayal et al. [2006] and Bhattacharya et al. [2010] suggested that Kopili fault zone is approximately 300 km long and 50km wide zone which is like a divider between Shillong plateau and Mikir massif and it is a north dipping strike slip fault.Nandy [2001] suggested that 1869 Cachar earthquake occurred on south eastern end of Kopili fault zone.Kopili fault is also responsible for 1943 Assam earthquake.

Naga and Disang thrust
The eastward subducting India plate under Eurasian plate caused formation of Indo-Burma ranges which comprises the Naga thrust and the Arakan-Yoma fold belt.The extension of the Naga thrust to the southwest is named as Disang thrust which extends as haflong fracture zone and then it joins the east-west trending Dauki fault [Chaudhury and Srivastava, 1976].

A.2
The area between Assam syntaxis and Mikir hills is seismically less active and named as "Assam Gap" [Khattri and Wyss, 1978] and "Aseismic Corridor" [Kayal, 1996].They suggested that a large earthquake may be expected in this area.

Churachandpur-mao fault
Churachandpur-Mao fault (CMF) is an N-S trending strike slip fault which is responsible for low magnitude earthquakes.Kumar et al. [2011] found 0.5 mm/year slip rate in southern part and 3.9 mm/year slip rate in northern part of the CMF by geodetic observations during year 2004-2005 and suggested that the southern part of CMF is less active while the northern part is deformed at higher rates.The 4 January 2016 Manipur earthquake occurred 15 km west to the CMF at depth of ~59 km, described by Gahalaut and Kundu [2016].

Sylhet fault
Sylhet fault is delineated in Figure 1.a and 1.b.NE-SW trending Sylhet fault is the active fault in the Bengal basin.The 1918 Srimangal earthquake (Mw=7.6)occurred beneath the Bengal basin due to the rupture along the Sylhet fault [Sarker et al., 2010].

Chittagong Coastal Fault (CCF), Kaladan, MAT, Eastern Boundary Thrust (EBT), Kabaw fault and Myanmar Central Basin (MCB)
Chittagong coastal fault (CCF), Kaladan, Mat, Eastern boundary thrust (EBT), Kabaw fault, Myanmar central basin (MCB) and other fractures are delineated in Figure 1.a and1.b.These are the major faults which affect the Indo Burmese wedge by strike slip faulting.Maurin and Rangin [2009] proposed that CCF is a new major fault that onset after the Kaladan fault due to the westward progression of the Indo Burmese wedge.The CCF is seismically less active but had experienced two major earthquakes in 1851 and 1865.While Sikder and Alam [2003] described that the Kaladan fault bounds the outer Indo Burmese wedge to the east.The Kabaw fault is in between the Indo Burmese wedge and Myanmar central basin as described by Sikder and Alam [2003].

Sagaing fault
About 60 % of relative motion of India plate and Sunda plate is taken up by Sagaing fault suggested by Maurin et al. [2010] and Vigny et al. [2003].Kundu and Gahalaut [2012] suggested that earthquakes are generated by strike slip faulting on this Sagaing fault.It is responsible for the May 23, 1912earthquake [Guzman-Speziale and Ni 1996, Maurin et al. 2010, Richter 1958, Tsutsumi and Sato 2009].

Selection of seismic sources in the study area
After removing the repeated events and declustering the foreshock-aftershock activities, only few earthquakes (Mw ≥ 6.0) are observed in some seismic source zone.The seismic source zones, where less than four Hence, in this study, we have found total 12 seismic source zones out of 29 seismic source zones where reasonable amount of earthquake events (Mw ≥ 6.0) have been found from the year 1737 to 2015 in the study area (see Table 1) and used for further study.These 12 seismic source zones have been delineated and highlighted in the Figure 3 (b).

Tectonics and seismicity of study region
The India plate in the north is subducting under Eurasian plate and in the east, it is subducting and colliding under Burmese plate [Molnar and Tapponnier 1975, 1977, Kayal 1996, 1998, Hall 1997, Olympa and Kumar 2015] as shown in Figure 4. Hall [1997] suggested that due to the collision of the India plate with Eurasian plate, the Burma plate consisting of IBW (Indo Burmese wedge), MCB (Myanmar Central Basin) along with Andaman Sumatra arc rotated clockwise.Actually Burma plate is sandwiched between the India plate to the west and Sunda plate to the east (see Figure 4).Vigny et al. [2003], Nielsen et al. [2004], Maurin et al. [2010] suggested that in the region of Indo Burmese wedge, India plate is moving with a rate of 35 mm/year towards the north with respect to the Sunda plate.Due to these complex tectonic settings various major earthquakes have been reported in the past.Gahalaut and Kundu [2016] suggested that 2 April 1762 and 24 August 1858 Arakan earthquake caused notable damages which were reported from Chittagong, Sylhet, Manipur valley and Cachar region.Oldham [1882] and Singh [1965] reported that January 10, 1869 Cachar earthquake (Mw=7.5)caused severe damages in Cachar valley, near Silchar and Manipur valley, near Imphal.Five deaths from Silchar and three deaths from Imphal were reported.
However, in 1897 Shillong earthquake [Mw=8.1,Olympa and Kumar 2015] the northern edge of the Shillong Plateau rose more than 11 meters and acceleration exceeded the gravitational acceleration (1 g) as observed by Bilham and England [2001].However, May 23, 1912 earthquake [M=8, Maurin et al. 2010 Myanmar earthquake [M 6.9, USGS] and August 24, 2016 Myanmar earthquake [M 6.8, USGS] were the some major earthquakes that had been experienced by the study area.Molnar [1990], Seeber and Armbruster [1981] suggested that due to ongoing subduction of India plate under Eurasian plate caused the 1950 earthquake.Kundu and Gahalaut [2012] suggested that the 1988 Manipur earthquake was the largest magnitude earthquake in Indo Burmese arc region in the past 50 years.

Probability distribution models
In probabilistic and statistical analysis, a probability distribution model provides a mathematical description of random events in terms of probabilities.Oc-currence of an earthquake is a random phenomenon in nature.Statistical analysis is the most appropriate and widely used method for forecasting of the future earthquakes.The statistical approach based on trends of earthquake such as seismicity patterns, seismic gaps, characteristics of earthquake is to be adopted for forecasting of earthquake in the northeast region of India.
(a) Gamma distribution The Gamma distribution is a two parameter family of continuous probability distribution used to model sums of exponentially distributed random variables.Gamma model is often used as a probability model for waiting times.Chi-square and Exponential distribution are children of the Gamma distribution because these two are special cases of the Gamma distribution [Dikko et al. 2013].
If T is the time in years between the two successive earthquake events then the probability distribution function (pdf ) is given by Gamma as; (4) Where 'a' is shape parameter and 'b' is scale parameter.The inverse of the scale parameter is also called a rate parameter [Stacy 1962, Box-Steffensmeier andJones 2004].
The cumulative probability Ø(t) and conditional probability p (t/T) for Gamma distribution are given as; (5) Here T denotes already elapsed time without any strong event and Γ(p,q) represents the incomplete Gamma function of the second kind [Parvez and Ram, 1999] i.e.
(7) This is the same as described below in the Weibull model.The expected recurrence interval of earthquake event can be estimated as;   (Maurin and Rangin, 2009).
which have a normally distributed logarithm.So this distribution is useful for representing quantities that can't have negative values, since log(x) exists only when x is positive.Actually x is the multiplicative product of many independent random variables, each of which is positive [Heyde 1963, Barakat 1976].Probability density function (pdf ), cumulative probability of the next earthquake event Ø(t) , and conditional probability p (t/T) of the two parameters lognormal distribution [Nishenko andBuland 1987, Yadav et al. 2008] are given as; (10) Where, T is the time between two successive events.On a logarithmic scale m and σ are termed as location and scale parameters of lognormal distribution.cumulative probability (11) cumulative probability (12) Here T denotes already elapsed time without any strong event and Ø(y) represents the error integral given as; (13) The expectation of the earthquake in period t is calculated with the expected value Weibull distribution is a continuous probability distribution.The Weibull distribution with two parameters (shape, scale) is the most widely used stochastic model for earthquake data analysis [Utsu 1984, 2002, Sil et al. 2015] If T is the time in years between the two successive earthquake events then the probability distribution function (pdf ) is given by Weibull [1951] as: Where, β is the shape parameter and θ is the scale parameter of Weibull distribution.The Weibull cumulative distribution function (cdf ) is expressed as; (17) The conditional probability p(t/T) that the next earthquake will occur during the time interval between T and T+t can be expressed as ( 18) Where, T denotes already elapsed time without any strong event.In quality control engineering this conditional probability is called hazard rate [Wesnousky et al. 1984 andRikitake 1991].
The expected recurrence interval of earthquake event can be estimated as [Hagiwara 1974, Yilmaz andCelik 2008]: (d) Log-logistic distribution [Yilmaz et al. 2004, Yilmaz andCelik 2008] Loglogistic distribution is a probability distribution of random variables, if logarithm of the random variable is logistically distributed.The loglogistic distribution is generally used to model events that experience an initial increase, followed by a rate decrease.If T is the time in years between the two successive earthquake events then the probability distribution function (pdf ), cumulative probability or cumulative distribution function (cdf ) are given by log-logistic [Ashkar and Mahdi 2006] Where m and σ are the model parameters of Loglogistic distribution.On a logarithmic scale m and σ are termed as location and scale parameters of Loglogistic distribution.In conditional probability T denotes already elapsed time without any strong event.
The expected or mean interval of earthquake occurrence time can be estimated as; (26) For σ ≥1.0,E(t) doesn't exist.

Logarithm likelihood function and maximum likelihood estimation method
The logarithm likelihood function of the parameter 'a' given the observed data 't', is presented by lnL(a ⁄ t).The logarithm likelihood function lnL(a ⁄ t) is not a probability density function.It is an important component of both frequentist and Bayesian analyses.If we compare the logarithm likelihood function at two parameter points and find that lnL(a 1 ⁄ t) > lnL(a 2 ⁄ t) then the sample we actually observed is more likely to have occurred if a=a 1 than a= a 2 .This could be interpreted as 'a 1 ' is more suitable value for 'a' than 'a 2 '.A higher value of this function suggests a better (or more suitable) model and lower shows worse [Sil et al. 2015].The loglikelihood function has the particularly simple form; (28) Thus, lnL(a ⁄ t) represents the log-likelihood of the parameter 'a' given the observed data 't', and as such is a function of 'a'.And f(t i ;a) represents the probability distribution function of the random variables't' given the parameter 'a'.
However, in this work maximum likelihood estimation has been done using the data statistics.In this case of the maximum likelihood estimation, let we have a probability distribution function (pdf ) f(t i ;a) of random variables t and we are interested in estimating 'a' parameter.a mle is the value of 'a' that maximizes value of likelihood function L(a ⁄ t).It is very difficult to maximize L(a ⁄ t) directly but it is much easier to maximize the log-likelihood function lnL(a ⁄ t).Since ln(•) is a monotonic function.The value of the 'a' that maximizes lnL(a ⁄ t) will also maximizeL(a ⁄ t).Therefore we may also define a mle as the value of 'a' that solvesMAX a lnL(a ⁄ t).We may find the Maximum likelihood estimate by differentiating lnL(a ⁄ t) with respect to parameter 'a' and equating zero [Myung 2003]. (29)

Artificial Neural Network (ANN)
A neural network is a computing system inspired by biological neural network.In an Artificial Neural Network (ANN), the processing elements called neurons are connected together to form a network, where every link that connects neurons has a numeric weight.This numeric weight refers to the strength or amplitude of a connection between neurons.The output from each processing element is calculated by applying an activation function to the sum of inputs multiple by the weight vector [Alarifi et al. 2012].The Artificial Neural Networks support number of activation functions, such as linear function, non-linear function, sign function, step function [Alarifi et al. 2012].
There are many kinds of neural networks depending on their structure, function, and their training method.In this study, we used a feedforward neural network with a backpropagation learning algorithm to train the network.A typical neural network in feedforward direction is given by [Bichkule 2014]: Where a i is the input vector, O j is the output vector, w i j is a weight factor between two neurons, b j is the bias weight vector.Among the different kinds of activation function, the Log-sigmoid function and the rectifier linear function have been adopted in this study.These activation functions were applied to the neurons in hidden and output layers, whereas Input layer neurons are linear.The Log-Sigmoid function is considered the most popular activation function used in multilayer networks that are trained using backpropagation algorithm because this function is differentiable [Alarifi et al. 2012].This function [Wong et al. 1992] is defined as: (31) The rectifier activation function [Maas et al. 2013] is defined as: (32) A neuron employing the rectifier function is also called rectified linear neuron, which is neither fully differentiable (not at zero) nor bounded [Glorot et al. 2011].However, its derivative can take only two values zero or one and it can be represented as; The backpropagation learning algorithm is based on a generalized delta-rule, accelerated by a momentum term.The momentum term is added to accelerate the convergence of error during the training, with good learning rate.The weight factors and bias are adjusted by using the following equations [Bichkule 2014]: ∂ln L(a mle /t) ∂a Where 'η' is the learning rate, 'α' is the momentum coefficient, 'w ij ' is the weight factor associated between two neurons 'b j ' is the bias weight, 'O' is the output, 'δ' is the gradient descent correction term and 'k' stands for number of pattern.Each pattern presentation represents an iteration and a presentation of the entire training set is called an epoch.The performance of the trained network was checked by estimating the sum of squared error (SSE) using eq.( 39 N j = the total number of output variables.

K-Fold cross validation
This test is performed to validate the estimated values in the neural networks.In k-fold cross validation, the dataset 'D' is randomly split into approximately equal size of 'k' subsets ('D1', 'D2', 'D3',…,'Dk').The training and testing iterations are performed 'k' times and in each iteration i ε (1, 2, 3,….., k) and the dataset is tested on D \ Di and tested on Di [Kohavi 1995].In this way, in k-fold cross validation, every data point would be part of training sample and part of validation sample.
The estimated results are validated by calculating the average error (in this study, 'SSE').To calcu-late the average SSE, the sum of all the SSE of the testing subsets is divided by 'k'.

Estimation of logarithm likelihood function and model parameters
To check the goodness of fit of models with actual data, logarithm likelihood function is estimated for all considered four models.A higher value of this function suggests a better (or more suitable) model and lower shows worse.Logarithm likelihood functions corresponding to all twelve fault/thrust zones for all distribution models are listed in Table 2.
In order to estimate model parameters three different approaches can be applied.Namely, maximum likelihood estimation (MLE), method of moments (MOM) and least square method (LSM).But Utsu [1984] concluded that there is no significant difference among the results of MLE and MOM methods.While Yilmaz and Celik [2008] described MLE appeared to be the best estimation method among these three methods.Sil et al. [2015] also used the MLE method to estimate the model parameters after observing its higher precision characteristics.As a result, we only used maximum likelihood estimation (MLE) method to estimate the model parameters of considered distribution models in this study.Model parameters corresponding to each fault/thrust zones and whole study area for their best fit distribution models are listed in Table 3.

Probability of recurrence
The expected or mean interval of earthquake occurrence time for all twelve seismic source zones and ZAROLA AND SIL 14 whole study area have been estimated using maximum likelihood estimation and is listed in Table 3.The conditional probabilities for all twelve seismic source zones and whole study area have been estimated for various combinations of elapsed time (T) and time interval (t) using earthquake catalogue (Mw ≥ 6.0) from year 1737 to 2015 with their best fit model and respective model parameters which are listed in Table 4 and plotted in Figure 5 (a-m).
Since, this research work is being done in the year 2016, so the conditional probabilities that an earthquake with magnitude Mw ≥ 6.0 will occur in the next t years (time interval 1-24 i.e. 2017-2040) for elapsed time up to 2016 from the occurrence of the last earthquake of magnitude Mw ≥ 6.0 also have been estimated for all twelve seismic source zones using their respective best fit model and model parameters.For example the last earthquake occurred in MBT and MCT region in the year 2011, therefore the conditional probabilities have been estimated for elapsed time T=2016-2011=5 years with different time intervals t (1-24 years i.e. 2017 to 2040).See Table 5 and Figure 6.

Training and testing of neural network
In this study we used the Artificial Neural Network (ANN) technique to estimate the magnitude of future earthquakes in the study region.Time of the future earthquakes has been estimated using different stochastic models for different seismic source zones in this research work.The proposed Artificial Neural Network (ANN) technique is based on feedforward backpropagation neural network model with single hidden layer and neural networks were trained for maximum 25000 numbers of training epochs to attain the sum of squared error between target outputs and estimated outputs as 0.009.
It is a best practice to pre-process the input data before using it in training the neural network.Normalization of the input and output parameters is required so that all the input and output parameters are at a comparable range.Without this normalization, training of neural networks is very slow [Jayalakshmi and Santhakumaran 2011]

Best fi t model
The conditional probability that an earthquake with Mw ≥ 6.0 will occur in the next t years, given that an elapsed time of T years have passed from the last earthquake event
To forecast the magnitude of future earthquake for any seismic source zone, first the neural network is to be trained for the past seismicity of that source zone.The Architecture of a typical neural network with three input parameters (longitude, latitude and elapsed time between two consecutive earthquake events) and one output parameter (magnitude of earthquake) is shown in Figure 7. Since, the study area has been divided into 29 identified seismic source zones and the dataset available for some zones is too small therefore, it is not appropriate to divide the dataset into training and testing datasets in such zones.The seismic source zones, where less than 10 patterns (data points) are available, such kind of seismic source zones are not possible to analyze using 5-fold cross validation method.Therefore, out of 12 seismic Estimation of conditional probability that an earthquake with Mw ≥ 6.0 will occur in the next t years, given that an elapsed time of T years have passed from the last earthquake event.(e) (f ) testing subsets (see Table 6.c) to perform 5-fold cross validation.
Firstly, the neural networks have been trained for different values of learning rate (η) and momentum coefficient (α).It is observed that a higher value of η leads to faster learning but the weights starts oscillating, while a lower value of η leads to slower learning process and long training time.To reduce the training time, the best option is to add momentum term into the weight update procedure, with a smaller value of η.It is observed that for η= 0.7 and α= 0.9 network yields stable solution and it reaches to desired best output value very fast.It has been also observed that for most of the seismic source zones the training SSE gets minimized to 0.009 before completing the 25000 numbers of training epochs while for some other seismic source zones the desired level of training SSE is not achieved even after completing the 25000 numbers of training epochs.
After training of the neural networks for training subsets, the trained networks have been used to estimate output "O" for testing subsets and testing SSE has been reported (see Table 6.c).Neural networks have been also trained considering all the data points for the whole study area.The dataset of the whole study area is divided into five subsets to perform the 5-fold cross validation and average SSE is estimated (Figure 8, a-e).

Magnitude estimation of the future earthquake
In this study, we have to answer the three questions: (1) Location of future earthquake (2) Time of future earthquake and (3) How big would it be.Prediction of exact location (longitude, latitude) of the future earthquake precisely is a challengefull task as our knowledge become limited still about the generation of earthquake process at the source and heterogeneity of the rock layers.The time of future earthquake has been estimated using conditional probabilities for various combination of elapsed time 'T' and time interval't' using four stochastic models (Gamma, Lognormal, Weibull and Log-logistic) in this research work.
To estimate the magnitude of future earthquake for any particular seismic source zone, first the neural network has been trained for the past seismicity of that source zone.To use this trained neural network for es-timation of magnitude of next earthquake, input parameters (longitude, latitude, and time) of next earthquake are required.The longitude and latitude of future earthquake have been taken within the seismic source zone where already maximum magnitude of earthquake has been recorded in the past years.For example the maximum observed magnitude of earthquake for MBT and MCT region is 7.7 Mw (see Table 1) therefore longitude 93.85°E and latitude 28.65°N have been taken as the location of future earthquake.The magnitude of the earthquake has been estimated for elapsed time T=0 year and time interval t in which seismic source zone shows conditional probability greater than 0.8.For example, for elapsed time zero year MBT and MCT shows conditional probability 0.78-0.93after about 5-8 years (see Table 4).For this time interval of 5 to 8 years, two values of earthquake Table 5. Estimation of conditional probability (using best fit model) that an earthquake with magnitude Mw ≥ 6.0 will occur in the next t (time interval) years i.e. 2017-2040 for elapsed time T equals to 2016 minus year of last earthquake event in the seismic source zone.magnitude were obtained; one is corresponding to 5 years and second is corresponding to 8 years.The average value of these two earthquake magnitudes were calculated and listed as magnitude of future earthquake.The magnitude of future earthquakes has been estimated using both the activation functions (Log-Sigmoid and Rectifier).In Table 7 it is listed for all six seismic source zones for the elapsed time T=0 year and time interval (t) in which conditional probability reaches above 0.8 and for the location within the seismic source zone where already maximum magnitude of earthquake has been recorded in the past years.

Results and discussion
The expected or mean interval of earthquake occurrence time for all twelve seismic source zones has been estimated using maximum likelihood estimation and is listed in Table 3.It is found that EBT and Kabaw zones has the smallest recurrence period (2.77 years) indicating highly potential and active seismic source zones among other sources, whereas Sylhet fault zone has highest recurrence period (38.01 years) for earthquake magnitude Mw ≥ 6.0.The mean recurrence interval for the whole northeast India is found to be 1.35 years.
It is observed that for elapsed time T=0 year, EBT and Kabaw zone shows highest conditional probability among all twelve seismic source zones showing chances of occurrences 0.7 to 0.91 (see Table 4) after about small time interval of 3-6 years (2014-2017; since last earthquake of Mw ≥ 6.0 occurred in the year 2011, see Table 1) for an earthquake magnitude Mw ≥ 6.0.Whereas, Sylhet zone shows lowest value of conditional probability among all twelve seismic source zones and it reaches 0.7 (see Table 4) after about a large time interval of 48 years (year 2045, since last event (Mw ≥ 6.0) occurred in the year 1999, see Table 1).For the same elapsed time (T=0), the conditional probability reaches above 0.8 for MCB, MBT and MCT and Sagaing seismic source zones after about small time interval of 3-6 years; for Kaladan, after about 6-9 years; for CCF, after about 12-15 years; for Mishmi and CMF, after about 15-18 years; for A.2, after about 18-21 years; for Dauki and Kopili, after about 33-36 years (see Table 4).And for elapsed time T=0 year, the whole study area shows chances of occurrences 0.85 to 0.92 after about very small time interval of 2-3 years (2014-2015; since last earthquake of Mw ≥ 6.0 occurred in the study area in 2012, see Table 1).It is also noticed that Lognormal and Log-Logistic models show higher conditional probabilities for smaller values of elapsed time (T) and lower conditional probabilities for larger values of elapsed time.While Weibull model shows just opposite results, means it shows lower conditional probabilities for smaller value of elapsed time and higher conditional probabilities for larger value of elapsed time and probability graph follows the same pattern (shown in Figure 5, a-m).
Since, this research work has been done in the year 2016, so the conditional probabilities that an earthquake with magnitude Mw ≥ 6.0 will occur in the next t years (time interval 1-24 i.e. 2017-2040) for elapsed time up to 2016 from the occurrence of the last earthquake of magnitude Mw ≥ 6.0 also have been estimated.See Table 5 and Figure 6, which shows that after

Table 6c.
Testing of an artificial neural network using 5-fold cross validation.
higher value of ln L function and for three seismic source zones (Sylhet, Kopili, CCF) Weibull model shows the higher value of ln L function.
Thereafter, the possible magnitude of future earthquakes has been forecasted using Artificial Neural Network (ANN) for six seismic source zones (where more than 10 data points are available) and the whole study area (see Table 7).When the probable magnitude of the future earthquake is estimated considering all the data for the whole study area, the average SSE in the 5-fold cross validation is found as 0.0364 using logsigmoid and 0.0346 using rectifier activation functions.(See Figure 8, a-e).
The results obtained in this study have been validated in terms of location, time and magnitude of earthquake events that occurred in the year 2016 in the study area and it is presented in Table 8.Gahalaut and Kundu [2016] concluded that the 4 January 2016 Manipur earthquake (Mw 6.7) occurred 20 km west of the Churachandpur-Mao fault (CMF).The USGS reported that the 13 April 2016 Myanmar earthquake (Mw 6.9) occurred at latitude 23.094° N and longitude 94.In this study, it is observed that for conditional probability 0.77-0.82 the neural network gives the magnitude of future earthquake as Mw 6.6 in CMF region for elapsed time T=0 year with time interval 15 to 18 years (2014-2017) and for conditional probability 0.76-0.92 it gives Mw 6.8 in MCB region for elapsed time T=0 year and time interval 4 to 7 years (2013-2016) and for conditional probability 0.8-0.94 it gives Mw 6.5 in EBT and Kabaw region for elapsed time T=0 year and time interval 4 to 7 years (2015)(2016)(2017)(2018).These observations confirm the reliability of the results and we believe that these results would be very useful for assessment of seismic hazard and mitigation of earthquake risk in the study area.

Conclusion
The present study describes the probabilistic assessment and forecasting of magnitude of future earthquakes in the northeast India using four models (Gamma, Lognormal, Weibull and Log-logistic) and artificial neural networks, respectively.There are many identified geologically active faults in the northeast region of India.In this study it has been noticed that out of 29 identified seismic source zones only 12 zones are seismically active and capable in producing earthquakes of Mw ≥ 6.0.Out of these 12 zones, 6 zones are seismically very active and produce earthquakes (Mw ≥ 6.0) frequently.Therefore, the probable magnitude of future earthquakes has been estimated for these six zones only and other six zones are kept to collect more earthquake events for future work.The artificial neural network model presented in this study yields higher forecasting accuracy for magnitude of earthquakes in the northeast region of India because of its capability of capturing non linear relationship.It is observed that, for most of the seismic source zones the rectifier function increases the training speed of the neural network and the sum of squared error gets minimized to 0.009 very fast.The results show that the seismic source zones CMF, MCB, EBT and Kabaw, Kaladan, MBT and MCT and Sagaing fault zone are highly potential seismic source zones having very small recurrence interval and capable in producing the significant earthquakes in northeast region of India.Therefore, earthquake forecasting and the proper seismic microzonation projects are required for such regions which would help in risk assessment, disaster mitigation and also in designing of safe structures to reduce loss of human lives.

FORECASTING
OF EARTHQUAKE IN NER, INDIA

Figure 2 .
Figure 2. .The figure shows a histogram of earthquakes (Mw ≥ 4.0) that occurred from year 1737-2015 in the study area.

Figure 3 .
Figure 3. (a).Map shows epicentral distribution of earthquakes (Mw ≥ 6.0) used for estimation of earthquake recurrence in the study area collected from the year 1737-2015.(b).Map shows the selected 12 seismic source zones for forecasting of future earthquakes.These selected zones were highlighted in the map.
Lognormal distributionThe Lognormal distribution is a statistical, continuous probabilistic distribution of random variables FORECASTING OF EARTHQUAKE IN NER, INDIA

Figure 4 .
Figure 4. .Map shows the motion of the India plate with respect to Burma, Sunda and Eurasian plates, respectively (Maurin and Rangin, 2009).
= the target value of output variables.O = the estimated value of output variables.N k = the total number of training data sets.

Figure 5
Figure5 (a-d).Conditional probability of earthquake magnitude Mw ≥ 6.0 estimated considering different values of elapsed time (T) and time interval (t) for seismic source zones Sylhet, Main Boundary Thrust (MBT) and Main Central Thrust (MCT), Sagaing and Eastern Boundary Thrust (EBT) and Kabaw.

Figure 5
Figure 5 (i-l).Conditional probability of earthquake magnitude Mw ≥ 6.0 estimated considering different values of elapsed time (T) and time interval (t) for seismic source zones Kaladan, Chittagong Coastal Fault (CCF), Churachandpur Mao Fault (CMF) and A2.

Figure 5
Figure 5 (m).Conditional probability of earthquake magnitude Mw ≥ 6.0 estimated considering different values of elapsed time (T) and time interval (t) for the whole study area.

Figure 6 .
Figure 6.Map shows conditional probability of earthquake magnitude Mw ≥ 6.0 estimated for the years 2017 to 2040 considering total elapsed time upto 2016 for all twelve seismic source zones.

Figure 8 .
Figure 8. Figures show the 5-fold cross validation of the testing datasets considering earthquakes of Mw ≥ 6.0 that occurred in the study area from 1737-2015.

Table 3 .
. Data normalization makes the training of the network faster and efficient.It yields FORECASTING OF EARTHQUAKE IN NER, INDIA Estimation of model parameters using maximum likelihood estimation method.

Table 6a .
Normalized dataset for different seismic source zones.

Table 6b .
Training of an artificial neural network using 5-fold cross validation.
.91 for MBT and MCT region, 0.85 to 0.88 for EBT and Kabaw zone, 0.81 to 0.85 for MCB region, 0.79 to 0.84 for CCF zone and 0.75 to 0.79 for Sagaing fault zone.Whereas, for the same time interval conditional probability is less than 0.45 for Kaladan, Mishmi, CMF, A3, Dauki, Kopili and Sylhet zones.It is also observed that, for different probability dis-tribution models ln L values are very close to each other.In general, all the four models can be used to estimate the conditional probability of earthquake occurrences in northeast India.Indeed, for five seismic source zones (MBT and MCT, EBT and Kabaw, Mishmi, MCB, A2) Lognormal model shows higher value of ln L, while for four seismic source zones (Sagaing, Dauki, Kaladan, CMF) and whole study area Log-Logistic model shows FORECASTING OF EARTHQUAKE IN NER, INDIA

Magnitude estimation of future earthquakes Sr. No.
FORECASTING OF EARTHQUAKE IN NER, INDIA

Table 7 .
Estimation of magnitude of future earthquakes.

Table 8 .
Validation of results.