We applied the HIST-ETAS model (e.g., Ogata et al. 2003; Ogata 2004; Bansal and Ogata 2013; Ueda et al. 2021) to the observed seismicity rate to investigate spatial variations in the seismicity characteristics and evaluate the probability of a given event being a background event. This model is a point-process model that includes the Omori-Utsu law (Utsu 1961; Utsu et al. 1995), which formulates the typical aftershock temporal decay, the Utsu-Seki law (Utsu and Seki 1954), which formulates the relationship between the aftershock area and mainshock magnitude, the decay in aftershock density with distance (e.g., Felzer and Brodsky 2006), and a branching process, such that each earthquake, regardless of magnitude, has the ability to increase the probability of triggering a future earthquake (Iwata 2009).
The earthquake occurrence rate \(\lambda\) at time and location \(\left(t,x,y\right)\) and occurrence history \({H}_{t}\) can be expressed by:
$$\begin{array}{c}\lambda \left(t,x,y|{H}_{t}\right)=\mu +{\sum }_{{t}_{i}<t}\frac{K}{{\left(t-{t}_{i}+c\right)}^{p}}{\left[\frac{\left(x-{x}_{i}, y-{y}_{i}\right){S}_{i}\left(\begin{array}{c}x-{x}_{i}\\ y-{y}_{i}\end{array}\right)}{\text{e}\text{x}\text{p}\left\{\alpha \left({\text{M}}_{\text{i}}-{\text{M}}_{\text{c}}\right)\right\}}+d\right]}^{-q} \#\left(1\right)\end{array}$$
,
where \(\mu\) is the background seismicity rate, and the second term expresses the rate of the earthquake occurrence triggered by a magnitude \({\text{M}}_{\text{i}}\) earthquake at time and location \(\left({t}_{i}, {x}_{i}, {y}_{i}\right)\), where \(K\), \(p\), and \(c\) are the parameters of the Omori-Utsu law, \({S}_{i}\) is a non-dimensional positive definite symmetric matrix for anisotropic clusters that is determined by identifying the aftershock cluster using a magnitude-based clustering algorithm (Ogata et al. 1995; Ogata 1998) and then choosing the best-fit ellipsoid that represents the cluster, \(\alpha\) is the aftershock magnitude sensitivity, \(q\) is the aftershock spatial decay rate, \(d\) is a constant, and \({\text{M}}_{\text{c}}\) is the cutoff magnitude, \(4.0\) in this analysis. We note that \({H}_{t}\) includes earthquakes that occurred before the target period (precursory period) because they are potentially influential to the seismicity in the target period (Ogata 2011).
The five seismicity parameters (\(\mu\), \(K\), \(\alpha\), \(p\), and \(q\)) are given as a function of space and are expressed as:
$$\begin{array}{c}\begin{array}{c}\mu \left(x,y\right)=\widehat{\mu }exp\left\{{\vartheta }_{\mu }\left(x,y\right)\right\}\\ K\left(x,y\right)=\widehat{K} exp\left\{{\vartheta }_{K}\left(x,y\right)\right\}\\ \alpha \left(x,y\right)=\widehat{\alpha }exp\left\{{\vartheta }_{\alpha }\left(x,y\right)\right\}\\ p\left(x,y\right)=\widehat{p}exp\left\{{\vartheta }_{p}\left(x,y\right)\right\}\\ q\left(x,y\right)=\widehat{q}exp\left\{{\vartheta }_{q}\left(x,y\right)\right\}\end{array}\#\left(2\right)\end{array}$$
,
where \(\widehat{\mu }\), \(\widehat{K}\), \(\widehat{\alpha }\), \(\widehat{p}\), and \(\widehat{q}\) correspond to the geometric mean value of each parameter averaged over the analysis region. We adapt each \(\vartheta \left(x,y\right)\) value to the data by expressing each function using many coefficients that are placed in the locations of each earthquake epicenter and some additional points on the boundary of the analysis region. Each \(\vartheta \left(x,y\right)\) value at an arbitrary location is linearly interpolated using the three values at the vertices of each Delaunay triangle.
The unknown parameters can be estimated via the maximum likelihood estimation method. The log-likelihood is expressed as:
$$\begin{array}{c}\text{log}L={\sum }_{k}\text{log}\lambda \left({t}_{k}, {x}_{k}, {y}_{k}|{H}_{{t}_{k}}\right)-{\int }_{0}^{T}{\iint }_{S}^{}\lambda \left({t}^{{\prime }}, x, y|{H}_{{t}^{{\prime }}}\right)dxdyd{t}^{{\prime }} \#\left(3\right)\end{array}$$
,
where \(k\) is the number of events in the analysis, \(S\) is the analysis region, and \(\left[0,T\right]\) is the analysis interval. However, it is hard to estimate the seismicity parameters stably because the number of unknown parameters is about five times larger than the number of data (events) used in the analysis. We obtained stable solutions by penalizing the spatial gradient of the parameter functions under the assumption of using smoothed functions (the facets of the piecewise linear function being as flat as possible):
$$\begin{array}{c}Q\left(\theta |\tau \right)=\sum _{k=1}^{5}{w}_{k}{\iint }_{A}^{}\left\{{\left(\frac{\partial {\vartheta }_{k}}{\partial x}\right)}^{2}+{\left(\frac{\partial {\vartheta }_{k}}{\partial y}\right)}^{2}\right\}dxdy \#\left(4\right)\end{array}$$
,
where \(\theta\) is a set of seismicity parameters, and \(\tau =\left({w}_{1}, \dots , {w}_{5}\right)\) is a set of weights that is used to tune the strength of the constraints. We then estimated the seismicity parameters that maximized the penalized log likelihood function (Good and Gaskins 1971):
$$\begin{array}{c}R\left(\theta \right)=\text{ln}L\left(\theta \right)-Q\left(\theta |\tau \right) \#\left(5\right)\end{array}$$
.
The weights \(\tau\) were objectively tuned by maximizing the integrated posterior distribution:
$$\begin{array}{c}\varLambda \left(\tau \right)={\int }_{{\Theta }}^{}L\left(\theta \right)\pi \left(\theta |\tau \right)d\theta \#\left(6\right)\end{array}$$
,
where \(\pi \left(\theta |\tau \right)\propto \text{e}\text{x}\text{p}\left\{-Q\left(\theta |\tau \right)\right\}\) is the probability distribution. The solution of \(\theta\) was then obtained for fixed weights \(\tau\) by maximizing the penalized log-likelihood in Eq. (5), which yields the optimal maximum a posteriori estimates. See Ogata et al. (2003) and Ogata (2004) for the details of the parameter estimation.
We estimated the uncertainties in the HIST-ETAS parameters at each location following the method outlined in Ueda et al. (2021). We first resampled the data 100 times by randomly extracting 90% of the earthquake data used in this analysis. Note that we did not exclude the M9 mainshock from the resampled data set. We then applied the HIST-ETAS model to each resampled dataset and estimated the spatial patterns of the five seismicity parameters. We note that the reference parameters (\(\widehat{\mu }\), \(\widehat{K}\), \(\widehat{\alpha }\), \(\widehat{p}\), and \(\widehat{q}\)) may vary significantly among the 100 resampled datasets owing to the trade-offs between each parameter. Herak et al. (2001) investigated the interdependence of the ETAS model parameters using the aftershock sequences of the 1996 Ston-Slano earthquake and highlighted that the degree of correlation can be large, especially for pairs of aftershock parameters. Therefore, we normalized the model parameters by dividing the reference parameters estimated from each resampled dataset, and then calculated the standard deviation of \(\vartheta \left(x,y\right)/\text{ln}10\) (common logarithm of the normalized parameter) at each location to evaluate the significance of the relative differences in each parameter, independent of the trade-offs between parameters.
We used the JMA unified hypocenter catalog (the preliminary Determination of Epicenters) as the earthquake data in our analysis. We applied the HIST-ETAS model to the \({\text{M}}_{\text{j}} \ge 4.0\) earthquakes that occurred between the timing of the 2011 Tohoku-Oki mainshock (11 March 2011) and 20 February 2021 and were located at ≤100 km depth beneath Japan (20°–50°N, 120°–150°E). The \({\text{M}}_{\text{j}}\ge 6.0\) earthquakes that occurred since 1922 and the \({\text{M}}_{\text{j}}\ge 4.0\) earthquakes that occurred since February 2011 were used as the precursory occurrence history of the HIST-ETAS model.