Estimation of date and magnitude of four major earthquakes using integration of precursors obtained from remote sensing data

A single precursor is not usually an accurate, precise and adequate measure to predict earthquake parameters. Therefore, it is more appropriate to exploit parameters extracted from several other single precursors, so that their simultaneous combinations may reduce the uncertainty of the prediction. In this study, remote sensing observations in different modalities acquired from several days before impending earthquakes have been investigated to extract earthquake parameters. They are observations in electron and ion density, electron temperature, Total Electron Content (TEC), Land 10 Surface Temperature (LST), Sea Surface Temperature (SST), Aerosol Optical Depth (AOD), and Surface Latent Heat Flux (SLHF). Regarding the ionospheric precursors, the geomagnetic indices Dst, Kp, Ap and F10.7 were used to detect preearthquake disturbances from frequent anomalies associated with geomagnetic activity. In this study, three methods of median, support vector regression (SVR) and random forest (RF) have been used to detect anomalies. When anomalies associated with impending earthquakes are detected, the number of prior days associated with the earthquake is estimated 15 based on the type of precursor. Then, by estimation of the amount of anomaly deviation from the normal state, the magnitude of the impending earthquake is estimated. The final earthquake parameters (such as date and magnitude) can be obtained by integrating the earthquake parameters extracted from different earthquake precursors using mean square error (MSE) method.


Introduction
When an earthquake is happening, energy transmission is generated due the destructive effects of the earthquakes to the environment. The occurrence of these changes before and/or after the earthquake may have various physical and chemical effects on lithosphere, atmosphere and ionosphere making the earthquake more accurately predictable. The abnormal variations in lithospheric, atmospheric and ionospheric parameters are taken as "earthquake precursors". They serve as 25 alarms for impending earthquakes. Many studies have been carried out on earthquake predictions using precursors in lithosphere, atmosphere and ionosphere. The problem arise when some of these major abnormalities may not appear on occasion an earthquake. There are several studies based on the observation of the seismic Lithosphere Atmosphere Ionosphere Coupling (LAIC) anomalies which confirm that the anomalies begin several days before the earthquake and remain a few days after. None of the earthquake precursors can be used alone as an accurate and independent parameter for 30 estimating earthquake parameters without generating some level of uncertainty. Hence, it is necessary to integrate different types of earthquake predictors or earthquake precursors. By integrating a variety of earthquake parameters extracted from different precursors, a more accurate and suitable estimation of final earthquake parameters may be obtained. ‫ﺑﺪوح‬ ‫ﺑﺪوح‬ 3 ‫ﺑﺪوح‬ ‫ﺑﺪوح‬ ultraviolet. Its spatial resolution is 0.25° Lat×0.25° Lon. The OMI products are available via: 75 https://giovanni.gsfc.nasa.gov/giovanni/.

AVHRR Data
Two products of AVHRR (Advanced Very High Resolution Radiometer)including Sea Surface Temperature (SST) and Surface Latent Heat Flux (SLHF) data have been used in this study. Sea surface temperature (SST) anomaly can be related to near coastal seismic activity (Ouzounov and Freund, 2004), but conditions and currents can strongly affect SST. Due to a 80 large thermal inertia of the seawater, its temperature changes more slowly; therefore, in the case of the SST anomalies, some mechanisms of LST anomalies are not applicable (Jiao et al. 2018). The NOAA 0.25° daily Optimum Interpolation Sea Surface Temperature (OISST) is an analysis constructed by combining observations from different platforms (satellites, ships, buoys) on a regular global grid. A spatially complete SST map is produced by interpolating to fill in gaps. The SLHF products are available via:https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.highres.html. 85 SLHF is (Wm -2 ) the heat flux absorbed or released by the phase transition (i.e., condensation, evaporation, and melting) of water from the Earth's surface to the atmosphere (Jiao et al. 2018). SLHF is one of the important components of Earth's surface energy budget, which is mainly affected by the atmospheric relative humidity, wind speed, surface temperature, and season (Jiao et al. 2018). Due to the underground fluid movement and the interaction among the underground, surface, and atmosphere, the SLHF anomaly that occurs prior to earthquakes is considered (Alvan et al. 2013). The SLHF products are 90 available via:https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html.

Geomagnetic Indices
The ionospheric parameters measured by satellite are mainly influenced by the geomagnetic storms and geomagnetic field disturbances, particularly in the equatorial and polar regions. However, in the case of an impending Earthquake, it may be affected further more in the form of anomaly. Such anomalies should also be removed. In order to distinguish anomalies 95 caused by seismic activity from anomalies created by geomagnetic and solar activities, the geomagnetic and solar indices i.e. Dst, Kp, Ap, and F10.7acquired from Space Physics Data Facility (SPDF)have been utilized in this study. In conditions where the quiet solar geomagnetic is established (i.e. Kp< 2.5, -20 nT<Dst< 20 nT and F10.7 < 120 SFU) the situation will be regarded as normal. However, if unusual changes are seen in the ionospheric data, it can be considered due to seismic activity.The geomagnetic and solar indices are available through NOAA (https://www.ngdc.noaa.gov/stp/stp.html). 100

Anomaly detection method
Anomaly detection is the process of identifying unexpected items or events in data sets, which differ from the norm. Anomaly detection has two basic assumptions: 1) anomalies only occur very rarely in the data, and 2) their features differ from the normal instances significantly. In order to detect anomalies in each set of time series data from remotely 105 sensed precursors, identification of the normal behaviour of the phenomenon is necessary.

Normal behaviour modelling
In order to model time series behaviour, two common machine learning methods namely the Support Vector Regression (SVR) and Random Forest (RF) have been used in this study.
In SVR-based regression solutions, the input vectors are mapped into a higher-dimensional feature space, then by employing 110 a linear regression in the feature space, the input vectors are separated as far apart as possible. In this study, which is operating in a large space, a kernel function is used (Khosravi, Jouybari-Moghaddam, and Sarajian, 2017

‫ﺑﺪوح‬ ‫ﺑﺪوح‬
Random forest algorithm which is a non-parametric machine learning ensemble method, is an extended version of CART (Classification and Regression Trees) model proposed by Breiman (2001). It is based on the information combination method in which a large number of decision trees are generated and then the results of all the trees are combined for 115 prediction (Cutler et al. 2007). The RF is the multitude of regression trees which performs based on the variance of the data (Liaw and Wiener, 2002).

Anomaly detection
In order to detect anomalies, a reasonable range for regular variations in the time series data should be specified. Signals with normal behaviour fluctuate inside both the upper and lower bounds. Signals outside the boundaries will be detected as 120 anomalies. In this study, the median and inter-quartile method is used to determine the upper and lower bounds (Liu et al. 2004).
The upper and lower bounds are determined using the following Eq. (1): where x low , x high , m and iqr are respectively the lower bound, upper bound, median value and inter-quartile range of x. If the signal lies within the range of the upper and lower bounds, the behaviour of the signal is considered as normal. To determine 125 the intense of abnormal behaviour of the signal, the parameter D x as deviation of x is calculated by Eq. (2): If the absolute value of the parameter D x is greater than k (i.e. |D x |>k), the behaviour of the parameter is regarded as anomalous. Also, the percentage of parameter deviation from the natural state can be calculated using the Eq. (3) (Saradjian and Akhoondzadeh 2011):

Preliminary parameters estimation 130
The Earthquake parameters are preliminarily estimated for each of individual precursors. The D x value obtained from the previous step is relatively suitable parameter for calculating Earthquake magnitude. Saradjian and Akhoondzadeh (2011) showed the relationship between this parameter and the Earthquake magnitude that can be extracted from Also, according to the day when the anomaly is observed, an impending Earthquake's approximate date can be estimated. 135 Based on observations so far, as an average, a 15-day interval in ionospheric and atmospheric precursors, and as an average, a 16-day interval in thermal precursors from the anomaly observation till earthquake day is considered (Saradjian and Akhoondzadeh 2011).

Finalizing parameters estimation method
After the earthquake parameters (i.e. date and magnitude) are estimated through various precursors using the preliminarily 140 Earthquake parameters estimation, the final value of parameters of the earthquake can be estimated by combining their https://doi.org/10.5194/angeo-2021-41 Preprint.
where V and M are variance and median of the predicted values of the earthquake parameters for all precursors, respectively; and x is earthquake parameter value for any precursor. Any parameter that has a minimum MSE value is considered as the 145 final parameter of the earthquake.

Case Studies and results
Four major earthquakes with Magnitude Mw>6have been investigated in this study. These earthquakes occurred in Samoa Islands, Sichuan (China), Kermanshah and Bam (Iran). The characteristics of these earthquakes have been presented in Table   2.The ionospheric parameters obtained from the DEMETER and GPS satellites have been studied and analyzed over the 150 relevant periods of time for each earthquake for areas selected according to Dobrovolsky Formula R = 10 0.43M which relates the size of affected area to the magnitude of the earthquake (Dobrovolsky et al. 1979).The rest of the time series data for other precursors have been acquired for the relevant periods of time for areas of about 2.5×2.5 degrees in size around each epicentre. [

Kermanshah Earthquake
In the case study of Kermanshah Earthquake, all time series data were provided for the period of 21 September to 21 November 2017.The TEC anomalies associated with Kermanshah earthquake were observed on November 3 and 4 (Table   3). By observing these anomalies, it can be concluded that an earthquake with magnitude ranging from 7 to 8Mw between November 4 and 19, 2017 would have happened ( Fig. 1(d)). 160

[Table (3) near here]
By investigating the Terra day-time LST (°C) anomalies, the minimum value of -3.2°C on October19indicates an earthquake between December 20 and November 4 with Mw>8 would have happened ( Fig. 2(a)). The anomaly observed in the AOD data obtained from the OMI sensor, with a maximum of 370.96% on October 26, 2017, indicates an impending earthquake with Mw>8 between October 27 and November 10 ( Fig. 2(b)). Also, the estimated earthquake magnitude for the observed anomaly on October 31, 2017 related to the AOD precursor obtained from the Terra and Aqua sensors would be greater than 8Mw (Fig. 2 By investigating the changes in SLHF, a sharp increase of 346.46%, 117.38% and 288.4% have been observed on 3, 5 and 6 November, respectively ( Fig. 2(e)). Due to these anomalies, it can be predicted that an earthquake with magnitude more than 8 Mw will occur in the region.
The results obtained by SVR and Random Forest methods can be found in tables 4 and 5, respectively. 175

[Table (5) near here]
In the case study of Kermanshah, the earthquake parameters deduced based on median from the different precursors using the MSE method indicate that an earthquake would occur between November 4 and 19 with a magnitude more than 8 Mw https://doi.org/10.5194/angeo-2021-41 Preprint.  (Table 15). Also by using MSE method for the obtained results from both SVR and RF methods, the predicted magnitude of earthquake will be from 7 to 8 Mw between November 4 and 18, 2017 (Table 15).

Samoa Earthquake
For the Samoa earthquake, the time series data of precursors were provided for the period of 21 August to 21 October 185 2009.TEC anomaly fluctuations began on September 18. The estimated magnitude of the Earthquake for this anomaly would have been between 7 and 8Mwhappening between September 19 and October 3. Anomalies have also been observed on September 25 and 28, 2009 implying that an Earthquake with magnitude 7<Mw<8 would be expected to happen (Fig. 3(d)).

[Table (6) near here]
The data on the changes in total ion density, electron and ion density and electron temperature recorded by DEMETER IAP 190 and ISL sensors for Samoa Earthquake are given in Table 6. The data were recorded when the satellite's orbits were near the epicentre of the Earthquake (i.e. less than 1500 km). With an increase in total ion density observed around 10:30 local time on 25 September, 2009, it was expected that an Earthquake with magnitude between 7 and 8Mw would occur between 26 December and 10 October, 2009 ( Fig. 4(a)). Also, the changes in electron density showed unusual behaviour. These changes By investigating changes in SST, there were sharp increases of 135.84% and 65.32% on September 6 and 7, 2009, respectively ( Fig. 4(g)). The irregularity seen on September 6 indicates that an Earthquake with magnitude greater than Mw = 8.0 would have occurred between September 7 and 27.

4(h)). 210
Unusual AOD changes from the Terra sensor, with a maximum anomaly of 2.07 and 2.08 on August30 and September 21, 2009, respectively indicates that an Earthquake with magnitude ranging from 7 to 8 Mw would have happened ( Fig. 4(i)).
The characteristics of detected anomalies using SVR and Random Forest methods are shown in Tables 7 and 8, respectively.
[ In order to finalize each Earthquake's date and magnitude parameters using the parameters obtained from different 215 precursors, the MSE fusion method has been applied. In the case study of Samoa, according to the MSE method, it was predicted that an Earthquake with magnitude between 7 and 8 Mw would occur between September 25 and October 9, 2009 indicated by both Median and SVR results (Table 15). Based on Random Forest method, the parameters would be the same except for one day earlier starting date (i.e. September 24).

Sichuan Earthquake
In the case study of Sichuan Earthquake, all time series data were provided for the period of 21 March to 21 May 2008.Some 225 intense anomalies related to TEC data have been observed on April 21 (6 UTC), May 3 (6 and 18 UTC), and May 9 (10 and 14 UTC) (Fig. 5(d)). These strong anomalies indicate an impending Earthquake with magnitude between 7 and 8Mw. The characteristics of other detected anomalies are seen in Table 9.
[ The variations of the various parameters extracted from the DEMETER experimental data over the Sichuan region have been 230 presented in Table 9. A sudden and unusual change in ion temperature has been observed three days prior to the Earthquake around 22:30 local time ( Fig. 6(a)). This means that an Earthquake with magnitude ranging from 7 to 8 Mw would have occurred between May 10 and 24, 2008. An anomaly has also been observed on May 21, 2008 in electron density around 10:30 local time which implies an impending earthquake as strong as 7 to 8Mw between 22 April and 6 May ( Fig. 6(b)). An  By using the MSE method and the combination of the earthquake parameters obtained from different precursors, it was predicted that an earthquake 7<Mw<8 would have happened between 4 and 18 May 2008 (Table 15). Also, the MSE method for Random Forest method indicates that an Earthquake with the magnitude between 7 to 8 Mw would have been occurred between 4 and 18 May, 2008 and for SVR method with the magnitude ranging from 7 to 8Mw between April 22 and 18, 245 2008 (Table 15).
[ Table ( Earthquake with the magnitude between 6 and 7Mw would have happened (Fig. 7(d)). 255 [ Also, a sudden decrease in temperature was observed in the day-time LST data on December 12 and 13, which indicates that an Earthquake with magnitude ranging from 7 to 8 Mw would have been occurred. By investigating the night-time LST fluctuations, extreme increases in temperature on December 6 and 7, 2003 indicate an Earthquake with magnitude ranging from 6 to 7Mw would have been happened. Moreover, by observing the maximum value of 2.14 on December 6, it can be 260 concluded that an Earthquake with magnitude between 7 and 8 Mw would have happened between December 7 and 22.
Unusual change in day time LST data (Terra-MODIS sensor) with the minimum value (D x = -3.32) occurred on November 30. Due to this anomaly, an Earthquake with Mw>8 would have been occurred (Fig. 8(a) to (d)).

270
The anomalies obtained by SVR and Random Forest methods, respectively, are shown in tables 13 and 14.
[ By combining the predicted parameters obtained from different predictors using the MSE method, it is predicted that an Earthquake would occur between December 13 and 28, 2003, for median anomaly detection method. The combination of anomalies obtained from the SVR and Random Forest by using MSE method predicts an Earthquake between December 16 275 and 31, 2003 and between December 16 and 30, 2003, respectively. The magnitude of this Earthquake is estimated to be 7<Mw<8 for both median and random forest methods and to be 6<Mw<7 for the SVR method (Table 15).

Conclusions
Assuming that estimation of Earthquake parameters using each predictor individually is accompanied by some uncertainties, this study considered integrating the capabilities of different earthquake parameters extracted from some of the same earthquake predictors to better estimate earthquake parameters. To identify the anomalous states that may be associated with 285 impending earthquakes, variations of different earthquake precursors have been analysed for four earthquakes by using Median, SVR and Random Forest methods. For each precursors, the date and magnitude were estimated according to the earthquake signals. By integrating the earthquake parameters obtained from all precursors, the final earthquake parameters were estimated more accurately.
Since different precursors have been used to analyse the final earthquake parameters, therefore, for all earthquakes, the 290 estimated earthquake parameters for each earthquake are close to the actually recorded parameters. This can lead to accurate estimation of earthquake parameters with respect to the number and variety of earthquake precursors. Based on the results, it seems that methods such as SVR and Random Forest dealing with nonlinear and complex behaviours of time series are more sensitive than the Median method. Therefore, it can be accounted that these methods are suitable tools for detecting anomalies in nonlinear time series related to changes in seismic precursors. 295 Since various factors can cause unusual behaviour in different ionospheric, atmospheric or lithospheric parameters, more careful studies should be conducted to distinguish the anomalies caused by the daily changes from anomalies due to seismic activities.