Monitoring of Gas Emissions in Light of an OEF Application

: This study analyzes the possibility to use geophysical and geochemical parameters in an OEF (Operational Earthquake Forecasting) application correlated with short-term changes in seismicity rates using a magnitude–frequency relationship. Tectonic stress over the limits of rock elasticity generates earthquakes, but it is possible that the emission of gases increases as a result of the breaking process. The question is how reliable is the emission of radon-222 and Carbon Dioxide (CO 2 ), with effects on air ionization and aerosol concentration, in an OEF application? The ﬁrst step is to select the seismic area (in our study this is the Vrancea area characterized by deep earthquakes at the bend of the Carpathian Mountains), then determine the daily and seasonal evolution of the forecast parameters, their deviations from the normal level, the short-term changes in seismicity rates using a magnitude–frequency relationship and ﬁnally to correlate the data with recorded seismic events. The results of anomaly detection, effect evaluation and data analysis alert the beneﬁciaries specialized in emergency situations (Inspectorate for Emergency Situations, organizations involved in managing special events). Standard methods such as the standard deviation from the mean value, time gradient, cross correlation, and linear regression are customized for the geological speciﬁcity of the area under investigation. For detection we use the short-time-average through long-time-average trigger (STA/LTA) method on time-integral data and the daily–seasonal variation of parameters is correlated with atmospheric conditions to avoid false decisions. The probability and epistemic uncertainty of the gas emissions resulting from this study, in addition to other precursor factors such as air ionization, time between earthquakes, temperature in the borehole, telluric currents, and Gutenberg Richter “a-b” parameters, act as inputs into a logical decision tree, indicating the possibility of implementing an OEF application for the Vrancea area. This study is novel in its analysis of the Vrancea area and performs a seismic forecasting procedure in a new form compared to the known ones.


Introduction
An OEF (Operational Earthquake Forecasting) application is based on the real-time monitoring of geophysical and geochemical parameters as well as the long-term data analysis. In this context, the monitoring of gas emissions due to tectonic stress could be an important component in making decisions to reduce losses caused by a seismic event. S. Pulinets Sergey and D. Ouzounov point out in [1] the meaning of short-term earthquake forecasts and the role of gas interaction between the lithosphere and the atmosphere. These gas emanations result in air ionization and ultimately to seismically induced aerosols. This study is a more deterministic approach to the forecast with respect to the presently available purely probabilistic OEF models because it is more heavily based on the analysis

Monitoring Network
The multidisciplinary monitoring network of the NIEP (National Institute for Earth Physics) includes five stations measuring radon, CO 2 , CO, air ionization, atmospheric electrostatic field, earth and magnetic field, temperature in the borehole, meteorological station, solar radiation, telluric field in addition to seismic speed and acceleration sensors. The monitored area is characterized by faults (an example in Discussion and Conclusions paragraph) and in the southern part by muddy volcanoes. Figure 1 presents the seismicity of the area where deep earthquakes occur (marked in dark blue) and the name and position of the stations used in this study.  Table 1 details the monitoring structure. Each location has a unique identification code that allows access to data on the http://geobs.infp.ro/platform. The radon detector Radon-Scout PLUS is produced by SARAD. The measuring range is up to 10 MBq/m 3 , sensitivity 1.8 cpm @ 1000 Bq/m 3 independent of humidity, and the working temperature is −20 to 40 • C. The measurement chamber, equipped with a semiconductor detector and high voltage collection, is immune to ambient humidity. The unit includes a barometer (800-1200 mbar), a temperature sensor, a humidity sensor and an inclinometer. The sample rate is 3 h. CO 2 and CO monitoring was conducted in real time. This information along with humidity, temperature and dew point were measured with DL 302-DL303 produced by ICPDAS. Its specifications are: range 0 to 9999 for CO 2 , resolution 1 ppm, accuracy +/− 30 ppm, response time 20 s, range −10 to 50 • C, and resolution 0.1 • C. The sample rate is 1 s. The borehole was performed for seismic measurements, was filled with water and includes an accelerometer and a temperature sensor.

Data and Detection Methods
Data monitoring was conducted through independent acquisition systems that perform the first level of detection based on the trigger and detrigger levels with real-time message transmission to the http://geobs.infp.ro/web platform. The radon-222 monitoring equipment is autonomous, having the ability to stay connected on batteries for 2 months. The data were then downloaded manually using a software application from the device manufacturer. The information was processed automatically and transformed into a format compatible with the other data, namely ASCII text with a TAB delimiter. The analysis program converted the sampling periods from Table 1 to 1 min for large amounts of information. The data were saved in 1 h files which were transmitted to a database and viewed/queried through the http://geobs.infp.ro/platform.
An example of STA/LTA events detection [28][29][30][31] is shown in the Figure 2. Three methods were used: PC-SUDS, Processing of Seismic Data Stored in the Seismic Unified Data System, [32,33], Allen [34], and Hilbert transforms, which correlated with the envelope function [35]. The PC-SUDS was used for the seismic detection of the P wave. This method is useful for faster signals (100 Hz sample rate). For very low signals the Hilbert method is better (days sample rate). For our case the equivalent sample rate used was 60 s and we acquired good results with the STA/LTA method. This value was also used to integrate signals. Even if the data acquisition was conducted at 1 s, we used sampling periods of 60 s for the analysis of longer periods of time. Figure 2 shows the result of the Allen STA/LTA method with two limits-threshold pick and threshold valley, which were applied to the integral radon signal. In this case, we used two time windows and calculated the sum of the samples power for the two cases. Finally, we made the ratio between them (STA/LTA) and the result is a time function that indicates the signal anomalies. Depending on the case (radon, CO2, temperature), the two time windows were dimensioned.
The authors of [36] describe the methods for detecting anomalies applied to radon time series in correlation with seismicity. The standard deviation from the mean value, and the relationship between radon and seismicity are presented in Figure 3. In Figure 3a the standard deviation and mean values are shown for the whole period of time. Normally, this is carried out over limited periods of time in cycles in which the analyzed parameter varies. Figure 3a shows that the maximum annual variations are in November for Bisoca station. The intersections between the radon signal and the +/− 2 STD are not always correlated with seismicity ( Figure 3). Seasonal variations introduce errors in the application of this method. These could be eliminated by applying this detection of anomalies for shorter periods of time or by filtering seasonal variations. In the procedure described in the next paragraph, the analysis is conducted daily and seasonally to reduce false detections. Figure 3b shows the correlation between seismicity and the integral of the radon signal. We chose this method to highlight the dependence of radon concentration on time.

OEF Procedure
The OEF procedure is general and involves the following steps: (1) select the seismic area ( Figure 1); (2) determine the a, b parameter from Gutenberg Richter law [37][38][39][40][41], and other forecast parameters such as air ionization, radon concentration, CO 2 , temperature in borehole, telluric field; (3) compare the seismicity with precursor factors (geophysics, geochemical) and determine the probability and epistemic uncertain forecast time; (4) evaluate the emergency state using a logic tree; (5) decide the event state and send the information.
For this purpose, we developed two software applications that allow the updating of the predicted seismicity according to the real one and analyze the precursor factors. Part of the seismic information comes from ISC web site: International Seismological Centre, Online Bulletin, http://www.isc.ac.uk, Internatl. Seismol. Cent., Thatcham, United Kingdom, 2016, Romplus Catalog, and from EMSC (European-Mediterranean Seismological Centre).
(1.) Select the seismic area: For the selection of the seismic zone, a large area was initially selected, after which we chose the area of interest. In this stage it is necessary to have geological information. The location of the sensors was established in the faulty areas (see Discussion and Conclusions paragraph) where it is estimated that the gas emission will be higher (radon concentration is one of precursor factor). A seismic analysis requires data from a catalog or from the ISC (International Seismological Centre). In Figure 1, the selected area is marked with a red polygon and the method can be applied to any area for which the ISC provides the required data. The preparation zone is also used to establish the area of interest. This is determined using an empiric of Dobrovolsky [42] formula.
(2.) Determine short-term changes in seismicity rates using the Gutenberg Richter magnitude-frequency relationship and the daily-seasonal evolution of forecast parameters: The Gutenberg Richter (GR) law [37][38][39][40][41] states that earthquake magnitudes are distributed exponentially as where N(m) is the number of earthquakes with a magnitude larger or equal to "m", "b" is a scaling parameter and "a" is a constant. We use the least square regression method for determining "a" and "b". The magnitude of completeness Mc is determined by applying a regression algorithm until the error starts to grow. Detection algorithms use both level and duration thresholds ( Figure 2, STA/LTA case) to highlight deviations from the norm. The choice of these parameters is also made according to the daily (night-day) and seasonal variations [43]. Even if the analysis is carried out over the same season, differences may occur. We chose the years 2017, marked with (a), and 2018, marked with (b), in order to highlight the annual variations of radon and CO 2 in Bisoca and Nehoiu, respectively. Figures 4-7 present radon variation in Bisoca station in spring, summer, autumn and winter and Figures 8-11 represent the CO 2 variations in Nehoiu. Radon registers higher values at midnight in the warm period of the year and becomes uniform before decreasing in winter. CO 2 has a similar behavior. Figure 6b shows an increase in radon preceding the earthquake with a magnitude of 5.8R. In the case of CO 2 , Figures 9a and 10b indicate a maximum of CO 2 before the seismic event at a time of day when it should have decreased. This type of analysis requires a large amount of information in order to eliminate false decisions.        A high variation of radon or CO 2 could be seasonal without having a connection with seismicity. If a variation of the GR parameters is superimposed in this case, it can lead to a false decision. A machine learning method could be a solution to this problem [36].
(3.) Compare the seismicity with precursors factors (geophysics, geochemical) and determine the probability and the epistemic uncertain forecast time: Figure 12 presents in the same plot the geophysical information next to GR "a-b" evolution for Vrancea area in period 2018-2020/03/31. The signals related to gas emission and air ionization are time integral ("*dt" notation on graphs is 60 s) to incorporate the effect of accumulations over time. The first trace, BISRd (see Table 1), represents the time series of radon but the detection was performed on the integral signal because this method is more efficient. The following signals represent the radon in Lopatari, Nehoiu, Vrancioaia and Plostina. The magnitude is Mdpvs_M and is higher when the GR a-b is low long time. Air ionization (NhCOb_Iop, NhCOb_Ion) and CO 2 (NhCOb_Co2) are measured from Nehoiu station. The temperature in a borehole (Table 1, Tf and Tforaj in PLOR station) is a precursor factor analyzed in Figure 12 along with the gas emission. In all stations we also have meteorological information. The measured values should be corrected with temperature, humidity and atmospheric pressure. The equipment that determines the concentration of radon and CO 2 have these sensors embedded in order to make these corrections. For a general analysis we consider that the monitored area is small and homogeneous from a meteorological point of view. Sometimes this is not always the case. Monitoring stations are located in isolated areas near villages where no polluting industrial activities take place that could affect our measurements. However, in two cases we were surprised to find the presence of CO, indicating pollution. These situations were partially explained (in one case a saw was propelled by a diesel engine near the monitoring station, and in another there were live fires due to natural well gas emissions in the Lopatari area). In the first situation it was decided to move the sensors (Vrancioaia) and for the second case (Lopatari) the research was at the beginning. For a more detailed analysis we have information about the movement of clouds and their electrostatic charge through a network of radars (Boltek equipment).
The correlation between GR and seismicity is presented in Figure 13. The time between the moments of earthquakes (DeltaT, Figure 13) is useful but it indicates an increase in seismicity not necessarily an event. So, when the GR-a decreases to an approximate value of 2.9 (a minimum for Vrancea) for a period longer than 41 days we will have a seismic with Mw> 4.5. GR a-b were calculated on a 10-week window and with a 1-week (7-day) step.
We calculate the probability and epistemic uncertain forecast time shown, for example, in Figure 12. The results can be found in Table 2, while the determination data are presented in Table 3. It is important to take into consideration the distance between the measuring location and the hypocenter of the earthquake. S. Süer et al. in [44] demonstrate that the concentration of radon-222 at an active fault is correlated with the total earthquake energy (TEE). This parameter is inversely proportional to the square of the hypocenter distance. For this reason there are differences between the values determined in Table  2. These differences can be minimized with a proper selection of the investigated area. The choice of the monitoring area was made according to the preparation zone determined by the Dobrovolsky [42] formula. The GR a-b parameters were included in the table because they represent the magnitude-frequency relationship used by OEF applications (Figure 13).    Table 2 is the result of the analysis of evolution of all parameters from Figure 12. We measured the time determined by the detection of precursor anomalies and the moments at which the earthquakes occurred. We decide with TRUE (1) the cases in which, after an anomaly, we had an earthquake and false otherwise, FALSE (0). The total number of 1 over the all determinations is the probability from Table 2. The standard deviation of forecast times (TRUE cases) is the epistemic uncertainty [45]. For example, for radon we have 20 of TRUE values from 48 determinations (Table 3) and 4 out of 7 cases for BISRd. The uncertainty of the time forecast related to radon in Bisoca (BISRd) is determined by 11.9, 20.4, 21.2, 12.5 forecast days corresponding to the 4 earthquakes produced in the analyzed time interval (Table 3). We select earthquakes greater than 4.6R in the area marked in red in Figure 1.
In Table 3, slope + means a rising curve, max a maximum value and "dt" a deviation visible using a derivative function.
The detection of events is realized in two steps: in real time directly in the monitoring area and second offline with specialized software that generates graphs such as the ones in Figures 12 and 13. The first detection aims to draw attention to an event. This is analogous to the difference between the first localization of an earthquake (which is normally carried out automatically) and its revision (what you can find at International Seismological Centre ISC or in a seismic catalog). The web platform http://geobs.infp.ro/ displays real-time warnings and parameter values over 5 min, such as those in Figure 14.     Another method of analysis and detection is cross correlation [46]. An example is in Figure 15, in which we analyze the CO 2 correlation in Vrancioaia and Nehoiu before producing a 5.2R earthquake. The results may be better in case of higher seismicity. V.D. Rusov et al. in [47] use a cross correlation between radon and magnetic field. The common element is the temperature on which the radon concentration depends, but which also has a local piezoelectric effect that modifies the magnetic and telluric fields. The gas emission depends on the soil temperature which is correlated with the solar radiation [48]. In Plostina station we have a net radiometer (Table 1) to determine the direct and reflected solar radiation that allows the determination of the soil temperature. A study [49] showed a correlation between seismicity and albedo. As a result, solar radiation is also a factor to consider in a decision logic tree.
The radon, CO 2 , +ions+ and -ions graphs from this study represent the filtered values of time series with 1 min-resolution. In Table 1 each parameter has a defined sampling rate ranging from 3 h to 0.1 s. Higher sampling frequencies are useful in spectral analysis but a 1 min sampling period is optimal in analyzing data over long periods of time without a significant loss of information affecting the detection of events.
(4.) Evaluate the emergency state using a logic tree: The probabilities and epistemic uncertainties from Table 2 form the foundation of a decision in a logical tree. Epistemic uncertainties were quantified for gas emissions, temperature in the borehole in correlation with GR a-b parameters that represent the short-term changes in seismicity rates. The answer to how a logical tree would work for geophysical and geochemical parameters can only be given by experimentation. A decisionmaker requires certain information in order to declare a seismic event and take actions to mitigate the effects in a cost-benefit approach. Seismic sources and the monitoring network involved in this study are presented in Figure 1. The result of Table 2 can be presented in the form of a logical tree with weights associated for each parameter. Bommer et al. [50] presents transparent weighting procedures for logical tree branches. Probabilities can be used in other vote-based methods, such as many seismic digitizers [33]. In the case of Table 2, the weights were calculated using the criterion that the sum of the branches is 1 in the logical tree ( Figure 16). The Figure 16 is a branch of a bigger logical tree that is a decision element in an OEF structure. This method allows the collaboration between the tasks of a complex application. A logical tree made for geophysical and geochemical parameters becomes a branch in a complex time-dependent seismic hazard assessment application. The collaboration can be carried out through a flexible structure, such as the one in Table 2, which can have new stations and other types of sensors. In this study we use only a part of our multidisciplinary network. A Bayesian probability assessment can be applied in this case using the data from Table 2 and adding new information to improve the time-dependent seismic hazard assessment. For example, new parameters can be added such as soil temperature, propagation of VLF-ULF-ELF (Very Low, Ultra Low, and Extremely Low frequency) radio waves, infrasound, the magnetic field and meteorological data. In addition to the measured information, the results of the analysis of other information, such as the delay time between earthquakes and the assessment of cumulative seismic energy, can also be used. A probabilistic seismic hazard assessment considering epistemic uncertainties, logical trees and a Bayesian approach is presented by Gottfried Grünthal et al. in [51].
SELENA is another tool with application in seismic risk and loss assessment using a logical tree [52]. In this case the inputs are ASCII plain text files that represent tables with rows and columns. Another example to convert a logic tree into a XML file is: <?xml version="1.0" encoding="UTF-8" standalone="yes"?> <data-set xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance\T1\textquotedblright> <record> <Type>Radon</Type> <Station>BISRd</Station> <Probability>0.57</Probability> <Uncertain>4.98</Uncertain> <Weight>0.18</Weight> </record> <record> <Type>Radon</Type> <Station>LOPRd</Station> <Probability>0.29</Probability> <Uncertain>3.93</Uncertain> <Weight>0.14</Weight> </record> . . . (5.) Decide the event state and send the information: A probabilistic seismic hazard assessment is based on information updated by automatic detection or an offline analysis that sends the information to a logic tree. Detection must be validated by the daily and seasonal variations of the forecast parameters. A large earthquake for the Vrancea area would be 7.2R or above, but we did not have such data in our analysis. We assume that in such a case both the radon and theCO 2 levels from all stations would increase. If GR a-b would indicate variations such as the ones in Figure 12 we can declare an emergency state. The first detection level automatically sends data to a web platform. The second level, the offline analysis (Figure 12), combined with a logical tree will decide the level of emergency. This information is retrieved by a server, Figure 17a, which sends alerts to the entire network using software (b), and displays the earthquake parameters (c). The software presented in Figure 17 works in an EEW (Earthquake Early Warning) system.  The traffic light symbol from Figure 17b has 5 states: green-no warning elements; green-yellow-attention; yellow-warning; yellow-red-alert and red-major alert.

Discussion and Conclusions
The presented procedure highlights five steps. The first one, selecting the seismic zone, depends on the geological information we have. Gas emissions are predominant in fault zones ( Figure 18). In [53], the crustal structure of the Vrancea area is shown and Figure 19 presents a resistivity tomography of a particular monitoring location, Bisoca, located on the Casin-Bisoca fault, and marked in the general plan of the faults; Figure 18. Knowing the structure of the area, an optimal positioning of the monitoring station can be achieved.  All monitoring stations have been installed in fault areas, but this is not enough. Improper installation of the sensors will affect the measurements. Figure 20 shows how the CO 2 detector was initially mounted in a room where the data was affected by staff activity (top) and in the second stage it was introduced in a box and fixed to the outside (bottom). In the second phase of the procedure the daily-seasonal evolution of forecast parameters (Figures 4-11) was presented for only two years. For a correct analysis, it would have been beneficial to process larger periods of time, especially since the annual variations are affected by the global warming phenomenon. The analysis in this study focused on the correlation between gas emission and seismicity but in an OEF application many parameters are used to reduce the number of false errors. For example, N. Kastelis and K. Kourtidis in [54] analyze the correlation of an atmospheric electric field with CO 2 and S. Abbad et al. in [55] present the influence of meteorological and geological parameter variables on the concentration of radon in soil. It is more difficult to determine the interdependence between parameters such as radon activity and CO 2 flux in soil, which implies a correction of measurements [56]. Another mentioned correlation is between radon variations and the magnetic field [48]. This correlation could not be demonstrated in the case of radon measurements in the Muntele Rosu tunnel (Romania), because the humidity exhibits large oscillations dependent on water infiltrations in the mountains and on air flow from the tunnel, as presented by A. Mihai et al. [57]. For this reason the radon detector was moved to another location (code PLRd2 in Table 1). This is an example that shows how a method cannot be generalized because it depends on the specifics of the location.
We used long-term data to test the hypotheses, methods and to analyze events. In reality we have to make decisions on real-time information. In Figure 3 we have an example in which the difference between a long (e) and limited (f) period of time analysis is observed.
Regarding the next phase (3) of the procedure, Table 2 shows the result of an offline anomaly detection analysis performed by an experienced operator. If we had applied the automatic detection method based on the STA/LTA algorithm, the result would have been faster but more inaccurate. The results also depend on the number of sensors. For CO 2 or the drilling temperature Tf, only one location was used, which determines the maximum weight. The uncertain forecast time was determined using the simplest method that is the standard deviation in our case ( Table 2). N. Ridler et al. in [45] combined all uncertainties using the Law of Propagation of Uncertainty (LPU) using the root-sum-squares (RSSs) method of standard deviations to give a global uncertainty (a global standard deviation).
In phase (4), the use of a logical tree can be simplified by using a voting system [33] or a decision tree that does not use weights [36]. The way in which the logical tree is made depends on the previous stages and only through experimentation can we adjust the weights of each branch.
In the final phase we must not forget that the methods are empirical in many cases and that decision-making requires responsibility if it involves departments specialized in intervention.
The main conclusion is that there is not a single precursor factor, not even a model, and that only a multidisciplinary network allows a complex analysis that reduces the number of false errors and increases the probability of a correct decision. In the case of a complex OEF system comprising a large monitoring area with different geological structures (e.g., Europe), the recommended solution is to decentralize the decision by applying the procedure described to independent zones that transmit alerts to a general information portal. In this context, the answer to the question J. Zechar wrote in his article [58], "Is Europe-wide operational earthquake forecasting (OEF) possible?", is affirmative.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.