Using geoelectric field skewness and kurtosis to forecast the 2016 / 2 / 6 , ML 6 . 6 Meinong , Taiwan Earthquake

The earthquake-alarm model developed by Chen and Chen (Nat. Hazards 2016) is investigated to validate its forecasting performance for the 2016/2/6, ML 6.6 Meinong, Taiwan earthquake. This alarm model is based on geoelectric field skewness and kurtosis anomalies. The model parameters, such as the detection range and predicted time window, allow us to estimate the empirical relationships between geoelectric anomalies and large earthquakes. As a result, the skewness and kurtosis anomalies are shown to appear before the Meinong earthquake on the four neighboring stations (LIOQ, WANL, KAOH, and CHCH). According to the model analysis a time lag exists between anomaly clusters and earthquakes, depending on local geological features, as well as the durations over which anomalies are continuously observed, which might also display time dependence. In conclusion, this alarm model is able to correlate earthquakes and geoelectrical anomalies, with promising usefulness in forecasting large earthquakes. Article history: Received 29 August 2016 Revised 31 October 2016 Accepted 1 November 2016


InTrodUcTIon
This work investigates an earthquake-forecasting algorithm based on geoelectric field skewness and kurtosis (Chen and Chen 2016).The M L 6.6 Meinong earthquake hit the southern part of Taiwan (120.54°E,22.92°N) at 03:57 AM, 2016/2/6 (UTC+8) at a depth of 14.64 km, causing the death of over one hundred people, the collapse of tens of buildings, and other significant infrastructure damage (cf. the listed Website links: GEER Association 2016; NCDR 2016; Wikipedia 2016).Earthquake forecasting, therefore, is a critical task.The effective earthquake forecasting should provide information that includes the time, location, magnitude and probability for an impending earthquake.Effective forecasting should also help scientists simulate shaking hazard maps, as well as warn the public in areas that are likely to be seriously damaged if a large earthquake strikes.
Evidence of the precursory nature of these indexes proposed in these previous works was provided via dependence measures, but without tests of their forecasting skills or within systematic analyses on full time series.To solve the ambiguity, Chen and Chen (2016) developed a model to examine the empirical relationship between earthquakes and anomalous geoelectric field statistics.They proposed the "Geoelectric Monitoring System's Time of Increased Probability" (GEMSTIP) model and applied this model to estimate the occurrence time of the anomalies, the predicted time window for the earthquakes and other parameters associated with the anomalies.The GEMSTIP model, as a ruler measuring the distance between two points, helps us understanding the empirical connection between earthquakes and geoelectric fields.This paper provides a validation step for this model using the available datasets from 2013/1/16 to 2015/12/31 as the training set, and the datasets from 2016/1/1 to 2016/3/31 as the forecasting set.We expect that the best parameters obtained from the training period allow us to forecast the M L 6.6 Meinong earthquake.Chen and Chen (2016) found time lags between two earthquakes with M L ≥ 6 and their preceding clustered anomalies.A time lag is the elapsed time between the end time of the clustered precursory anomalies and an earthquake occurrence.We modify their model in this paper by taking the time-lag effect into account.Other questions addressed in this paper are as follows: (1) What is the best model, the original one or the modified one?(2) Why does a time-lag effect exist? (3) What should be the optimal training window size for better forecasting performances?Answering those questions may help to obtain a deeper understanding of the potential links between earthquakes and the geoelectric field.

dATA
The earthquake catalog is routinely processed and released by the Central Weather Bureau (CWB) of Taiwan.The geoelectric field is registered by the Geoelectric Monitoring System (GEMS), which is maintained by Prof. Chien-Chih Chen and his team.We used all M L ≥ 5 earthquakes and calculated the skewness and kurtosis of the daily geoelectric field distribution, using a sampling rate of 1 Hz (cf.Chen and Chen 2016).Figure 1 shows the geoelectric station spatial distribution and the earthquakes with M L ≥ 5.The solid magenta star is the 2016/2/6 M L 6.6 Meinong earthquake, and LIOQ is the nearest station to this event.We used both the earthquake catalog and geoelectric data from 2013/1/16 to 2015/12/31 as the training set, and data from 2016/1/1 to 2016/3/31 as the forecasting set.We also selected differ-

METhods
This section is composed of (i) the "Geoelectric Monitoring System's Time of Increased Probability" (GEMSTIP) analysis algorithm workflow, (ii) the original GEMSTIP model proposed in Chen and Chen (2016), (iii) the modified GEMSTIP model considering the time-lag effect, (iv) the priority-search area concept (PSA), and (iv) model ranking used to select the best model in similar performance cases.

GEMsTIP Analysis Algorithm Workflow
Chen and Chen (2016) developed the "Geoelectric Monitoring System's Time of Increased Probability" (GEM-STIP) model and tested the relationship between earthquakes with M L ≥ 5 and the anomalous skewness and kurtosis of geoelectric fields.Figure 2 shows the methodology workflow proposed in Chen and Chen (2016), which we shall also follow in this paper.The skewness and kurtosis of the daily geoelectric field distributions for North-South (NS) and East-West (EW) components are first calculated.There are four indexes per day.The start and end time of a training period is selected next.Chen and Chen (2016) used the training period from 2012/1/1 to 2014/12/31, whereas we use the training period from 2013/1/16 to 2015/12/31 and also discuss the effect of different training windows between the start and end times.The upper and lower thresholds used to detect anomalies are defined for each index as the median ±3 times the interquartile range in the selected time window (shown in the upper panel of Fig. 3).An anomalous index is detected if the observed index value is beyond one of the thresholds.The anomalous index number (AIN) is then counted and summed according to the anomalous indexes of the skewness and kurtosis per component (shown in the lower panel of Fig. 3).Note that AIN is a discrete value.The GEMSTIP model is finally applied to the AIN series and the target earthquakes.The best model parameters are then determined for the considered training window.

original GEMsTIP Model
We consider the Chen and Chen model ( 2016) as the original GEMSTIP model.In the original model the parameter vector G o is: G o = (Rad, Dep, t cg , N thr , t thr , t obs , t pred ), where the superscript O means the original GEMSTIP model.The meaning of each parameter is shown in Fig. 3.The parameter Rad is the radius for selecting earthquakes within the region centered at a given station.The parameter Dep is the depth above which we select earthquakes.The parameter t cg is the coarse-grained time of an earthquake.The parameter N thr is a threshold number to label a day as anomalous if AIN ≥ N thr for that day.The parameters t thr , t obs , and t pred are time windows, which mean that if the number of anomalous days within the observation time window t obs is greater than or equal to the threshold duration t thr , then the immediate future time window t pred is alarmed as a time of increased probability (TIP) (cf.Chen and Chen 2016).Note that the parameter t cg is equal to t pred hereafter.
We then define the coarse-grained earthquake (CGEQ) time function as follows: , , where Mag stands for the magnitude of the earthquake.We herein use Mag ≥ 5. We also define the prediction index time function, i.e., time of increased probability (TIP) as follows: , , greater than or equal to the value t thr in the period t et obs ≤ t ≤ t e , the TIP is declared to be 1 in the period t e ≤ t ≤ t e + t pred .Note that it is not necessary that the anomalous days are consecutive.In Fig. 3, we see a continuous TIP segment and a continuous CGEQ segment.We denote them as t TIP and t CGEQ , respectively.We then define a continuously anomalous time t ANO as follows: In Eq. ( 4), TP means the number of true positives, FP the number of false positives, FN the number of false negatives, and TN the number of true negatives.For example, CGEQ = (1, 0, 0, 0, 1) and TIP = (1, 1, 0, 0, 0).In this case, TP = 1 based on the first elements of CGEQ and TIP, FP = 1 based on the second, FN = 1 based on the fifth, and TN = 2 based on the third and the fourth.We then define the model performances as follows: 1( ) In Eq. ( 5), C1 is the TN performance, F1 is the TP performance, and R is the ensemble performance.We set the parameters Rad = (30:5:60, 70:10:100) (km), Dep = (10:10:50, 100:50:300) (km), N thr = (1:1:4), t thr = ( : : t 1 1 3 2 obs ) (day), t obs = (5:5:30, 40:10:100) (day), t pred = (5:5:30, 40:10:100) (day).Note that the data format is (initial value : increment : end value) which means that the values of one parameter are from an initial value to an end value with an incremental step.The number of original tested models thus amounts to 1361360.

Modified GEMsTIP Model
Chen and Chen (2016) found that there is a time lag between a large earthquake and the clusters of geoelectric anomalies.We take the time-lag effect into consideration in this paper.We, hereafter, called the GEMSTIP model considering the time-lag effect as the modified GEMSTIP model.
The parameter vector G M is: G M = (Rad, Dep, t cg , N thr , t thr , t obs , t pred , t lag ), where the superscript M means the modified GEMSTIP model.The time lag parameter t lag is the elapsed time between the end of the observation time window t obs and the beginning of the predicted time window t pred (shown in the lower panel of Fig. 3), whereas the other parameters have the same meanings as in the original model.The parameter t cg is also equal to t pred hereafter.

Priority-searching Area
The priority-search area concept (PSA) is now introduced.Chen and Chen (2016) plotted the C1-F1 diagrams to observe the model performance patterns and the location of the best model.However, it is not convenient to compare the C1-F1 patterns when the number of the models with different tested conditions is strongly increased.For the sake of simplifying the C1-F1 diagram comparisons, we divided the C1-F1 area into 7 PSAs, as shown in Fig. 4. PSA 1 corresponds to R > 1, C1 > 0.5, and F1 > 0.5; PSA 2 to R ≤ 1, C1 > 0.5, and F1 > 0.5; PSA 3 to R > 1 and F1 ≤ 0.5; PSA 4 to R > 1 and C1 ≤ 0.5; PSA 5 to R ≤ 1, C1 > 0.5, and F1 ≤ 0.5; PSA 6 to R ≤ 1, C1 ≤ 0.5, and F1 > 0.5; and finally PSA 7 to C1 ≤ 0.5 and F1 ≤ 0.5.Intuitively, PSA 1 is the best region for the model performance, while PSA 7 is the worst region.Hence, PSAs 1 to 7 correspond to the model performances in descending order.We looked for the corresponding PSA of each model and eventually determined the best model parameters for each PSA.Based on the PSA number, we are able to image the C1-F1 patterns without plotting the C1-F1 diagrams (cf.Figs. 4 and S3 of Chen and Chen 2016), which is rather helpful when comparing the enormous tested model sets.Moreover, when the best model occurred in PSA 5, 6, or 7, we did not use this as the best model.

Model ranking
Consider the original GEMSTIP model as an example and assume a M L 6 earthquake occurred 40 km from a given station, with a depth of 39 km. 4 is overestimated.One can imagine similar considerations for the other parameters.This amounts to a kind of LASSO (least absolute shrinkage and selection operator) penalization of the parameters (Tibshirani 1996).We solve this problem using the spatial LASSO and temporal LASSO in order to rank the models.The spatial LASSO is the detected volume: The temporal LASSO is the 'volume' of anomalies: If there are models with the same best ensemble performance R, we first choose the model with the minimum of the spatial LASSO L S .If models with the same R value still occur, we then choose the model with the minimum temporal LASSO L T .In practice, we observed that the spatial LASSO is sufficient to rank models with similar R values most of the time.

rEsULTs
Figure 5 shows the skewness and kurtosis series at 4 stations, which are nearest to the epicenter of the 2016/2/6 M L 6.6 Meinong earthquake.LIOQ, the closest station (at 17.30 km), exhibits anomalous skewness and kurtosis clusters within 2015/6/18 to 2015/11/13, corresponding apparently to the Meinong earthquake.There is thus a time lag of approximately 86 days between the end time of the anomalous clusters and the earthquake occurrence.At WANL, the second nearest station (at 37.17 km), the skewness and kurtosis start to deviate from their thresholds from 2015/7/10.Yet, before, during and after the Meinong earthquake, the anomalies persist.Note that WANL started recording on 2012/2/2, and thus features a gap in the data at the beginning of its time series (Fig. 5b).A similar situation holds at KAOH, the third nearest station (at 39.30 km).KAOH's skewness and kurtosis start to deviate from its thresholds from 2015/5/31, and keeps on up to the end of the dataset.CHCH, the fourth nearest station (at 51.2 km), displays kurtosis deviating from its thresholds from 2015/9/16, but its skewness behaves relatively quiescently.For WANL and KAOH, we also observe different continuous anomaly durations from 2015/10/14 to 2016/3/31 and from 2015/6/1 to 2016/3/31, respectively, before and after the Meinong earthquake occurrence.For other periods, at those 4 stations, there are other anomalies corresponding apparently to other large earthquakes.Because it is hard to determine the relationship using the naked eye, we applied the "Geoelectric Monitoring System's Time of Increased Probability" (GEMSTIP) model to quantitatively measure the relationship between those anomalies and earthquakes.
After the GEMSTIP analysis, Fig. 6 shows the CGEQ-TIP matching diagrams for the best original and modified models.in Table 1.We then used the determined parameters to forecast the earthquakes with M L ≥ 5 occurring from 2016/1/1 to 2016/3/31.We first compared the CGEQ-TIP matching results in the training set.Overall, the original and modified models have a similar pattern for the CGEQ-TIP matchings.However, comparing their performances, the modified model ( C1 t = 0.99, F1 t = 0.94, R t = 1.29) slightly outperforms the original model ( C1 t = 0.97, F1 t = 0.93, R t = 1.26).Note that, hereinafter, the operator .means the average of a certain parameter over the 20 stations.The subscript t means the training set, while the subscript f means the forecasting set.For the original model, we obtained Rad = 64 (km), Dep = 44.5 (km), N thr = 2, t thr = 10.85 (days), t obs = 43.5 (days), t pred = 48.75(days) (averaged from Table 1a).For the modified model, we obtained Rad = 69 (km), Dep = 49 (km), N thr = 2, t thr = 10.65 (days), t obs = 49.5 (days), t pred = 34 (days), t lag = 22.5 (days) (averaged from Table 1b).Those averages are summarized in Table 2.The average predicted time window t pred in the modified model is smaller than that in the original one, which allows us to point out more precisely when an earthquake will occur.Focusing on the forecasting set, the performances of the original model are C1 f = 0.64, F1 f = 0.34, and R f = 0.83, while those of the modified model are slightly better with C1 f = 0.70, F1 f = 0.49, and R f = 0.90.WANL forecasted the M L 6.6 Meinong earthquake and the other two earthquakes with M L ≥ 5 in both the original and modified models.On the other hand, LIOQ and YULI do not forecast the Meinong earthquake in the original model, but forecast the earthquake in the modified one.In the supplementary materials, Figs.S1 and S2 show the CGEQ-TIP matching results of the original and the modified models with different training windows, and Tables S1 and S2 list their best model parameters.Figure 7 shows plots of averages of the best parameters versus the size of the training windows.For the parameters Rad , Dep , N thr , t thr , t obs , t pred in both the original and modified models, and the parameter t lag in the modified model, we see that changing the training window size changes slightly the average parameters.These results suggest that the inverted models are rather robust.The predicted time windows t pred in the modified models are smaller than in the original models for all training window sizes, except for a training window of 180 days.The other average parameters have similar values in both the original and modified models.The original model proposed in Chen and Chen (2016) gives a predicted time window for an earthquake that is overestimated because of the time-lag effect between earthquake occurrence and geoelectric field anomalies.We truly improved the original model by considering the time-lag effect.Due to the smaller predicted time windows, we could forecast the earthquake occurrence time more accurately.
We next compared the average performances of the original models to those of the modified models when using different training windows.Figure 8 shows the plots of the average performances ( C1 , F1 , R ) versus the training windows (180:90:1080) (days) for both the original and modified models.We observe in Fig. 8a

dIscUssIon
In this work, we found that of the four nearest stations to the M L 6.6 Meinong earthquake, the clustered anomalous skewness and kurtosis in LIOQ became relatively quiescent prior to the Meinong earthquake.We consider this as the time-lag effect.Chen and Chen (2016) also observed time lags for two earthquakes with M L ≥ 6 in 2013 near PULI.Varotsos and his team also found time lags between SESs and earthquakes in Greece (Varotsos et al. 1993(Varotsos et al. , 2009) ) and Japan (Varotsos et al. 2013).In rock experiments, Triantis et al. (2008) found the pressure stimulated current recordings (PSC) featured a peak before failure.They speculated that the time lags were produced by the formation of the fracture plane.This plane forms at high stresses and the extensive cracking process that is associated with it limits the available conductive path in the sample bulk, inducing a consequent obstacle on the emitted PSC.We applied the modified GEMSTIP model to quantify the time lag between the defined anomalies and earthquakes.However, not all stations show such a time-lag effect.When observing the column corresponding to the time lag t lag in Table 1b, almost all stations in the coastal region display no time-lag, such as SHCH, HERM, FENL, SIHU, CHCH, RUEY, KAOH, and WANL, while stations in the central mountain region possess one.We suspect that the time lag depends on the local geology.The observed low resistivity under the central   mountain region suggests the presence of abundant ground water (Chen andChen 1998, 2000;Chen et al. 1998;Bertrand et al. 2009Bertrand et al. , 2012)).In the initial stage, small-scale irreversible deformations and fractures occur that generate electromagnetic radiation.Hence, the skewness and kurtosis of the geoelectric field vary largely.However, the abundant ground water starts to saturate those newly formed (or reactivated) cracks.It may also interfere with currents generated by the irreversible deformations.The generated electromagnetic radiation can then decay significantly when passing through the crack-saturated area.This scenario is depicted in Fig. 9.This might explain why there is a relative quiescence in skewness and kurtosis in the geoelectric field just prior to large earthquakes at stations deployed in the saturated area.
Besides the time-lag effect, we found that the duration of anomalous skewness and kurtosis periods vary.Consider SHRL, KUOL, TOCH, and HUAL for instance (in Fig. S2e of the supplementary materials): the continuous TIP times t TIP have different durations.Based on Eq. (3), it is suggested that the durations of lasting anomalies t ANO are also different.The varied lasting anomalous times seem to rely on the preparation process of different earthquakes.This phenomenon suggests that the earthquake preparation processes are different even within the same region.Therefore, when analyzing geoelectric fields in order to forecast an earthquake, the varied lasting anomalous times largely affect the forecasting performance.The different durations of lasting anomalies not only have a spatial dependence but also a temporal dependence within the same region.This issue might be solved by relating it to the crustal parameters that mainly influence the electrical structure, such as permeability, porosity and so on.
The forecasting TP performance F1 f m = 0.49 with a training window of 360 days is larger than F1 f o = 0.34.For training windows larger 360 days, F1 f o and F1 f m have similar values.The better forecasting performance is produced by smaller predicted time window values t pred in the modified model.YULI and LIOQ in the modified model can identify earthquakes in the forecasting set, whereas they cannot in the original model.On the other hand, the larger t pred of the original model easily leads to longer TIP times t TIP , and to longer coarse-grained earthquake times t CGEQ .This situation leads to unreasonable forecasts.Therefore, it Focusing on the KAOH and CHCH stations in Fig. 6, the CGEQ-TIP matchings are all empty.They also do not give any alarm for the Meinong earthquake.This is because their best models are located within PSA 5, and their performances C1 = 1 and F1 = 0. Hence, their best models have minimum parameters in the given ranges.For models located within PSA 5, 6, and 7, there is less relationship between earthquakes and anomalies skewness and kurtosis in the training periods.When meeting this situation, we moved the training window back to the previous period in which there are large earthquakes and anomalies.We used the reselected training datasets to re-optimize the models, and obtain the best model, called the historic best model.This historic model is substituted for the original best model.Based on Figs.5c and d, we selected the retraining periods from 2013/4/1 to 2014/3/26 for KAOH and from 2012/4/26 to 2013/4/20 for CHCH.The historic best parameters for KAOH and CHCH are listed in Table 3. Figure 10 shows the new GEMSTIP matching diagrams for the original and modified models.We found that KAOH and CHCH can give alarms and identify an earthquake with M L ≥ 5, but they cannot identify the Meinong earthquake within their detection ranges.They share the common parameter Dep = 10, but the depth of the Meinong earthquake is 14.64 km.In Fig. 8, we found that the retrained TP performances F1 f o = 0.39 and F1 f m = 0.53 (marked as red symbols) are all increased compared to their original values, which suggests that substituting the historic best parameters solved this no TIPs and no CGEQs problem in the training set.Although the retrained TN performances C1 f o = 0.56 and C1 f m = 0.62 both decrease, the retrained ensemble performances R f o = 0.81 and R f m = 0.88 remained almost the same.The modified GEMSTIP model may help to forecast earthquakes, but some problems remain.For instance, optimizing the large number of models is extremely timeconsuming and increasing the execution efficiency will be one important task.The abovementioned issue of the lasting anomalous time t ANO also should be carefully investigated to increase the forecasting accuracy of large earthquakes.The selection of ways to quantify the success rate of a given model should be considered.We choose a C1-F1 diagram to optimize the GEMSTIP models.However, there are other performance scores in a binary classification, such as a ROC diagram, a Molchan diagram, etc.We have to understand the advantages and disadvantages of those performance scores and select the suitable ones for the GEMSTIP models.
In summary, the GEMSTIP model is useful in estimating the empirical relationship between earthquakes and anomalous statistical geoelectric field indexes.We first showed that the anomalous skewness and kurtosis of the geoelectric field appearing before large earthquakes.We next showed that the time-lag effect can exist between the end of clustered anomalies and a large earthquake, and that this effect depends on the local geological features.Third, we found that the durations of continuous anomalies are different, which means that the preparation process of each earthquake is different and shows spatial and temporal dependences.Finally, the GEMSTIP model identified the M L 6.6 Meinong earthquake and other large earthquakes.The GEMSTIP model helps us to understand which factor affects the connection between earthquakes and geoelectric   fields (which even could be extended to any geophysical data), and also forecast large earthquakes.Although some problems remain to be resolved, this model is promising for building an earthquake-forecasting system and make it practical in the future.
ent training windows in the training period (2013/1/16 to 2015/12/31) to discuss the effect of the training window size on the forecasting performances.The training windows thus span from 180 days to 1080 days, with a step of 90 days, simply denoted as (180:90:1080) (days).The first training set spans from 2015/7/5 to 2015/12/31 with a training window of 180 days.The second training set spans from 2015/4/6 to 2015/12/31 with a training window of 270 days.The third training set spans from 2015/1/6 to 2015/12/31 with a training window of 360 days, and so forth.The end time for all training sets is 2015/12/31.In the forecasting period spans from 2016/1/1 to 2016/3/31.Seven target earthquakes are denoted as colorful stars in Fig. 1.Four of them are inland and around LIOQ, the others are located on the eastern coast.

Fig. 1 .
Fig. 1.Spatial distribution of the geoelectric stations and the earthquakes with M L ≥ 5 in Taiwan.Open stars are the earthquakes with M L e [5, 6), and solid stars are the earthquakes with M L ≥ 6. Gray and dark stars are the earthquakes before 2015/12/31 within the training set, while pink and magenta stars are the earthquakes from 2016/1/1 to 2016/3/31 within the forecasting set.(Color online only) If the models G O 1 = (50, 40, 2, 4, 10, 15) and G O 2 = (60, 40, 2, 4, 10, 15) have the same ensemble performance R, we prefer model G O 1 because the parameter Rad = 60 in G O 2 is overestimated in this case.Or if models G O 3 = (50, 40, 2, 4, 10, 15) and G O 4 = (50, 50, 2, 4, 10, 15) have the same ensemble performance R, we choose model G O 3 because the parameter Dep = 50 in the model G O Fig. 4. Schematic diagram of priority-searching areas (PSA) of the C1-F1 diagram.The vertical line is the TN performance C1 = 0.5, the horizontal line is the TP performance F1 = 0.5, and the quarter circle curve is the ensemble performance R = 1.PSAs 1 through 7 are ranked in descending order of performance.If the best model occurs in PSAs 5, 6, or 7, it should not be used for forecasting.(Color online only) Fig. 6.CGEQ-TIP matching diagrams of (a) the best original model and (b) the best modified model.The training set is from 2015/1/6 to 2015/12/31 with a training window of 360 days, and the forecasting set is from 2016/1/1 to 2016/3/31.Black and blue horizontal lines are TIPs.Gray and green horizontals are CGEQs, the time expansions of earthquake occurrence.Open stars are the earthquakes with M L e [5, 6), and solid stars are the earthquakes with M L ≥ 6.The gray-level part on the left concerns the training set, while the colorful part on the right deals with the forecasting set.Gray dotted lines are auxiliary guides for the eye.(Color online only) Figure7shows plots of averages of the best parameters versus the size of the training windows.For the parameters Rad , Dep , N thr , t thr , t obs , t pred in both the original and modified models, and the parameter t lag in the modified model, we see that changing the training window size changes slightly the average parameters.These results suggest that the inverted models are rather robust.The predicted time windows t pred in the modified models are smaller than in the original models for all training window sizes, except for a training window of 180 days.The other average parameters have similar values in both the original and modified models.The original model proposed inChen and Chen (2016) gives a predicted time window for an earthquake that is overestimated because of the time-lag effect between earthquake occurrence and geoelectric field anomalies.We truly improved the original model by considering the time-lag effect.Due to the smaller predicted time windows, we could forecast the earthquake occurrence time more accurately.We next compared the average performances of the original models to those of the modified models when using different training windows.Figure8shows the plots of the average performances ( C1 , F1 , R ) versus the training windows (180:90:1080) (days) for both the original and modified models.We observe in Fig.8athat the TN performances C1 t o and C1 t m in the training sets (subscript t) of the original model (superscript o) and the modified model (superscript m) are mostly above 0.9, and slightly decrease with the increasing training windows size beyond 540 days.The TN performances C1 f o and C1 f m in the forecasting sets (subscript f) are all decreased by approximately 0.25 when compared to the training TN performances C1 t o and C1 t m .Additionally, the forecasting TN performances C1 f m are mostly larger than C1 f o except for training windows of 180 and 450 days.Focusing on Fig. 8b, the training TP performances F1 t o and F1 t m decrease approximately

Fig. 7 .
Fig. 7. Line charts of averages of best parameters versus training windows of the original and modified models for (a) the average Rad, (b) the average Dep, (c) the average N thr , (d) the average t thr , (e) the average t obs , (f) the average t pred , and (g) the average t lag .(Color online only)

Fig. 8 .
Fig. 8. Line charts of average model performances versus training windows, (a) the TN performance C1, (b) the TP performance F1, and (c) the ensemble performance R. Black lines with circles are the training performances of the original models, and blue lines with circles are for the forecasting performances.Black lines with squares are the training performances of the modified models, and blue lines with squares are for the forecasting performances.The red circle is for the forecasting performances of the retrained original model, and the red square is the forecasting performances of the retrained modified model.The operator .means the average of parameters for the 20 stations.The superscript o means the original model, the superscript m means the modified model, the subscript t means the training set, and the subscript f means the forecasting set.(Color online only)
Imagine that the window set of t obs and t pred find anomalies in the t obs only one time, so that t TIP = t pred and t ANO = t obs .Hence, the start time point of t ANO is that of t TIP minus t obs , and the end time point of t ANO is that of t TIP minus t pred .
The relationship between t ANO and t TIP in the modified model is the same as in Eq. (3).However, imagine that the window set of t obs , t pred , and t lag find anomalies in the t obs only one time, so that t TIP = t pred and t ANO = t obs , but there is a time gap t lag between t ANO and t TIP .Hence, the start time point of t ANO is that of t TIP minus (t obs + t lag ), and the end time point of t ANO is that of t TIP minus (t pred + t lag ).

Table 1 .
Best GEMSTIP parameters of (a) the original model and (b) the modified model.The training set is from 2015/1/6 to 2015/12/31 with a training window of 360 days.

Table 2 .
Averages of (a) best parameters and (b) performances of the original and modified models.The training set is from 2015/1/6 to 2015/12/31 with a training window of 360 days.

Table 3 .
Best GESMTIP parameters of (a) the original model and (b) the modified model with retrained CHCH and KAOH.The (re)training set of CHCH is from 2012/4/26 to 2013/4/20 with a training window of 360 days, while the (re)training set of KAOH is from 2013/4/1 to 2014/3/26 with a training window of 360 days.

Table S2 .
Best parameters of the modified GEMSTIP models with different training windows (T w ).