SHAHRISVAND_finalonline_Layout 6 ANNALS OF GEOPHYSICS, 57, 5, 2014, A0543; doi:10.4401/ag-6612 A0543 Detection of gravity changes before powerful earthquakes in GRACE satellite observations Mohsen Shahrisvand1, Mehdi Akhoondzadeh1,*, Mohammad Ali Sharifi2 1 University of Tehran, College of Engineering, Surveying and Geomatics Engineering Dept., Remote Sensing Division, Iran 2 University of Tehran, College of Engineering, Surveying and Geomatics Engineering Dept., Geodesy Division, Iran ABSTRACT 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 inde- pendent 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 asso- ciated with impending earthquakes. Afterwards, Nonlinear Auto-Re- gressive 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 oc- currences 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. 1. Introduction Earthquakes are one of the most destructive dis- asters that affect normal life of millions of people every year. Therefore, prediction of earthquake before its oc- currence can save the life of so many victims. During an earthquake preparation, a dynamic process involves an energy transfer due to crust dis- placement 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 lith- osphere, 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 re- searches in domain of earthquake prediction have re- sulted in recognition of many earthquake precursors in the lithosphere, atmosphere and ionosphere. The ef- fects of preseismic activity on the lithosphere can be investigated using gravity change. During the primary period of an earthquake, the gravity field of the seis- mic 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 in- vestigated 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 promis- ing 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 space- crafts co-orbiting at low altitude (~450 km), separated by ~220 km and linked with K-band micro-wave rang- ing (KBR) satellite-to-satellite tracking (SST) system. Article history Received June 24, 2014; accepted September 4, 2014. Subject classification: Earthquake, Precursor, Gravity, Grace, Ionosphere, Gravity anomalies. Little changes in the distance between GRACE satel- lites (measured by KBR) reflect mass anomalies within the Earth system and indicate changes in Earth’s grav- ity 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 man- tle, thus including a long term postseismic relaxation process [Einarsson et al. 2010]. Previous studies demon- strated 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 measure- ments) due to great Sumatra-Andaman (SA) earth- quake (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 post- seismic 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 vari- ations 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 earth- quakes using wavelet analysis. Han and Simons [2008] used Slepian’s spatiospectral localizing basis functions, as an alternative to the usual spherical harmonic func- tions. 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 expo- nential-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 hydrologi- cal models. They are able to detect both co- and post- seismic signals of the SA earthquake. Other studies using GRACE to detect or constrain coseismic or post- seismic 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 con- firmed that gravity perturbations are a real and sys- tematic occurrence, then it is possible to consider them as short-term precursors that happen before earth- quakes. 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 there- fore it is required to exploit and integrate different pre- cursors 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 calcu- lation 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. 2. 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 normal- ized 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 (C20) 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 gravi- tational gradient tensor (second derivative of gravita- tional 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- SHAHRISVAND ET AL. 2 3 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 de- scribed as striping patterns. Several authors have shown that non-isotropic filters [e.g., Swenson and Wahr 2006] account for the GRACE stripes more favourably com- pared with traditional isotropic Gaussian [Jeckeli 1981]. Therefore, we make use of non-isotropic DDK3-filter [Kusche 2007, Kusche et al. 2009] that is specifically de- signed for this purpose. Among five independent components of gravita- tional gradient tensor, some components are less sen- sitive to GRACE striping errors. As mentioned earlier, GRACE striping errors distributed in north-south di- rection, thus these strips generate fluctuations in the east-west direction [Swenson and Wahr 2006]. The north-south component (x-component) of the gravita- tional 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 ∞Vxz 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 apply- ing any destriping and spatial smoothing filtering pro- cedure and the achieved results are shown in Figure 1. As it can be seen in Figure 1, tacking derivative re- lated to x-axis caused dramatically suppress stripes er- rors. Similar results can be obtained for ∞Vxx component. Due to the less sensitivity behaviour of GRACE strip- ing errors and more signal to noise ratio rather than other components. By considering this fact, we have used only ∞Vxx and ∞Vxz 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 fre- quency components of the Earth gravity field and re- veal more details in temporal domain. The changes in GRACE residues due to hydrolog- ical effect can be equal or even greater than the seismic effect of earthquake. These effects are larger in conti- nent area compared to offshore area. In the final solu- tions 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 [Ray and Luthchke 2006]. Model errors of the S2 tidal wave cre- ate an alias at 161-d period that is clearly visible in the GRACE signal. Another source of variability in GRACE residues are seasonal and inter-annual changes in ocean circulation. These variations are smaller than other hy- GRAVITY CHANGES BEFORE POWERFUL EARTHQUAKES Figure 1. GRACE global ∞Vxz 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). drologic effects and their amplitude is often at noise level [de Linage et al. 2009]. From these considerations, to separate the above- mentioned effects in the GRACE data, we adopted the following strategy. At each point of a 1*1 grid, we si- multaneously fitted to the selected gravitational gradi- ent components time series the following time-function [Chen et al. 2007]: (1) Where t is time and t0 is time of earthquake occur and model parameters are: 1) a1, {1, a2, {2 are amplitude and phase of the an- nual and semi-annual waves to model the seasonal and annual variations of hydrology and long period oceanic circulation before earthquake 2) a3, {3 are amplitude and phase of a 161-d sine curve to correct the errors on the S2 tidal wave before earthquake 3) a4, {4, a5, {5 are amplitude and phase of the an- nual and semi-annual waves to model the seasonal and annual variations of hydrology and long period oceanic circulation after earthquake 4) a6, {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 earth- quake 6) A and A are mean of signal before and after earthquake. We computed these parameters by nonlinear least- squares 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 earth- quake. 3. Methodology At each case study, using the reported geographic latitude and longitude concerning the earthquake epi- center, time series of gravitational gradient components (∞Vxx and ∞Vxz) obtained from GRACE measurements during some years before earthquake has been ana- lyzed. By considering this fact that other preseismic anomalies (e.g. ionosphere precursors) does not occur in the vertical projection of earthquake epicenter, out- skirt of each epicenter has been tested in order to de- tect anomaly. In order to detect anomalies in a time series, a reasonable range of gravitational gradient vari- ations 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. 3.1. Inter-quartile range method The Inter-quartile range of data is used as the pri- mary method for anomaly detection. The upper and lower bounds of time series is calculated using the fol- lowing equations: (2) (3) Where x, xhigh, xlow, M, IQR and Dx are data value, upper bound, lower bound, median value, inter-quar- tile range of data and differential of x, respectively. Ac- cording 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. Ac- cording to Equation (4), indicates the percentage of parameter change from the undistributed state. Therefore in the anomaly detection from ∞Vxx and ∞Vxz time series, M and IQR values were computed over total time interval of interest for each case study. 3.2. 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, dy 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-Av- erage (MA) or AutoRegressive Moving-Mverage (ARMA) introduced by Box and Jenkis [Pollock et al. 1999]. These .x M k IQRhigh = + .x M k IQRlow = - p k Dx k 100!= -c m y t p+t^ h , ...,y t p f y t y t dy+ = -t^ ^ ^^h h hh ;x x x k IQR x M k Dx IQR x M low high &1 1 1 1- - = - SHAHRISVAND ET AL. 4 , cos cos y t t A Bt a t A B t a t if t t if t t i i i i i i i i 0 1 3 3 3 1 3 0 0 # # ~ { ~ { = + + + + + + = + + = l l ^ ^ ^ h h h Z [ \ ]] ]] / / (4) 5 models are linear and are unable to cope with certain non-stationary signals, and signals whose mathemati- cal model is not linear. Artificial neural networks are a type of flexible nonlinear models that can discover pat- terns 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 fore- casting non-linear time series [Akhoondzadeh 2013]. Among equations which require input/output re- lationships, ANNs are trained by sample data to build the input/output vector maps in an implicit way. There- fore, ANNs can solve highly nonlinear problems with- out 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 net- work leads to over fitting and poor generalization per- formance. 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 net- works and it has been proven that they are well suited for modeling nonlinear systems and specially time se- ries [Diaconescu 2008]. The NARX model is commonly used in time series modeling, where the next value of the dependent out- put signal y(t) is regressed using the previous values of the output signal and previous values of an independ- ent 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 configu- ration, the effective parameters, which influence the value of predictive error, including the number of pat- tern of input; lags of input and output; the number of hidden layers and their number of neurons; the acti- vation functions and the learning algorithm have been obtained via iteration process to assess the minimum predictive error when the training process was imple- mented. The best results indicate the linear and sigmoid ac- tivation functions for output and hidden layers, respec- tively. 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 algo- rithms were evaluated and the results were poorer. To implement prediction process, N observations y1,y2,...,yN are selected as training set and the remain ones yN+1,yN+2,...,yN+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 out- put layer has been proposed. The training patterns in the proposed network are: , , ..., , , , ..., y t f y t y t y t n x t x t x t n 1 2 1 2 y x = - - - - - - ^ ^ ^ ^^ ^ ^ ^ h h h h h h hh GRAVITY CHANGES BEFORE POWERFUL EARTHQUAKES Figure 3. NARX neural network for time series prediction.Figure 2. A Static feed forward neural network. 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. 4. Implementation In order to ensure reliable detection of gravity field variation using GRACE data before the earthquake, re- sults for three most recent earthquakes with magnitude greater than 8.0 in Maule (Chile), Tohoku-Oki ( Japan), and Indian Ocean have been studied. 4.1. 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 (Fig- ure 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 co- seismic signature and slip history of this devastating event. In addition, spaceborne gravimetry data from GRACE have been used to observe the coseismic de- formation 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 ∞Vxx and ∞Vxz components, obtained by GRACE measure- ments 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 ∞Vxx and ∞Vxz compo- nents by scale of 1013 (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 deriv- ative operator shift phase of mass anomaly in the spa- tial domain. Therefore we should not expect that preseismic anomaly for different components of gravi- tational gradient components occur in the same loca- tion. Therefore, the locations for time series calculation for ∞Vxx and ∞Vxz components are different. When implementing the interquartile method, Dx is calculated using Equation (4) for both components. The ∞Vxx and ∞Vxz values exceed the lower bounds (M- 2.6*IQR for ΔVxx and M-3.1*IQR for ΔVxz) two weeks before earthquake with values of 29% and 30% of the lower bounds, respectively. Figures 6a and 6b represent the observed and predicted ∞Vxx and ∞Vxz 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 ∞Vxx and ∞Vxz val- ues 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 train- ing data. Therefore, the observed values at mentioned week could be considered as earthquake anomaly. Thus, as indicated previously, the two applied meth- ods on gravitational gradient time series of Chilean PE y t k t ky k N 0 = - -- - t^ ^^ h hh/ yt , , , , , , , , , , , , , , , y f y y y t t t y f y y y t t t y f y y y t t t N N N N N N N N N N N N N N N m N m N m N m N m N m N m 4 1 2 3 1 2 3 5 2 3 4 2 3 4 3 2 1 3 2 1 = = = + + + + + + + + + + + + + + + + - + - + - + - + - + - ^ ^ ^ h h h , , , , , , , , , , , , , , , y f y y y y f y y y t t t y f y y y t t t t t t 2 2N N N N N N N 4 1 2 3 1 2 3 5 2 3 4 2 3 4 3 1 3 1 = = = - - - - - - ^ ^ ^ h h h SHAHRISVAND ET AL. 6 Figure 4. The geographic location of Maule earthquake epicenter. A vacant star represents epicenter. (Figure taken from: http://earth quake.usgs.gov). 7 earthquake detect a noticeable anomalies two weeks prior to the earthquake. 4.2. Tohoku-Oki earthquake (March 11, 2011) The Tohoku-Oki earthquake, Mw 9.0, which oc- curred 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]. There are several publications indicating the preseismic TEC anomalies concerning Tohoku-Oki earthquake [e.g., Heki and Enomoto 2011, Akhoondzadeh 2013]. In addition Tsuboi and Naka- mura in 2013 demonstrate sea surface gravity change GRAVITY CHANGES BEFORE POWERFUL EARTHQUAKES 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 ∞Vxx and ∞Vxz. Figure 5. Weekly time series of ∞Vxx and ∞Vxz 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. Plate configurations of the Japanese islands. The focal mech- anism of the Tohoku-Oki earthquake is taken from the Global Cen- troid-Moment-Tensor Project. The red arrows indicate relative motion between the two plates at a plate boundary [Ozawa et al. 2011]. three months before this earthquake has been meas- ured by shipboard geophysical observations. But there are no reports showing preseismic gravity changes in GRACE measurements. By implementing the Interquartile method (Equa- tion 4), the Dx values for ∞Vxx and ∞Vxz 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 dis- tribution 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 ana- lyze signals of gravity field after the earthquake for anomaly detection, the median and Interquartile re- lated to this period should be used. Figures 9a and 9b represent the observed and pre- dicted ∞Vxx and ∞Vxz values using NARX neural net- work, respectively, during the time period of 2008 to 2013. Figures 9c and 9d represent the differences be- tween observed and predicted values during the testing data. These figures clearly show unusual differences three weeks prior earthquake. The ∞Vxx and ∞Vxz val- ues 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 train- ing 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 To- hoku-Oki earthquake detect a noticeable anomaly on three weeks prior to earthquake. 4.3. 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 SHAHRISVAND ET AL. 8 Figure 9. Same as Figure 6 for Tohoku-Oki earthquake. Figure 8. Same as Figure 5 for Tohoku-Oki earthquake event. 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 ∞Vxx and ∞Vxz 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 earth- quake, therefore we computed time series during 2005 (one month after Sumatra earthquake) to 2013. By implementing the interquartile method the Dx values for ∞Vxx exceed the upper bound five weeks be- fore earthquake with values of 18% of higher bound. GRAVITY CHANGES BEFORE POWERFUL EARTHQUAKES Figure 10. The geographic location of Indian earthquake epicen- ter. A vacant star represents epicenter. (Figure taken from: http:// earthquake.usgs.gov). Figure 12. Same as Figure 6 but for Indian earthquake. Figure 11. Same as Figure 5 for Indian earthquake. On three weeks prior the earthquake, the ∞Vxz value passes the upper bound by 24% . Figures 12a and 12b represent the observed and pre- dicted ∞Vxx and ∞Vxz values using NARX neural network, respectively, during the time period of 2008 to 2013. Figures 12c and 12d represent the differences be- tween observed predicted values during the testing data. These figures clearly show differences five weeks prior earthquake. The ΔVxx 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 ΔVxz value exceeds the upper bound with value of 68% of higher bound five before earthquake. There is another anom- aly 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. There- fore, the observed values at mentioned weeks could be considered as earthquake anomaly. Thus, as indicated previously on gravitational gradient time series of In- dian Ocean earthquake, a noticeable anomaly was de- tected in three and five weeks prior to the earthquake. 5. Conclusion It has been shown in this study that weekly time variable gravity solutions obtained from GRACE satel- lite data are able to detect striking anomalies in the gravity field, in the vicinity of recent great earthquake epicenters several weeks before their happenings. Tra- ditionally, to examine preseismic gravity field activity in the seismic region, gravitational changes have been analyzed. Due to GRACE stripy errors, analysis of grav- itational changes does not show any anomalies before earthquakes. In addition, the efficiency of integrating in- terquartile and NARX neural network methods to de- tect anomalies in gravitational gradient time series has been shown in this study. It may be worthwhile to men- tion 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 earth- quakes during several years data over the region of im- pending earthquakes. References Akhoondzadeh, M., M. Parrot and M.R. Saradjian (2010). Electron and ion density variations before strong earthquakes (M>6.0) using DEMETER and GPS data, Natural Hazards and Earth System Sciences, 10, 7-18; doi:10.5194/nhess-10-7-2010. Akhoondzadeh, M. (2013). A MLP neural network as an investigator of TEC time series to detect seismo- ionospheric anomalies, Advances in Space Research, 51, 2048-2057; http://dx.doi.org/10.1016/j.asr.2013. 01.012. Broerse, D.B.T., L.L.A. Vermeersen, R.E.M. Riva and W. van der Wal (2011). Ocean contribution to co- seismic crustal deformation and geoid anomalies: application to the 2004 December 26 Sumatran-An- daman earthquake, Earth and Planetary Science Letters, 305, 341-349. Cambiotti, G., A. Bordoni, R. Sabadini and L. Colli (2011). GRACE gravity data help constraining seis- mic models of the 2004 Sumatran earthquake, Jour- nal of Geophysical Research, 116, B10403; doi:10.10 29/2010JB007848. Chen, J.L., C.R. Wilson, B.D. Tapley and S. Grand (2007). GRACE detects coseismic and postseismic deformation from the Sumatra-Andaman earth- quake, Geophysical Research Letters, 34, L13302; doi:10.1029/2007gl030356. Cheng, M., and B.D. Tapley (2004). Variations in the Earth’s oblateness during the past 28 years, Journal of Geophysical Research, 109, B09402; doi:10.1029/ 2004JB003028. de Linage, C., L. Rivera, J. Hinderer, J.P. Boy, Y. Rogis- ter, S. Lambotte and R. Biancale (2009). Separation of coseismic and postseismic gravity changes for the 2004 Sumatran earthquake from 4.6 yr of GRACE observations and modelling of the coseismic change by normal mode summation, Geophysical Journal International, 176, 695-714. Diaconescu, E. (2008). The use of NARX neural net- works to predict chaotic time series, WSEAS Trans- actions on Computer Research, 3, 182-191. Dorffner, G. (1996). Neural networks for time series processing, Neural Network World, 4/96, 447-468. Einarsson, I., A. Hoechner, R. Wang and J. Kusche (2010). Gravity changes due to the Sumatra–Andaman and Nias earthquakes as detected by the GRACE satel- lites: a reexamination, Geophysical Journal Interna- tional, 183, 733-747; doi:10.1111/j.1365-246X.2010. SHAHRISVAND ET AL. 10 11 04756.x. Han, S.-C., C.K. Shum, C. Jekili, C.-Y. Kuo and C. Wilson (2005). Nonisotropic filtering of GRACE temporal gravity for geophysical signal enhancement, Geo- physical Journal International, 163, 18-25. Han, S.-C., C.K. Shum, M. Bevis, C. Ji and C.-Y. Kuo (2006). Crustal dilatation observed by GRACE after the 2004 Sumatra Andaman earthquake, Science, 313 (5787), 658-662; doi:10.1126/science.1128661. Han, S.-C., and F.J. Simons (2008). Spatiospectral local- ization of global geopotential fields from the Grav- ity Recovery and Climate Experiment (GRACE) reveals the coseismic gravity change owing to the 2004 Sumatra–Andaman earthquake, Journal of Geophysical Research, 113, B01405; doi:10.1029/20 07JB004927. Han, S.-C., J. Sauber and S. Luthcke (2010). Regional gravity decrease after the 2010 Maule (Chile) earth- quake indicates large-scale mass redistribution, Geo- physical Research Letter, 37, L23307; doi:10.1029/ 2010GL045449. Han, S.-C., J. Sauber and R. Riva (2011). Contribution of satellite gravimetry to understanding seismic source processes of the 2011 Tohoku-Oki earthquake, Geophysical Research Letter, 38, L24312; doi:10.10 29/2011GL049975. Han, S.-C, R. Riva, J. Sauber and E. Okal (2013). Source parameter inversion for recent great earthquakes from a decade-long observation of global gravity fields, Journal of Geophysical Research, 118, 1240- 1267; doi:10.1002/jgrb.50116. Heki, K., and Y. Enomoto (2013). Preseismic ionospheric electron enhancements revisited, Journal of Geo- physical Research: Space Physics, 118, 6618-6626. Jekeli, C. (1981). Alternative methods to smooth the Earth’s gravity field, Technical Report 327, Geodetic Science, Ohio State Univ., Columbus, OH. Kusche, J. (2007). Approximate decorrelation and non- isotropic smoothing of time-variable GRACE-type gravity field models, Journal of Geodesy, 81, 733- 749; http://dx.doi.org/10.1007/s00190-007-01433. Kusche, J., R. Schmidt, S. Petrovic and R. Rietbroek (2009). Decorrelated GRACE time-variable gravity solutions by GFZ, and their validation using a hy- drological model, Journal of geodesy, 83, 903-913; http://dx.doi.org/10.1007/s00190-009-0308-3. Lan, S.-C., T.-T. Yu, H. Cheinway and K. Ricky (2011). An Analysis of Mechanical Constraints when Using Superconducting Gravimeters for Far-Field Pre-Seis- mic Anomaly Detection, Terrestrial, Atmospheric and Oceanic Sciences, 22, 271-282. Matsuo, K., and K. Heki (2011). Coseismic gravity changes of the 2011 Tohoku-Oki earthquake from satellite gravimetry, Geophysical Research Letter, 38, L00G12; doi:10.1029/2011GL049018. Menezes Jr, J.M.P., and G.A. Barreto (2006). A New Look at Nonlinear Time Series Prediction with NARX Recurrent Neural Network, In: SBRN 06: Proceedings of the Ninth Brazilian Symposium on Neural Networks, 28-33. Ogawa, R., and K. Heki (2007). Slow postseismic re- covery of geoid depression formed by the 2004 Sumatra-Andaman Earthquake by mantle water dif- fusion, Geophysical Research Letters, 34, L06313; doi:10.1029/2007gl029340. Ozawa, S., T. Nishimura, H. Suito, T. Kobayashi, M. To- bita and T. Imakiire (2011). Coseismic and postseis- mic slip of the 2011 magnitude-9 Tohoku-Oki earthquake, Nature, 475, 373-376; doi:10.1038/natu re10227. Panet, I., V. Mikhailov, M. Diament, F. Pollitz, G. King, O. de Viron, M. Holschneider, R. Biancale and J.-M. Lemoine (2007). Coseismic and post-seismic signa- tures of the Sumatra 2004 December and 2005 March earthquakes in GRACE satellite gravity, Geophysical Journal International, 171, 177-190; doi:10.1111/j.13 65-246X.2007.03525.x. Pollock, D.S.G, R.C. Green and T. Nguyen (1999). Handbook of time series analysis, signal processing, and dynamics: Access Online via Elsevier. Ray, R.D., and S.B. Luthcke (2006). Tide model errors and GRACE gravimetry: towards a more realistic as- sessment, Geophysical Journal International, 167, 1055-1059. Seiji, T., and T. Nakamura (2013). Sea surface gravity changes observed prior to March 11, 2011 Tohoku earthquake, Physics of the Earth and Planetary In- teriors, 221, 60-65. Shen, W.B., D.J. Wang and C.W. Hwang (2011). Anom- alous signals prior to Wenchuan earthquake detected by superconducting gravimeter and broadband seis- mometers records, Journal of Earth Science, 22, 640-651. Shrivastava, M.N., and C.D. Reddy (2013). The Mw 8.6 Indian Ocean earthquake on 11 April 2012: Coseis- mic displacement, Coulomb stress change and af- tershocks pattern, Journal of the Geological Society of India, 81, 813-820. Swenson, S., and J. Wahr (2006). Post-processing re- moval of correlated errors in GRACE data, Geo- physical Research Letters, 33; doi:10.1029/2005gl02 5285. Tapley, B.D., S. Bettadpur, J. Ries, P.F. Thompson and M. Watkins (2004). GRACE measurements of mass variability in the Earth system, Science, 305, 503-505. Wang, L., C.K. Shum and C. Jekeli (2012a). Gravita- GRAVITY CHANGES BEFORE POWERFUL EARTHQUAKES tional gradient changes following the 2004 Decem- ber 26 Sumatra-Andaman Earthquake inferred from GRACE, Geophysical Journal International, 191 (3), 1109-1118; doi:10.1111/j.1365-246X.2012.05674.x. Wang, L., C.K. Shum, F.J. Simons, A. Tassara, K. Erkan, C. Jekeli, A. Braun, C. Kao, H. Lee and D-N. Yuan (2012b). Coseismic slip of the 2010 Mw 8.8 Great Maule, Chile, earthquake quantified by the inver- sion of GRACE observations, Earth and Planetary Science Letters, 335-336, 167-179; doi:10.1016/j.epsl. 2012.04.044. Wang, L., C.K. Shum, F.J. Simons, B. Tapley and C. Dai (2012c). Coseismic and postseismic deformation of the 2011 Tohoku-Oki earthquake constrained by GRACE gravimetry, Geophysical Research Letter., 39, L07301; doi:10.1029/2012GL051104. Yang, C.-C, S.O. Prasher and G.R. Mehuys (1997). An artificial neural network to estimate soil tempera- ture, Canadian Journal of Soil Science, 77, 421-429. *Corresponding author: Mehdi Akhoondzadeh, University of Tehran, University college of Engineering, Surveying and Geomatics Engineering Department, Remote Sensing Division, Tehran, Iran; email: makhonz@ut.ac.ir. © 2014 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. SHAHRISVAND ET AL. 12 << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles false /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.3 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.1000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams true /MaxSubsetPct 100 /Optimize false /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile (None) /AlwaysEmbed [ true /AndaleMono /Apple-Chancery /Arial-Black /Arial-BoldItalicMT /Arial-BoldMT /Arial-ItalicMT /ArialMT /CapitalsRegular /Charcoal /Chicago /ComicSansMS /ComicSansMS-Bold /Courier /Courier-Bold /CourierNewPS-BoldItalicMT /CourierNewPS-BoldMT /CourierNewPS-ItalicMT /CourierNewPSMT /GadgetRegular /Geneva /Georgia /Georgia-Bold /Georgia-BoldItalic /Georgia-Italic /Helvetica /Helvetica-Bold /HelveticaInserat-Roman /HoeflerText-Black /HoeflerText-BlackItalic /HoeflerText-Italic /HoeflerText-Ornaments /HoeflerText-Regular /Impact /Monaco /NewYork /Palatino-Bold /Palatino-BoldItalic /Palatino-Italic /Palatino-Roman /SandRegular /Skia-Regular /Symbol /TechnoRegular /TextileRegular /Times-Bold /Times-BoldItalic /Times-Italic /Times-Roman /TimesNewRomanPS-BoldItalicMT /TimesNewRomanPS-BoldMT /TimesNewRomanPS-ItalicMT /TimesNewRomanPSMT /Trebuchet-BoldItalic /TrebuchetMS /TrebuchetMS-Bold /TrebuchetMS-Italic /Verdana /Verdana-Bold /Verdana-BoldItalic /Verdana-Italic /Webdings ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 150 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.10000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.10000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.08250 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName (http://www.color.org) /PDFXTrapped /Unknown /CreateJDFFile false /SyntheticBoldness 1.000000 /Description << /ENU (Use these settings to create PDF documents with higher image resolution for high quality pre-press printing. The PDF documents can be opened with Acrobat and Reader 5.0 and later. These settings require font embedding.) /JPN /FRA /DEU /PTB /DAN /NLD /ESP /SUO /NOR /SVE /KOR /CHS /CHT /ITA >> >> setdistillerparams << /HWResolution [2400 2400] /PageSize [595.000 842.000] >> setpagedevice