Detection of gravity changes before powerful earthquakes in GRACE satellite observations

The variations of gravity field have been analyzed in this article, in order to find disturbances in the vicinity of recent great earthquakes epicenters including Chile (February 27, 2010), Tohoku-Oki (March 11, 2011) and Indian Ocean (April 11, 2012) prior to the events. For this purpose, the 10 years long time series of gravitational gradient components obtained from weekly Gravity Recovery And Climate Experiment (GRACE) solutions have been used. Some of gravitational gradient components are independent to GRACE stripy error and amplify high frequency components of gravity field. Therefore that preseismic activity can be better illustrated. The interquartile method has been used to construct the higher and lower bounds in time series to detect outlying solution outside the bounds associated with impending earthquakes. Afterwards, Nonlinear Auto-Regressive models with eXogenous inputs (NARX) neural network has been proposed in the detection process of prominent gravity field anomalies prior the earthquakes. Both methods detect considerable anomalous occurrences during 2-5 weeks prior to the earthquakes. Our results statistics show that the anomalous deviations before the earthquakes have different signs and amplitudes at different cases.


Introduction
Earthquakes are one of the most destructive disasters that affect normal life of millions of people every year.Therefore, prediction of earthquake before its occurrence can save the life of so many victims.
During an earthquake preparation, a dynamic process involves an energy transfer due to crust displacement and, at the time of the shock a breakdown between the source and the environment occurs.The change prior to the earthquake or along with it, may have different physical and chemical effects on the lithosphere, atmosphere and ionosphere, and accordingly makes it possible to be detected.These variations of lithosphere, atmosphere and ionosphere parameters before the main earthquakes are considered as a hint of impending earthquakes (earthquake precursors) [Akhoondzadeh et al. 2010].
Earthquakes prediction is mainly based on the observation of precursory phenomena.Extensive researches in domain of earthquake prediction have resulted in recognition of many earthquake precursors in the lithosphere, atmosphere and ionosphere.The effects of preseismic activity on the lithosphere can be investigated using gravity change.During the primary period of an earthquake, the gravity field of the seismic region may change due to the crust deformation or medium-property variation.Therefore, the effects of the pre-seismic activity on the gravity field can be investigated by measuring the variations of functional of gravity field.
Recently, several ground based observations have shown the possibility of gravity disturbance generation by earthquake preparation processes [e.g., Lan et al. 2011, Shen et al. 2011].There is lack of extensive ground experiments to monitor geophysical parameters in most areas.But due to the vast coverage of the seismic zones, satellite experiments are regarded as a promising mean to study earthquakes.
The Gravity Recovery and Climate Experiment (GRACE) has been the first space mission which can measure temporal variations of the Earth's gravity field.Therefore, it can be a useful tool for assessing changes in the gravity field related to earthquakes.The GRACE mission is a space borne gravity mission launched in March 2002 jointly by the National Aeronautics and space Administration (NASA) and German Aerospace Center (DLR).GRACE consists of two identical spacecrafts co-orbiting at low altitude (~450 km), separated by ~220 km and linked with K-band micro-wave ranging (KBR) satellite-to-satellite tracking (SST) system.
Little changes in the distance between GRACE satellites (measured by KBR) reflect mass anomalies within the Earth system and indicate changes in Earth's gravity field [Tapley et al. 2004].
Displacement taking place immediately during the earthquake are called coseismic.The dislocated masses disturb the isotactic equilibrium of the crust and mantle, thus including a long term postseismic relaxation process [Einarsson et al. 2010].Previous studies demonstrated coseismic and postseismic gravity field changes related to recent great earthquakes.Han et al. [2006] for the first time revealed coseismic gravity changes by Level-1 GRACE data (range and range-rate measurements) due to great Sumatra-Andaman (SA) earthquake (Mw = 9.1-9.3).Ogawa and Heki [2007] processed the GRACE 36 monthly gravitational field data (Level 2 or spherical harmonics) to derive coseismic and postseismic geoid variations due to Sumatra earthquake.Chen et al. [2007] showed that improved GRACE monthly gravity fields are able to reveal the coseismic changes with greater spatial resolution.They used variations of equivalent water height with P3M6 and 300 km Gaussian filter.Panet et al. [2007] showed coseismic and postseismic signal for both Sumatra and Nias earthquakes using wavelet analysis.Han and Simons [2008] used Slepian's spatiospectral localizing basis functions, as an alternative to the usual spherical harmonic functions.They estimated the coseismic jump at the time of the Sumatra earthquake.Han and Simons [2008] also use Slepian basis functions estimated from in-situ data.They determine a coseismic jump as well as an exponential-like post-seismic signal.de Linage et al. [2009] determined the earthquake signal of the Sumatra earthquake.They use the regularized CNES solution as well as the monthly CSR release four solution.They examine the possible effect of hydrological signals on the estimated earthquake signals by using hydrological models.They are able to detect both co-and postseismic signals of the SA earthquake.Other studies using GRACE to detect or constrain coseismic or postseismic deformations of great earthquake include: i) Einarsson et al. [2010], Broerse et al. [2011], Cambiotti et al. [2011], Wang et al. [2012a], Han et al. [2013] for the 2004 Sumatra-Andaman and Nias earthquake; ii) Han et al. [2010], Wang et al. [2012b], Han et al. [2013] for the 2010 Chile Maule event; iii) Han et al. [2011], Matsuo and Heki [2011], Wang et al. [2012c], Han et al. [2013] for the 2011 Tohoku-Oki earthquake; iv) Han et al. [2013] for 2012 Indian Ocean earthquake.
Despite all the efforts to detect gravity anomalies before large earthquakes [e.g., Seiji and Nakamura 2013], previous studies only detected coseismic and postseismic gravity field changes.However, there are no reports associated with detection of preseismic gravity changes using GRACE measurements.In this study, we analyzed GRACE observations before recent megathrust earthquake in order to detect gravity anomalies as an earthquake precursor.If it can be confirmed that gravity perturbations are a real and systematic occurrence, then it is possible to consider them as short-term precursors that happen before earthquakes.Since any particular earthquake precursor can not be considered as an accurate stand-alone tool for earthquake prediction and it should not be expected that every precursor appear in every earthquake therefore it is required to exploit and integrate different precursors from different experiments to achieve more accurate results.
In order to determine whether some gravity changes before earthquake have unique characteristics related to seismic activity, this study examines 10 years of gravitational gradient components obtained from GRACE satellite data in the seismic region.Section 2 will explain the method of gravitational gradient calculation by GRACE data, applied methods for anomaly detection are discussed in Section 3 and then seismo gravitational gradient anomalies observed for three last earthquakes with magnitude greater than 8.0 have been studied in Section 4.

GRACE data processing
In this study, we used weekly Level 2 (L2) global gravity solutions of GeoForschungsZentrum (GFZ, 30 Potsdam RL05) which are composed of fully normalized spherical harmonic coefficients complete to degree and order 30.A total 448 gravitational field solution covering from first week of January 2004 to third week of February 2013 are used.
The main reason for using weekly data instead of monthly solutions is that weekly data reveals short-time (shorter that one month) variations of gravity field in the seismic region.The Earth's oblateness values (C 20 ) has been replaced by those of the Satellite Laser Ranging (SLR) because of its accuracy [Cheng and Tapley 2004].Then in order to compute variations of gravity field the spherical harmonic coefficients of each individual weekly solutions have been subtracted from mean of solutions between January 2004 and December 2012.
For the study period, we computed the full gravitational gradient tensor (second derivative of gravitational potential) in spherical coordinates in 1 degree by 1 degree spherical regular grid.To better demonstrate the local mass anomaly, the local north-east-down frame (NED) at each point with spherical coordinate is introduced: x axis is directed to the north, the y-axis to the east, and the z-axis 1 downward.Based on the prin-cipal of converting between coordinate systems, the full gravitational 2 gradient tensor in this local NED frame is obtained [Wang et al. 2012a].
To obtain reasonable estimates of time-varying gravity filed signals, spherical harmonic coefficients should be filtered because high frequency components (high SH degree) are poorly determined [Han et al. 2005].This problem leads to creation of features that elongated in the north-south direction, usually described as striping patterns.Several authors have shown that non-isotropic filters [e.g., Swenson and Wahr 2006] account for the GRACE stripes more favourably compared with traditional isotropic Gaussian [Jeckeli 1981].Therefore, we make use of non-isotropic DDK3-filter [Kusche 2007, Kusche et al. 2009] that is specifically designed for this purpose.
Among five independent components of gravitational gradient tensor, some components are less sensitive to GRACE striping errors.As mentioned earlier, GRACE striping errors distributed in north-south direction, thus these strips generate fluctuations in the east-west direction [Swenson and Wahr 2006].The north-south component (x-component) of the gravitational gradient tensor has less sensitivity to GRACE striping errors and have more signal to noise ratio than other components.
As an example, we compute coseismic ∞V xz changes induced by the March 11, 2011, Mw 9.0 Tohoku-Oki earthquake using CSR Release (RL) 05 GRACE Level 2 monthly gravitational field data products without applying any destriping and spatial smoothing filtering procedure and the achieved results are shown in Figure 1.
As it can be seen in Figure 1, tacking derivative related to x-axis caused dramatically suppress stripes errors.Similar results can be obtained for ∞V xx component.Due to the less sensitivity behaviour of GRACE striping errors and more signal to noise ratio rather than other components.By considering this fact, we have used only ∞V xx and ∞V xz components for time series calculation.
In addition for geophysical phenomena of regional scale which are dominated by high frequency content (such as earthquake), gravitational gradient changes are more useful.These components enhance high frequency components of the Earth gravity field and reveal more details in temporal domain.
The changes in GRACE residues due to hydrological effect can be equal or even greater than the seismic effect of earthquake.These effects are larger in continent area compared to offshore area.In the final solutions models that were used for de-aliasing GRACE raw data, some errors were introduced at long periods which is the case for ocean tidal models drologic effects and their amplitude is often at noise level [de Linage et al. 2009].
From these considerations, to separate the abovementioned effects in the GRACE data, we adopted the following strategy.At each point of a 1*1 grid, we simultaneously fitted to the selected gravitational gradient components time series the following time-function [Chen et al. 2007]: (1) Where t is time and t 0 is time of earthquake occur and model parameters are: 1) a 1 , { 1 , a 2 , { 2 are amplitude and phase of the annual and semi-annual waves to model the seasonal and annual variations of hydrology and long period oceanic circulation before earthquake 2) a 3 , { 3 are amplitude and phase of a 161-d sine curve to correct the errors on the S2 tidal wave before earthquake 3) a 4 , { 4 , a 5 , { 5 are amplitude and phase of the annual and semi-annual waves to model the seasonal and annual variations of hydrology and long period oceanic circulation after earthquake 4) a 6 , { 6 are amplitude and phase of a 161-d sine curve to correct the errors on the S2 tidal wave after earthquake 5) B and B are linear trends before and after earthquake 6) A and A are mean of signal before and after earthquake.
We computed these parameters by nonlinear leastsquares adjustment.Then we removed effects of linear trend, annual, semi-annual and S2 tidal wave from each time series of each grid point before and after earthquake.

Methodology
At each case study, using the reported geographic latitude and longitude concerning the earthquake epicenter, time series of gravitational gradient components (∞V xx and ∞V xz ) obtained from GRACE measurements during some years before earthquake has been analyzed.By considering this fact that other preseismic anomalies (e.g.ionosphere precursors) does not occur in the vertical projection of earthquake epicenter, outskirt of each epicenter has been tested in order to detect anomaly.In order to detect anomalies in a time series, a reasonable range of gravitational gradient variations are required to be determined.Therefore Inter-Quartile Range (IQR) and Artificial Neural Network (ANN) methods have been implemented to observe the discord patterns in the gravity variations.

Inter-quartile range method
The Inter-quartile range of data is used as the primary method for anomaly detection.The upper and lower bounds of time series is calculated using the following equations: (2) (3) Where x, x high , x low , M, IQR and Dx are data value, upper bound, lower bound, median value, inter-quartile range of data and differential of x, respectively.According to this, if the absolute value of the data would be greater than k.IQR (|x|> k.IQR), the behaviour of the relevant parameter (x) is regarded as anomaly.According to Equation (4), indicates the percentage of parameter change from the undistributed state.Therefore in the anomaly detection from ∞V xx and ∞V xz time series, M and IQR values were computed over total time interval of interest for each case study.

Artificial Neural Network method
A time series is a sequence of vectors that depend on time.In order to forecast time series of gravitational gradient variations before and after occurrence of each earthquake, ANN technique has been implemented.
The problem of prediction can be formulated as finding a function f to obtain an estimated of the vector y at time t + p (p = 1,2,...), given the values of y up to time t: Where y(t) represents output of model as time t, d y is the lag of the system and f a nonlinear function.Typically, p=1, measure one-step ahead, but can take any value larger than 1 at multi-step ahead [Dorffner 1996].
Before 1980s, for time series prediction, methods of linear parametric AutoRegressive (AR), Moving-Average (MA) or AutoRegressive Moving-Mverage (ARMA) introduced by Box and Jenkis [Pollock et al. 1999 models are linear and are unable to cope with certain non-stationary signals, and signals whose mathematical model is not linear.Artificial neural networks are a type of flexible nonlinear models that can discover patterns adoptively from the data.
In many researches, ANN has been proposed as a promising approach for time series forecasting.Several successful applications of this method demonstrated that ANN can be a very effective tool in modeling and forecasting non-linear time series [Akhoondzadeh 2013].
Among equations which require input/output relationships, ANNs are trained by sample data to build the input/output vector maps in an implicit way.Therefore, ANNs can solve highly nonlinear problems without the need to define the relationship between inputs and outputs [Yang et al. 1997].
A neural network can be classified into static and dynamic categories.The static feed forward network has no feedback elements and contains no delays (Figure 2).
In these networks output is calculated directly from the input through feed forward connections like Back Propagation (BP) neural networks.
In general, BP is effective only when the network architecture is chosen correctly.Too small network cannot learn the problem well, but too large sized network leads to over fitting and poor generalization performance.In dynamic networks, the input depends not only on the current input to the network, but also on the current or previous, outputs or states of the network.
In this paper, the architectural approach proposed to predict gravitational gradient changes time series is one based on "Nonlinear Autoregressive models with eXogenous inputs (NARX model)", which are therefore called NARX recurrent neural networks [Menezes et al. 2006].This is a powerful class of dynamic neural networks and it has been proven that they are well suited for modeling nonlinear systems and specially time series [Diaconescu 2008].
The NARX model is commonly used in time series modeling, where the next value of the dependent output signal y(t) is regressed using the previous values of the output signal and previous values of an independent input signal.Figure 3 shows the basic architecture of a NARX model capable in time series forecasting.The defining equation for the NARX model is: In order to determine the best network configuration, the effective parameters, which influence the value of predictive error, including the number of pattern of input; lags of input and output; the number of hidden layers and their number of neurons; the activation functions and the learning algorithm have been obtained via iteration process to assess the minimum predictive error when the training process was implemented.
The best results indicate the linear and sigmoid activation functions for output and hidden layers, respectively.In general, in function approximation problems, for networks that contain up to a few hundred weights, the Levenberg-Marquardt algorithm will have the faster convergences [Diaconescu 2008].Therefore the Levenberg-Marquardt optimization has been selected as training algorithm.However other training algorithms were evaluated and the results were poorer.
To implement prediction process, N observations y 1 ,y 2 ,...,y N are selected as training set and the remain ones y N+1 ,y N+2 ,...,y N+m are considered as test set.The number of input nodes corresponds to the number of lagged observations to discover the underlying pattern in a time series.In this study a network with six nodes in the input layer, 15 nodes in the first hidden layer, 10 nodes in the second hidden layer and one node in output layer has been proposed.Usually, the evaluation of prediction performance is done by finding optimized weights such the prediction error (PE) is minimized.PE equation can be written as: (5) Where is the output of network.The testing patterns are, Finally, if the difference value between the actual value and predicted values is outside the predefined bounds, the anomaly is detected.

Implementation
In order to ensure reliable detection of gravity field variation using GRACE data before the earthquake, results for three most recent earthquakes with magnitude greater than 8.0 in Maule (Chile), Tohoku-Oki ( Japan), and Indian Ocean have been studied.

Chile Maule earthquake (February 27, 2010)
The February 27, 2010, Maule Chilean earthquake (Mw = 8.8), which was caused by the subduction of the Nazla plate underneath the overlying South America plate, is the sixth largest event in the seismic record (Figure 4).Data from teleseismic networks, coastal/river markers, tsunami sensors, Global Positioning System (GPS), and Interfero-metric Synthetic Aperture Radar (InSAR) have been used to observe and model the coseismic signature and slip history of this devastating event.In addition, spaceborne gravimetry data from GRACE have been used to observe the coseismic deformation signal of 2010 Maule earthquake [Han et al. 2010, Wang et al. 2012b, Han et al. 2013].
Figures 5a and 5b illustrate the variations of ∞V xx and ∞V xz components, obtained by GRACE measurements from first week of January 2004 to third week of February 2013.A vertical dotted line indicates the earthquake time.The x-axis represents actual time and y-axis represents variations of ∞V xx and ∞V xz components by scale of 10 13 (unit is s -2 ).
By visual inspection in both time series, we can see a strange decreasing in both signals two weeks before the earthquake date.
It should be noted that there is a horizontal derivative operator shift phase of mass anomaly in the spatial domain.Therefore we should not expect that preseismic anomaly for different components of gravitational gradient components occur in the same location.Therefore, the locations for time series calculation for ∞V xx and ∞V xz components are different.
When implementing the interquartile method, Dx is calculated using Equation (4) for both components.The ∞V xx and ∞V xz values exceed the lower bounds (M-2.6*IQR for ΔV xx and M-3.1*IQR for ΔV xz ) two weeks before earthquake with values of 29% and 30% of the lower bounds, respectively.
Figures 6a and 6b represent the observed and predicted ∞V xx and ∞V xz values using NARX neural network, respectively, during the time period of 2008 to 2013.
Figures 6c and 6d represent the differences between observed and predicted values during the testing data.
These figures clearly show unusual differences two weeks prior to the earthquake.The ∞V xx and ∞V xz values exceed the lower bounds with values of 14% and 49%.In other words, on these weeks, the NARX neural network has not been able to predict the gravitational gradient variation based on model deduced from training data.Therefore, the observed values at mentioned week could be considered as earthquake anomaly.
Thus, as indicated previously, the two applied methods on gravitational gradient time series of Chilean   earthquake detect a noticeable anomalies two weeks prior to the earthquake.

Tohoku-Oki earthquake (March 11, 2011)
The Tohoku-Oki earthquake, Mw 9.0, which occurred at 05:46 UT, on March 11, 2011, at the Japan Trench east of NE Japan (Figure 7), ruptured the fault as large as 500~200 km [Ozawa et al. 2011].Earthquake of this size and magnitude cause mass-dislocation on the scale that is measurable by GRACE gravity satellite mission.Some articles demonstrated coseismic and postseismic gravity changes due to this earthquake [Han et al. 2011, Matsuo and Heki 2011, Wang et al. 2012c, Han et al. 2013   three months before this earthquake has been measured by shipboard geophysical observations.But there are no reports showing preseismic gravity changes in GRACE measurements.
By implementing the Interquartile method (Equation 4), the Dx values for ∞V xx and ∞V xz exceed the lower bounds, three weeks before earthquake with the values of 18% and 26% , respectively (Figure 8).
After the event there are some outlier points of predefined bounds.Due to the changes of the mass distribution in the studied region, the same behavior should not be expected from the signal before and after the earthquake in the seismic region.In order to analyze signals of gravity field after the earthquake for anomaly detection, the median and Interquartile related to this period should be used.
Figures 9a and 9b represent the observed and predicted ∞V xx and ∞V xz values using NARX neural network, respectively, during the time period of 2008 to 2013.Figures 9c and 9d represent the differences between observed and predicted values during the testing data.These figures clearly show unusual differences three weeks prior earthquake.The ∞V xx and ∞V xz values exceed the lower bounds with values of 30% and 38%.In other words, on these weeks, the NARX neural network has not been able to predict the gravitational gradient variations based on model deduced from training data.Therefore, the observed values at mentioned weeks could be considered as earthquake anomaly.
Thus, as indicated previously the two applied methods on gravitational gradient time series of Tohoku-Oki earthquake detect a noticeable anomaly on three weeks prior to earthquake.

Indian Ocean earthquake (April 11, 2012)
The Mw 8.6 earthquake occurred on April 11, 2012, at the depth of 40 km (Global CMT), located at NW Indian, Australian and Sunda plates (~500km, west  of December 26, 2004, earthquake) (Figure 10).It is the largest interplate strike-slip earthquake in the known history.On the same day, with delay of ~2 hours, ~120 km away from the main event, another earthquake of Mw 8.1 occurred [Shrivastava and Reddy 2013].Some authors compute coseismic displacement [Shrivastava and Reddy 2013] and gravity changes observed by GRACE measurements [Han et al. 2013].
Figures 11a and 11b show the time series of ∞V xx and ∞V xz variations during 2005 to 2013 obtained from GRACE measurements.The sudden jump in earthquake date represents coseismic jump.In this case study, since the December 26, 2004, Sumatra-Andaman earthquake (Mw 9.1-9.3)occurred near the epicentre of this earthquake, therefore we computed time series during 2005 (one month after Sumatra earthquake) to 2013.
By implementing the interquartile method the Dx values for ∞V xx exceed the upper bound five weeks before earthquake with values of 18% of higher bound.On three weeks prior the earthquake, the ∞V xz value passes the upper bound by 24% .
Figures 12a and 12b represent the observed and predicted ∞V xx and ∞V xz values using NARX neural network, respectively, during the time period of 2008 to 2013.
Figures 12c and 12d represent the differences between observed predicted values during the testing data.These figures clearly show differences five weeks prior earthquake.The ΔV xx value exceeds the upper bound with values of 56%.There is another anomaly that has been detected by NARX neural network with values of 30% three weeks before earthquake that IQR method has not been able to detect it.The ΔV xz value exceeds the upper bound with value of 68% of higher bound five before earthquake.There is another anomaly has been detected by NARX neural network three weeks before earthquake with value of 60% that IQR method has not been able to detect it.In other words, on these weeks, NARX neural network was unable to predict the gravitational gradient variations based on the model that was deduced from training data.Therefore, the observed values at mentioned weeks could be considered as earthquake anomaly.Thus, as indicated previously on gravitational gradient time series of Indian Ocean earthquake, a noticeable anomaly was detected in three and five weeks prior to the earthquake.

Conclusion
It has been shown in this study that weekly time variable gravity solutions obtained from GRACE satellite data are able to detect striking anomalies in the gravity field, in the vicinity of recent great earthquake epicenters several weeks before their happenings.Traditionally, to examine preseismic gravity field activity in the seismic region, gravitational changes have been analyzed.Due to GRACE stripy errors, analysis of gravitational changes does not show any anomalies before earthquakes.
In addition, the efficiency of integrating interquartile and NARX neural network methods to detect anomalies in gravitational gradient time series has been shown in this study.It may be worthwhile to mention that ANN has been used for the first time in this study to successfully detect anomalies in time series of gravity field variations.
The results indicate that highest deviations from normal state that were regarded as anomaly appeared within the time interval 2-5 weeks before earthquakes.It should be pointed out that preseismic gravity anomalies before earthquake can be positive as well as negative.
In this study, it has been shown that detected anomalies using applied methods can be related to the three recent great earthquake since variations of grav-ity field represented by gravitational gradient changes shows the anomalies a few weeks prior these earthquakes during several years data over the region of impending earthquakes.
Figure 1.GRACE global ∞V xz changes between of 2011 and 2012 and mean of 2009 and 2010 without applying any destriping and spatial smoothing filtering procedure (in unit of mili-Eotvos).

Figure 3 .
Figure 3. NARX neural network for time series prediction.Figure 2. A Static feed forward neural network.

Figure 2 .
Figure 3. NARX neural network for time series prediction.Figure 2. A Static feed forward neural network.
Figure 6.(a) and (b) variations of the observed (green curve) and predicted (blue curve) gravitational gradient obtained from NARX method.(c) and (d) Variations of the differences between the observed and the predicted values of ∞V xx and ∞V xz .

Figure 5 .
Figure 5. Weekly time series of ∞V xx and ∞V xz components obtained from GRACE Level 2 data after DDK3 filtering.Seasonal (annual and semiannual and S2 tidal wave) signals and linear trend are removed from time series using least squares.The black vertical lines indicate the date of the earthquake.The red horizontal lines indicate the upper and lower bounds.The green horizontal lines indicate median values.

Figure 7 .
Figure 7. Plate configurations of the Japanese islands.The focal mechanism of the Tohoku-Oki earthquake is taken from the Global Centroid-Moment-Tensor Project.The red arrows indicate relative motion between the two plates at a plate boundary [Ozawa et al. 2011].

Figure 12 .
Figure 12.Same as Figure 6 but for Indian earthquake.

Figure 11 .
Figure 11.Same as Figure 5 for Indian earthquake.