Earthquake detection is carried out following the coherence-based approach proposed by Adinolfi et al. (2020). The detection module performs a time-shifting and stacking of characteristic functions computed on seismic traces recorded at different stations, aiming to identify and preliminarily locate earthquakes inside a pre-defined spatial grid of potential source locations. In this module, only the S-wave arrivals are used for the detection. Then, the module classifies detected seismic events as real or false based on the presence of a minimum number of triggered stations, their distances from the epicenter, and the coherence of their arrival times (see Heimann et al., 2017; Adinolfi et al., 2019).

The picking of P- and S-wave arrival times is based on the algorithm of Ross and Ben-Zion (2014) using the triggered stations from the Earthquake Detection Module. Specifically, the P-pick is performed by requiring the exceedance of a threshold of the STA/LTA characteristic function around the theoretical P-wave arrival time. The S-pick is performed by searching the maximum point of the kurtosis derivative function.

This module implements the procedure proposed by Zollo et al. (2021) for event locations following a probabilistic approach based on the P- and S-phase picks and the travel-time grids generated by the NLLoc software (Lomax et al., 2000; http://alomax.free.fr/nlloc). The module gives as output the earthquake location and a map of the PDF distribution with the maximum likelihood hypocenter. The location errors are defined by considering the 31% and 68% significance levels of the probability density function (PDF), which corresponds to ±1σ.

This module exploits an approach proposed by Picozzi et al. (2017) for earthquake early warning purposes and extended to monitoring applications by Picozzi et al. (2018, 2019), and Spallarossa et al. (2021b). For each record, it obtains an estimate of the apparent energy by integrating the S-transform instantaneous spectrum, after correcting it for the geometrical spreading. Then, the event radiated energy is computed as the average among single station estimates. The corrective factor for geometrical spreading is set in the configuration file. In addition, this module measures peak ground velocity (PGV) and acceleration (PGA) for the three component recordings.

With the aim of determining the seismic moment (M₀) and the moment magnitude (M_w) of an earthquake, this module exploits the time evolution of the P-wave amplitude as proposed by Colombelli et al. (2014, 2015), and Nazeri et al. (2019). The logarithm of the P-wave peak (acceleration and velocity) amplitude is measured on the vertical component of recorded waveforms, in expanding time windows after the P-wave arrival, and after correction of the amplitude for the geometrical spreading. Corrective factors are set in a configuration file.

This module inverts the observed, absolute initial P-wave peak amplitudes, corrected for the geometrical attenuation effect, and the P- polarities to estimate focal mechanisms in a Bayesian framework (Tarantino et al., 2019). Data for the focal mechanism estimation are selected using an SNR threshold, which is set in the configuration file along with the minimum number of readings to perform the estimation.

This module evaluates the average V_P/V_S ratio along the ray paths to each seismic station from observations of the difference between P- and S- arrival times divided by the P-wave travel time from the source to the station (Wadati and Oki, 1933; Kisslinger and Engdahl, 1973; Lucente et al., 2010; Chiarabba et al., 2009). The temporal evolution of the V_P/V_S ratio for a set of stations can enable identification of both spatial and temporal changes in the medium properties (Amoroso et al., 2018; De Landro et al., 2022).

We validated the TREMOR platform by applying it to a dataset of continuous waveforms from the Northern Nagano region (Japan) that was recorded by the Japanese High Sensitivity Seismograph Network (Hi-net, https://www.hinet.bosai.go.jp/), operated by the National Research Institute for Earth Science and Disaster Prevention (National Research Institute for Earth Science and Disaster Resilience, 2019). We selected the Nagano area in Honshu, Japan, for two reasons: 1) its high seismogenic potential and seismicity due to the presence of the Itoigawa–Shizuoka Tectonic Line (ISTL), which cuts the Island of Honshu and extends for ∼150 km from Itoigawa City on the Sea of Japan to Shizuoka City on the Pacific Ocean (Figure 2A); and 2) the high station density in this area, with an average inter-station distance of about 10 km over a maximum aperture of 70 km.

The ISTL strikes NNE–SSW to NNW–SSE and is one of the most active faults in Japan, since it forms one arm of the triple junction of the Eurasian, North American, and Philippine Sea plates. The northernmost segment of the ISTL, the Kamishiro fault, has been shown to be an active thrust fault (Sato et al., 2004; Takeda et al., 2004). On the other hand, the Otari–Nakayama fault, which runs parallel to the Kamishiro and East Matsumoto Basin faults, shows no geomorphological evidence of late Quaternary activity (Matsuta et al., 2004; Ueki, 2008).

We selected a study area of about 40 km × 40 km (blue rectangle in Figure 2B) internal to the larger area covered by the network and collected continuous three-component velocity recordings from 25 borehole stations of Hi-net network in and near this area, (yellow triangles in Figure 2C) with a sampling rate of 100 Hz. Continuous waveform data were downloaded from NIED Hi-net website (https://www.hinet.bosai.go.jp/). We considered the period from the 1st to the 31st of May 2011, for which the Japan Meteorological Agency (JMA) reported 101 earthquakes with magnitudes ranging between −0.6 and 2.6 and depths up to 15 km (Figure 2C, circles). The JMA catalog magnitude of completeness estimated for this time period is 0.2.

Before running the code, we formatted the data in mseed format, organized it in files of one-hour length, and wrote input files with all necessary information (e.g., seismic network, velocity model). Moreover, we modified each control file for TREMOR modules with a set of parameters properly tuned for the characteristics of Nagano seismicity to be analyzed (for more details see Supplementary Tables S1–S6).

For the earthquake detection, we defined a 3D grid of potential source locations with a size of 70 km × 75 km × 20 km, with spacing equal to 1.5 km and centered on the barycenter of the seismic network used (blue rectangle in Figure 2B). For the parameter tuning, a specific analysis on a 10-day training dataset was carried out by modifying the parameters and optimizing them based on the number of real and false detections. After a visual inspection of the waveforms, we compared the preliminary results of this analysis with seismicity data from the JMA catalog. The chosen detection parameters are: coherence greater than 700, a signal-to-noise ratio greater than or equal to 3, which is used for the definition of triggered stations (with 10 s of signal, 3 s of noise), and three triggered stations among the 10 closest stations to the epicenter. Figure 3 shows an example of the detection module application for an M_V 2.5 earthquake (velocity magnitude as calculated by JMA), which occurred on 16/05/2011 at 13:25:04 (JST) in the study area. Waveforms (Figure 3A) and characteristic functions (Figure 3B, normalized amplitudes) are calculated for each station to evaluate the vertical stack of the coherence map for the search region (Figure 3C) and global image functions corresponding to the best fit source time (Figure 3D).

Despite this positive performance, it is worth noting that the TREMOR event detection algorithm missed 5 earthquakes which are in the JMA catalog, with Mv < 0, while other six events that were external to the seismic network and target area which were mislocated at the border of the location grid. These cases can occur when very small magnitude events (M_v < 0) are not detected or are detected with low SNR, resulting in mis-location inside the target area or at the border of the grid due to too few and/or incorrect phase readings.

To properly set the parameters that guide the picking procedure, we performed an optimization analysis based on comparisons between manual and automatic P-wave arrival time identifications on a sub-set of Nagano events. We explored a wide range of possible parameters and selected as the “best” configuration the one that minimizes the average differences between manual and automatic picks, for all source-station pairs and for both the vertical and horizontal components (average on the order of 0.05 s). In the optimization procedure, we also tried to minimize the number of missed picks (less than 15% for the “best” parameters). The picking parameters that we selected are listed in the Supplementary Table S4.

For the earthquake locations, we used the same 1D velocity models that JMA uses for locating Japanese seismicity (Ueno, 2002). We parameterized the crustal structure using a 3D grid with spacing 0.5 km and size of 70 km × 75 km × 20 km. The theoretical travel times from the grid nodes to the seismic stations were calculated using the Eikonal finite-difference scheme of (Podvin and Lecomte, 1991), as implemented by the Non LinLoc software (Lomax et al., 2000) for absolute earthquake locations in 3D velocity media. We used the same 3D grid for both location steps during the analysis. We located earthquakes with at least 2 P- and 2 S-phase picks, after discarding arrival times with residuals for P- and S-waves higher than 0.5 s and 1 s, respectively. Although it is common practice to use a higher minimum number of phases, we have chosen to use only four phases to maximize the number of micro-earthquake locations.

Figure 4 shows an example of the picking module application for an M_V 2.5 earthquake which occurred on 16/05/2011 at 13:25:04 (JST) in the study area. The characteristic functions (Figure 4A) and the phase picking on the vertical (Figure 4B) and horizontal (Figure 4C) velocity records are shown.

Figure 5 shows the probabilistic location for the same event. After the first run of data processing with the defined spatial grid, 406 earthquakes were detected by the platform of which 339 were effectively located. Based on their location, pick number, and location quality, 157 earthquakes were identified as located within the smaller target grid (red rectangle in Figure 2B). Of the 157 events, 96 are also in the JMA catalog. Five events in the JMA catalog (with M ≤ 0) were missed by TREMOR, while it provided 61 new detections. After visual inspection, we confirmed that all of the new detected events are real earthquakes, with no false detections.

For the 157 detected and located events in the study area, TREMOR provided 1286 P-wave picks and 1216 S-wave picks in the two steps of phase picking. To evaluate the quality and robustness of the automatic picking, we manually picked a subset of 100 events among the final catalog, selected in order to be representative in terms of location and magnitude of the entire seismicity, and we compared the automatic and manual P and S picks (see Supplementary Material). The mean values of the time differences between manual and automatic P- and S-wave onsets are equal to −0.06 and 0.08 s, respectively, with standard deviations of 0.29 and 0.31 (Supplementary Figure S1). These mean values, along with the histograms in Supplementary Figure S1, indicate that the automatic arrival time picks tend to be slightly later than the manual picks.

In the first step of the locations, the events were located with mean errors of 1.2, 2.0, and 2.0 km for the NS, EW and vertical directions, respectively, with an average root-mean-square (RMS) of 0.5 s and a mean of 5 P- and 5 S- phases. The final earthquake locations, shown in Figure 6, have mean location errors equal to 1.1, 1.9, 2.0 km along the NS, and EW horizontal and vertical directions, respectively, a mean RMS (weighted) of 0.2 s, and a mean of 8 P- and 8 S-phases per event (Figure 6A).

The seismicity pattern delineated by the 157 located events correlates well with the main seismicity features of the area (Figure 6B; Supplementary Figure S2). The seismicity appears mainly aligned along the two segments of the ISTL, the Kamishiro and the Otari–Nakayama faults, down to a depth of 15 km, under the currently dominant east–west compressional stress field (Ando et al., 2017; Panayotopoulos et al., 2014; Panayotopoulos et al., 2016). The hypocenters of the 157 events show the same features of the earthquakes clustered mainly along Otari–Nakayama fault and its NE branch fault, in the southern and eastern part of Matsumoto basin, defining in depth an E-dipping plane (Supplementary Figure S2).

The JMA assigns earthquakes a velocity magnitude (M_JMA), determined from the maximum amplitudes of velocity seismograms (Funasaki and Earthquake Prediction Information Division, 2004). Unfortunately, we could not use the M_v relation directly because the distance and depth corrections are not accessible from JMA web service. So, to better compare our results with those of the JMA catalog, we calculated a velocity magnitude (M_v) based on a linear regression between the maximum horizontal velocities at selected seismic stations and JMA catalog magnitudes. Looking at Figure 7, a good agreement exists between the magnitude estimates of TREMOR and the JMA magnitude values. Approximately linear scaling relations with M_JMA are evident for M_v, M_e and M_w, with mean differences of 0.09, −0.34, −0.55 and standard deviations of 0.17, 0.17, 0.39, respectively.

Fault plane solutions range from reverse to strike-slip (Figures 8, 9A), in accordance with the complex stress regime of this region. We found 57 focal mechanism solutions with an average RMS equal to 0.18 and a mean mismatch on polarity equal to 0.21 (see Supplementary Material). Matsumoto basin is characterized by the interaction of the Kamishiro fault with dip 30°–45° SE and the Otari-Nakayama fault that dips more steeply, at 50°–65° SE. According to Panayotopoulos et al. (2016), focal mechanism solutions for aftershocks of the 2014, M 6.7 Northern Nagano earthquake exhibit both reverse and left-lateral components, with a large strike-slip component distributed along a branch fault with a steep dip-angle (i.e., Otari-Nakayama), consistent with the fault plane solutions estimated by TREMOR.

In order to validate the V_P/V_S values obtained by the automatic module application, we compared our results with the tomographic results from the work of Panayotopoulos et al. (2014) (Figures 6C, 8). First, we calculated a V_P/V_S ratio value for each station by averaging the values obtained for the entire month (Supplementary Figure S3). Then, we interpolated the V_P/V_S ratios at the stations to obtain a map that can be compared to the tomographic image. The V_P/V_S map shown in Figure 9B agrees with the results by Panayotopoulos et al. (2014)). Indeed, they found a range of V_P/V_S ratio between 1.5 and 1.9, high V_P/V_S anomalies along the ISTL fault and a mean V_P/V_S value of about 1.73 for the other areas.

Figure 10 compares the frequency-magnitude distributions (Gutenberg and Richter, 1942; Aki, 1965; Wiemer, 2001) for the TREMOR and JMA catalogs for the same area and time period. We observe that the TREMOR catalog has more small-magnitude events than the JMA catalog and a lower magnitude of completeness, Mc, of M_v −0.12. In comparison, the JMA catalogue is complete for M_v 0.20 and larger.