Spatio-temporal evolution of low-magnitude seismicity before the May 24, 2013, Sea of Okhotsk earthquake recovered by waveform cross correlation. Is it an earthquake prediction case?

 

Abstract

According to the International Data Centre (IDC), the Sea of Okhotsk earthquake occurred at 05:44:49.7 on May 24, 2013, had coordinates 54.89°N,153.31°E, m_b=6.27, and depth of 604 km. The USGS moment magnitude is 8.3. The previous event detected by the IDC in the surrounding volume 53°N-57°N, 151°E-155°E, depth from 400 to 700 km occurred on February 6, 2012. Using the same seismic data from the stations of the International Monitoring System together with detection and phase association methods based on waveform cross correlation, a series of low-magnitude earthquakes was recovered immediately before this major earthquake. More than 200 events obeying the Event Definition Criteria adapted by the IDC were found  between May 13 and the mainshock, with a sudden increase in their occurrence rate starting on the afternoon May 19. The evolution of the numbers of these low-magnitude earthquakes in various ranges of statistical significance within the source volume demonstrates some features, which can be related to the approaching initiation of the Sea of Okhotsk earthquake.

 

Key words: Sea of Okhotsk earthquake 2013, earthquake prediction, IDC, IMS, waveform cross correlation

 

Introduction

The Sea of Okhotsk is a well-known area of deep earthquakes. The subducting Pacific plate reaches here the boundary between the upper and lower mantle with the deepest earthquakes located at ~700 km as reported by the International Data Centre (IDC) using data from seismic stations of the International Monitoring System (IMS). The deepest earthquakes are of higher interest not only because of their seismic features but also they provide valuable data for other geophysical fields. The evolution of seismicity at the depths of 400 to 700 km may differ from that near the surface as related to temperature, lithostatic pressure, rheology,  chemical and mineralogical structure of the plate relative to the surrounding mantle, etc.

            The IMS provides continuous seismic data to the IDC since 2001. It includes seismic arrays which are characterized by much higher sensitivity and resolution than collocated three-component (3-C) stations. The usage of arrays allows to reduce the detection threshold worldwide and to detect earthquakes with much lower magnitudes than with the 3-C networks. Dense regional seismic networks of 3-C stations can detect events with lower magnitudes, but within the remote regions such as the Sea of Okhotsk there are no possibility to install enough bottom seismic stations just for the deep earthquakes study. With the array stations, the IMS is the best instrument to observe the evolution of seismicity in remote regions of the earth and the IDC provides a quality list of reliable event hypotheses as defined by the adapted version of Event Definition Criteria [Coyne et al., 2012]. The interactive analysis serves as a quality filter for the Reviewed Event Bulletin (REB) - the IDC final internal product of data processing. The REB contains up to 40% of events added by analysts, which have no counterpart in the automatically generated standard event lists (SELs) of the IDC.

            The natural seismic process is characterized by repeatability of events in location and size. The magnitude distribution of repeating earthquakes is well described by recurrence curves introduced in [Gutenberg,  Richter, 1936,1945]. As the time sequence of earthquakes in a given area is characterized by close locations one can effectively use not only the arrival times of the signals from them to find new events, but also the full waveform [Schaff, Richards, 2025]. This allows for finding many more events with much lower magnitudes in the vicinity of well-defined historical sources [Israelsson, 1990; Joswig, 1990; Schaff, Richards, 2004;  Gibbons, Ringdal, 2004, 2006; Gibbons et al., 2007, 2011, 2017; Waldhauser, Schaff, 2008; Adushkin et al., 2015, 2017, Kitov, Sanina, 2022]. The matched filter is an optimal detector for signals similar in shape [Turin, 1960] and it uses waveform cross correlation (WCC) to maximize the signal-to-noise ratio. When applied to deep earthquakes, the WCC-based detection and phase association techniques may recover the seismic events missed by standard IDC processing of the IMS data.

            The main objective of this study is to recover the low-magnitude seismic events during eleven days before the May 24, 2013 earthquake in the Sea of Okhotsk and three days after it. According to the IDC, this earthquake occurred at 05:44:49.7 on May 24, 2013, had coordinates 54.89°N,153.31°E, m_b=6.27, and depth of 604 km. The USGS moment magnitude is 8.3. There was no events detected by the IDC in the surrounding area 53°N-57°N, 151°E-155°E since February 6, 2012. The same seismic data from the stations of the IMS is used together with detection and phase association methods based on WCC. The recovered events are studied additionally in various magnitude ranges in order to estimate the evolution of the size distribution just before the mainshock. The sequence of mechanical effects before the mainshock may be revealed as the evolution from micro-cracking to the final rupture runs over increasing sizes of fracturing.  

 

Data and method

Data sources and general features

The IMS seismic network consists of arrays and 3-C stations. There are IMS primary stations defining the statistical significance of the created event hypotheses matching the EDC and auxiliary stations assisting in location and magnitude estimation of the valid hypotheses [Coyne et al., 2012]. This division is prescribed by the Comprehensive Nuclear-Test-Ban Treaty (CTBT) as providing uniform distribution of detection threshold worldwide. The civilian applications of the IMS data are not constrained by this rule and all stations are treated according to their performance in the  statistical sense. The IMS data is available via virtual Data Exploitation Centre (vDEC) of the CTBT Organization (CTBTO) [vDEC, 2026].

            As only the IMS data is used in this study, the REB is a natural source of events and associated detections used for comparison with the result of WCC processing. For the purposes of WCC processing, the REB provides event locations (epicenters and depths), magnitudes m_b and M_s, origin times as well as arrival times at IMS stations with a 0.001 s precision, and signal-to-noise ratio (SNR) for the associated signals. Selected properties of the deep seismicity in the Sea of Okhotsk region are also obtained from the REB. Figure 1 shows the recurrence curve for the region 42°N-65°N, 140°E-165°E and depth>410 km obtained from the REB. This curve uses the data from 2001 to the end of 2025 with the total number of approximately 1000 events. The corner seismic magnitude is in the range between 3.0 and 3.4. The regression line is extended beyond the corner magnitude to demonstrate the increasing deficit of low-magnitude earthquakes, which can be partially recovered by WCC.

            For the purposes of WCC processing, the signals associated with a set of REB events are used as waveform templates. At 3-C stations, the templates may have one- or three components depending on seismic phase and the station-event distance. For the IMS arrays, multichannel templates are used, with only vertical channels involved. Cross correlation coefficient, CC,  is a common measure of the angle between two data vectors representing digital time series of a template and a sought signal.  As a measure of similarity, this angle is calculated by a textbook dot product of two data vectors. The CC is not the best variable  for detection purposes, however, as it depends on signal spectrum and duration. In the presence of ambient seismic noise, both the template and the sought signal may differ significantly even being identical in their unique source.

            The search for a similar signal with a high-quality template, i.e. using a signal with large SNR, meets the noise as an important obstacle. By definition, the matched filter is an optimal detector maximizing SNR only in presence of stochastic and additive noise [Turin, 1960]. On average, such theoretical noise has the minimum CC value with any regular signal.  In practice, there is no such actual noise realization, which could be considered as stochastic and allowing for perfect detection conditions. Any ambient noise realization can have a component that is coherent with any real signal, from very low similarities up to a fully identical case. The latter case is related to the search of a signal from an aftershock of a catastrophic earthquake, as a multitude of other aftershocks in the same area simultaneously generate almost identical signals. As a result, the immediate aftershocks of large-magnitude earthquakes are difficult to find.

           

Figure 1. Recurrence curve for the region 42°N-65°N, 140°E-165°E and depth>410 km. The trend line extended below the corner magnitude of 3.0 to 3.4 indicates the number of potentially missed events. A portion of them can be recovered from the IMS data using WCC.

 

            Deep earthquakes are different in the post-seismic activity form shallow events. There are only 33 aftershocks of the M8.3 Sea of Okhotsk earthquake observed during the first 60 hours as Figure 2 shows. The decay of event magnitude with time during the initial 4 hours is a common feature of the major earthquakes. There is an effect of elevated noise as related to the mainshock wavefield, which likely prohibits to find aftershocks with lower magnitudes. The time gap between the mainshock and the first aftershock is approximately 18 minutes. It is longer than for the Tohoku March 11, 2011, and the Kamchatka July 29, 2025 earthquakes with moment magnitudes closer to 9.0. The second aftershock was found ~11 minutes after the first one. As well known from practice, deep earthquakes generate short impulsive signals not usual for near-surface sources except explosions. This effect may impact the performance of WCC detector. The signals from deep events are very similar to each other and are characterized by an almost vertical incidence angle at the closest station. In that sense, they are not coherent to signals from the shallow sources, which are by far the main contributors to the ambient noise. 

 

Figure 2. Magnitudes (black dots) and depths (red dots) of the mainshock and 33 aftershocks as a function of the elapsed time.

 

            In this retrospective analysis, the aftershocks of the May 24, 2013 earthquake can serve as master events (MEs) with corresponding templates at IMS stations together with other REB events found before and after the mainshock. Figure 3 demonstrates the distribution of 70 selected MEs and their catalog is listed in Appendix 1. The biggest earthquakes are less effective in detection of the smallest ones as their signals differ in shapes as related to their spectral content. However, after band-pass filtering the templates from the m_b>5.5 events are useful as they have very large SNR compared to the low-magnitude MEs.

            Altogether, the IMS data set and the IDC bulletin provide a set of waveforms, events and detections/templates for an analysis of deep seismicity in the remote Sea of Okhotsk region. The IMS array stations allow to reduce the detection threshold in standard IDC processing by a factor of 3 to 5 [Schweitzer et al., 2012]. The WCC-based processing allows for another 3 to 5 times reduction in detection threshold and finding of 60% to 100% of new events matching the REB requirements for a valid hypothesis [Bobrov et al., 2014; Bobrov et al., 2016ab; Bobrov et al., 2017; Kitov, Rozhkov, 2024; Kitov, Sanina, 2025b].

 

WCC-based detection

In accordance with the matched filter approach, the WCC-based detection procedure uses the signal-to-noise ratio technique  applied to the CC time series. Not to be confused with the standard SNR, this ratio is known as SNRcc. A set of CC time series is calculated using various template lengths or cross correlation window widths (CCWW). Waveform filtering is a standard procedure, and at the IDC it is realized in routine processing as a comb of band-pass filters [Coyne et al., 2012] covering the range of the expected signals related to the CTBTO monitoring mandate. The same procedure is applied to the IMS data before the CC  is calculated in a given CCWW [Bobrov et al., 2014]. The number of CC time series is equal to the product of the numbers of CCW and filters.

            The CCWWs and the comb of band filters are both varied to make the WCC procedure reliable as there is no a priory information on the sought signal length and spectrum above the ambient noise. It is  mandatory to cover the entire spectrum of interest for the first P-phases used in the current WCC version adapted to low-magnitude events: from 0.75 Hz to 8 Hz. The potential signal lengths are between a few seconds for impulsive P-phases and up to 100-120 seconds for regional sequence of signals from P_g/P_n to L_g. The length step depends on source type (shallow and deep earthquakes, volcano eruption, landslide, explosions, quarry blasts, other artificial sources), station-source distance, and properties of the propagation path.

 

Figure 3. Distribution of 70 MEs. For more information see Appendix 1.

 

            For the Sea of Okhotsk mainshock and its aftershocks, Figure 4 presents examples of templates at IMS stations PETK (epicentral distance ~3°) and MKAR (~45°). At PETK, the mainshock (the upper trace) creates a long wave-train (the original record is clipped) related to the source size and propagation path. Several aftershocks of various magnitude (see scales for relative sizes) also demonstrate long signals filtered between 0.6 Hz and 4.5 Hz. The CCWW set is adapted to the observed wave-pattern. At MKAR, the signals from the mainshock (and likely its immediate aftershocks hidden in the self-generated noise) are also lengthy, but the aftershocks, which are different from those at PETK, are characterized by impulsive and high-frequency first arrivals (3.0-6.0 Hz band filter is applied). Since our goal is to find the weakest possible signals like the aftershocks in Figure 4, the best at the level of the ambient noise, the low-magnitude aftershocks of the Sea of Okhotsk earthquake could be the optimal choice to be used as MEs. The initial segment of the mainshock signal is not similar to the sought signals and would be appropriate to find only the biggest events in the area.

 

a)

b)

Figure 4. Examples of waveform templates at stations a) PETK (filtered in 0.6-4.5 Hz band), b) MKAR (filtered in 3.0-6.0 Hz band). The scales on the left show their relative amplitudes.

 

            The detection thresholds at IMS stations have to be tuned to optimal fluxes of valid and false arrivals [Saragiotis, Kitov, 2020]. Statistical significance of the created event hypotheses depends on the probability of associated phases not to be randomly generated but rather to be related to a set of seismic sources to find. Therefore, the frequency distribution of SNRcc values determines the reliability of the created event hypotheses, their statistical significance and power. A detection list consists of the signals having SNRcc above the predefined threshold, but also matching a few additional quality parameters, with the main being the standard SNR threshold. By varying the SNRcc threshold one can control the detection flux. A reasonable value is 60 detections per hour with possible scattering ±30 for different MEs and time intervals. For quiet days, as defined by the absence of REB events in the studied area, the most of these detections are likely to be false in terms of relation to the used MEs. In some cases, such detection can be the result of the side sensitivity of the WCC detector to very intensive sources like the Tohoku earthquake. A spontaneous above the threshold SNRcc is also possible due to data problems like zeroes or spikes. The quality check is able to suppress such cases by the requirement of a tight sequence of SNRcc values above the threshold. There is also a possibility to use adaptive thresholds in order to adjust to the current valid events flux and noise level. Such adaptive thresholds are useful in tune WCC routine processing to low and high seismicity periods.

            A quasi-random SNRcc distribution for a given station can be generated within the WCC processing when the threshold is very high, say 500. In this case, detection fails and all SNRcc values represent noise. For an hour processing interval, the number of SNRcc values should be 3600 times the sampling rate. The IMS stations have sampling rates from 20 Hz (e.g., ASAR, BRTR, CMAR, ILAR, KSRS) to 100 Hz  (FINES, but it was 40 Hz in 2013) with a majority of IMS seismic stations working at 40 samples per second.  For a random SNRcc generation case during quiet time periods without detected signals similar to a given template, the frequency distribution is expected to be quasi-exponential. Any deviation from this basic case would reveal the presence of valid detections as the template reacts to the presence of similar signals in data. The SNRcc level where this deviation starts may be used as the detection threshold.  

            Figure 5 shows several SNRcc frequency distributions for IMS stations PETK and MKAR. Station PETK is the closest to the mainshock, and station MKAR is the second most sensitive to deep events in the Sea of Okhotsk as per the REB. Four MEs were selected from the top of template quality list as defined by SNR values of corresponding template signals. Two detection thresholds are used: 3.1 and 500. There are three cases appropriate for the comparison. Firstly, a 6-hour interval on jdate=2013133 (May 13, 2013) with a zero total number of events in the cross-correlation standard event list (XSEL) - an analog of the final automatic standard event list of the IDC (SEL3), as obtained in this study, represents the most quiet time period. Secondly, a 6-hour interval after the start of visible increase in seismic activity found by WCC on 2013139 (May 19, 2013). There were no REB events found in this interval, but 45 new XSEL events were detected. The third time interval starts with the mainshock and also lasts 6 hours. There were 31 aftershock with only 20 (65%) of them having seeds in the SEL3. Other (pure) REB events are added by IDC analysts using different sources of information and procedures. The XSEL includes 51 new events, and matches 26 REB events, including 18 with seeds from SEL3 and 8 from the pure REB subset.

 

Figure 5. Frequency distribution (0.1 bin width) of SNRcc at PETK and MKAR for three days. Four MEs are shown with two different detection thresholds: 3.1 and 500. The difference between the quiet (133) day and the days with only hidden (139) and IDC detected (144) activity is significant. The patterns also differ between stations.

 

            The earthquake preparation process has to include the growth in the sizes of the sources emitting the observed signals. At the same time, the noise at the recording stations from the increasing seismic emission includes the very same signals from all sources in the zone of the mainshock preparation. When beamforming is used as a detection technique [Schweitzer et al., 2012], this noise may mask the smaller events at array stations as the suppression is based on the assumption of destructive interference. The same effect is observed for the WCC-based methods,  but there is a possibility to improve the array resolution by adding numerically generated random time series to real data [Adushkin et al., 2025]. In this study, we applied this method and found several aftershocks missed by the IDC in the first 30 minutes after the mainshock.  

            When interpreting the curves in Figure 5, it is important to understand what does mean a valid signal with SNRcc above the predefined detection threshold. A sought signal can  generate a long sequence of relatively large SNRcc values even below the threshold (e.g., 500). When correlated, a template and a sought signal create a sequence of CC values higher than those obtained by correlation with the pre-signal noise with a length compared to the CCWW. This sequence of CC values transforms into a sequence of SNRcc values above those created by the noise. This is similar to standard STA/LTA detector where a signal above noise generates a sequence of elevated SNR values.

            When a signal is detected, the LTA has to be frozen in order not to add the signal energy to the SNR or SNRcc. This makes a longer sequence of higher SNR values. For the WCC detector, the effect of frozen LTA is approximately the CCWW. In any case, the presence of sought signals create longer sequences of higher SNRcc values, as one can observe in Figure 5. However, there is only one of all the SNRcc in a given sequence to be used to characterize the detected signal, which is called "arrival" after it passes quality checks. For the  arrivals, many parameters are calculated such as arrival time, amplitude, period, SNRcc, SNR, and CC. When detection failed, the sequence of higher SNR and SNRcc is shorter and their values are lower as the LTA is not frozen. To tune a detection threshold a special procedure has to be implemented [Saragiotis, Kitov, 2020], which balances the detection rate and the portion of the generated arrivals associated with valid event hypotheses in the final bulletin (SEL or REB).

            In this study, the detection thresholds at individual IMS stations were estimated to retain the total number of arrivals below the predefined tolerance level. They were tuned in a series of test calculations and then are fixed thought the calculations. This makes the following phase association procedure to be unbiased in terms of the arrival lists quality over the studied period. The MEs are selected from the REB and have different numbers of associated stations and quality of corresponding templates. There are places where only relatively poor REB events with just a few associated stations are available.  Therefore, the detection rate at individual stations and the total number of detections produced by a given ME vary in a wide range. This effect may create some disproportions in the MEs input to the final XSEL. This is a normal situation for an automatic bulletin. For example, SEL3 has many seed events at ranges of 90° to 180° from the final REB solutions.  An automatic bulletin has to be approved in an interactive analysis or compared to a high quality bulletin for the same area. Fortunately, the REB can be used as a reference for the XSEL for the studied region. A direct comparison of arrivals associated with the events in both bulletins is possible at the same IMS stations by their arrival times. Such a comparison needs only 1 common arrival within the allowed retiming window (±10s) to match an XSEL and REB event. Such a procedure is often used by the IDC for comparison with alternative solutions.

 

 

Local association

            The WCC arrival list consists of individual lists for each of the 70 MEs (see Figure 3). These individual lists do not intersect but may contain arrivals corresponding to the same physical signals. By design, all templates of a given ME are able to find only signals from other events within a small footprint around its epicenter. For the sought events within 1/4 of the signal wavelength, CC can be between 0.95 and 0.99 if the noise is negligible [Arrowsmith  et al., 2006; Baisch et al., 2008]. Such a situation was observed for the signals from the DPRK announced underground tests [Selby, 2010;  Kitov et al., 2025]. For very high CC and SNRcc values, the arrival time precision for the sought signal can be as low as 0.001 to 0.005 s. As a result, the RMS travel time residuals for the sought event located using the WCC arrivals can be of the same magnitude of 0.005 to 0.01 s. In terms of spatial location accuracy, it corresponds to 40-80 meters for the apparent wave velocity of 8 km/s, e.g., P_n-phase. This is the accuracy of relative location and the absolute location accuracy of the sought event is practically the same as for the corresponding ME. The relative location using all the associated phases (stations) and confined to a small area around a ME is called Local Association (LA). The LA version used in this study is well described in [Bobrov et al., 2014; Bobrov et al., 2016ab; Kitov, Sanina, 2025b].

            The CC and SNRcc values rapidly fall with distance from a given ME and with the amplitude of the sought signal relative to the noise level. For an array station, a virtual event location moved over the LA grid may change the scalar slowness and the azimuth in the arriving plane wave, which can improve or suppress the multi-channel cross correlation with a template. The change in actual event location may also affect the source function of the sought event and the difference in propagation paths may disturb the sought signal shape.

            The search area for a given ME is defined by the degradation of correlation with distance and by the presence of different MEs. For the MEs in Figure 3, the search area of 40 to 90 km in radius is appropriate as the MEs also have confidence ellipses of from one digit to tens of km and some intersection of the search areas would be helpful. For such an area, a grid search is a natural choice as any position within the search radius can be tested for the sought event location without extensive calculations. For example, a grid with 35 nodes and 2 km step in both directions from a ME can be used. The depth of the sought event is fixed to that of the ME as the depth sensitivity is very low and thus increases uncertainty. To accommodate the depth dependence into the WCC processing, MEs in 40 km depth intervals are independently selected where available. 

            The grid search for valid event hypotheses is based on the origin times for virtual solutions related to the nodes of the grid. A virtual event origin time in the (0,0) node can be calculated as the WCC arrival times less the empirical ME-station travel times as obtained from the REB. For the other nodes, the travel times have to be corrected for the dot product of the vector to the node and the vector slowness of the signal at the related station. This is a very accurate time correction as the node vector is a predefined value and the slowness vector for a given station does not change much between the nodes of a small grid search area compared to the ME-station distances.

            For a valid event hypothesis, the origin times for three or more associated arrivals (stations) have to be within the predefined tolerance. At the IDC, the allowed travel time residual relative to the theoretical time is 2 s for the first P-phases, and  2.2 s for the P_g- and P_n- waves. These residuals correspond to the desired statistical significance of the REB events as adapted in the EDC. At least three first P phases at different stations have to be associated with a valid event. These rules are used in the WCC processing to match the IDC approach [Coyne et al., 2012]. As the XSEL arrivals potentially have higher relative arrival time accuracy, the origin time tolerance can be much lower than for the IDC without a significant loss in performance efficiency expressed in the number of valid (i.e., matching REB) events. For automatic IDC bulletins (e.g., SEL3), the allowed travel time tolerance is much larger than 2 s and can reach 5 to 7 s. The WCC-processing, despite being automatic, strictly follows the entire list of the REB requirements, except the division into the primary and auxiliary stations.

            Three or more arrivals with origin times within ±2 s interval from their average value are not a guarantee of a reliable event hypothesis. There are many arrivals at many stations and the event hypotheses have an additional degree of freedom - varying location. The REB provides valuable information on the IMS station performance for the historical events in any given seismic region. There are stations which contribute to almost all REB events is a specific area as Australian stations WRA and ASAR to the events within Australia. Thus, an event hypothesis within Australia without these stations, when they are operational, would be rather false despite other parameters may match the EDC. There are also IMS stations contributing only to the largest events in this region or not contributing at all. The difference between the most sensitive and the ineffective stations can be described by their participation rates in the REB events. Such rates are calculated for each IMS station for all regions with repeating seismic events including the Sea of Okhotsk. For example, the first P phases detected by PETK have the highest probability to be associated with the events in the studied area and can be characterized by a station weight (or probability) of 1.0.  Stations MKAR has a weight of ~0.9.  Other stations have weighs less than 0.85.        To create a statistically significant event hypothesis having three associated stations, with PETK and MKAR included, one can define the minimum event weight, i.e. the sum of all associated station weights, of 2.7.  Then there will be just a few valid three-station configurations.

            With the increasing number of associated stations, the requirement of total event weight plays lower and lower role as the probability of many random arrivals to be within a narrow origin time interval rapidly decreases. An important parameter for the statistical significance of  an event hypothesis is the requirement of large SNRcc at one or two best stations. According to Figure 5, random arrivals with SNRcc>5 are almost excluded on quiet days. A similar approach to the event quality using the arrival probabilities is the sum of individual SNRcc values. This parameter introduces a series of thresholds for 3- to 7-station event hypotheses. For example, for a 3-station event the sum of SNRcc is set to the level of 12.0. One of these stations has to have SNRcc above 5.0 and the other two have to add at least 7.0 in any combination. A 4-station event has a threshold of 15.5. With the increasing number of stations above 5, the threshold increases by a constant predefined step.

            All in all, the most important LA parameters are as follows: origin time tolerance, total event weight, weight of the best station, number of associated stations, sum of SNRcc values of all associated arrivals, and the minimum SNRcc for one of the best stations. The grid radius and spacing between nodes are also important for the final XSEL. There are less influential parameters, which are retained unchanged through all calculations after tuned in the test calculations.

 

Conflict resolution

            There can be dozens or even thousands of MEs in WCC processing. The global coverage in routine processing includes more than 42,000 MEs. Each of the MEs generates an individual arrival list and an XSEL. Two or more neighboring MEs can find the same physical event, i.e. they can detect the same physical signals and use them to create similar hypotheses. But a unique physical signal must belong to one and only one XSEL event. A similar conflict of interests exists in any seismological agency including the IDC. To resolve such a conflict a special procedure is developed and applied [Coyne et al., 2012]. In automatic processing, Conflict Resolution between the event hypotheses competing for one arrival is based on an assumption that the event with a larger number of associated phases wins.  If this number is equal between the competing events, the event weight comes into play. After the CR procedure is finished, the arrival is associated with a stronger event hypothesis and all competing event hypotheses lose this arrival. Then, the losing hypotheses are checked for the EDC match and can be rejected if not meeting the requirements.

            The IDC event weight includes the input of associated phases. Each phase can add a specific value related to the residuals of arrival time, azimuth and slowness. These values depend on phase and station type (array or 3-C).  For a secondary phase to add to the event weight, it has to be the first P-wave associated at the same station. The sum of individual values related to all associated phases at primary IMS stations defines the event weight and the threshold weight defines the statistical significance of the event. This threshold is tuned to the level than only the event hypotheses with 3 or more primary stations with associated first P-phases can be promoted to the REB.

            The number of conflicts in the final XSEL creation is much higher as several MEs usually conflict for the same arrivals at many stations. Following the IDC approach and the event weight definition used in the WCC processing, all conflicts are resolved in favor of the event with the highest weight, then follows the number of associated phases. If these two parameters do not resolve the conflict, the RMS residual origin time is used. The final XSEL contains only one arrival per a unique physical signal.

            An event hypothesis losing one or more arrivals is checked against the EDC related to the WCC processing. If it does not match any of them it is rejected. The loss of an associated phase does not preclude this hypotheses to exist if it still matches the EDC. There are a few  cases of earthquake doublets close in time and space which are difficult to distinguish and potential mis-association of phases from different source with one out of the two.  The SEL3  gives some examples of  "joint" (arrivals from two different events are associated with a single event) and "split" (one event is split into two different events) events. Such events are mainly belong to the aftershock sequences of the largest earthquakes when secondary phases can be mis-associated and the arrival times are not accurate. These are the problems for the automatic WCC processing as well.

 

Measuring the growing activity with increasing event magnitudes before the mainshock

            For routine WCC processing, a permanent set of LA parameters is determined in a tuning procedure. This set has to generate the XSEL with the largest number of matched REB events and the lowest number of false event hypotheses. A few experiments of interactive review of the automatic XSEL according to the IDC rules made it possible to conclude that the REB misses 50% to 100% of REB-ready events. And these newly found events are weak, but fully comply with the EDC.  These smallest events are the principal objective of the CTBT compliance. The biggest events tell for themselves.

            The study of the Sea of Okhotsk earthquake is different from the principles of the IDC operation with well established thresholds, as the seismicity before the mainshock has to be found and characterized according to the mechanics of earthquake preparation process. The IDC has to limit the event magnitude by 2.0 as one can observe from Figure 1, where no events are observed below this threshold. There is no limit to the WCC processing resolution except the requirement that the XSEL should include valid events according to the  adapted EDC.

            The mantle convection drives the plate tectonics with the solid plates subducting under the continents to a depth of 700 km and likely deeper. The deformations and stresses in the plate play defining role in the near surface tectonics and seismicity. With increasing depth, the earthquake mechanism may change as the impulsive character of the observed signals indicates. There are various scenarios proposed and the correct one has to rely on actual observations of the seismic process it describes. One of the possibilities to look into the details of this seismic process is opened by the use of WCC-based processing.

            The evolution of seismic process before major seismic events should be the same as in laboratory experiments from the point of view of increasing sizes of smaller events before the mainshock.  For shallow  earthquakes, small cracks in the very beginning evolve through small faults to the level of main rupture. For  deep earthquakes, there are three major mechanisms proposed: transformational faulting, dehydration  embrittlement, and thermal runaway. Their combination in a unique event is also possible. Physical models of the preparation process need observations before such deep events in order to tune their defining parameters.

            In order to recover the seismic process before the Sea of Okhotsk earthquake, the WCC-based processing was split into several cases with different resolution and sensitivity. In the most sensitive case, the XSEL has to include the weakest event hypotheses allowing to find the earlier signs of the start of the earthquake preparation process. In the most strict case, the XSEL has to be tuned up to the level where it can find and match the most of the REB events in the aftershock process. This allows to prove the statistical power of the entire procedure. In this case, the number of XSEL events extra to the REB should be within the range of 50% to 100%, as observed in the previous interactive exercises. Below this strict XSEL level on the way to the weakest set of XSEL defining parameters, the same number or more REB events can be matched and the number of XSEL events may grow. The increasing number of XSEL events may include those which cannot be approved in an interactive review as the analysts have to have visual confirmation of the signals. Nevertheless, these XSEL events can be valid as radar applications without human interference prove. The signals deep below the noise level can be safely detected and used.

            Before the Sea of Okhotsk earthquake, seismic activity may progress up from the weakest events around magnitude 2.0 in Figure 1 to the events which can be potentially detected by the IDC. Therefore, various cases of WCC-processing have to focus on specific parts of the recurrence curve from below 2.0 to 3.0-3.4 in order to detect and describe the pre-seismic process evolution. There is no reliable quantitative information on the uncertainty of the magnitude estimates for the WCC detections near the noise. One can only judge on the sensitivity of a given WCC processing scenario by the number of XSEL events built in the predefined time intervals. The full set of the WCC scenarios is based on two limit versions - strict and weak. The strict version has the search radius of 43 km, the best station minimum weight of 0.855, total event weight of 2.5, the SNRcc sum thresholds of 15.0, 18.5, 22.0 for 3-, 4-, and 5-station events, respectively, and the SNRcc threshold of 5.0 for one of the best (with weight > 0.855) associated stations. The weak version has the search radius of 80 km, the best station minimum weight of 0.8, total event weight of 1.8, the SNRcc sum thresholds of 14.0, 17.5, 21.0, and the best station SNRcc threshold of 4.5. All other parameters are the same. The statistical significance of XSEL events in each version is additionally controlled by six values of the origin time tolerance, Δt: 0.25, 0.5, 1.0, 2.0, 3.0, and 5.0 seconds. There are 12 version/case pairs: v₁c₁, ...v₁c₆, v₂c₁,...,v₂c₆. The narrowest origin time window allows practically only actual event hypotheses to be created as the probability of  3 or more random detections at the best stations to occur in a 0.5 s time window is extremely low considering the limit on the detection rate of 60±30 per hour.  On the other hand, the narrowest window may potentially miss many valid XSEL events consisting of weaker signals with poorer arrival time estimates as related by the noise influence. The v₁c₁ pair should not be able to detect too many new XSEL events extra to the REB. Overall, the XSEL events matching the REB events have to be of a higher quality related to the SNR values of the associated signals.

            The weak WCC version with the shortest origin time window, v₂c₁,  has to be focused on the new XSEL events with the best quality just below the corner magnitude of the recurrence curve. These events are of high statistical significance and may include weak WCC arrivals. The v₂c₆ pair has the highest resolution below magnitude 2.0. It may contain many valid events and likely some false events, with both types formally matching the EDC. The strict and the weak versions with the same Δt are focused on a specific event quality range, which depends on event magnitude and mechanism as well as on the noise level at stations around the arrival times of all associated phases. Formally, for pairs v_ic_j, we introduce quality boundaries Bⁱ_j, where i=1,2 is the version - weak or strict, and j=1,...,6 is the tolerance case. By design, B¹_j ≤ B²_j and Bⁱ_j+1 ≤ Bⁱ_j. The quantitative meaning of these boundaries is revealed by WCC calculations with an increasing number of event hypotheses generated with increasing j for a given i. When the sought seismic activity from the quiet state reaches the lower boundary Bⁱ_j  the number of weak events grow and the number of strict events, which must have larger magnitudes and statistical significance, stays still. The change in the ratio of the event numbers in the weak and strict XSEL versions, B²_j/B¹_j, expresses the change of seismic process in this range. When the activity starts to grow from a quiet state, the B²_j/B¹_j ratio also increases with the growing number of XSEL events above B²_j. The activity may not reach its upper boundary B¹_j and start to fall to the initial level. In some cases, the activity may reach the upper boundary and move above it as, for example,  happens after larger earthquakes. Then the ratio may fall because of large number of events above the upper boundary and constant or decreasing number above the lower boundary. The larger earthquakes may generate noise masking the smaller events. In an aftershock sequence of the largest earthquakes, the  B²_j/B¹_j ratios in the magnitude ranges much below the corner magnitude may reach the lowermost levels and express the distribution of earthquakes in the recurrence curve. With the decaying post-seismic activity,  the ratio grows again as the number of events above B¹_j has to decrease fast and the number above B²_j decreases much slower. All in all, there are six quality, potentially related to magnitude,  ranges to follow the seismic process before the earthquake.

 

 

Results

Statistical power. Aftershocks and new XSEL events

            The mainshock at 5:44:50 on May 24, 2013 is a reference for all calculations. Initially, the recovered seismicity was analyzed at a daily rate. The processing period was 25 hours per day, with 1 hour in the next day serving to find the events in the processed data day but with associated arrivals in the next day. Such an overlap is used to make processing continuous. The revealed evolution of seismic activity made it more appropriate to present data in 6-hour intervals in order to better determine the start of the change in the rate of the event occurrence and the evolution of the aftershock sequence since its principal phase lasted only 12 hours.

            There was no aftershock before 6:02:45, i.e. ~18 minutes after the mainshock. Such a still period is observed after the major earthquakes and can be related to the absence of any large magnitude post-seismic activity or to a much higher noise level hiding the immediate aftershocks. The latter problem can be resolved with the coherent noise suppression based on the added stochastic noise [Kitov et al., 2026]. In total, there were 40 aftershocks between 06:00 and 18:00 on May 24. Then, the postseismic activity stopped. From these 40 aftershocks, 27 had SEL3 seeds and 13 (33%) were built by IDC analysts.  Table 1 lists the numbers of events in the XSELs for all 12 versions/cases used in WCC processing for two 6-hour intervals within the 12-hour period between 6:00 and 18:00. The best strict WCC version (Δt=3.0 s ) finds almost all (19 out of 20 for the interval 6:00-12:00, and 7 out of 7 for the 12:00-18:00) REB events with SEL3 seeds. For the pure REB events, the best result is in the Δt=3.0 s case - 8 out of 11 for the first interval and 1 of 2  in the second. In total, there were 27 out of 31 REB events matched in the  Δt=3.0 s case of the strict version.

            For the v₁c₁, the XSEL has matched 25 REB events, i.e. practically the same proportion of the total as in the best case v₁c₅.  The number of new XSEL events, not matching the REB, in the 12-hour segment increases with Δt from 23 for 0.25 s to 74 for 5.0 s. For the v₁c₁, the number of new XSEL event hypotheses corresponds to the experience with the XSEL interactive analysis. These are the events that match the IDC rules and could be confirmed in interactive analysis. The v₁c₆ contains all the new events found by the v₁c₁, but they have more associated stations as the origin time residuals between 0.25 s and 5.0 s cannot be included in the latter case. Statistical significance of the v₁c₆ may suffer the probability fall related to the larger tolerance, but this fall is partially or fully compensated by the number of associated arrivals. The number of added stations is from 0 to 4. When no stations are added, the XSEL solutions for the v₁c₁ and v₁c₆ are identical. The importance of additional arrivals is in the change of event location and magnitude.

            There are many additional to the v₁c₁ events created by v₁c₆. Formally, they all match the EDC and the quality criteria of WCC processing. Most of them have to be valid as their statistical significance similar or even larger (e.g., larger event weight, more associated stations, lower RMS origin time residual) to the XSEL events matching the REB events on May 24. The overall WCC performance for the Sea of Okhotsk aftershock sequence demonstrates high statistical power of the XSEL events confirmed by the REB match. This conclusion on the statistical power of the XSEL events is important for the analysis of the events missed in the REB. The XSEL events have the same statistical significance through the whole study, and thus, the same statistical power as that during the aftershock period.

            A special question is the quality of the REB events not matched in the XSEL. The REB may contain some poor solutions often related to the higher noise level after the mainshock forcing phase mis-timing and mis-interpretation/association [Kitov, Sanina, 2025b]. Appendix 2 lists bulletins (only P-phases) for those six REB events which were not matched in the v₁c₁. The SNR values of associated phases are predominantly below the IDC detection threshold. These arrivals were added manually. Low SNR values indicate poor quality of these arrivals - the reason they were not detected in automatic processing. The arrivals with SNR between 1.0 and 2.0 in Appendix 2 are rather a challenge for visual identification. Two out of these six not-matched events consist only of such poor arrivals, one event has just one arrival above the automatic detection threshold. Five out of the six not-matched events occurred during the hour of the most intensive aftershock activity when the detection rate was the highest due to many aftershocks and very high number of secondary phases coherent to valid first P phases. These events are poor and are likely unreliable.  Therefore, the XSEL may fail to create similar hypotheses. An extended study is needed with more powerful detection method based on WCC with noise stochatization [Adushkin et al., 2025].

             

 Table 1. Statistics of WCC performance between 6:00 and 18:00, May 24, 2013.

Δt - origin time tolerance; Start, End - start and end hour of 6-hour intervals; SEL - number of REB events with SEL3 seeds; REB - number of pure REB events, SEL/Total - share of REB events with SEL3 seeds in the REB. For Strict and Weak, SEL and REB - number of matched REB(SEL3) and pure REB. Total REB - total number of matched REB.  Total XSEL - total number of XSEL events. NEW - number of XSEL events new to the REB.

 

 

Finding the aftershocks in the still using noise suppression

            The still period between 5:45 and 6:03 was processed together with a noise suppression technique based on a stochastic noise component of varying amplitude added to the filtered waveforms. This technique has shown high efficiency and allowed to find several aftershocks in the still period after the Kamchatka July 29, 2025, earthquake [Kitov et al., 2026]. The WCC processing is essentially the same except the stochastic component destroying the template coherent noise to a larger extent than the sought signal. The ambient noise after the major earthquakes mainly consists of the signals from the aftershocks, i.e. the events with lower than the mainshock magnitudes but similar source functions and propagation paths to IMS stations. This component of the ambient noise is well correlated with any template from the aftershock area and thus generates higher CC values than the usual seismic noise consisting of signals from a multitude of sources with various magnitudes and at different distances, not mentioning the wind and human activity as noise sources.

            These higher CC values generate higher LTA at the CC-traces. Therefore, the sought signal has to have higher amplitudes in order to generate elevated CC values that can be detected at the same threshold level. An SNRcc value of 2.5 means that the STA is 2.5 times larger than LTA. By adding a stochastic noise component equal to the actual LTA one can reduce the CC level and thus the  new LTA by a factor larger than 2.0, as numerous test calculations show [Adushkin et al., 2025; Kitov et al., 2026]. For the sought signal, the STA for the new CC-trace may lose just a small portion of around 20%. Then a conservative estimate of the new SNRcc will be (STA*0.8)/(LTA*0.5) = 4.0. The SNRcc gain is 60%. As the amplitude of the coherent noise component is not known in advance, by varying the stochastic noise amplitude in a wide range one can reach the almost maximum possible gain in SNRcc. This procedure needs extensive calculations and is most effective in finding the weakest signals in the noise coherent to the template and the sought signal. The basic version of WCC processing is efficient enough for different cases when hidden seismicity has to be recovered in a usual noise realizations.  

            There are many ways to add stochastic noise to real data. The simplest is to scale it to the maximum amplitude of the signal in a given time interval, which can be from several minutes to a day. In this study, one hour intervals are used with a 20-minute overlap with the next interval, which covers the travel times for all first P-wave phases. As the sought signal amplitudes can span over a few orders of magnitude, the number of various stochastic noise factors, StN, covering the whole range can be very large. The stochastic noise amplitude scaling might be more flexible with standard STA and LTA playing the role of local original noise amplitude to suppress. The standard STA and LTA time series are calculated first and then used stochastic noise scaling. For the LTA case, the product of a LTA estimate, a random number between 0.0 and 1.0 and StN is calculated for each time count and then added to the original trace for further WCC processing. For the STA case, a similar procedure is applied to each time count where standard STA/LTA is below a predefined threshold, which is not the detection  threshold but a tuned parameter. When the STA/LTA ratio reaches the threshold, the LTA value for this count is used  instead of STA for generation of the stochastic noise component. The length of fixed LTA window is defined by the STA/(fixed LTA) above the threshold rule.

            The simplest approach with scaling of the maximum amplitude in the interval is used. The range of  StN is from 0 to 20. The first value corresponds to the basic WCC progressing as an alternative to the noise suppression technique. The WCC processing is applied to each StN value in the range separately and the maximum SNRcc value is taken for each count. This is similar to standard detection procedure [Coyne et al., 2012]. Two 1-hour time intervals between 5:00 and 7:00 were processed at all involved IMS stations and the results for the pre-mainshock and the after-mainshock still periods are presented in Table 2.  No events are in the REB except the mainshock which is excluded from the list. All 12 version/case pairs were used, but Table 2 presents only the limit pairs StN_v₁c₁, StN_v₁c₆, StN_v₂c₁, StN_v₂c₆. As a reference, the basic pair 0_v₂c₂ is listed. There were 13 XSEL events found by the StN_v₂c₆ pair before the mainshock with magnitudes from 2.64 (±0.46) to 3.45 (±0.47).  The number of associated stations, Nass, varies from 5 to 12. The other three version/case pairs found from zero (StN_v₁c₁) to 4 (StN_v₁c₆) out of the 13 events. The 0_v₂c₆ pair detected 5 out of 13 events, but also added one event hypothesis not found by the StN_v₂c₆. This 0_v₂c₆ event (highlighted bold) at 5:10:43 was very close to the 5:11:04 event in the v₂c₆ and could be missed due to the CR rules. These numbers demonstrate that the technique using stochastic time series to suppress the ambient noise can find 2 to 3 times more events than the basic WCC processing despite the latter is able to find many events not detected in standard processing.

           

Table 2. Event hypotheses for the period 05:00-06:02, May 24, 2013. The still period is between 5:44:50 and 6.02:45.

 

In the 18-minute interval between the mainshock and the first aftershock, the StN_v₂c₆ pair found 3 events and the 0_v₂c₆ pair found an event different from these 3 events at 5:59:01 (highlighted bold), which is again close to the 5:58:41 event. The first aftershocks at 5:56:20 is confirmed by all pairs but with different number of associated stations and m_b values. There is still a gap between 5:45 and 5:56 without aftershocks. Unlike the Kamchatka July 29, 2025, earthquake, the stochastic noise technique did not demonstrate a better performance than the basic WCC processing during the high-noise period. The source functions of the deep aftershocks might be too impulsive and the noise created by the such events in a relatively small volume is rather a sequence of impulses with varying amplitude at array stations where all secondary phases with higher slowness are effectively suppressed if also impulsive.  In any case, further analysis is needed with various scaling of stochastic noise.

 

The evolution of seismic activity since May 19 recovered by WCC

            The main results of the WCC processing for two weeks between May 13 and May 26, 2013 are shown in Figure 6. The weak (a)  and strict (b) versions are presented separately, but each version includes all six origin time tolerance cases. There was revealed a significant rise in the number of found XSEL events for all twelve v_ic_j pairs. The ratio of the number of new XSEL strict events during the fourth quarter of 2013139 and the average number in the twenty six 6-hour intervals between the Q1/2013133 and the Q4/2013139 varies from 6 (52/9) for v₁c₆ to 37 (37/1) for v₁c₁. The latter value is especially important as these 37 new XSEL events in a narrow 6-hour interval are likely the most statistically significant among all cases. For the weak version, the same ratios decreased to 4 and 25 for v₁c₆ and v₁c₁, respectively. As a function Δt, these ratios practically follow a power law with the exponent of -0.67 (R²=0.97) for the strict version and -0.59 (R²=0.97) for the weak version. By all means, the spike from 1 to 37 per 6 hours is outstanding.  No of these events were found by the IDC, but one can consider the Q4/2013139 as a start of the mainshock preparation process. Then the behavior of the seismic activity in various event quality bins Bⁱ_j may reflect the evolution of the earthquake preparation process. The number of events in the 6-hour interval after the mainshock was increased to the level of 45 for v₁c₁ to 83 for  v₁c₆. For the weak version, ratio is 47/108. These events include the REB matched by the XSEL and new XSEL events.

 

a)

b)

Figure 6. Comparison of the number of detected events in 6-hour intervals using six different origin time tolerances and two versions of phase association parameters: a) strict; b) weak. Vertical red line - Q4/2013139, black line - mainshock.

 

            In Figure 7, the evolution of seismic activity between 2013133 and 2013146 is presented by a) depth and b) event weight for the v₁c₄ pair. The origin time residual of 2 s corresponds to the threshold adapted by the IDC for the REB and the strict version likely gives events more reliable from statistical point of view. The depths of XSEL events are equal to the depths of the ME winning the potential CR. Therefore,  the best MEs generate XSEL sequences with constant depths. Before Q4/2013139, the XSEL events In Figure 7a have lower depths on average and after the start of activity the events move to the depths where the most of aftershocks reside together with the mainshock. The event weights are lower in the low seismicity interval and grow significantly with the start of increased activity. On average, the XSEL events matching the REB events have larger weights but there are many events, mainly pure REB, which have the weights between 4.0 and 6.0 as the events in the pre-mainshock interval.

 

a)      

  

b)

c)

Figure 7. Depth a), event weight b), and the number of associated stations c) as a function of time for the v₁c₄ XSEL events. The start of seismic activity and the mainshock are presented by red and black vertical lines,  respectively. The XSEL event matching aftershocks in the REB are highlighted by red circles.

 

            Relative location of the XSEL events for v₁c₄ and v₂c₄ is shown in Figure 8. For the strict version the footprint is smaller than for the weak one as defined by the search radii. There are three different periods: Q1/2013133-Q3/2013139, Q4/2013149-Q1/2013144, and Q2/2013144-Q4/2013146.  Overall, the ME are to the south of the mainshock and there are not too many XSEL events above 55° N despite the search area covers this part of the observed seismicity and some immediate aftershocks also were located in this part of the plate. The events before the mainshock do not reveal any specific gravity points with very high density, but their density within the area occupied by the northern MEs is relatively low during the pre-mainshock period for both strict and weak versions. The pure XSEL events after the mainshock also dispersed over the search area, but those matching the REB are close to them and concentrate near the northern MEs and the mainshock.

           

 

Figure 8. Locations of the XSEL events for v₁c₄ and v₂c₄ pairs.

 

            Accurate location of the XSEL events is not the main objective because the low-magnitude events usually have large confidence ellipses and their absolute locations are not reliable. Statistical significance and power are essential for a reliable estimate of the seismicity evolution before the mainshock. For a more strict set of the LA defining parameters, including a smaller search area, one could obtain a tighter XSEL distribution with a lower number of events. The statistics of the event flux would change as a result, however. At the same time, even if to admit that a part of the XSEL events are false according to some strict rule, their presence is definitely the result of an elevated flux of signals similar to templates as Figure 5 demonstrates. In that sense, even the potentially false events, e.g., generated by valid first P-phases and some secondary phases created by reflections from the plate boundaries,   support the observation of the increasing seismic energy generation started in the third or fourth quarter of May 19.

            After the sharp peak in Q4/2013139 there were four peaks in the XSEL events number in 6h intervals with decreasing amplitude before the mainshock. The observed variations in relatively short time intervals partially mask the pattern of increasing/decreasing activity before the mainshock. Figure 9 shows the curves in Figure 6 smoothed by a simple moving sum of the four 6-hour intervals, MS(4), to retain the causality principle with no data from the future. A downward trend starts in the first quarter of 2013142, but for the weak version it is less significant, especially for larger Δt values.              Seismic activity in the strict version after 2013133 is characterized by 3 groups with close figures: for Δ of 3.0 and 5.0 s,  1.0 and 2.0 s, and 0.25 and 0.5 s. The gap between these groups is larger than the distance between them in their pairs. On 2013138 the curve Δ=0.5 s diverts to the second group and the  Δ=0.25 s  curve stays alone despite it also starts to grow. For the weak version, all curves almost all time are separate and do not group for long. This is an important observation since the set of Bⁱ_j  was constructed to reveal variations in the ratio of the events in the strict and weak versions. The B²_j/B¹_j ratio should reveal the upcoming activity in various quality ranges defined by the origin time tolerance. The larger the case j, the earlier the signs of activity should be observed.

           

a)

b)

Figure 9. Same as in Figure 6 smoothes by a moving sum MS(4) applied. The start time of the occurrence rate growth in shown by vertical red line.

 

            Figure 10 presents the weak-to-strict number ratio for all the six Δ values. Before Q3/2013139, these ratios are varying in a wide range and do not reveal any coherent behavior since the numbers of new XSEL event in the strict version are relatively low and demonstrate strong variation (see Figure 6).  In Q4/2013139, the curves start to evolve almost in sync and the numbers of new events in the strict and weal versions rise significantly. After approximately two days of low variations, the curves start to grow, with those with the larger Δ starting earlier than for  the Δ=0.25 s and Δ=0.5 s cases. The Δ=2.0 s and Δ=5.0 s curves reach their peak values first among all cases, with the curves Δ=3.0 s and Δ=1.0 s reaching their peaks immediately after. The Δ=0.5 s curve grows slowly to Q3/2013143 and then jumps to its peak on Q1/2013144. The  Δ=0.5 s curve is approximately constant to Q3/2013143 and then also increases by almost  a factor of 2 to its peak on Q1/2013144. Figure 11 shows this behavior in details with the only difference that not only new XSEL events are counted in, but the XSEL events matching the REB are also included. Therefore, these curves present "total" numbers, which are identical with "new" before the mainshock. The curves with Δ>1.0 s slowly fall to their minimum level in Q2/2013144 practically ignoring the mainshock. The other three curves also fall from Q1/2013144 to their respective minima in Q2/2013144. For the Δ=0.25 s case, this minimum is close to 1.0. When the IDC aftershock sequence stops, the curve start to diverge again.

 

Figure 10. The ratio of the number of new XSEL events in weak and strict versions in Figure 9. 

 

            The described behavior of the weak/strict ratios follows the expected scenario of the low-magnitude seismicity related to the earthquake preparation progressing from the smallest to the largest sizes or emitted seismic energies. The whole process of the Sea of Okhotsk earthquake preparation spans over 4.5 days. In takes 2.5 days of synchronized higher seismicity in all six event quality bins between Q4/2013139 and Q1/2013142 (see Figure 9). Then the seismic activity  reaches its peaks in the lower quality bins without significant changes in the higher quality bins with Δ=0.25 s and Δ=0.5 s. In the final preparation stage, seismic activity reaches the top quality bin signaling the readiness to the mainshock initiation.

           

Figure 11. The ratio of the total (new and REB matching) number of XSEL events in weak and strict versions for the period between 2013139 and 2013146.

 

            This qualitative description of the observed seismicity evolution is one of possible explanations. A comprehensive WCC processing over the entire dataset since 2001 has to be fulfilled to prove that this example is a case specific and there are no other similar seismic activity sessions without catastrophic earthquakes. Alternatively, other major earthquakes can be analyzed by the same WCC procedure in order to evaluate the evolution of pre-event seismicity in the same or different quality bins. The Sea of Okhotsk earthquake is a unique one. The source volume is isolated with no seismic events detected by the IDC foe more than one year. This makes any seismic activity within this volume to be directly related to the major event. Such earthquakes as the Tohoku and the July 29, 2025 Kamchatka occur in the extremely tectonically and seismically active environment. In some sense, this makes them random cases with unpredictable rupture trajectory and generated seismic energy. For example, the July 20, 2025 Kamchatka earthquake was very close to the July 29 event in time and space but failed to progress  along the same rapture path. These and other major events have to be studied in a procedure similar that developed in this study.

 

 

Discussion

The WCC-based detection and phase association methods allow to reduce the detection threshold by almost an order of magnitude, create statistically significant event hypotheses of very weak sources, and to improve the estimates of relative and absolute location of these low-magnitude seismic events. The usage of array stations provides additional advantages to the WCC detection and estimation of azimuth and scalar slowness of valid arrivals. The primary IMS seismic network includes almost 30 arrays together with 6 arrays of the auxiliary network. The quality of the REB used in the CTBT monitoring and many civilian applications is mainly defined by the global distribution of array stations.             When applied to the IMS seismic data, the WCC processing found additional 50% to 100% of new events in the studied areas, all confirmed by IDC analysts and promoted to lebview database account of the IDC as REB-ready events.

            No events were detected by the IDC in the studied area of the May 24, 2013, Sea of Okhotsk earthquake for more than 14 months before it occurred. More than 200 events extra to the REB were found  between May 13, and the mainshock of the Sea of Okhotsk earthquake using data from the IMS seismic stations and the WCC-based methods of detection and phase association. A sudden increase in the occurrence rate of the low-magnitude events was detected on the afternoon May 19th. The absence of events found by standard detection and phase association methods is an evidence of the problem in the study of the earthquake preparation processes. The relevant events occur below the detection threshold even when they are not hidden in the surrounding seismicity. More sensitive methods are needed.

            The transition to the elevated seismicity at an ultra-low-magnitude level is likely related to the change in mechanical state of the volume of the following major earthquake as well as the change in other geophysical fields relevant to the deep earthquakes. The evolution of the numbers of these events in various magnitude ranges within the source volume may demonstrate some features, which can be related to the physical processes leading to the Okhotsk Sea earthquake. The observed pattern corresponds to the propagation of seismic activity from the lowest magnitudes and sizes to those compatible with the smallest detected aftershocks. This study is mostly focused on the solution of general methodological problems related to detection and phase association. A trial-and-error method was thoroughly used to obtain the first results demonstrating the potential of the WCC-based methods. Further studies are needed to adjust the set of defining parameters to the optimal resolution of the ultra-low-magnitude seismicity.   

            The final stage of the preparation started from 18 to 12 hours before the mainshock when concentration of the found events started to increase near the lower boundary of the strict LA version with the origin time tolerance of 0.25 s. This was a signal of the readiness to start. If to presume that the preparation processes before major earthquakes follow the same sequence a special attention has to be focused on the  search  for similar features of seismic activity below standard detection thresholds of global seismic networks. The study of a few major earthquakes has been started and the first results are promising.

            The Sea of Okhotsk earthquake is a specific, almost laboratory, case since the source volume is almost fully isolated and its seismic activity is very low between the larger events. These features allow for very low-amplitude signals to be detected using the WCC-based methods. There is a possibility even to reduce the detection threshold by noise stochastization also called noise "whitening" as was demonstrated for the aftershocks of the Kamchatka July 29, 2025 earthquake and in this study. This method can be effectively used in noisy environment in seismically active areas. 

 

           

References

            Adushkin V.V., Kitov I.O., Konstantinovskaya N.L., Nepeina K.S., Nesterkina M.A., and Sanina I.A. Detection of ultraweak signals on the Mikhnevo small-aperture seismic array by using cross-correlation of waveforms // Dokl. Earth Sci. 2015. V. 460. № 2. P. 189–191.

Adushkin V.V., Bobrov D.I., Kitov I.O., Rozhkov M.V., Sanina I.A. Remote detection of aftershock activity as a new method of seismic monitoring // Dokl. Earth Sci. 2017. V. 473. № 1. P. 303–307.

            Adushkin V. V., Kitov I. O., and Sanina I. A. Further Development of the Matched Filter Method for Solving Seismological Problems // Doklady Earth Sciences. 2025. V. 523:13. https://doi.org/10.1134/S1028334X25606182

Arrowsmith S. J., Eisner L. A technique for identifying microseismic multiplets and application to the Valhall field, North Sea // Geophysics. 2006. V. 71. P. 31‒40.

Baisch S., Ceranna L., Harjes H.-P. Earthquake Cluster: What Can We Learn from Waveform Similarity? // Bull. Seismol. Soc. Am. 2008. V. 98. P. 2806‒2814.

Bobrov D., Kitov I., Zerbo L. Perspectives of cross correlation in seismic monitoring at the International Data Centre // Pure Appl. Geophys. 2014. V. 171. № 3–5. P. 439–468.

Bobrov D.I., Kitov I.O., Rozhkov M.V., Friberg P. Towards Global Seismic Monitoring of Underground Nuclear Explosions Using Waveform Cross Correlation. Part I: Grand Master Events // Seismic Instruments. 2016a. V. 52. № 1. P. 43–59.

Bobrov D.I., Kitov I.O., Rozhkov M.V., Friberg P. Towards Global Seismic Monitoring of Underground Nuclear Explosions Using Waveform Cross Correlation. Part II: Synthetic Master Events // Seismic Instruments. 2016b. V. 52. № 3. P. 207–223.      

Bobrov D., Kitov I., Rozhkov M. Studying seismicity of the Atlantic Ocean using waveform cross-correlation // NNC RK Bulletin. 2017. 2(70). P. 5-19.

Comprehensive Nuclear-Test-Ban Treaty. 1996. Protocol to the Comprehensive Nuclear-Test-Ban Treaty. https:www.ctbto.org.sites.default.files.Documents.treatytext.tt.html

Coyne J., Bobrov D., Bormann P., Duran E., Grenard P., Haralabus G., Kitov I., Starovoit Yu. Chapter 15: CTBTO: Goals, Networks, Data Analysis and Data Availability / Ed. P. Bormann // New Manual of Seismological Practice Observatory. 2012. P. 1–41. DOI: 10.2312/GFZ.NMSOP-2_ch15  

Gibbons S., Ringdal F. A waveform correlation procedure for detecting decoupled chemical explosions // NORSAR Scientific Report: Semiannual Technical Summary № 2. NORSAR, Kjeller, Norway. 2004. P. 41–50.

Gibbons S.J., Ringdal F. The detection of low magnitude seismic events using array based waveform correlation // Geophys. J. Int. 2006. V. 165. P. 149–166.

            Gibbons S.J., Sorensen M.B., Harris D.B. and Ringdal F. The detection and location of low magnitude earthquakes in northern Norway using multichannel waveform correlation at regional distances // Phys. Earth Planet. Inter. 2007. Vol. 160. N. 3–4. P. 285–309.

            Gibbons S.J., Schweitzer J., Ringdal F., Kværna T., Mykkeltveit S. and Paulsen B. Improvements to seismic monitoring of the European Arctic using three component array processing at SPITS // Bull. Seismol. Soc. Am. 2011. Vol. 101. N. 6. P. 2737–2754.                  

Gibbons S.J., Pabian F., Näsholm S.P., Kværna T., Mykkeltveit S. Accurate relative location estimates for the North Korean nuclear tests using empirical slowness corrections // Geophys. J. Int. 2017. V. 208. № 1. P. 101–117. 

            Gutenberg, B., & Richter, C. F. (1936). Magnitude and energy of earthquakes. Science, 83, 183–185. 

            Gutenberg, B., & Richter, C. F. (1945). Earthquakes magnitude, intensity, energy and acceleration. Bulletin of Seismological Society of American, 46(3), 105–145.

Israelsson H. Correlation of waveforms from closely spaced regional events // Bull. Seismol. Soc. Am. 1990. V. 80. № 6. P. 2177–2193.

            Joswig M. Pattern recognition for earthquake detection // Bull. Seismol. Soc. Am. 1990. V. 80. P. 170–186.

            Joswig M., Schulte-Theis H. Master-event correlations of weak local earthquakes by dynamic waveform matching // Geophys. J. Int. 1993. V. 113. Issue 3. P. 562‒574. DOI: https://doi.org/10.1111/j.1365-246X.1993.tb04652.x

Kim W.-Y., Richards P., Schaff D. et al. Identification of Seismic Events on and near the North Korean Test Site after the Underground Test Explosion of 3 September 2017 // Seismol. Res. Lett. 2018. V. 89. № 6. P. 2120–2130.

Kitov I.O., Rozhkov M.V. New Applications at the International Data Centre for Seismic, Hydroacoustic and Infrasound Expert Technical Analysis. Vienna. Austria / Eds M. Kalinowski et al. // Twenty-five Years Progress of the Comprehensive Nuclear-Test-Ban Treaty Verification System. PTS Preparatory Commission for the CTBTO. 2024. P. 233‒254.        

Kitov I.O., Sanina I.A. Analysis of Sequences of Aftershocks Initiated by Underground Nuclear Tests Conducted in North Korea on September 9, 2016 and September 3, 2017 // Seism. Instr. 2022. V. 58. P. 567–580.

Kitov I.O., Sanina I.A. Recovery of the Aftershock Sequence of the North Atlantic Earthquake Using Waveform Cross-Correlation // Geodynamics & Tectonophysics. 2025a.  16 (4), 0838. doi:10.5800/GT-2025-16-4-0838

Kitov I.O., Sanina I.A. Quality assurance of seismological bulletins using waveform cross-correlation // Russian Seismological Journal. 2025b. V. 7, N1. P. 26–41.  (in Russian) DOI: 10.35540/2686-7907.2025.1.02.

Kitov I.O, Sanina I.A., Sokolova I.N., Vinogradov Yu.A. Study of Seismic Activity between the Mainshock and the First Aftershock of the Kamchatka Earthquake on July 29 2025 // 2026. Russian Journal of Earth Sciences, (in press)

            Kitov I.O., Volosov S.G., Kishkina S.B., Konstantinovskaya N.L., Nepeina K.S., Nesterkina M.A., and Sanina I.A.  Detection of Regional Phases of Seismic Body Waves Using an Array of Three Component Sensors // Seismic Instruments. 2016. Vol. 52. N. 1. P. 19–31.

            Kitov I. O., Sanina I. A.,  Volosov S. G., Konstantinovskaya N. L.  The 20th Anniversary of the Installation of the Small-Aperture  Mikhnevo Array for Monitoring Induced Seismicity // 2025. Izvestiya, Physics of the Solid Earth. V. 61. N. 2. P. 288–304. DOI: 10.1134/S1069351325700181

            Saragiotis C., Kitov I. Tuning IMS station processing parameters and detection thresholds to increase detection precision and decrease detection miss rate // EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-8949. 2020.  DOI: 10.5194/egusphere-egu2020-8949  

            Schaff D.P. and Richards P.G. Repeating seismic events in China // Science. 2004. Vol. 303. P. 1176–1178.

            Schaff D. P. and Richards P.G. Improvements in magnitude precision, using the statistics of relative amplitudes measured by cross correlation // Geophys. J. Int. Seismology. 2014. 197(1). P. 335–350. DOI: 10.1093/gji/ggt433

            Schaff, D.P., Kim, W.-Y. and Richards P.G. Background Seismicity for parts of the northern Korean peninsula // 2025. CTBT Science & Technology Conference. P1.2-714. Vienna, Austria.

            Schweitzer J., Fyen J., Mykkeltveit S., Gibbons S.J., Pirli M., Kühn D. and Kværna T. // Seismic arrays, in New Manual of Seismological Practice Observatory. Bormann. P., Ed. 2012. Ch. 9. doi: 10.2312/GFZ.NMSOP2_ch9.

Selby, N. Relative location of the October 2006 and May 2009 DPRK announced nuclear tests using International Monitoring System Seismometer arrays // Bull. Seismol. Soc. Am. 2010. V. 100. P. 1779-1784.

            Turin, G. L. An introduction to matched filters //  IRE Transactions on Information Theory. 1960.V. 6. PP. 311–329. doi:10.1109/TIT.1960.1057571. S2CID 5128742

vDEC. 2026. /https://www.ctbto.org/resources/for-researchers-experts/vdec

Waldhauser, F., Schaff, D. Large-scale cross correlation based relocation of two decades of northern California seismicity // J. Geophys. Res. 2008. V. 113. B08311. doi:10.1029/2007JB005479.

 

 

Appendix 1.

Table 1A. List of 70 master events

ORID (IDC)	epoch time, s	lat, deg	lon, deg	Nass	mb	depth, km
10135295	1380598702.792	53.168	152.861	222	5.78	582.31
9813664	1369374289.657	54.886	153.312	163	6.27	604.236
5761403	1260412253.013	53.448	152.746	163	5.35	660.977
9853071	1370948198.355	54.175	152.882	122	4.24	636.182
9824827	1369731519.192	54.258	153.314	89	3.84	629.503
9813875	1369380725.728	54.942	153.55	72	3.58	614.955
17205602	1554592969.038	53.327	153.648	63	3.4	514.569
5066131	1228431005.988	53.078	152.769	61	3.65	572.169
7135786	1297181407.808	53.181	153.345	58	3.78	518.525
4826051	1215223924.014	53.905	153.028	58	5.85	627.028
9813345	1369518064.658	53.61	153.426	56	3.56	539.297
23854429	1680745370.046	53.204	153.505	50	3.2	510.653
9884947	1372297849.439	54.151	153.47	50	3.58	635.21
21500256	1639353003.391	53.308	153.183	48	3.49	543.782
9813788	1369379233.433	54.641	153.26	47	3.28	616.858
9810781	1369441175.714	54.567	153.129	45	3.66	620.434
10168383	1381576364.328	54.425	153.153	38	3.28	613.307
9841029	1370474984.910	54.212	153.31	38	3.19	631.664
3230677	1115676914.851	53.31	153.481	35	3.59	501.793
14489716	1496701706.135	53.083	152.718	34	3.22	547.991
9820186	1369572061.070	54.869	153.089	33	3.21	620.501
9813674	1369406005.282	54.164	153.452	32	3.27	637.238
10721708	1396757164.906	53.153	153.228	31	3.24	515.628
9810966	1369381262.379	54.56	153.149	29	3.07	617.183
8369312	1328497880.060	53.702	153.743	29	3.11	500.868
4055723	1170136677.811	53.098	153.366	29	3.35	511.446
9811200	1369377323.512	54.345	153.366	28	3.25	627.911
9813886	1369403385.967	54.219	153.175	26	2.97	626.165
9813792	1369376719.598	54.382	153.041	26	3.66	628.772
9813791	1369379662.599	54.176	153.509	25	2.86	633.172
9811242	1369398901.028	54.524	153.163	24	3.05	618.703
9811209	1369377661.761	54.27	153.028	24	2.97	624.959
9809264	1369377742.847	53.161	153.169	22	2.98	511.99
2699433	1092728752.311	53.639	153.856	22	3.23	522.065
9811104	1369421983.175	54.558	153.103	20	2.83	616.142
9810705	1369377079.908	55.131	153.765	19	3.23	603.144
24428519	1690114315.736	53.204	153.673	18	2.8	503.885
10886687	1400716754.619	54.463	153.322	17	2.88	613.972
7572355	1304827534.681	53.368	153.5	17	2.8	506.299
3088584	1109401670.402	53.491	153.694	17	3.02	501.396
1972498	1067433219.035	53.344	153.554	17	3.27	514.242
9810611	1369390865.824	54.827	153.532	16	2.63	620.625
9809372	1369382410.885	54.184	153.365	16	2.75	633.731
29423659	1764193387.507	53.475	153.836	15	3.04	510.526
28634412	1753832718.971	55.927	152.869	15	3.8	585.023
12375814	1441527964.074	53.186	153.509	15	2.79	508.303
9822022	1369557714.879	53.543	153.734	15	2.88	508.593
9813031	1369385141.779	54.437	153.511	15	2.84	616.076
9809262	1369382861.729	54.649	153.85	15	2.78	567.235
10234531	1383401230.684	54.276	153.045	14	2.7	622.781
4429214	1193146082.919	54.006	153.95	14	3.13	509.094
9813147	1369387000.640	54.05	153.364	13	2.62	627.19
9812995	1369385016.762	54.129	153.365	13	2.68	633.355
9811489	1369377821.856	54.531	153.235	13	2.81	609.585
9811181	1369377279.646	54.18	153.506	13	3.16	641.86
9821981	1369545820.709	54.901	153.316	12	2.65	565.503
9811284	1369375365.180	55.282	154.02	12	3.93	606.697
24970171	1698470745.129	53.003	152.641	11	2.84	540.441
9887080	1372297878.828	54.934	152.816	11	2.73	699.052
9857445	1371270785.773	54.847	153.268	11	2.72	614.875
9853538	1371015355.260	54.752	153.587	11	2.85	599.315
9811275	1369411887.652	54.863	153.404	11	2.71	620.76
9811250	1369400392.260	54.109	152.842	11	2.63	628.53
9811178	1369377204.923	54.262	153.127	11	3.14	638.519
3175041	1112732892.358	53.267	153.544	11	2.76	508.515
22386776	1655810590.838	53.833	154.078	10	2.98	507.946
9826638	1369989297.098	54.425	153.371	10	2.53	612.616
9820010	1369507278.262	54.707	152.895	10	2.43	627.294
9813299	1369389498.103	54.794	153.197	10	2.5	608.356
9813000	1369385056.226	54.176	153.144	10	2.57	626.751

Appendix 2.

Table. The not-found REB events for the strict version and Δt=0.25 s