Determination of possible responses of Radon-222, magnetic effects, and total electron content to earthquakes on the North Anatolian Fault Zone, Turkiye: an ARIMA and Monte Carlo Simulation

Around the world, earthquake forecasting studies have become very important nowadays due to the increase in number of fatal earthquakes annually. This paper proposes to achieve a possible relationship between soil radon gas concentration and atmospheric total electron content (TEC) during earthquakes taking into account magnetic effects on the North Anatolian Fault Zone (NAFZ) in Turkiye. The ARIMA and Monte Carlo simulation (MCS) are employed for determining radon gas concentrations by taking into account magnetic effects as an innovative approach. In the study area, relatively small and medium-scale earthquakes have taken place during the observation period. As a result of the investigations, the relationships between each of the parameters and earthquakes are determined; hence, a good relationship is obtained between Rn gas anomaly and micro-seismic activity. In the ionosphere, geomagnetic activity has a primary impact and long duration on TEC distribution, but due to micro-seismic events it has rather small influence. The proposed ARIMA and MCS simulations to detect changes in soil Rn gas concentrations have significant results for detecting micro-seismic activity anomalies.


Introduction
Earthquake forecasting studies have become very popular recently, especially in the last five years, huge earthquakes have started to occur worldwide with magnitudes more than 7. Scientists hope to be able to predict earthquakes to keep society away from death and injury. The relationship of 222 Rn with earthquakes has been studied individually by scientists (Kamislioglu and Külahcı 2014;Külahcı and Şen 2014;Külahcı and Çiçek 2015;Kamislioğlu and Külahcı 2016;Külahcı 2020) but research on the relation between earthquakes, 222 Rn, and atmospheric TEC is an extensive and very complicated research area. Requirement to combine joint experience in different fields is necessary, such as seismology, atmosphere with ionosphere physics, and geomagnetic activity (Le et al. 2013). In the earth's crust, most of the rocks and grains contain uranium element and since radon gas is a part of the uranium ( 238 U) decay chain, they are continuously produced in the soil. Most of the 222 Rn emissions can be detected in the earth's crust, and they can move large distance in the earth and atmosphere (Nazaroff and Nero 1988;Scholz 2019). It is well-known that one of the important and most often used precursory for predicting earthquakes is 222 Rn gas detection (Rikitake 1968;Ghosh et al. 2009;Külahcı and Çiçek 2015;Muhammad et al. 2020;Papastefanou 2007). First research about pre-earthquake 222 Rn anomaly was observed in 1966 for the Tashkent City of Uzbekistan earthquake (M = 5.2) capital, where the study area was monitored over 6 years before the earthquake, and these observations indicated increasing 222 Rn concentration until the earthquake occurrence time. In China (1976), before some severe earthquakes occurrences, unusual changes were observed in radon gas Ulomov and Mavashev 1967;Kuo et al. 2011;Friedmann 2012). Rn concentrations were determined in fault lines higher than normal soils (Singh et al. 1999).
Ionospheric anomalous due to seismic events is another precursory acceptance by scientists. Research on this subject is rapidly increasing in the last decade. Prior to high magnitude earthquake occurrences, mechanical energy is accumulated in the plates, which causes to positive and negative ions generations through compression of the rocks in the ground, and acoustic gravity waves (AGW) through earth crust motion and hotbed gases emission of into the atmosphere (Miyaki et al. 2002;Liperovsky et al. 2008;Md Yusoff and Hwee San 2016). These processes from the ground to air atmosphere can cause for temperature rise, pressure rise in the troposphere. TEC perturbation in ionosphere occurs with variation in an earthquake precursory Tariq et al. 2019;Xia et al. 2011). One of the significant in investigation about ionosphere is observable TEC parameter. TEC monitoring provides a direct link to connect to atmospheric activity and causes for understanding the complex phenomena that occurs in the ionosphere and plasmasphere. TEC is defined as an electron density profile amount between two points in a cylinder tube of (1 m 2 ) cross-sectional area, which corresponds to the total number of electrons. The TEC unit is defined as (TECU) and 1 TECU = 10 16 electron/m 2 (Nayir et al. 2007;Simha et al. 2014;Xia et al. 2011). Determination of the reasonable effects on the ionospheric parameters variation is important. The influences of geomagnetic activity are very essential from above and seismic activity below the ionosphere (Bhattarai et al. 2016). Solar winds carry the sun's energy and transfer it into the earth's magnetosphere and as a result geomagnetic storms generation occurs (Spruit and Roberts 1983). During geomagnetic field disturbance, a huge amount of energy enters into the earth magnetosphere and ionosphere, which causes to ionospheric parameter changes, such as composition, temperature, and height changes (Bhattarai et al. 2016). The first seismo-ionospheric copulation study was recorded in 1964 about the Alaska earthquake (Davies and Baker 1965;Leonard and Barnes Jr 1965). After about 20 years later satellite data are used for seismo-ionospheric investigation in the early 1980s (Gokhberg 1983). On Jan 17, 1994, for the first time (GPS) technology was used during the earthquake in M w = 6.7 Northridge (Calais and Minster 1995).
The ARIMA method or the Box-Jenkins method (1960-1970 was found after the two statisticians George Box and Gwilym Jenkins, who applied a model to find the best fit of a time series to past values of a time series data. ARIMA method is derived from Autoregressive Integrated Moving Average (Box et al. 2011;Ho and Ting 2015). Based on the past values of a data time series information, ARIMA method simulation can be made about future values. Simulation models produce an effective estimating to find unknown quantities so that this technique is valuable to estimate and forecast soil radon concentration . The MCS is a probabilistic method; it is useful for predicting future values of a variable. The basic idea behind the MCS method is to run simulations to get the probability distribution of all possible outputs with the best and worst scenarios (Farrance and Frenkel 2014). The first work in the literature on the Monte Carlo was published by applying Monte Carlo techniques on the measurement of the spectra Bremsstrahlung photon beam of radiotherapy by using the Sodium Iodide (NaI) detector from Bentley et al (1967) (Bentley et al. 1967). The use of the Autoregressive Integrated Moving Average (ARIMA) method together with Monte Carlo Simulations (MCSs) technique is a rather new technique and it is the first time used for prediction in this field with these parameters.
The objective of this study is to determine a possible relationship between 222 Rn gas concentration and TEC variations with earthquakes by considering magnetic effects, and studying all variables jointly. Another innovation is that ARIMA and Monte Carlo simulations help to identify a model variation and 222 Rn concentration prediction on the North Anatolian Fault Zone (NAFZ), Turkiye.

Data sources and method of analysis
The North Anatolian Fault Zone (NAFZ) is one of the most famous dextral strike-slip faults, especially in Turkiye, and generally, in the Eastern Mediterranean region because of its significant seismic activity (Akin et al. 2013). The NAFZ stretched (1200 km) from Saros Gulf in the Northern Aegean Sea to the Karliova triple junction in Eastern Anatolia and have a broad arc-shape. It forms a border between Anatolian and Eurasian plates almost parallel to the Black Sea coast maintaining a relatively regular distance of 100 km from the coast (Akin et al. 2016;Bayrak et al. 2011;Bozkurt 2001). The study area, Yolkonak station is shown in Fig. 1.
The study area, Yolkonak station is located 37 km north-east of the Tokat city center at 40.53932N-36.89443E on the eastern part of the seismically active NAFZ segment and the left embankment of the Kelkit River, the altitude from sea level is 355 m. NAFZ has high seismic activity for posing a high risk to areas that have a high population around the fault. It is very important for the understanding of large strike-slip fault systems (Bayrak et al. 2011). Study area elevation 3D graphs are shown in Fig. 2.

Earthquake measurement
The seismic activity in the Tokat-Yolkonak region during the observation period between 2007 and 2010 can be obtained from the Bogazici Kandilli Observatory (http:// www. koeri. boun. edu. tr/ scrip ts/ lasteq. asp). A total of 54 earthquakes in a circle with 150 km radius around the station were recorded and 29 earthquakes with a magnitude between M L (2.6-4.7) at prediction range of data. The largest one is M L = 4.7, which occurred inside the observation area.

Radon measurement
Rn gas concentration data for period 2007-2010 in the study area is donated by the AFAD-Ministry of Interior Disaster and Emergency Management Presidency (https:// en. afad. gov.  tr/). Data transferred to the center from all online stations. In the present study, the online 222 Rn measurement system consists of the detector (Alpha Meter 611, manufactured by Alpha Nuclear Co., Canada), power source, data logger unit, and global system mobile communication (GSM) modem (Thomas et al. 1992). Soil 222 Rn gas was deliberated continuously every 15 min interval by using alpha meter detectors. The radon gas emanation variation and earthquakes for the study period are showed in Fig. 3. There is one Rn station established on the fault line in Tokat Yolkonak. Continuous Rn concentration measurements were taken from this station.
The nearest and available TEC stations for the period (01 Mar 2007 to 12 Feb 2010) from the study area are three stations (Ank, Tubi, and Ista). The TEC daily data are recorded at every 2.5 min for a single and multiple-station computation in the range of consecutive days. The upper and lower bounds are obtained through the statistical equation (μ ± 2σ), where μ is the mean and σ is the standard deviation for 10 days before and 5 days after the event ).

Solar geomagnetic activity
Three important parameters are global geomagnetic activity (K p -index), disturbance time storms (D st ), solar flux (F10.7) are used in this study to explain the effects of the geomagnetic activity on the ionospheric TEC. The parameters are obtained from OMNI internet-based data through explorer system project developed by (NASA) (www. omniw eb. gsfc. nasa. gov). These parameters are global geomagnetic activity (K p -index), which defines the earth's magnetic field disturbance caused by the solar Fig. 3 Soil 222 Rn concentration levels and seismic activity in the study period time wind, and it values are between 10 and 90. When its value is below 40 then the condition is considered as quiet, K p -index value multiplied by 10 in this study (Menvielle and Berthelier 1991). Disturbance time storms (D st ) show the geomagnetic activity in one-hour time, when D st value between is in the range −20 and 20 nT, and this condition is considered as quiet geomagnetic (Cander and Mihajlovic 1998). The solar flux (F10.7) has unit sfu (solar flux unit = 10 -22 W m −2 Hz −1 ), and it defines the density of a microwave flux having 10.7 cm wavelength when the F10.7 value is below (70 sfu), and the condition is considered as a quiet (Guo et al. 2015).

Results and discussion
Rn gas concentration is monitored continuously for about three years from 2007 to 2010 as a capable precursor of seismic activity at NAFZ. All data are used in this study equal to (103,235) radon data for every 15 min interval. Prior to the ARIMA and MCS methods application, the data are divided into train and test parts as 80% and 20%, respectively. The test prediction part occupies 20% in the tail of data, which is used to evaluate accuracy of the model in the prediction.
To determine the best model for data time series, Box-Jenkins recommended assurance of stationarity in time series. For this reason, more than one test it taken into consideration to assure the time series stationarity (Kwiatkowski et al. 1992), such as the KPSS test, augmented Dickey-Fuller (ADF) test, skewness test, (one-sample and paired-sample) t-test, and Leybourne-McCabe test (Cheung and Lai 1995;Hobijn et al. 2004;Leybourne and McCabe 1994). After checking the results of (KPSS test, Leybourne-McCabe, skewness test, and ttest), it is observed that all support that the data is non-stationary despite ADF test. Hence, it is necessary to take the differentiation until the series becomes stationary. The same is valid for above-mentioned tests for radon data differentiation in a second time. After proving the time series data as stationary, the next step is the ARIMA model usage to find the best model fit to the Rn time series. Model identification needs to determine values of p, and q from lags crossing the confidence significant lags bounds of the PACF (Partial Autocorrelation Function) and ACF (Autocorrelation Function) plot, respectively. The ARIMA (8,1,13) model is the most suitable one for the prediction of the future value of 222 Rn concentrations and the model has a satisfactory estimate, because simulation values are in good agreement with the real radon values, and the smallest standard deviation belongs to ARIMA (8,1,13) model, which is defined in Eq. 1 in the below.
where coefficients of autoregressive is ϕ, and coefficients of moving average is θ, and c is a constant. Different ways are examined to test the validity of the recommended ARIMA model. For this reason, the ARIMA mode fit, residual plot, histogram plot, and Q-Q plot are obtained. To evaluate the power of the model in Fig. 4a, b, it shows that the model has good fitting with radon data and the residual distribution plot fluctuates around the mean with some anomalies as in Fig. 4c, d one can see the histogram (number of repetitions of similar values for residual values regularity changes. The Q-Q graph produces a scatter plot between standard normal quantities with quantities of the input sample. In the sample, the data mostly fit the model with a little diffraction in both tails, which conforms that data fit with same anomalies. In Fig. 5, the model data prediction average is shown. After finding the best model fit to the time series, the ARIMA and the MCS methods can be employed for simulations to get probability distribution of future radon concentration values and identification of all possible future probability values. MCS is also used for time series data prediction testing part, which consists of 20% of data in the tail of the time series data as in Fig. 6. It can be seen that the MCS prediction result has a good fit with real radon data in (Panel A) and four different possible outputs for MCS prediction are shown

Analysis of the relation between radon gas and TEC with earthquakes
The prediction part of radon data equal, which includes 20% of all data and MCS results to soil radon data are used in this study based on the earthquake concentration. It is divided into four different groups, where each group includes 15 days for discussion about the relation and influence between soil radon gas and total electron content (TEC) with seismic activity by considering the geomagnetic field activity effect.

First group
The first group starts from 12 to 27 July 2009 with six micro-earthquake magnitude occurrences as 2.6-3.1. The TEC and radon concentration for the same period are shown in Fig. 7.
In the first group, a positive anomaly can be seen in the soil radon concentration and increase over prediction MCS redline value by 40 Bq/m 3 , which is related to the first two earthquakes magnitudes, 2.7 and 3 in two days and 18 h before the event, respectively. The TEC variations for both events are not suitable to count as an anomaly in the range between the upper and lower boundaries (mean ± two standard deviation) in all stations. Three earthquakes with 3.1, 2.8, and 2.6 magnitudes occurred on July 20, 2009, during 18 h before earthquakes increasing (50 Bq/m 3 ) of 222 Rn and over the MCS prediction there is an anomaly is related to earthquakes and a short negative anomaly in the TEC can be seen in all stations especially in the nearest one to Ankara TEC station, where one day before is related to events, because of geomagnetic activities are quite well in this period. Last earthquake in this group with M L = 2.7 during three days before the earthquake, there is an anomaly in the radon sharp increase over the MCS prediction due to the earthquake. In contrast, total electron content uprising on 22 July is not reliable as earthquake anomaly due to variations in the (K p and D st ) magnitudes as shown in Fig. 8.
From Fig. 8, on July 22, 2009, the disturbance in the geomagnetic can be seen as a result of disturbance in the earth's magnetic field, which is caused by the solar wind (K p ) equal to 58, which is more than 40 and it counted as a quiet condition disturbance time storms (D st ) is equal normally to −80 nT, and it should be in the range from 20 to −20 nT with the source for TEC anomaly in all stations. In the same graph, one can see the density of a solar microwave flux having as 10.7 cm wavelength value, which is below 70 sfu; hence, the condition is considered as quiet.

Second group
The second group evaluates four earthquakes with other parameters starting from August 16, 2009, to October 1, 2009. The results for TEC variations for three nearest international stations are shown in Fig. 9 with radon gas concentration and seismic activities.
In Fig. 9 graph shows that there is no clear and significant pre-seismic or co-seismic anomaly in the 222 Rn gas concentration. This may be due to the large distance epicenter with radon station and low magnitude, but after the earthquake (M L = 2.7) there is a post-seismic increase in soil 222 Rn gas concentration. The small uprising can be seen in the TEC peaks on 19, 20, and 22 of August middays due to earthquake M L = 3.1 in 26 August., The increase is clearly that more in the (Ankara city) station compared to other stations. The reason for this small variation in the TEC is the earthquake with low magnitude, pressure production in the ground, which causes to stress generation on rocks and heavily charged ion clusters, distribution in the electric field between ground and atmosphere could be lead to this TEC variation in the ionosphere a few days before the earthquake. Another anomaly in the TEC on 30 August can be seen, but this variation is related to geomagnetic variations as shown in Fig. 10.
In Fig. 10, all three parameters at (19, 20, and 22 of August) have normal values. It approves that the TEC variation source in earthquake M L = 3.1 on 26 August causes to geomagnetic activities on 30 to 31 August with an increase in the (K p index) over normal value (40) increasing toward it to (58). The variation in disturbance storm time (D st ) is under −20 nT), with these anomalies in the geomagnetic activity, they are the main factors for ionospheric TEC anomalies on 30 August over the upper bound, which is equal to mean plus two standard deviations (μ + 2σ) as can be seen in all three stations in the midday.

Third group
The third group includes observation radon and TEC data with earthquakes from August 30, 2009, to October 16, 2009, which leads to discussion on the relation between parameters as shown in Fig. 11. This group consists of five micro-earthquakes magnitudes between 2.7 and 3.9. According to observations, there is a positive Rn anomaly on 4, 8, and 10 October with approximate increase (40 Bq m −3 ) over the prediction magnitude estimation by MCS. This variation is counted as a pre-seismic anomaly with 3.9 and 2.8 magnitude earthquakes that occurred on 10 October. For the second series of micro-seismic events that occurred on 12 October, there is a co-seismic increase in the Rn gas emanated in the soil due to produced pressure on soil pores and the number of micro-cracks increases is helpful to allow radon gas migration. TEC response for seismic events and soil Rn variations under quiet geomagnetic at the period of increasing Rn gas on 4 October, the TEC increased one day after variation Rn gas. There is another anomaly in the TEC Ankara, Tubi, and Istanbul stations on 11 October exactly one day before the second series of seismic events. The effect of this positive anomaly in the TEC in all stations is a seismic activity because the geomagnetic activity for the same time is approximately quiet, and the variations of geomagnetic activity parameters (K p , D st , and F10.7 flux) are shown in Fig. 12.
The geomagnetic activity or ground base due to seismic activity are considered during observations for clarifying the TEC anomalies sources. Variation in the TEC on 4 October can see a small negative anomaly under the lower bound as equal to mean minus two standard deviations (μ-2σ) in all stations. It is due to increasing soil radon gas at the same time because all three (K p , D st , F10.7) parameters and are approximately quiet. For the second TEC anomaly on 11 October, the seismic events are the main source for these positive anomalies over the upper bound in all stations, because during observations, K p index, F10.7 flux is in the quiet range with a small increase in the D st .

Fourth group
The fourth group includes the observation data from Jan 22, 2010 to Feb 7, 2010, and the most powerful seismic events during this study belongs to this group and the effects on the soil radon gas concentration and TEC are presented in Fig. 13.
The Alpha meter 611 detector detects an increase 70 Bq m −3 in a short time in the Rn gas concentration 18 h before the earthquake magnitude M L = 3.4 on 24 January and sharply decreases after the earthquake to the previous level of concentration. Two days later radon increases 15 h before the M L = 3.3 magnitude earthquake on 27 January. Both earthquakes create an effect on soil Rn emanation, and TEC positively responds with record anomaly above the upper bound (μ + 2σ), at the same time for Rn variations, but cannot be counted as earthquake precursory due to high solar F10.7 flux. The most powerful earthquake among all earthquakes is observed in this study magnitude equal to 4.7, which has some aftershocks with magnitudes 3.1, 2.9, 2.8, and 2.7, respectively. The average 222 Rn gas increases due to these seismic events at approximately 50 Bq m −3 in three days, which is not satisfactory for the measure as an anomaly, but, the epicenter distance to Yolkonak radon station is the main reason for this kind response of radon concentration to seismic events.
The abnormal variation in the 222 Rn concentration on 1 February nearly for about five days could be due to earthquake blasts for false recording, high rainfall and snow for the same period, and water was accumulated on the detector as a consequence of flooding at the measurement site; surface sensor severely affected by humidity for short-circuits, by continuation of irregular and unpredictable records up 34,000 Bq/m 3 just for 15 min. This kind of records for Alpha Meter detector could only be a possibility in the case of shortcircuit occurrence inside the Alpha meter detector. This effect is confirmed in the literature by this article (Thomas et al. 1992).
The variation of TEC in the all-stations path through the upper and lower boundaries on 3 to 4 February is a response for negative variation in the disturbance time storm (D st ) as in Fig. 14.
Results in Fig. 14. show that the abnormal increase in the solar flux F10.7 to the range 77 sfu on 27 January is the primary source for variation in the TEC in all stations, nonetheless other two parameters, namely, D st and K p index record the normal quite magnitude at the same time.
It can be said that the main source for both positive and negative anomalies in TEC in all stations is due to the disturbance time storms D st variation, which is equal to −24 nT, and other parameters have the quite normal magnitude in the period.

Relation between lithosphere and ionosphere
Relation between soil radon concentration, TEC and seismic activity for all stations are shown in Fig. 15. This study discusses all three parameters together for the first time in the literature.
The peaks are time lag between the 222 Rn gas, and the TEC with seismic events. It is important that the anomaly peaks are bigger and wider in stations far from the epicenter, which means that after emanation radon gas by earthquake produced pressure, and it takes some time to affect the ionosphere and TEC. The reason is that the wider peaks are functioned for exchange in the energy and effect between the lithosphere and ionosphere through the seismic events in the preparation period before the main-shock time (Ryu et al. 2014;Shah et al. 2019).

Conclusions
During the observation period in the study area, only approximately small and mediumsize earthquakes have occurred. Results show that a possible good correlation exists between soil 222 Rn gas anomalies and micro-seismic activities observations suggest that the soil 222 Rn data variations in the study area in NAFZ can be relied on as an earthquake precursor, and they can be useful for predicting future earthquakes.
ARIMA and MCS usage together with detect radon gas anomalies is a rather new method in this field. The proposed ARIMA model and Monte Carlo forecasts of the soil radon gas data related to micro-earthquake has a significant result and effective finding and estimating other parameters. They are also suitable in forecasting soil radon gas related to macro-earthquake prediction, and other radioactive substances in the environment.
The soil radon gas anomaly amplitude and peak sizes depend on the magnitude of impending earthquakes and the distance to the epicenter. Earthquake magnitude has a positive influence on the radon concentration but the earthquake epicenter depth has a negative impact.
The main source for ionospheric TEC distributions related to the magnetic effects and the distribution can be seen in all the TEC stations, while the distributions due to the micro-seismic events have a lesser extent with smaller impacts.
A direct relationship between radon and ionosphere in seismic periods would be mathematically complicated to establish. To achieve this, a number of multidisciplinary results needed to be presented, which would give insight on physical mechanisms involved in the process, as well as other atmospheric processes. Hence, the relationship is better explained when perceived as indirect. There are several factors accompanying radon injection which can affect the ionosphere. Among them are acoustic-gravity waves generated by heating or movement of the lower atmosphere, as well as perturbation of electric current in the global circuit caused by charged aerosols injection. These factors can occur simultaneously alongside radon emanation (Sorokin et al. 2020). The study hypothesis is centered around the notion that; "Since the TEC-Earthquake and Rn-Earthquake pairs show distinct physical changes separately, both changes can together result to meaningful interpretations." Of course, there are presently many proposed mechanisms for the Lithosphere Atmosphere Ionospheric (LAIC) phenomenon. Here, we adopt some possible explanations of seismically induced events, which might relate Rn and TEC in seismic periods. As stress accumulates in rocks due to ground movements, electron migration from unstressed part of the rocks (especially the surface) induces an electric field along the stress gradient (Liu et al. 2000;Freund et al. 2009;Petraki et al. 2015). Geochemical gasses such as carbon dioxide, radon, and methane, etc. get exhaled. This exhalation is also accompanied with charged aerosols (Sorokin et al. 2007). Ionization of the lower atmospheric molecules occur in the process. Moreover, Rn progeny attach to these aerosols, continuously ionizing them. This evolution of ion-molecular interactions leads to cluster formation, and the region gets electrified. According to Namgaladze et al. (2018), the charged aerosols play a great role in the generation of tropospheric electric field, which we suppose is enhanced by ionization caused from Rn and its progeny. Governed by atmospheric convective forces, gravitational sedimentation, turbulent diffusion and other forces over the seismic region, the vertical motion of charged aerosols and particles having different masses and polarities generate an EMF, which induces an upward electric current density. This induced current is termed "extraneous current," it adds to the conductive current (fair weather 1 3 current) of the global electric circuit, which might result in causing noticeable ionospheric disturbance. The generation of this extraneous current strongly depends on the nature of atmospheric processes involved. More explanation on this can be found in (Liu et al. 2016;Namgaladze et al. 2018;Sorokin et al. 2020).