Forecasting seismicity rate in the north-west Himalaya using rate and state dependent friction law

ABSTRACT In this study, rate and state Coulomb stress transfer model is adopted to forecast the seismicity rate of earthquakes (MW ≥ 5) in the north-west Himalaya region within the testing period 2011–2013. Coulomb stress changes (ΔCFF), considered to be the most critical parameter in the model, exhibit stress increase in the whole study region, excluding the Chaman fault of the Kirthar range where significant stress shadow has been observed. The estimated background seismicity rate varies in the range 0.0–0.7 in the region, which is preoccupied by low aftershock duration of <50 years. Furthermore, a low b-value that varies between 0.54 and 0.83 is observed in Kirthar ranges, Karakoram fault and Pamir-Hindukush region. However, areas like Hazara syntaxis of the northern Pakistan and northern Pamir of the Eurasian plate exhibit higher b-values in the range 1.23–1.74. Considering constant constitutive properties of the faults (i.e. Aσ = 0.05 MPa), our forecast model for variable ΔCFF and heterogeneous b-value successfully captures the observed seismicity rate of earthquakes. Results have been verified using statistical S-test. However, the model fails to capture the observed seismicity rate during the period when reconstructed for average b-value to be 0.86 and no change in ΔCFF (ΔCFF = 0).


Introduction
In last few decades, studies and interests in earthquake forecasting have increased tremendously. The rate and state friction law, as introduced by Dieterich (1979Dieterich ( , 1994 on the basis of laboratory friction experiments, has been frequently used to model seismic activity. The model can be categorized into aseismic or seismic (earthquakes) based on chosen parameters (physical constitutive properties of faults, the stressing rate, and the background seismicity rate) incorporated into this rate and state friction law (Helmstetter and Shaw 2006). In unstable condition, the rate and state law provides a relation between stress history and seismicity (Dieterich 1994).Therefore, this model can be used to forecast seismicity rate changes caused by Coulomb stress changes (DCFF) after the mainshocks. Moreover, several recent studies have discussed this model and applied it in different seismically active regions for understanding the temporal and spatial response of seismicity due to successive stress changes, including, toggling, decay, and aftershock migration for seismic hazard assessment (King et al. 1994;Nalbant et al. 1998;Stein 1999Stein , 2003. CONTACT  Several worldwide evidences have been reported that the small changes in DCFF in rate and state Coulomb stress transfer model can enhance the seismicity rate in areas of high background seismicity and can also suppress the same in areas of stress shadows with previously high seismicity rates . Parsons (2002) has also successfully demonstrated that the highest probability rate for triggered seismicity is associated with the DCFF after the mainshocks, corroborating the model proposed by Dieterich (1994). Gomberg et al. (2000Gomberg et al. ( , 2005 also elucidated the model's ability in explaining the temporal behaviour of seismicity and in justifying the variation of time delays between successive events. Several investigators have considered the DCFF and background seismicity rate (r) for investigating seismicity, earthquake migration and seismic hazard assessment Nalbant et al. 2005). Toda et al. (2005) calculated the expected seismicity rate of next ten years in southern California. Further, Toda and Enescu (2011) modified the methodology which was previously adopted by several authors (Hainzl et al. 2009;Toda et al. 2005) for forecasting M W 5 shallow crustal earthquakes for Japan's mainland.
The seismically active north-west (NW) Himalayan region is formed due to oblique compression between the Indian and Eurasian plates, thereby resulting in large scale thrusting and rotation of faulted blocks (Valdiya 1991). The folds and thrust belts that exist in the study region are also considered to be seismically most active zones in the world (Khwaja and Monalisa 2005) (Figure 1). The seismological investigation of the area reveals that the occurrence of earthquakes (M W 4) is closely associated with both the surface and blind faults (Khwaja and Monalisa 2005). Seismic hazard evaluation in this seismically active Himalaya is therefore important because earthquakes cause several types of threats to the people living adjacent to this gigantic mountain chain. In addition, it is widely believed that numerous active faults in the geological past have potential and may be Dashed lines indicate the political boundary. Red lines indicate five broad seismogenic source zones on the basis of seismicity, tectonics and focal mechanism of earthquakes as suggested by Yadav et al. (2011). All zones are indicating by number i.e. 1, 2, 3, 4 and 5. (To view this figure in colour, see the online version of the journal). reactivated in future (Philip et al. 2014). They also elucidated that these active faults may trigger seismic activity (>80% seismic activity) in the study region.
In present study, we forecast the seismicity rate during the period 2011-2013 in the NW Himalaya and its adjoining region using modified rate and state model as proposed by Toda and Enescu (2011). We also determine the probability of M W 5.0 shallow crustal mainshock occurrences using homogenized M W earthquake catalogue from International Seismological Centre (ISC) and Global Centroid moment tensor (GCMT) accessible at http://www.isc.ac.uk/bulland http://www.globalcmt. org/CMTsearch.html, respectively, in the five broad seismogenic zones of NW Himalaya and its surrounding regions. These broad source zones have been delineated by Yadav et al. (2012) by considering the information of historical and instrumental seismicity, tectonics, geology, paleoseismology and other neo-tectonic properties of the studied region (Tandon and Srivastava 1974;Chandra 1978;Bapat et al. 1983;Oldham 1883). The present forecasting study is commonly known as pseudo-prospective testing and the forecasting model largely considers the constrained input parameters, based on the physical mechanism of earthquake occurrences, for estimating seismicity rates for the region in the testing period (Cocco et al. 2010;Strader and Jackson 2015).

Earthquake catalogue and fault plane solution
It is well accepted that going back in time leads to lower quality and inadequacy of data in comparison with the more recently obtained ones (Leptokaropoulos et al. 2012). Hence, we have analysed seismicity data between 1975 and 2010 from the earthquake catalogues of ISC and GCMT for forecasting the seismicity rate using rate and state friction law. On the basis of recording wave types of earthquakes, the catalogues comprise of different magnitude scales, i.e., local magnitude (M L ), body wave magnitude (m b ), surface wave magnitude (M S ) and moment magnitude M W . However, the magnitude scales such as M L , m b and M S saturate at different higher magnitudes and therefore need to be converted into a single non-saturated magnitude scale for avoiding incorrect estimation of any statistical analyses (Chingtham et al. 2014). For this purpose, a homogenized earthquake catalogue in M W is compiled from ISC and GCMT databases, for the period 1975-2010, by utilizing the regression technique as adopted by Yadav et al. (2011) for the study region bounded by latitude, 25 N -40 N and longitude, 65 E -85 E (Figure 2(a)). The spatial and temporal windowing method of Knopoff (2000) is employed for declustering the catalogue (Figure 2(b)). We have also checked the completeness of magnitude (M C ) for each sub-region by fitting the power law to the frequency magnitude distribution of earthquakes at 90% confidence level (Figure 3(a)-(e)) Here, the Maximum Curvature method is adopted for estimating M C in each source zone. The calculated M C for source zones 1, 2, 3, 4 and 5 are found to be 4.1, 4.3, 4.0, 4.0 and 4.3, respectively. During the period 1975-2010, the annual rates of declustered seismicity are then calculated in the square cells of 0.5 £ 0.5 by considering only the earthquake magnitudes above the minimum magnitude in each source region. Then Gaussian function of Frankel (1995) with correlation distance of 100 km is applied to smooth the annual rates of declustered seismicity. For our forecasting purpose, 20% of the minimum rate is assigned to the particular cells having zero seismicity rates. Furthermore, we obtained the focal mechanism solution from GCMT catalogue for computing the DCFF. The key problem for computing DCFF computation is precisely defined by receiver's fault mechanisms in each source zone. In present study, we computed coulomb stress changes based on specified fault for each zone except Zone 3. For Zone 3, we assumed the multiple receiver fault orientations as suggested by Hainzl et al. (2010) and Steacy et al. (2005). Details of the earthquake sources used in the analyses for the study region are mentioned in Table 1.

Seismic sources
We consider mainshocks of magnitude M W 5 that occurred after 1975 as the seismic sources for computing DCFF (Figure 4). Some moderate to large earthquakes that occurred in this study region are listed in Table 1. A variable slip model has been used for 8 October 2005 Kashmir earthquake, while a simple uniform slip model obtained from the established empirical relations of Wells and Coppersmith (1994) is considered for other earthquakes for estimating DCFF.

Receiver faults
The basic information on the fault geometry and slip direction of the receiver faults is essential and necessary for computing DCFF. Two different approaches are primarily adopted while taking into account the types of receiver faults for estimating DCFF (King et al. 1994;Toda and Enescu 2011). The first approach commonly known as specified fault considers that the regional dominant faulting mechanisms of the receiver source faults have similar focal mechanisms with those of the mainshock source faults due to the complex fault systems (Nostro et al. 2005). The other approach comprises of resolving the optimally oriented planes for Coulomb failure (King et al. 1994;King and Cocco 2001), and thereby considering the maximized total stress tensor obtained from the optimized values of strike, dip and rake of the source faults (Toda et al. 2005). As such, the predicted focal mechanisms associated with optimally oriented planes strongly depend on the orientation and magnitude of the regional stress fields.
In this study, the specified fault approach is adopted because the regional stress tensor is not well described for the study region. The Main Himalayan Thrust (MHT), Karakoram fault and Chaman fault (CF) are some of the major regional active faults and geological structures that control the seismotectonics of the study region. Central Himalayan region is well dominated by thrust faults with low dip angles (Shanker et al. 2011) whereas Pamir-Hindukush region exhibits reverse as well as strike-slip faults with active structures (Fan et al. 1994;Shanker et al. 2011). Mixtures of both the normal and right lateral strike-slip faulting earthquakes and also active geological structures are quite prevalent in the Karakoram Himalaya ( Figure 5(a)). Both the focal mechanisms and active fault information are utilized to select dominant regional mechanisms for each zone and subsequently, the grid map of these typical focal mechanisms at every 0.5 o interval is plotted by   that reverse faults generated during NE-SW shortening were reactivated as normal faults during NE-SW extension. Such NE-SW extension in the upper crust seems to be related to rampÀlike geometry of the underthrusting zone along the MHT and the motion of Himalayan rocks over such ramp as suggested by Hintersberger et al., 2010. Furthermore, Negi et al. (2017 has suggested that the weak evidence of normal fault mechanism at shallow depth is non-seismogenic feature with no persistent record at depth. Therefore in present forecast modelling we did not consider the fault mechanism other than regional fault mechanism for Zone 4.

Coulomb stress changes and its associated parameters
In an elastic half-space, the static Coulomb stress changes (DCFF) caused by each mainshock (Okada 1992) is calculated by considering 0.25 and 32 GPa as Poisson's ratio and shear modulus respectively. The simplifying assumptions of pore pressure effects (King et al. 1994) are primarily taken into account to compute DCFF by using the following equation: where Dt is the shear stress in the direction of slip on the assumed causative fault plane. Ds n is the normal stress changes (positive for unclamping or extension) and m 0 is the effective coefficient of friction for all the fault.
We have used uniform slip model for computing DCFF for all the events except 8 October 2005 Kashmir earthquake, where finite slip model is given by Parsons et al. (2006). Then the maximum DCFF has been computed from the seismogenic depth range for each zone i.e., for Zone 1, 0-15 km, Zone 2, 0-15 km, Zone 3, 0-35 km, Zone 4, 0-15 km, and Zone 5, 0-20 km (depth range for each zone is detailed in Supplementary Figure S1). As discussed earlier, DCFF for each zone is computed on specified fault except for Zone 3 where the characteristics of the mainshock occurrences are totally different due to complex fault systems. Figure 6 depicts that the study area is mostly dominated with high DCFF, except at few patches in Zone 1 and Zone 5 with low DCFF.
As we assumed the uniform slip model to compute the DCFF, sensitivity analysis assuming the uniform slip model with 10%, 20% and 30% variable in slip has been done to test the sensitivity of  Figure S2). We have also compared our results with available slip model for Kashmir earthquake (Parsons et al. 2006) independently; variations in DCFF are shown in Figure S2. In case of the uniform slip model, despite of changing the slip amount, the overall stress patterns are found to be consistent. Parsons et al. (2006) have done similar test for 2005 Kashmir earthquake with changing friction coefficient and target fault depth, and found the overall stress change pattern to be consistent.

Background seismicity rate
A fair understanding of r is a prerequisite for studying the effects of DCFF on the seismicity of a region. Although, the definition and the measure of r is still controversial (Hainzl and Ogata 2005;Lombardi et al. 2006;Lombardi and Marzocchi 2007), there are two different approaches generally reported in several literatures for retrospective seismicity forecasting (Toda and Enescu 2011). Toda et al. (2005) defined the first approach as the seismicity rate r before the occurrence of large earthquakes while the other approach considers r as a time independent smoothed seismicity rate estimated from the declustered catalogue in a prescribed time window (Catalli et al. 2008). In this paper, the approach provided by Catalli et al (2008) is used to estimate r from the declustered catalogue as shown in Figure 7.

Time-dependent rate and state formulation
Seismicity of a region for a particular time period can be examined through the history of stress perturbations and its corresponding state variable and seismicity rate changes associated with multiple mainshocks. The spatial heterogeneities of DCFF at each node of grids is estimated for every moderate to large mainshocks. The seismicity rate R of a region is related to the tectonic secular shear stressing rate, _ t , and the state variable, g, (Dieterich 1994) by the following equation: where r is the background seismicity rate of the considered area and g is the state variable that evolves with time and stressing history. For a constant_ t, g at each node takes the value given by the following equation: which accordingly gives R = r. This implies that the seismicity rate is the background seismicity rate without the presence of stress perturbation. Therefore, the initial g o changes to a new g due to increased/decreased DCFF after the mainshocks, following the equation: where g nÀ1 and g n are the values of the state variable just before and after the applied DCFF, respectively. The product, As is an important controlling parameter for fault friction and here its value is taken to be 0.05 MPa (Toda and Stein 2003). Hence, seismicity can be easily obtained by substituting this new state variable in Equation (4) during the time of stress perturbations. An increased DCFF on a fault finally increases the seismicity rate by producing large slip due to decreased g while a sudden stress drop increases g and lowers the seismicity rate (Toda and Enescu 2011). Then, g can be approximated for any next time step Dt by using the relation: In rate-and -state dependent friction model, aftershock duration, t a , can be expressed with As and _ t (Dieterich 1994) through the given equation: where t a is the duration of transient effects or the aftershock duration that defines the characteristic time for aftershocks to return to the background seismicity rate. Toda and Enescu (2011) suggested that the stress release takes several time periods (e.g. decades to centuries), mostly on the inactive faults. Along the plate interface, the stressing rate is found maximum but diminishes from the interface as a function of time. As such, the duration of such observed t a of the collision/subduction event is found to be much smaller in comparison to the mainshocks located in the stable continent. The stressing rate for a constant As taken throughout time and space in the Dieterich model is found to be inversely proportion to t a (equation 4). By considering the plate interface's distance and spatially clustered aftershock distributions of earlier large earthquakes, the spatially variable t a is arbitrarily assigned to each node of the grids (Figure 8). Thus, it is practically observed that the spatial distribution of t a strongly influences the initial condition of g. We start our calculations by assuming spatially non-homogeneous steady state variable ð g o Þ as the collision region generally exhibits lower g, whereas high g is mostly observed in the stable region. Furthermore, DCFF at discrete time step is computed that involves g changing with respect to time and space (equations (4) and (5)).

b-value calculation and extrapolation of large earthquakes
A maximum likelihood procedure (Aki 1965) is adopted to compute b-value at each node from at least 50 events encapsulated by 0.5 £ 0.5 grid. However, the b-value was not calculated due to incompleteness issue for those nodes where the total number of earthquakes having magnitudes greater than or equal to M C was less than 25. Hence, the estimation process tends to introduce spatial gaps due to constraints on the estimation of parameters by the number of earthquakes (Chingtham et al. 2014;Thingbaijam et al. 2008Thingbaijam et al. , 2009). In such cases, an average b-value is arbitrarily assigned to those particular nodes for forecasting purpose. Finally, the constant b-value estimation at each node is utilized for extrapolating the occurrence rate of small earthquakes to large earthquakes as suggested by Toda and Enescu (2011) (Figure 9).

Spatial likelihood test (S-test)
Statistical testing of earthquake forecast is an important task for introspecting the consistency of the forecast model and for this several procedures/methodologies exist. Kagan and Jackson (1995)   introduced three different tests namely the number test (N-test), the likelihood test (L-test), and the ratio test (R-test) to quantify the earthquake predictability for improving the seismic hazard estimation. In our present analysis, we tested our forecast model by using the spatial likelihood test, i.e. Stest as shown in Figure 12.
In this test, the spatial information is firstly extracted by summing the forecasted expected rates over the magnitude bins and later normalized its sum so that the sum gives the total number of observed earthquakes. Following the common notations used by Zechar et al. (2010), we can define the expected rates and the observed earthquakes as given below: where N obs and N fore indicates the total observed earthquakes and forecasted expected rates in the total magnitude-space bins, i.e. M and S defined in the study region. Here, w i; j ð Þ and λ i; j ð Þ also gives the number of earthquakes observed and forecasted expected rates in the particular magnitude-space bin (i, j). The observed joint log-likelihood which is estimated by considering these values is given by In order to account for forecast uncertainty in the model, a set of simulated catalogues ðVÞ, is obtained by simulating 1000 times, where each catalogue can be represented aŝ wherev x ði; jÞ gives the number of simulated earthquakes in bin (i, j). Then, the setŜ È É is constructed with the xth member equal to the joint log-likelihood of the xth simulated catalogue from each simulated catalogue and later compared with the observed joint log-likelihood. If this observed joint log-likelihood falls in the higher tail of simulated joint log-likelihood, then we can conclude that the observed and forecasted spatial seismicity patterns are consistent. However, this will not be true if we consider statistical S-test as a two-tailed test. In this study, we have considered S-test as one-tailed test by taking into consideration of those previous Collaboratory for the Study of Earthquake Prediction related studies (Schorlemmer et al. 2007;Zechar et al. 2010).
Finally, the quantile score z of S-test is estimated by using the following relation: WhereŜ x is the observed joint log-likelihood. If the z value is very small then, we can say that the observed spatial distribution is found inconsistent with the forecast. However, more and more details about this test can be obtained from the literature of Zechar et al. (2010).

Results and discussion
In order to forecast the seismicity rate, the input model parameters have to be constrained a priori value based on the available data and information of the target study area. These input parameters are then incorporated into the modified rate and state model for forecasting the seismicity rate of earthquakes (M W 5) during the period 2011-2013 ( Figure 10).

Coulomb (DCFF) stress analysis
First, we compute DCFF for mainshocks (M W 5) at seismogenic depth based on specified fault approach (King et al. 1994). There are 41 M W 5 earthquakes spread over the entire zone that are listed in Table 1. Stress increase and/or decrease have been observed for significant earthquakes in the study region. The noticeable increase in DCFF are primarily associated with 17 June 1990 Pakistan earthquake (Zone 1), 19 October 1991 Uttarkashi earthquake, 28 March 1999 Chamoli earthquake and 8 October 2005 Kashmir earthquake (Zone 4). The 17 June 1990 strike-slip Pakistan earthquake may also have important contribution to the stress shadow observed in the Chaman fault (Kirthar Range). In Zone 5, most of the earthquakes that are of strike-slip and normal faulting type are responsible for deformation in the Tibetan plateau. Ogata (2005) showed that even small changes in DCFF i.e. increment/decrement can trigger respective activation/lowering of microseismicity. Toda et al. (2005) and Cocco et al. (2010) have also suggested that DCFF may not alone induce extraordinary seismic activity, but the combination of all other parameters, i.e. As, the stressing rate and the background seismicity rate incorporating to the model, may result in significant triggering or rate decrease. Furthermore, Bhloscaidh et al. (2014) affirmed that the forecast ability of Coulomb-based models can also be considerably improved by employing spatial variations of r:

Seismicity rate changes
Considering the importance of r; we have calculated the spatial variations of annual background seismicity rate using declustered catalogue for the period 1975-2010. The value of r varies from 0-0.7 in the study region (Figure 7). The maximum seismicity rates are found in Zone 3 where the complex faults pattern and high seismicity is exhibited. The moderate background seismicity rate is found in Zones 1 and 5, where mostly strike-slip fault systems exist. In these zones, few large earthquakes have been documented in recent past (Table 1). Zone 2 is characterized by the diffused seismicity and moderate hazard. The high rate of seismicity in study region is primarily associated with the complex Himalayan thrust fault system formed during the continuous collision between the Indian and Eurasian plates (Avouac et al. 2001;Yadav et al. 2011)..

Analysis of aftershocks duration
From equation (6), t a is found to be inversely proportional to _ t when the parameter As is kept constant (0.05 MPa) throughout the region. In our study region, t a lies within the range 0-650 years. Figure 8 shows that the study region is predominantly occupied by low t a values except some high values seen in Zone 5, which might be caused by longer stress relaxation times of the earthquakes occurred in this zone. Dieterich (1994) pointed out that t a values range from 0.2 to 12 years for different regions worldwide; however, Toda et al. (2005) later calculated the values of t a ranging between 5 and 500 years for constant As (0.05 MPa). Catalli et al. (2008) estimated t a values in between 5714.357 and 57143 years by considering constant As (0.04 MPa) in the Umbria-Marche sector of Northern Apennines (Italy). Recently, Toda and Enescu (2011) have computed the t a values within the range 0-100 years by assuming constant of As (0.05 MPa) for Japan. Figure 9 shows the spatial distribution of b-value in the study regime of Himalaya. The low b-value observed in the Kirthar ranges, nearby the eastern part of CF (Zone 1), can be correlated with the occurrences of large earthquakes, indicating high differential stress build-up in the zone. Higher bvalues are found in the Hazara syntaxis of Northern Pakistan (Zone 2) suggesting an active ongoing creeping tectonics (Thingbaijam et al. 2009). However, low b-values are observed in Zones 3 and 4 that can be primarily associated with the presence of high differential stress accumulation within the rock mass. Moreover, high-stress pockets associated with low b-value arising out of the clustering events is observed along the southern part of Karakorum fault in the Tibetan Plateau region (Zone 5), which corroborates with findings of other researchers (Chingtham et al. 2015;Loannis et al. 2011;Shanker and Sharma 1998;Thingbaijam et al. 2009;Yadav et al. 2012).

Forecast model
The forecasted and the observed number of mainshocks having magnitude M W 5 during the years 2011-2013 are shown in Figure 10. The reference forecast model for the expected number of M W 5 mainshocks is prepared from the assumed DCFF = 0 and average b-value (0.86) during the period 2011-2013 (Figure 10(a)). Figure 10(b) shows the expected number of M W 5 mainshocks estimated from the assumed DCFF = 0 and spatially variable b-value during the period 2011-2013. However, Figure 10(c) depicts the forecasted frequency of M W 5 mainshocks, reproduced by considering average b-value (0.86) in our model, while stress perturbations and heterogeneous b-value are incorporated to predict the number of mainshocks with magnitude M W 5 in Figure 10(d). The expected number of M W 5 mainshocks estimated from the reference forecast model (Figure 10(a)) is more or less similar to the calculated seismicity rate r (Figure 7) due to the absence of stress perturbation (equations (2) and (3)). Figure 10(d) demonstrates that the forecasted seismicity rate is substantially well correlated with the observed seismicity during the time period 2011-2013, whereas the forecast model incorporated with the constant DCFF and average b-value could not infer the observed seismicity, as shown in Figure 10(a)-(c). Figure 10(d) clearly suggests that the forecasted seismicity is strongly influenced by DCFF and amplified background seismicity rate. However, the absence of observed mainshocks in high background seismicity rates can be primarily contributed to the shadow effect of DCFF. On the contrary, the forecast model also exhibits rather high seismicity in few areas of stress shadows. There are also some cases where model fails to explain observed seismicity rate increases, for example, in Zone 1. Such increase/decrease might be explained by dynamic Coulomb stress changes, which are sensitive to the direction of rupture propagation, local pore fluid effects and geothermal effect (Leptokaropoulos et al. 2014) or perhaps to locally anomalous rate/state parameters. Our forecasted seismicity rate could not explain the smaller changes due to small scale heterogeneities in DCFF. It is also noticed that the heterogeneity of spatial b-value increases the spatial heterogeneity of expected number of mainshock with magnitudes M W 5 (Figure 10(d)). Also, low b-value found in particular Kirthar ranges, Himalayan Frontal Thrusts and Tibetan Plateau region highlights the impact of the stress perturbation thereby increasing the production of expected number of lager mainshocks in the study region. More in-depth information regarding the effects of stress perturbation due to heterogeneous b-value can be obtained by simply examining Figures 9(d) and 10(b). However, the observed low b-value is solely responsible for high forecasted seismicity rate in southern Kirthar region of Zone1 and Tibetan plateau of Zone 5. Toda and Stein (2003) affirmed that the high background seismicity along with long observational period is necessary to examine the stress shadows in the region. However, it is quite complicated to detect the small number of earthquakes occurred in the stress shadows from the maps depicting the forecasted number of mainshocks. This may be due to strong dependence of computed forecasted seismicity rate on the assumed stress history and input parameters that can vary over several orders of magnitude. For proper understanding of the absence of mainshocks in stress shadows, Figure 11 is presented by subtracting the expected numbers of M W 5 earthquakes prepared with the assumed DCFF = 0 from the ones with spatially variable stress perturbations. Furthermore, a careful observation indicates that the contribution of older events (29 July 1980 Nepal earthquake) to the seismicity is very less while the recent events, i.e. 7 April 2005 Xizang earthquake, 8 October 2005 Kashmir earthquake and 28 October 2008 Pakistan earthquake have remarkably increased the  Based on this static stress triggering, we find that the elastic stress transfer along with high r and low b-value promote seismicity rate changes consistent with the occurrence time and the location of the mainshocks. Our results agree with the main findings of Nostro et al. (2005) and Catalli et al. (2008) who calculated static stress changes for the 1997 Umbria-Marche sequence (central Italy) using rate/state friction model. Hence, the approach we adopted provides satisfactory results in forecasting the number of mainshocks with magnitudes M W 5 in the NW frontier province of Himalaya during the period 2011-2013 (Figure 10(d)). In summary, the forecast model applied in present study provides promising results, as it was able to forecast future seismicity rates, in spite of the aforementioned assumptions and uncertainties.

Consistency of the forecast model
In this study, we also attempt to assess the goodness-of-fit statistics for the forecast and the observed data ( Figure 12). A histogram of 1000 simulated joint log-likelihoods along with observed data (dashed blue line) are plotted (Figure 12(a)). The intersection of the observed joint log-likelihood (dashed black line) with the corresponding empirical cumulative distribution functions of the simulated joint log-likelihoods (solid blue line) is also presented in (Figure 12(b)). It is apparent from these two figures that the quantile score z of S-test is found to be 72%, i.e. more than the critical region (shaded box) given by a = 2.5%. This implies that z lies outside the shaded critical region, thereby indicating the consistency of the spatial forecast and observed spatial distribution.

Conclusions
Coulomb stress transfer model incorporating rate and state dependent friction law is used in the study to forecast the seismicity rate of mainshocks with magnitudes M W 5 in the NW frontier of Himalaya during the period 2011-2013. For this purpose, we have computed DCFF, the most dominant parameter of the model, for each mainshock (M W 5) to examine the increase/decrease of seismicity rate due to stress perturbations. Also, other constituent parameters of this model such as smoothened background seismicity rate, aftershock duration,, and spatially variable b-value have been calculated for each source zone using compiled earthquake catalogue for the period 1975-2010. In this model, the heterogeneous b-value distribution tremendously increases the spatial heterogeneities of the forecasted mainshocks (M W 5) in the region. Although long observational periods and high background seismicity rate before perturbing mainshocks are essentially required, the stress shadows associated with substantial seismicity rate are difficult to be detected because of sustained stress changes in the study region. The DCFF distributions that exhibit significant stress increase nearby the mainshock sources are found to be diminished with the inverse of the distance from the fault rupture plane, associated with the large mainshock. Despite of the presence of uncertainties in the model parameters and inadequate knowledge of dynamic stress triggering (i.e. propagation of seismic waves, postseismic slip and creeping, pore fluid diffusion, viscoelastic relaxation, and rheological properties), the rate and state transfer model captures the observed number of M W 5 mainshocks fairly well. The technique is found much better than that estimated using the stress changes alone. In addition, S-test has been performed to check the stability of the results, which confirm that the spatial distribution of forecasted and observed earthquakes match well. However, some mismatch between the observed and modelled seismicity rates may be associated because of additional uncertainties arising from the rate and state parameters, such as, aftershock duration, spatially uniform values of As and the background seismicity rate.
However, the key aspects of the failure of forecasted seismicity rates map often depend on the poorly constrained parameters, whose values are chosen based on the preconceptions. In such case, forecasting of seismicity rates in any seismic regime will be poor. We can improve our model, if we choose the actual parameters based on the laboratory experimental value. Furthermore, recognizing the uncertainties would enable us to decide how much credence to place in the seismicity rates map and the forecast map should undergo rigorous and objective testing to compare their predictions to those of null hypotheses, including ones based on uniform regional seismicity.