Our strategy for quasi-real-time precise earthquake relocation is to divide the target earthquakes into fine time windows and perform relocation each time new time window data is added. However, as it is, this approach does not constrain the relative locations of earthquakes in different time windows, and the shorter the time window, the larger the error in the relative locations in the overall distribution. We overcome this problem by incorporating a traditional simple relative relocation method; we constrain the relative locations of events in the new window from reference events already relocated in the previous time windows. To implement this, we first modify the equation used in the DD method (Waldhauser and Ellsworth, 2000) (Subsection 2.1) and then describe the quasi-real-time relocation algorithm (Subsection 2.2).

## 2.1. DD method with reference events

The arrival time \({t}_{k}^{i}\) of the seismic wave of a given event \(i\) at a given station \(k\) can be expressed as the integral along the path \(s\) as

$$\begin{array}{c}{t}_{k}^{i}={\tau }^{i}+{\int }_{s}^{}u\left(s\right)ds,\#\left(1\right)\end{array}$$

where \({\tau }^{i}\) is the origine time, and \(u\left(s\right)\) is the slowness. We first set the initial hypocentral parameters \({\varvec{m}}^{i}=({x}^{i},{y}^{i},{z}^{i},{\tau }^{i})\), Taylor expand \({t}_{k}^{i}\) around \({\varvec{m}}^{i}\) at the first order to linearize it and obtain the time difference \({r}_{k}^{i}\) for a slight change in the hypocentral parameters \(\varDelta {\varvec{m}}^{i}={\left(\varDelta {x}^{i},\varDelta {y}^{i},\varDelta {z}^{i},\varDelta {\tau }^{i}\right)}^{T}\).

$$\begin{array}{c}{r}_{k}^{i}=\frac{\partial {t}_{k}^{i}}{\partial x}\varDelta {x}^{i}+\frac{\partial {t}_{k}^{i}}{\partial y}\varDelta {y}^{i}+\frac{\partial {t}_{k}^{i}}{\partial z}\varDelta {z}^{i}+\varDelta {\tau }^{i}=\frac{\partial {t}_{k}^{i}}{\partial \varvec{m}}\bullet \varDelta {\varvec{m}}^{i},\#\left(2\right)\end{array}$$

The DD method considers earthquake pairs \(i\) and \(j\). By using \({dr}_{k}^{ij}={r}_{k}^{i}-{r}_{k}^{j}\) at the same station \(k\),

$$\begin{array}{c}{dr}_{k}^{ij}=\frac{\partial {t}_{k}^{i}}{\partial \varvec{m}}\bullet \varDelta {\varvec{m}}^{i}-\frac{\partial {t}_{k}^{j}}{\partial \varvec{m}}\bullet \varDelta {\varvec{m}}^{j}.\#\left(4\right)\end{array}$$

This equation is matrixed for all stations and earthquake pairs as follows:

$$\begin{array}{c}\left(\begin{array}{c}d{r}_{{S}_{1}}^{{E}_{1}{E}_{2}}\\ d{r}_{{S}_{1}}^{{E}_{1}{E}_{3}}\\ ⋮\\ d{r}_{{S}_{L}}^{{E}_{N-1}{E}_{N}}\end{array}\right)=\left(\begin{array}{ccccc}\frac{\partial {t}_{{S}_{1}}^{{E}_{1}}}{\partial \varvec{m}}& -\frac{\partial {t}_{{S}_{1}}^{{E}_{2}}}{\partial \varvec{m}}& 0& \cdots & 0\\ \frac{\partial {t}_{{S}_{1}}^{{E}_{1}}}{\partial \varvec{m}}& 0& -\frac{\partial {t}_{{S}_{1}}^{{E}_{3}}}{\partial \varvec{m}}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & -\frac{\partial {t}_{{S}_{L}}^{{E}_{N}}}{\partial \varvec{m}}\end{array}\right)\left(\begin{array}{c}\varDelta {\varvec{m}}^{{E}_{1}}\\ \varDelta {\varvec{m}}^{{E}_{2}}\\ \varDelta {\varvec{m}}^{{E}_{3}}\\ ⋮\\ \varDelta {\varvec{m}}^{{E}_{N}}\end{array}\right), \#\left(5\right)\end{array}$$

where \({E}_{i} (i=\text{1,2},\dots ,N)\) shows the event number and \({S}_{i} (i=\text{1,2},\dots ,L)\) shows the station number. We denote the first matrix on the right-hand side by \(\mathbf{G}\) and the second vector by \(\mathbf{m}\).

The data of the standard hypocenter relocation method is the time residuals \({r}_{\text{o}\text{b}\text{s}}={t}^{obs}-{t}^{cal}\) between observed arrival time, \({t}^{obs}\), and calculated arrival time, \({t}^{cal}\), based on the initial hypocenter parameters as data. Unlike that, the DD method uses the difference between earthquake pairs:\(\begin{array}{c}{dr}_{{k}_{obs}}^{ij}={r}_{{k}_{obs}}^{i}-{r}_{{k}_{obs}}^{j}={\left({t}^{obs}-{t}^{cal}\right)}_{k}^{i}-{\left({t}^{obs}-{t}^{cal}\right)}_{k}^{j}={\left({t}_{k}^{i}-{t}_{k}^{j}\right)}^{obs}-{\left({t}_{k}^{i}-{t}_{k}^{j}\right)}^{cal},\#\left(6\right)\end{array}\)

which is called DD (Double Difference). The following matrix equation is solved to obtain hypocentral parameter \(\mathbf{m}\).

$$\begin{array}{c}{\mathbf{d}}_{\text{o}\text{b}\text{s}}=Gm,\#\left(6\right)\end{array}$$

where \({\mathbf{d}}_{\text{o}\text{b}\text{s}}\) is a single-column vector including \({dr}_{{k}_{obs}}^{ij}\). \(\mathbf{m}\) is a single-column vector with \(\varDelta {\varvec{m}}^{i}\), and \(\mathbf{G}\) is a matrix with the difference of partial differential coefficients, corresponding to the right-hand side of Eq. (5). Actual calculations are performed by applying different weights to each piece of data. The LSQR method (Paige & Saunders, 1982) is used when dealing with extensive data.

In Eq. (6), \(\mathbf{G}\) is a massive matrix of [the number of components of DD] × [four times the number of events], which makes the method computationally expensive. As described above, we split the earthquake data into multiple time windows and perform sequential DD relocation for each new time window. This reduces the matrix size considerably, but, as it is, the information on the relative earthquake locations in different windows is also lost. To avoid losing this information, we use events whose hypocentral parameters have already been determined in past time windows as reference events. We fix the hypocenter parameters of these reference events but use the differential arrival time data between them and new events to relocate the new events. The residuals for the new and reference event pairs are minimized by relocating the hypocentral parameters of the new events.

To do so, we change the equation to be solved from Eq. (5), as shown in Fig. 2. We first remove the rows of \(\mathbf{m}\) corresponding to the reference events and include only the hypocentral parameters of the new events in \(\mathbf{m}\). From \({\mathbf{d}}_{\text{o}\text{b}\text{s}}\), we removed rows containing differential arrival time data only between reference events, but retained the differential data only between new events as well as those between new and reference events. Thus, the size of \(\mathbf{G}\) is significantly reduced.

## 2.2. Relocation Algorithm

In preparation for quasi-real-time relocation, our algorithm first relocates the hypocenters of previous events and creates a reference database. When a new time window containing new events is added, these events are relocated using the method described in Subsection 2.1, with reference events from the database. The newly relocated hypocenters are then added to the reference database, and the above procedure is repeated each time a new time window is added.

The number of reference events increases for later time windows, which becomes much larger than that of new events. Our modification removes the unknown parameters of reference events, which significantly reduces the size of the matrix \(\mathbf{G}\) compared to the standard method of simultaneously solving for all the hypocentral parameters. As a result, computational efficiency is improved, and memory requirements are reduced.