\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Fault-valve behavior estimated from intensive foreshocks and aftershocks
of the 2017 M 5.3 Kagoshima Bay earthquake sequence, Kyushu, southern
Japan}



\author[1,1]{Yoshiaki Matsumoto}%
\author[1,1]{Keisuke Yoshida}%
\author[1,1]{Akira Hasegawa}%
\author[1,1]{Toru Matsuzawa}%
\affil[1]{Tohoku University}%


\vspace{-1em}



  
  \date{May 18, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Determining fluid migration and pore pressure changes within the Earth
is key to understanding earthquake occurrences. We investigated the
spatiotemporal characteristics of intense fore- and aftershocks of the
2017 ML 5.3 earthquake in Kagoshima Bay, Kyushu, southern Japan, to
examine the physical processes governing this earthquake sequence. The
results show that the foreshock hypocenters moved upward on a sharply
defined plane with steep dip. The mainshock hypocenter was located at
the edge of a seismic gap formed by foreshocks along the plane. This
spatial relationship suggests that the mainshock ruptured this seismic
gap. The corner frequency of the mainshock supports this hypothesis. The
aftershock hypocenters migrated upward along several steeply dipped
planes. The aftershock activity slightly differs from the simple
mainshock--aftershock type, suggesting that aseismic processes
controlled this earthquake sequence. We established the following
hypothesis: First, fluids originating from the subducting slab migrated
upward and intruded into the fault plane, reducing the fault strength
and causing a foreshock sequence and potentially aseismic slip. The
continuous decrease in the fault strength associated with an increase in
the pore pressure and the increase in shear stress associated with
aseismic slip and foreshocks caused the mainshock in an area with
relatively high fault strength. The change in the pore pressure
associated with post-failure fluid discharge contributed to aftershocks,
causing the upward migration of the earthquake. These observations
demonstrate the importance of considering fluid movement at depth not
only earthquake swarms but also foreshock---mainshock--aftershock
sequences.%
\end{abstract}%



\sloppy


\textbf{Fault-valve behavior estimated from intensive foreshocks and
aftershocks of the 2017 M 5.3 Kagoshima Bay earthquake sequence, Kyushu,
southern Japan}

\textbf{Yoshiaki Matsumoto\textsuperscript{1}} \textsuperscript{+}

\textbf{, Keisuke Yoshida\textsuperscript{1*}, Toru
Matsuzawa\textsuperscript{1}, and Akira Hasegawa\textsuperscript{1}}

\textsuperscript{1} Research Center for Prediction of Earthquakes and
Volcanic Eruptions, Graduate School of Science, Tohoku University,
Sendai 980, Japan

\textsuperscript{+} Now at Japan Meteorological Agency, Japan

*Corresponding author: Keisuke Yoshida
(\emph{keisuke.yoshida.d7@tohoku.ac.jp)}

Key Points:

\begin{itemize}
\tightlist
\item
  Intensive foreshocks migrate via one plane.
\item
  Aftershock hypocenters migrate toward shallower levels via several
  planes.
\item
  Upward pore pressure migration explains the occurrence of the
  foreshock--mainshock--aftershock sequence.
\end{itemize}

Abstract

Determining fluid migration and pore pressure changes within the Earth
is key to understanding earthquake occurrences. We investigated the
spatiotemporal characteristics of intense fore- and aftershocks of the
2017 M\textsubscript{L} 5.3 earthquake in Kagoshima Bay, Kyushu,
southern Japan, to examine the physical processes governing this
earthquake sequence. The results show that the foreshock hypocenters
moved upward on a sharply defined plane with steep dip. The mainshock
hypocenter was located at the edge of a seismic gap formed by foreshocks
along the plane. This spatial relationship suggests that the mainshock
ruptured this seismic gap. The corner frequency of the mainshock
supports this hypothesis. The aftershock hypocenters migrated upward
along several steeply dipped planes. The aftershock activity slightly
differs from the simple mainshock--aftershock type, suggesting that
aseismic processes controlled this earthquake sequence. We established
the following hypothesis: First, fluids originating from the subducting
slab migrated upward and intruded into the fault plane, reducing the
fault strength and causing a foreshock sequence and potentially aseismic
slip. The continuous decrease in the fault strength associated with an
increase in the pore pressure and the increase in shear stress
associated with aseismic slip and foreshocks caused the mainshock in an
area with relatively high fault strength. The change in the pore
pressure associated with post-failure fluid discharge contributed to
aftershocks, causing the upward migration of the earthquake. These
observations demonstrate the importance of considering fluid movement at
depth not only earthquake swarms but also
foreshock---mainshock--aftershock sequences.

\textbf{1 Introduction}

An earthquake is a natural phenomenon during which a high-speed rupture
propagates along a fault. Two factors control the occurrence of an
earthquake: an increase in the shear stress acting on the fault and a
decrease in the fault strength. The results of previous studies
suggested that the increase in the pore pressure plays an important role
in the earthquake occurrence (e.g., Hasegawa, 2017; Hubbert \& Rubey,
1959; Nur \& Booker, 1972; Sibson, 1992; Rice, 1992) because it reduces
the fault strength.

A well-known example of fluid-driven seismicity is the seismicity
induced by fluid injection for engineering purposes (e.g., Ellsworth,
2013). There is also growing evidence that natural earthquake swarms are
closely related to fluid movement at depth. In fact, the characteristics
of many natural seismic swarms are similar to those of fluid
injection-induced seismicity including the migration behavior of the
earthquake hypocenter (e.g., Fischer and Hor\selectlanguage{ngerman}álek, 2003; Parotidis et
al., 2003; Bianco et al., 2004; Yukutake et al., 2011; Shelly et al.,
2016; Yoshida et al., 2016a; Ruhl et al., 2016; De Barros et al., 2019).
Based on the determination of the hypocenters and focal mechanisms of
earthquake swarm at the 2009 Hakone volcano, the diffusion of
high-pressure fluid triggered the swarm (Yukutake et al., 2011). The
spatiotemporal evolution of seismic activity in the Long Valley Caldera,
California, indicates that a pore pressure transient with a
low-viscosity fluid initiated and sustained the swarm in 2014 (Shelly et
al., 2016). It has been hypothesized that several earthquake swarms that
occurred after the 2011 Tohoku-Oki earthquake were triggered by a
decrease in the fault strength due to upward pore pressure migration
(Terakawa et al., 2013; Okada et al., 2016; Yoshida et al., 2016a,
2019a).

Not only earthquake swarms but also foreshock--mainshock--aftershock
sequences may be closely related to the fluid behavior in the Earth
interior. Sibson (1992) established the fault-valve model in which the
pore pressure cycle controls the earthquake cycle due to overpressurized
fluids that rise from the deeper portion of the fault. In this model,
fault ruptures create a transient fracture permeability within the fault
zone, which acts as a valve, promoting the upward discharge of fluids
from deeper portions of the crust. This model is supported by various
geological and geophysical observations (Sibson, 2020). Hasegawa et al.
(2005) proposed a model for the deformation process in a subduction zone
based on various geophysical observations including seismic tomography
data obtained for northeastern Japan. In this model, fluids expelled
from the subducting slab migrate upward, reach the crust, and cause
anelastic crustal deformation including earthquakes.

The migration characteristics of earthquake hypocenters can be used to
infer the origin of the seismicity (e.g., Yukutake et al., 2011; Ruhl et
al., 2016; Yoshida \& Hasegawa, 2018a,b; De Barros et al., 2019). Pore
pressure migration and aseismic slip propagation are typical mechanisms
attributed to the migration of earthquakes. In the former mechanism, the
hypocenter migration is presumed to reflect the migration of fluids
(e.g., Shapiro et al., 1997; Talwani et al., 2007). In the latter
mechanism, the hypocenter migration is presumed to be a result of
aseismic slip propagation (e.g., Lohman \& McGuire, 2007; Roland \&
McGuire, 2009). The spatiotemporal distribution of earthquake
hypocenters can be more precisely estimated than other seismological
characteristics such as the three-dimensional seismic velocity
structure. By examining relocated hypocenters, we may extract
information on aseismic physical processes controlling earthquakes,
which is crucial to understanding the earthquake generation. The results
of previous studies showed that the seismic activity caused by aseismic
processes differs from that of the mainshock--aftershock sequence (e.g.,
Hainzl \& Ogata, 2005; Roland \& McGuire, 2009; Kumazawa \& Ogata, 2013;
Yoshida \& Hasegawa, 2018b). This suggests that investigations of the
seismicity may provide clues about aseismic processes governing
earthquakes.

The volcanic front on Kyushu Island in southern Japan formed due to the
subduction of the Philippine Sea Plate. Several of the most active
volcanoes in Japan are distributed along this volcanic front (e.g.,
Sakurajima and Aso). Kagoshima Bay is located at this volcanic front
(Fig. S1), which is characterized by a low-gravity anomaly that extends
from north to south. On July 11, 2017, an M\textsubscript{L} 5.3
strike-slip earthquake occurred at a depth of \textasciitilde{}10 km in
Kagoshima Bay (Fig. 1). Seismicity activity had been recorded near the
mainshock hypocenter since December 2016 (Fig. 1c). In total, 1.843
foreshock events were recorded and listed in the Japan Meteorological
Agency (JMA) unified seismic catalogue. The seismicity increased after
the mainshock; 12.595 events are listed in the JMA catalogue. Based on
the focal mechanisms of earthquakes in this region, these events were of
strike-slip type with a NW--SE P-axis (Fig. 1b). Only a small coseismic
step was detected by the national GNSS (Global Navigation Satellite
System) network (Fig. S2). Based on the spatiotemporal variation in the
b-value and the migration of the hypocenters, Nanjo et al. (2018)
suggested that fluid movement caused the earthquake sequence in
Kagoshima Bay, but the detailed physical process controlling this
foreshock--mainshock--aftershock sequence remains unclear.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 1} . (\textbf{a} ) Map showing southern Kyushu. Inverted
triangles indicate the seismic stations. We used arrival time data
obtained at both blue and green stations. We analyzed waveform data
recorded at the green stations. The black square shows the study area.
(\textbf{b} ) Hypocenter distribution of earthquakes that occurred in
Kagoshima Bay from January 1, 2003, to April 8, 2018, and their focal
mechanisms. The hypocenters and focal mechanisms were extracted from the
JMA unified catalog. The red square is defined as ``the area surrounding
the mainshock hypocenter'' in this study. The numbers above the focal
mechanisms indicate the JMA magnitude of each earthquake. (\textbf{c} )
M--T diagram and cumulative number of \(M_{\text{JMA}}\ \geq\ 1.0\)earthquakes that
occurred in the area surrounding the mainshock hypocenter (i.e., red
square in Fig. 1b) prior to the mainshock. The vertical red line denotes
the mainshock. (\textbf{d} ) Aftershock occurrence rate of events with a
magnitude\(M_{\text{JMA}}\ \geq 1.0\)(blue) and corresponding M--T diagram (gray).
The inset shows the correlation between the aftershock occurrence rate
and time on a log-log scale. The occurrence rate was estimated by
calculating the reciprocal of the time during which 10 events with
\(M_{\text{JMA}}\ \geq 1.0\)occurred.

In this study, we examined the physical processes that controlled the
M\textsubscript{L} 5.3 Kagoshima Bay earthquake sequence in Kyushu,
southern Japan. First, we determined the hypocenters and focal
mechanisms of this earthquake sequence and delineated the fault
structure. We also estimated the size of the mainshock and examined its
relationship with the fore- and aftershocks to obtain a comprehensive
view of this foreshock--mainshock--aftershock sequence. We then examined
the spatiotemporal characteristics of the intense fore- and aftershocks
to extract information about the aseismic phenomena governing this
earthquake sequence. Finally, by integrating the observations, we
established a model that can be used to explain the occurrence and
characteristics of the foreshock--mainshock--aftershock sequence
associated with the 2017 M5.3 Kagoshima Bay earthquake.



\textbf{2 Methods}

\textbf{2.1 Hypocenter relocations}

We relocated 21.102 events listed in the JMA unified catalogue for the
southern Kagoshima Bay region for the period from March 1, 2003 to April
8, 2018 using the Double-Difference (DD) method (Waldhauser \&
Ellsworth, 2000). This relative relocation method minimizes the
residuals between the observed and theoretical travel time differences
for adjacent earthquake pairs at each station. We applied the DD method
to differential arrival time data, which were estimated from the
waveform cross-correlation, and those listed in the JMA unified catalog.
The procedure was identical to that reported in Yoshida and Hasegawa
(2018a, b), which can be briefly described as follows.

First, we obtained precise differential arrival time data using waveform
cross-correlations. We used the waveform data observed at 20 permanent
seismic stations surrounding the focal area (Fig. 1a; green stations).
At each station, the ground velocity was measured using three-component
short-period seismometers (natural period of 1s) and a sampling rate of
100 Hz. We applied a 5--12 Hz Butterworth filter to the waveforms of
each target event. We used 2.8 and 4.3 s time windows for the P- and
S-waves, respectively, starting 0.3 s before their arrival. The arrival
times were obtained from the JMA unified catalogue. If arrival times
were not available, they were estimated using the one-dimensional
JMA2001velocity model (Ueno et al., 2002) and the hypocenters, and
origin times listed in the JMA unified catalogue. We calculated the
waveform cross-correlations of event pairs with hypocenters within 3 km
from each other and obtained differential arrival times when the
cross-correlation coefficients were greater than 0.8. In total, we
acquired 17.332.318 P-wave differential arrival time data points and
27.738.043 S-wave data points. We also derived the differential arrival
data from the arrival time data listed in the JMA unified catalog:
474.670 data for P-waves and 543.226 data for S-waves. For the
mainshock, only data derived from the JMA unified catalog were used
because of its long duration.

Second, we applied the hypo-DD algorithm (Waldhauser \& Ellsworth, 2000)
to the differential arrival time data. We used a spherical shell
two-layer model (Aki, 1965) for the hypocenter relocation. In this
model, the seismic velocities in each layer are proportional to the
power of the distance from the center of the Earth (Figure S3). The
medium parameters were determined for consistency with the seismic
tomography results obtained in the Kyushu region (Saiga et al., 2010).
We used the hypocenters listed in the JMA unified catalogue for the
initial locations for the relocation. Figures 2a and 3a show the
distribution of these initial hypocenters. Differential arrival time
data were weighted with respect to the square root of the
cross-correlation coefficient. The hypocenters were updated after 50
iterations of the relocation procedure. During the first ten iterations,
a higher weight was assigned to the catalogue data to constrain the
relative locations of large-scale features. In the latter 40 iterations,
a higher weight was assigned to the data derived by the
cross-correlations to delineate shorter-scale features. We evaluated the
uncertainty in the relative hypocenter locations by recalculating the
hypocenters 200 times based on bootstrap resampling of differential
arrival time data.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2.} Maps showing the distribution of the (\textbf{a} )
initial hypocenters listed in the JMA unified catalog and (\textbf{b} )
relocated hypocenters based on the DD method. Blue dots indicate the
locations of the hypocenters. The red lines labeled A to I indicate the
locations of the vertical sections shown in Figure 3.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 3.} Cross-sectional views showing the distribution of the
(\textbf{a} ) initial hypocenters listed in the JMA unified catalog and
(\textbf{b} ) relocated hypocenters based on the DD method. Blue dots
indicate the locations of the hypocenters. The nine figures (A--I)
represent the cross-sectional views along the vertical sections
indicated by the red lines in Fig. 2.



\textbf{2.2 Estimation of focal mechanisms}

We estimated the focal mechanisms based on the amplitude ratios of the
waveforms using the method of Yoshida et al. (2019b), which is similar
to that of Dahm (1996). We used six focal mechanisms determined by the
JMA (Fig. 1b) to represent effects of the path and site on the waveform.
We determined the focal mechanisms of 161 earthquakes with
M\textsubscript{L} [?] 2. We used displacement waveforms obtained by
integrating the velocity waveforms recorded at the 20 stations (green
triangles in Fig. 1a) surrounding the hypocenters. The vertical
component was used for the analysis of the P-wave, whereas radial and
transverse components were used for that of the S-wave. We applied a
2--5 Hz bandpass filter to the waveforms, cutting them out with time
windows of 2.8 s for P-waves and 4.3 s for S-waves starting 0.3 s before
their arrival.

We used waveform cross-correlations to measure the amplitude ratios
between target and reference events. The amplitude ratios were obtained
for pairs with absolute correlation coefficients above 0.75. We used
principal component analysis (PCA) to measure the amplitude ratios.

We only estimated the mechanism solutions when amplitude ratios were
obtained for more than 20 channels. We eliminated the results when the
Variance Reduction (VR) was below 80:

\(\text{VR}=\ \left(1-\frac{\sum_{k=1}^{n}\left(d_{k}-s_{k}\right)^{2}}{\sum_{k=1}^{n}d_{k}^{2}}\right)\ \bullet 100,\)(1)

where \(d_{k}\) and \(s_{k}\) are the observed and
calculated displacement amplitude ratios, respectively, at channel
\(k\).

\textbf{2.3 Estimation of the size of the mainshock source}

We estimated the size of the mainshock source based on the circular
crack source model (e.g., Sato \& Hirasawa, 1973; Madariaga, 1976). In
this source model, the source radius is related to the S-wave corner
frequency, \(f_{c}\), as follows:

\(r=\selectlanguage{greek}\selectlanguage{english}\frac{\text{kβ}}{f_{c}},\) (2)

where \(r\) is the source radius, \(k\) is a
constant, and \(\beta\) is the S-wave velocity close to the
source. Based on a rupture velocity of 0.9\(\beta\),
\(k\) is 0.44 in the model of Sato and Hirasawa (1973) and
0.32 in the model of Madariaga (1976) for P-waves (Kaneko \& Shearer,
2014). Because the estimated source size depends on the source model, we
computed the fault size using both models. We assumed a
\(\beta\) value of 3.4 km/s.

We used the spectral ratio method (e.g., Imanishi \& Ellsworth, 2006) to
estimate the corner frequency of the mainshock. In this method, effects
of the propagation and location on the seismic wave are empirically
removed using the waveforms of an adjacent small earthquake (empirical
Green's function, EGF, event). Based on the assumption that the source
spectrum, that is, \(S_{j}\left(f\right)\), follows the
\(\omega^{2}\)model (Aki, 1967; Brune, 1970), the theoretical ratio
between the velocity spectra of the mainshock, \(v_{i}(f)\), and
the EGF event,\(v_{i}^{\text{egf}}(f)\), at station \(i\) can be
calculated as follows:

\(\text{SSR}_{\text{ij}}\left(f\right)=\frac{v_{i}\left(f\right)}{v_{i}^{\text{egf}}\left(f\right)}=\frac{M_{0}}{M_{0}^{\text{egf}}}\selectlanguage{greek}\frac{R_{\text{θφi}}}\selectlanguage{greek}\selectlanguage{english}{R_{\text{θφi}}^{\text{egf}}}\selectlanguage{english}\frac{1+\left(\frac{f}{f_{c}^{\text{egf}}}\right)^{2}}{1+\left(\frac{f}{f_{c}}\right)^{2}},\)(3)

where \(M_{0}\) and \(M_{0}^{\text{egf}}\) are the seismic
moments of the target earthquake and EGF event, respectively;
\selectlanguage{greek}\(R_{\text{θφij}}\selectlanguage{english}\)and \selectlanguage{greek}\(R_{\text{θφi}}\selectlanguage{english}^{\text{egf}}\) are their radiation patterns
at station\selectlanguage{ngerman} i, respectively; and \(f_{c}^{\text{egf}}\) is
the corner frequency of the EGF event. Based on Eq. (3),
\(f_{c}\) can be estimated from the spectral ratios.

We calculated the spectral ratios by using P-wave velocity waveforms
observed at the 20 stations surrounding the source area (green inverted
triangles in Fig. 1b). The EGF events were earthquakes with M \selectlanguage{english}[?] 2 and a
distance from the mainshock below 1.0 km based on the relocated
hypocenters. The following procedure was performed (Yoshida et al.,
2017):

(1) For the target mainshock and EGF events, the waveforms of the three
components were extracted from a 2.0 s time window starting 0.3 s before
the arrival of the P-wave at each station. The multitaper method
(Thomson, 1982; Prieto et al., 2009) was applied to calculate the
spectra. (2) For channels with EGF observation spectra with a
signal-to-noise ratio \textgreater{} 4 at all frequencies from 0.5 to
30.0 Hz, the spectral ratio between the mainshock and EGF event was
calculated. We used waveforms up to 0.3 s before the arrival of the
P-waves for the noise window. (3) We calculated the geometric mean of
the spectral ratios\(\text{GSR}\left(f\right)\) of all channels at each frequency
for the EGF events, which satisfied the above-mentioned criterion at
five or more stations:\(\text{GSR}\left(f\right)=\prod_{i=1}^{N}{\left(\text{SR}_{i}\left(f\right)\right)^{\frac{1}{N}}\ ,}\)(4) where \(\text{SR}_{i}\left(f\right)\)
is the observed spectral ratio obtained at station i
and \(N\) is the number of stations. (4) By using the grid
search and minimizing the evaluation function\(J\), the
corner frequencies of the mainshock, \(f_{c}\), and EGF event,
\(f_{c}^{\text{egf}}\), were determined:\(J=\sum_{k=1}^{n_{\text{freq}}}\left|\log{\left(\text{GSR}\left(f_{k}\right)\right)-\operatorname{Alog}\left(\text{NSR}\left(f_{k};f_{c},f_{c}^{\text{egf}}\right)\right)}\right|,\)(5)
where\(\text{NSR}\left(f;f_{c},f_{c}^{\text{egf}}\right)=\frac{1+{(f/f_{c}^{\text{egf}})}^{2}}{1+{(f/f_{c})}^{2}}\),\(n_{\text{freq}}\) is the number of
frequencies, and \(f_{k}\) is frequency (at 0.5 Hz intervals
from 0.5 to 30 Hz).\(\ \)The grid search was performed for
\(f_{c}\) and \(f_{c}^{\text{egf}}\) by assuming a range from 0.1
to 100 Hz at 0.1 Hz steps. The amplitude ratio, \(A\), was
estimated using the least squares method for each grid search step.

We applied the spectral ratio method to 33 EGF candidates. We obtained
spectral ratios for 13 EGF events, which satisfy our S/N ratio and data
criteria. Figure S4 shows the spectral ratios of the 13 EGF events.

\textbf{2.4 Detection of aseismic processes from seismicity}

The Epidemic Type Aftershock Sequence (ETAS) model (Ogata, 1988), which
is based on the superposition of the modified Omori law (Ustu, 1961),
can be used to explain mainshock--aftershock seismicity. The ETAS model
assumes that the seismicity rate is the sum of the background rate of
independent events, \(\lambda_{0}\), and aftershocks triggered by
each event, \(\lambda_{i}(t)\):

\(\lambda\left(t\right)=\lambda_{0}+\sum_{i:t_{i}<t}{\lambda_{i}\left(t\right).}\)(6)

Based on the modified Omori law, each earthquake can trigger its own
aftershock sequence (Utsu et al., 1995):

\(\Lambda_{i}\left(t\right)=\frac{K_{0}}{\left(c+t-t_{i}\right)^{p}}e^{\alpha\left(M_{i}-M_{\min}\right)},\)(7)

where \(t_{i}\) is the occurrence time; \(M_{i}\) is
the magnitude of each event, \(i\), that occurred prior to
time \(t\); \(M_{\min}\) is the magnitude of
completeness of the earthquake catalogue;
\(K_{0}\),\(c\), and \(p\) are
constants; and \(t\) is the time that has elapsed since
the main event.

We applied the ETAS model to the seismicity observed after the mainshock
in Kagoshima Bay and investigated the difference between the simulated
and observed seismicity. The results show that the foreshock activity
cannot be explained by the ETAS model, likely because aseismic processes
mainly controlled the foreshock activity. We used the timings and
magnitudes of the earthquakes listed in the JMA catalogue. The lower
limit of the magnitude, M\textsubscript{C}, was set to 1.0. Figure S4
shows the magnitude--frequency distribution. The distribution follows
the Gutenberg--Richter law (Gutenberg \& Richter, 1944) when
M\textsubscript{JMA} \selectlanguage{english}[?] 1.0. The SASeis2006 algorithm by Ogata (2006) was
used to estimate the model parameters and calculate the residuals of the
ETAS model.

\textbf{3 Results}

\textbf{3.1 Fault structure and seismic gap}

We obtained the relocated hypocenters of 20.347 events and the focal
mechanisms of 61 events. Almost all events in the Kagoshima Bay
earthquake sequence can be accurately relocated with the DD algorithm.
The location data for 755 earthquakes were removed because their
hypocenters were located above the ground surface or they contained
outliers in the differential arrival time data. We computed the
differences between the maximum and minimum values in the 95\%
confidence interval of the hypocenter locations (Fig. S6) estimated from
the bootstrap resampling and obtained the medians as a measure of the
estimation error of the relative location: 0.00\(13\) in
longitude, 0.0011 in latitude, and 0.42 km in depth.

Figures 2b and 3b show the distribution of the relocated hypocenters.
Movie 1 shows the animation of the cross-sectional views of the
hypocenters along various lines. Most hypocenters are located within
\textasciitilde{}5 km from the mainshock hypocenter and are distributed
along several planes. These characteristics are in contrast to the
distribution of the initial hypocenters (Figs 2a and 3a), which were
scattered three-dimensionally, similar to a cloud. This significant
change in the hypocenter distribution is due to the improvements of the
relative locations of the hypocenters in this study based on the use of
many accurate differential arrival time data. Similar improvements of
the relative hypocenters, from cloud-like distribution to planar
structures, were previously reported for shallow earthquakes in Japan
based on a similar method and data (e.g., Yoshida \& Hasegawa, 2018a,
b). The cloud-like distribution of the initial hypocenters reflects the
errors in the hypocenter locations in the JMA unified catalog, which are
due to errors in the manual selection.

Figure 4 shows the spatial distribution of the focal mechanisms. Because
the reference focal mechanisms are in the northern part of the source
region (Fig. 1b), newly estimated focal mechanisms are mainly located in
the northern part. The figure shows that the nodal planes of most focal
mechanisms are parallel to the planar structures of the hypocenters,
suggesting that individual small earthquakes occurred on several
macroscopic planes.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure 4.} Estimated focal mechanisms plotted on the hypocenter
distribution. The left figure is a map view and the nine figures (A--I)
on the right are cross-sectional views along vertical sections indicated
by the red lines in the left figure.

Based on Figs 2b, 3b and 4, the fault structures of the 2017 Kagoshima
Bay earthquake sequence are complex, consisting of several subparallel
planes. However, the distribution of the hypocenters was relatively
simple before the mainshock. Figure 5 shows an enlarged view of the
spatial distribution of the hypocenters of the foreshocks (red dots).
Most hypocenters are evenly distributed in one plane, with a strike
parallel to those of the nodal planes of the focal mechanisms of the
mainshock and individual small earthquakes, suggesting that the
mainshock and most of the foreshocks occurred on this plane.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 5.} Hypocenter distribution of the foreshocks. Red and
blue circles represent the hypocenters of foreshocks and aftershocks,
respectively, Blue circles represent the hypocenters of the precursory
activity. (\textbf{a} ) Map showing the hypocenters of the precursor
activity and aftershocks. (\textbf{b} ) Map showing only the hypocenters
of the precursory activity. The broken ellipse indicates the location of
the seismic gap. (\textbf{c} ) Cross-sectional views along vertical
sections A to I shown in (\textbf{a} ). The yellow star indicates the
hypocenter location of the mainshock.

The hypocenters of the foreshocks are not uniformly distributed in the
plane, but they are distributed in form of a doughnut, that is, a
seismic gap forms in the center of the plane (broken ellipse in Fig.
5b). To demonstrate this distribution, we estimated the lateral
distribution of the moment release on the fault (Fig. 6a) during the
foreshock sequence following Yoshida et al. (2020a). We computed the
seismic moment release of each earthquake by assuming that its magnitude
is equal to the moment magnitude. Subsequently, we summed the moment
release values of the points that were evenly spaced every 0.04 km by
using the earthquakes within the nearest grid cell. The result shown in
Fig. 6b indicates that the moment release of the foreshock sequence is
smaller (\({<10}^{11}\text{\ Nm}\)) in the region corresponding to the seismic
gap than in the surrounding region (\({>10}^{11}\text{\ Nm}\)). The hypocenter
of the mainshock is located at the edge of this seismic gap. Figure 5
shows a comparison of the hypocenters of the foreshocks and aftershocks.
Although aftershocks occur inside the seismic gap based on the map (Fig.
5a), they actually occur in shallower areas than the foreshocks (Fig.
5c), that is, not within the seismic gap of the foreshocks.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{Figure 6.} Seismic gap of earthquakes in the foreshock period
projected on the dominant plane. (\textbf{a} ) Comparison of the size of
the seismic gap with the estimated sizes of the faults associated with
the mainshock. Blue circles indicate the hypocenters of the precursory
activity corresponding to the fault sizes by assuming stress drops of 10
MPa. The yellow star indicates the location of the mainshock hypocenter.
The red circle represents the size of the fault corresponding to the
mainshock estimated based on the model of Sato and Hirasawa (1973). The
green circle represents the size of the fault corresponding to the
mainshock estimated based on the model of Madariaga (1976). (\textbf{b}
) Moment release amount (color scale) computed for each 0.04 km grid
cell. The broken ellipse represents the seismic gap.

The median value of the estimated corner frequencies of the mainshock is
1.9 Hz (Fig. S4). The first and third quartiles are 1.8 and 2.5 Hz,
respectively. Based on the median corner frequency and the models
proposed by Sato and Hirasawa (1973) and Madariaga (1976), the source
radius of the mainshock is 787 and 572 m, respectively. It is much
smaller than the foreshock and aftershock regions but comparable to the
seismic gap (Fig. 6a).

\textbf{3.2 Foreshock and aftershock migration}

Figures 7a--c show the occurrences of the foreshocks on a color scale.
In Figs 8a--c, the occurrence of each earthquake is compared with the
longitude, latitude, and depth, respectively, to illustrate their
migration behavior. In the longitudinal direction (Fig. 8a), the
hypocenters expand nearly symmetrically in the first 230 days of
foreshock activity and concentrate in the east close to the hypocenter
of the mainshock during the last \textasciitilde{}70 days. In the
latitudinal direction (Fig. 8b), the hypocenters migrate from north to
south. In the depth direction (Fig. 8c), the hypocenters migrate both in
the shallow and deep directions, indicating that most earthquakes
occurred in the deeper part surrounding the mainshock hypocenter during
the last \textasciitilde{}70 days of activity.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{Figure 7.} Spatiotemporal evolution of the hypocenters
(\textbf{a} )--(\textbf{c} ) before and (\textbf{d)--(f} ) after the
mainshock. (\textbf{a} ) Projection of the hypocenters, (\textbf{d} )
map view, (\textbf{b} ) and (\textbf{e} ) east--west cross section, and
(\textbf{c} ) and (\textbf{f} ) north--south cross section. The symbol
sizes corresponds to the JMA magnitude. The hypocenters are colored
according to their occurrence time measured relative to that of the
mainshock, that is, the mainshock occurred at time 0 and negative and
positive numbers denote the days before and after the mainshock,
respectively.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}

\textbf{Figure 8.} Temporal evolution of hypocenters in the (\textbf{a}
) latitude, (\textbf{b)} longitude, and (\textbf{c} ) depth directions.
(\textbf{d} ) Temporal evolution of the aftershock hypocenters in the
depth direction. The red crosses and yellow curve in (\textbf{d} )
indicate the depth above which the shallowest 10\% of the hypocenters
are located (D10) for every bin with 400 events based on the occurrence
time. The circle size corresponds to the JMA magnitude. The yellow star
indicates the hypocenter of the mainshock.

Figures 7d--f show the distribution of aftershock hypocenters colored
based on the occurrence time of each event. Figure 8d shows the temporal
evolution of the aftershock hypocenters as a function of the depth. The
temporal evolution of the aftershock hypocenters in both the latitudinal
and longitudinal directions is shown in Fig. S7. Because the spatial
distribution of the aftershocks is complex, the spatiotemporal features
of the aftershocks are more difficult to determine than those of the
foreshocks. Overall, the aftershock hypocenters move upward with time,
as shown in Fig. 8d, which depicts the depths above which the shallowest
10\% of the hypocenters are located (D10) for each bin containing 400
events, as denoted by the red curve. Although earthquakes occur in a
relatively deep region immediately after the mainshock, the upper limit
of the seismic depth (D10) gradually moves in the shallow direction,
that is, the hypocenters gradually move to the shallower part with time
after the mainshock.

\textbf{3.3 Deviation of the seismicity from Omori's law}

We investigated the seismicity rate of the Kagoshima Bay earthquake
sequence after the mainshock. Figure 1d shows the seismicity rates of
the M\textsubscript{JMA} [?]1.0 events in the area surrounding the
hypocenter of the mainshock (red frame in Fig. 1b). The seismic rate was
obtained by calculating the reciprocal of the time required to generate
ten earthquakes that were arranged in chronological order. Based on Fig.
1d, the seismicity rate decreases by the power of the elapsed time
immediately after the mainshock, as described by the modified Omori law
(Utsu, 1961). The seismicity rate abruptly increases \textasciitilde{}44
days after the mainshock, which corresponds to the occurrence of the
largest aftershock (M\textsubscript{L} 4.4), suggesting that the
increase is due to secondary aftershocks. A period with a high
seismicity rate started approximately 20 to 40 days after the mainshock;
the seismic activity was temporarily strong despite the absence of large
aftershocks.

Based on maximum likelihood estimation, we obtained the following
parameters for the ETAS model: \(K_{0}=\ 34.205\),\(c=1.3163\times 10^{-2}\),
\(p=1.0685\), \(\alpha=1.5078\), and\(\mu=2.9603\times 10^{-2}\). Based on
Ogata (1992), the range of \selectlanguage{greek}α-\selectlanguage{english}values is 0.35--0.85 for swarm seismicity
and 1.2-- 3.1 for non-swarm seismicity. The \selectlanguage{greek}α \selectlanguage{english}value estimated for the
seismic activity in Kagoshima Bay is within the latter range.

In Fig. 9, the cumulative number of earthquakes simulated using the
estimated model parameters is compared with the observations. Overall,
the number predicted based on the ETAS model matches the observations.
However, the simulated number of earthquakes is lower than the observed
number 20--40 days after the mainshock. To quantitatively examine the
magnitude of the discrepancy between the model and observations, we
performed residual analysis using the transformed time, similar to Ogata
(1988). Figure 9c shows that the discrepancy between the model and
observations is high at a transformed time between 1.000 and 1.500,
corresponding to the period of 20--40 days after the mainshock. This
deviation is significant at the 95\% significance level based on the
assumption of a uniform distribution.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}

\textbf{Figure 9.} (\textbf{a} ) M--T diagram. (\textbf{b} ) Observed
cumulative number of aftershocks with \(M_{L}\ \geq 1.0\) (red solid
line) and predicted number based on the estimated ETAS parameters (blue
solid line). Each curve represents the cumulative numbers starting 0.1
days after the mainshock. (\textbf{c} ) Results of the residual
analysis, where the blue solid line shows the observed events with
respect to the transformed time on the horizontal axis and cumulative
number of observed \(M_{L}\ \geq 1.0\) earthquakes on the vertical axis.
The black dotted line represents the transformed time at which the
assumed model fully matches the observation. The red solid and red
broken lines indicate the two-sided 95\% and 99\% error bounds of the
Kolmogorov--Smirnov statistic, respectively. The gray zone in
(\textbf{c} ) shows the range of the transformed time corresponding to
the period of 20--44 days after the mainshock highlighted in gray in
(\textbf{b} ).

The large discrepancy between the predicted and observed seismicity
rates 20--40 days (\textasciitilde{}1.000--1.500 in Fig. 9c) after the
mainshock can be explained by a temporary increase in the background
seismicity, which was assumed to be constant over the entire period of
this analysis in the model. The transient increase in the background
seismicity rate suggests that the Kagoshima Bay earthquake sequence may
have been affected by physical processes other than
earthquake-to-earthquake interactions, especially during this period
(20--40 days) and that an aseismic process may have led to the largest
aftershock (M\textsubscript{L} 4.4) that occurred 44 days after the
mainshock. Contrarily, most aftershocks can be explained as general
mainshock--aftershock seismic activity, suggesting that stress changes
caused by the mainshock resulted in numerous aftershocks.



\textbf{4 Discussion}

Our results show that: (1) the foreshocks of the 2017 M5.3 Kagoshima Bay
earthquake sequence occurred on a single plane with a steep dip to the
east, whereas aftershocks occurred on several more complex planar
structures, (2) the foreshock hypocenters formed a seismic gap, and (3)
the foreshock and aftershock hypocenters exhibit clear migration
behaviors. In this section, we integrated these observations and propose
a simple model that can explain the occurrence of the
foreshock--mainshock--aftershock sequence of the 2017 M5.3 Kagoshima Bay
earthquake based on the upward fluid movement, which is similar to the
fault-valve model proposed by Sibson (1992).

\textbf{4.1 Migration of foreshock activity along a plane}

The clear hypocenter migration observed for the foreshock sequence
suggests that aseismic physical processes controlled this sequence. In
fact, the seismicity rate of the foreshock sequence could not be
reproduced by the ETAS model, suggesting that the
earthquake-to-earthquake interaction cannot explain this sequence. Thus,
the foreshock sequence must be understood as temporary increase in the
background seismicity rate, similar to that of the earthquake swarm.

In Fig. 10, the distances of the foreshock hypocenters from the mean
location of the first three events are plotted against time. The
expansion front of the pore pressure diffusion model reported in Shapiro
et al. (1997) is also shown, which can be expressed by the following
equation including various diffusion coefficients \(D_{h}\):

\begin{equation}
\begin{matrix}r=\sqrt{4\pi D_{h}t},\#\left(9\right)\\
\end{matrix}\nonumber \\
 \end{equation}

where \(r\) is the distance from the point pressure source
and \(t\) is the time. In this study, we set the initiation
time to 220 days before the mainshock because the seismicity rate
significantly increased at this time (Fig. 1c). We also show the
propagation fronts of the linear spread model that has been used for
aseismic slip propagation in the past (e.g., Vidale \& Shearer, 2006).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}

\textbf{Figure 10.} Temporal evolution of the distances between the
foreshocks and initial hypocenter. Blue circles represent the
hypocenters expressed by the size corresponding to the JMA magnitude.
The black curves show the fluid diffusion models with
D\textsubscript{h}= 0.01, 0.03, and 0.05 m\textsuperscript{2}/s. Gray
straight lines show the linear spread model with migration speeds of d =
0.001, 0.003, and 0.005 km/h.

The pore pressure diffusion model with a hydraulic diffusion coefficient
of \textasciitilde{}0.05 m\textsuperscript{2}/s matches the observations
better than the linear spread model. In previous studies, it has been
estimated that the hydraulic diffusion coefficient in the crust ranges
from \textasciitilde{}0.01--10 m\textsuperscript{2}/s (e.g., Talwani et
al., 2007; Shelly et al., 2016; Yoshida \& Hasegawa, 2018a), which is
similar to the foreshock migration speed of the M5.3 Kagoshima Bay
earthquake sequence. Based on the linear spread model, the propagation
velocity is \textasciitilde{}0.001--0.005 km/h. Based on previous
studies, the migration speed of aseismic slip propagation ranges from
0.1--1.0 km/h (e.g., Lohman \& McGuire, 2007; Kato et al., 2016), which
is significantly higher than the migration speed of the present
foreshock activity. If we advance the initiation timing of propagation,
the propagation speed decreases. Thus, according to the migration speed
and spatiotemporal pattern of the foreshocks, the pore pressure
diffusion model better explains the overall migration of the foreshock
hypocenters.

Aseismic creep related to the nucleation process of the mainshock might
be involved in the migration of the foreshocks. In fact, physical
simulations indicate interseismic creep in seismogenic patches from
external stable-slip regions before the occurrence of unstable slip (Tse
\& Rice, 1986). Such an expansion of quasi-static slip prior to the
mainshock may explain the current migration of the foreshocks (e.g.,
Dodge et al., 1996; Yabe \& Ide, 2018). However, the source of the
mainshock is smaller than the foreshock area (Fig. 6a), which
contradicts the hypothesis because the mainshock rupture zone should
include the nucleation area. Note that the size of the source of the
mainshock was estimated based on source models assuming a constant
rupture velocity (subshear rupture propagation). If the assumptions
differ from reality (e.g., supershear rupture propagation), the source
size may differ from our estimation, explaining this contradiction.
However, the aftershocks migrate upward on multiple planes (Fig. 8d),
which can be explained with the foreshock sequence if the pore pressure
migration model is adopted. Thus, we prefer the hypothesis that pore
pressure migration is primarily responsible for the generation of the
2017 M5.3 Kagoshima Bay earthquake sequence. The heterogeneity in the
permeability and/or pore pressure along the fault may explain the up-
and downward movement of the hypocenter along the plane (Fig. 8c).

However, recent observations of fluid injection-induced seismicity and
natural earthquake swarms suggest that an increase in the pore pressure
can cause aseismic slip (Cornet et al., 1997; Guglielmi et al., 2015;
Ruhl et al., 2016; Yoshida \& Hasegawa, 2018a; De Barros et al., 2020).
In the presence of fluids, the effective normal stress decreases and the
critical nucleation size increases; thus, the occurrence of aseismic
slip is likely (e.g., Scholz, 1998). The increase in the pore pressure
also accelerates creep in the stable-slip segment of the fault. Both
aseismic slip and fluid movement may have contributed to the occurrence
of foreshocks. Furthermore, the poroelastic effects associated with pore
pressure migration (Segall, 1989; Goebel et al., 2018) and the
earthquake-to-earthquake interaction (Helmstetter, 2002) may contribute
to the occurrence of earthquakes.

\textbf{4.2 Seismic gap of the foreshock and aftershock sequence in the
mainshock fault plane}

The doughnut-like pattern of the foreshocks (Fig. 6) is similar to the
``Mogi doughnut'' (Mogi, 1969). It has been reported that aftershocks do
not occur in the mainshock region (e.g., Mendoza \& Hartzell, 1988; Das
\& Henry, 2003; Woessner et al., 2006; Asano et al., 2011; Ebel \&
Chambers, 2016; Yoshida et al., 2016b and 2020a; Ross, 2017b, 2018;
Wetzler et al., 2018), which is likely because the shear stress was
released during the mainshock. In Figure 6a, the size of the seismic gap
is compared with the estimated size of the fault related to the
mainshock. Because the centroid location of the mainshock was not
determined, we assumed the centroid is located by a few hundred meters
into the shallower region such that the mainshock centroid is located in
the center of the seismic gap shown in the figure. The fault size of the
mainshock is similar to that of the seismic gap. This is consistent with
the hypothesis that the mainshock rupture occurred in the seismic gap of
the foreshock and aftershock activities. A similar spatial separation of
the mainshock and fore- and aftershocks in the rupture area was also
reported for a recent M5.2 intraplate earthquake that occurred in Akita,
NE Japan (Yoshida et al., 2020).

The seismic gap in the foreshock activity may originate from the spatial
heterogeneities in the frictional and material properties along the
fault plane. The fault strength of the mainshock rupture area may have
been higher than that of the surrounding area, as proposed in the
asperity model of Lay \& Kanamori (1981). Foreshocks can be understood
as failures of small seismogenic patches in the surrounding stable area.
Alternatively, the area may have been covered by an impermeable medium,
hindering fluid intrusion. The occurrences of foreshocks and aseismic
slip increased the shear stress in the future source region of the
mainshock rupture. The mainshock occurred in this region due to the
gradually increasing pore pressure and shear stress.

\textbf{4.3 Upward migration of the aftershocks along several planes}

The aftershock sequence of the Kagoshima Bay earthquake sequence follows
Omori's law (Fig. 1d), suggesting that this sequence was triggered by
the M5.3 mainshock. However, the aftershock sequence slightly deviates
from the prediction based on the ETAS model (Figs 9b--c). The transient
increase in the background seismicity rate suggests that the Kagoshima
Bay earthquake sequence may have been affected by physical processes
other than earthquake-to-earthquake interactions, especially during this
period (20 to 40 days) and that these aseismic processes may have led to
the largest aftershock (M\textsubscript{L} 4.4) that occurred 44 days
after the mainshock. Based on the model simulations of fluid
injection-induced seismicity, Hainzl and Ogata (2005) pointed out that
the background seismicity rate of the ETAS model is sensitive to the
amount of injected water. In previous studies, similar observations were
made for fluid injection-induced seismicity and natural earthquake
sequences (Llenos \& Michael, 2013; Yoshida \& Hasegawa, 2018b; Kumazawa
et al., 2019).

Our results indicate that the aftershock hypocenters migrated toward the
shallower portion on multiple planes. Such upward movements of
hypocenters have been previously reported for earthquake swarms
following nearby large earthquakes and it has been concluded that they
reflect the upward pore pressure migration associated with the
fault-valve behavior (Shelly et al., 2015; Ruhl et al., 2016; Yoshida \&
Hasegawa, 2018a, b). Examples are the earthquake swarms that occurred in
northeastern Japan following the 2011 Tohoku-Oki earthquake (Yoshida et
al., 2016a; Yoshida \& Hasegawa, 2018a, b). The earthquake swarms might
originate from the pore pressure increase because (1) they occurred in
the stress shadow of the 2011 Tohoku-Oki earthquake with a time delay of
a few weeks despite the reduction in the shear stress, (2) they are
located beneath the caldera structures that are believed to host shallow
igneous bodies, with hydrothermal fluids immediately below, (3) they are
located a few kilometers above S-wave reflectors and the low-velocity
zone including fluids, and (4) their hypocenters migrate upward (Yoshida
\& Hasegawa, 2018a; Yoshida et al., 2019a). The Kagoshima Bay swarm was
also located beneath an ancient caldera and involved the upward
migration of aftershocks, which can be explained by an increase in the
pore pressure. Fluid paths in the crust may have expanded due to the
deformation and shaking associated with the mainshock. Pore pressure
migration may explain deviations in the seismicity rate from Omori's
law. These observations are consistent with the prediction based on the
fault-valve model proposed in Sibson (1992), that is, upward fluid
discharge after the mainshock. In recent geodetic studies, a porosity
wave associated with the fault-valve action was detected (Rossi et al.,
2016 and 2018).

We presume that the subducting Philippine Sea Plate is the source of
fluids, similar to the model reported in Hasegawa et al. (2005), which
is based on the geophysical and geological observations in northeastern
Japan. This hypothesis is supported by seismic data obtained in Kyushu
using tomography, which indicate that the existence of an inclined
low-velocity layer continuously distributed in the mantle wedge and
reaching right below the volcanic front as northeastern Japan (Zhao et
al., 2012). The low-velocity zone is considered to represent the
ascending flow portion of the secondary convection within the mantle
wedge and therefore contains fluids from the slab and resultant melts
(Hasegawa et al., 2005). The buoyancy facilitated the upward migration
of the fluids, as shown in simulations (e.g., Iwamori, 1998; Wada et
al., 2015; Horiuchi et al., 2016), and the fluids reached the source
region of the Kagoshima Bay sequence.

\textbf{4.4 Comprehensive interpretation of the seismic activity in
Kagoshima Bay}

Here, we summarize our simple model that comprehensively explains the
observed results of the foreshock-mainshock-aftershock sequence of the
2017 M5.3 Kagoshima Bay earthquake.

First, fluids that have infiltrated the mainshock fault plane caused the
foreshock activity. The hypocenter migration of the foreshock activity
can be interpreted considered to be a reflection of fluid movement and
possibly triggered aseismic slip on the plane. Second, the occurrences
of foreshocks and aseismic slip increased the shear stress in the future
source region of the mainshock (seismic gap in Fig. 6). The mainshock
finally occurred in this region due to the gradually increasing pore
pressure and shear stress. Third, the change in the stress associated
with the occurrence of the mainshock primarily triggered aftershocks in
the area surrounding the mainshock including regions outside the
mainshock fault plane. Fluids started to move upward due to the
deformation and shaking associated with the mainshock. Together with the
fluids, the aftershock hypocenters moved to shallower regions. Thus, the
overall sequence of the 2017 M5.3 Kagoshima Bay earthquake can be
explained by consistent upward fluid movement.

\textbf{4.5 Implications to the foreshock--mainshock--aftershock
sequence}

The results of previous studies suggested that many earthquake swarms
are caused by the movements of crustal fluids (e.g., Mogi, 1989;
Italiano et al., 2001; Fischer and Hor\selectlanguage{ngerman}álek, 2003; Parotidis et al.,
2003; Bianco et al., 2004; Yukutake et al., 2011; Shelly et al., 2016;
Yoshida et al., 2016a; Ruhl et al., 2016; De Barros et al., 2019). The
results of the present study suggest that the generation mechanism of
the foreshock activity is the same as that of earthquake swarms, that
is, a temporary increase in background seismicity rate due to increasing
pore pressure and aseismic slip. The whole sequence of the Kagoshima Bay
seismicity can be understood as the transition from swarm activity to
the mainshock--aftershock sequence.

The 2008 Mogul earthquake swarm, Nevada, may be a similar example. This
sequence was also initiated by swarm activity but shifted to a
mainshock--aftershock sequence after the occurrence of the M4.9
mainshock. The upward migration of the earthquakes suggests that
fault-valve behavior is involved in the occurrence of this earthquake
sequence (Ruhl et al., 2016). The aftershock activities of the 2014
M\textsubscript{L} 4.8, Ubaye earthquake, France (De Barros et al.,
2019), and the foreshock and aftershock activities of the 2017 M5.2
Akita-Daisen event can also be understood as transitions from swarm
activity to mainshock--aftershock sequences (Yoshida et al., 2020b).
Similarly, aseismic slip may have caused the foreshocks and mainshock of
the 2011 M9 Tohoku-Oki earthquake (Kato et al., 2012); 2014 Iquique Mw
8.1 earthquake, Chile (Kato \& Nakagawa, 2014); and 2009 M6.3 L'Aquila
earthquake (Borghi et al., 2016). It is likely that pore pressure
migration and aseismic slip propagation occasionally coexist (Waite \&
Smith, 2002; Ross et al., 2017a; Yoshida \& Hasegawa, 2018; De Barros et
al., 2020) and contribute to the increase in the background seismicity
rate. Such aseismic processes may also cause mainshock--aftershock
activity without notable foreshocks. The 2019 M6.7 Yamagata-Oki
earthquake, NE Japan, may be an example. The earthquake occurred in the
stress shadow of the 2011 Tohoku-Oki earthquake and exhibited an upward
aftershock migration (Yoshida et al., 2020b). These observations suggest
that the monitoring of aseismic processes is crucial to understanding
the seismic activity.

\textbf{5 Conclusions}

The results of previous studies suggested that many earthquake swarms
have been caused by the movement of crustal fluids (e.g., Mogi, 1989;
Fischer and Horálek, 2003; Parotidis et al., 2003; Bianco et al., 2004;
Yukutake et al., 2011; Chen et al., 2012; Shelly et al., 2016; Yoshida
et al., 2016a; Ruhl et al., 2016; De Barros et al., 2019). In the
present study, the intense foreshock--mainshock--aftershock sequence of
the 2017 M5.3 Kagoshima Bay earthquake was examined. The results show
that the whole sequence can be explained by upward fluid movement: (1)
most foreshocks were located on a single plane with a steep dip to the
east and migrated along the plane. This foreshock migration can be
interpreted as a reflection of fluid movement and possibly triggered
aseismic slip on the plane; (2) The hypocenter of the mainshock was
located at the edge of a seismic gap with a size comparable to that of
the source of the mainshock rupture. This suggests that the mainshock
rupture was due to the slip of this seismic gap and the seismic gap was
a large seismogenic patch with higher fault strength, which finally
ruptured due to the increase in the pore pressure and aseismic slip in
the surrounding areas; and (3) Aftershocks occurred on several planes
with a steep dip to the east and moved from deeper to shallower regions.
The upward migration can be interpreted as a reflection of post-failure
fluid discharge. Thus, the overall sequence of the 2017 M5.3 Kagoshima
Bay earthquake can be explained by upward fluid movement, as presumed by
the fault-valve model (Sibson, 1992).

The whole sequence of the Kagoshima Bay seismicity can be understood as
the transition from swarm activity to a mainshock--aftershock sequence.
The results of the present study suggest that the generation mechanism
of the foreshock activity is the same as that of the earthquake swarms,
that is, a temporary increase in the background seismicity rate due to
increasing pore pressure and aseismic slip. Aseismic processes sometimes
cause a large earthquakes that is followed by numerous aftershocks; the
foreshock--mainshock--aftershock sequence may such a sequence.

\textbf{Acknowledgments}

We thank the editor Rachel Abercrombie, Giuliana Rossi, and an anonymous
reviewer for their constructive comments, which improved the manuscript.
KY thanks Dapeng Zhao for the information of seismic velocity structure
in Kyushu. In this study, hypocenter and P- and S-wave data reported in
the unified catalog of the JMA
(\emph{https://www.data.jma.go.jp/svd/eqev/data/bulletin/index\_e.html})
were used. The seismograms were collected and stored by the JMA,
National Universities, and National Research Institute for Earth Science
and Disaster Resilience (http://www.hinet.bosai.go.jp/?LANG=en). The
figures in this paper were created using GMT (Wessel \& Smith, 1998).
This research was supported by JSPS KAKENHI Grant Number JP 17K1437. The
results obtained for the hypocenters, focal mechanisms, and coseismic
slip distribution are available at
http://www.aob.gp.tohoku.ac.jp/\textasciitilde{}yoshida/pub/JGR2020b/.

\textbf{References}

Aki, K. (1965). A computer program for precise determination of focal
mechanism of local earthquakes by revising focal depths and crust-mantle
structure. \emph{Bulletin of the Earthquake Research Institute}
,\emph{43} , 15--22.

Aki, K. (1967). Scaling law of seismic spectrum. \emph{Journal of
Geophysical Research} , \emph{72} (4),
1217--1231.\emph{https://doi.org/10.1029/JZ072i004p01217}

Asano, Y., Saito, T., Ito, Y., Shiomi, K., Hirose, H., Matsumoto, T.,
.\selectlanguage{english} .\selectlanguage{english} . Sekiguchi, S. (2011). Spatial distribution and focal mechanisms of
aftershocks of the 2011 off the Pacific coast of Tohoku
Earthquake.\emph{Earth, Planets and Space} , \emph{63} (7),
669--673.\emph{https://doi.org/10.5047/eps.2011.06.016}

Bianco, F., Del Pezzo, E., Saccorotti, G., \& Ventura, G. (2004). The
role of hydrothermal fluids in triggering the July-August 2000 seismic
swarm at Campi Flegrei, Italy: Evidence from seismological and
mesostructural data. \emph{Journal of Volcanology and Geothermal
Research} , \emph{133} (1--4), 229--246.
https://doi.org/10.1016/S0377-0273(03)00400-1

Borghi, A., Aoudia, A., Javed, F., \& Barzaghi, R. (2016). Precursory
slow-slip loaded the 2009 L'Aquila earthquake sequence.\emph{Geophysical
Journal International} , \emph{205} (2), 776--784.
https://doi.org/10.1093/gji/ggw046

Brune, J. N. (1970). Tectonic stress and the spectra of seismic shear
waves from earthquakes. \emph{Journal of Geophysical Research}
,\emph{75} (26),
4997--5009.\emph{https://doi.org/10.1029/JB075i026p04997}

Cornet, F. H., Helm, J., Poitrenaud, H., \& Etchecopar, A. (1997).
Seismic and aseismic slips induced by large-scale fluid injections.
In.\emph{Pure and Applied Geophysics} . Berlin: Springer, 150(3--4),
(563--583).\emph{https://doi.org/10.1007/s000240050093}

Cox, S. F. (2016). Injection-driven swarm seismicity and permeability
enhancement: Implications for the dynamics of hydrothermal ore systems
in high fluid-flux, overpressured faulting regimes--an invited
paper.\emph{Economic Geology} , \emph{111} (3),
559--587.\emph{https://doi.org/10.2113/econgeo.111.3.559}

Dahm, T. (1996). Relative moment tensor inversion based on ray theory:
Theory and synthetic tests. \emph{Geophysical Journal International}
,\emph{124} (1),
245--257.\emph{https://doi.org/10.1111/j.1365-246X.1996.tb06368.x}

Das, S., \& Henry, C. (2003). Spatial relation between main earthquake
slip and its aftershock distribution. \emph{Reviews of Geophysics}
,\emph{41} (3).\emph{https://doi.org/10.1029/2002RG000119}

De Barros, L., Baques, M., Godano, M., Helmstetter, A., Deschamps, A.,
Larroque, C., \& Courboulex, F. (2019). Fluid-Induced Swarms and
Coseismic Stress Transfer: A Dual Process Highlighted in the Aftershock
Sequence of the 7 April 2014 Earthquake (Ml 4.8, Ubaye,
France).\emph{Journal of Geophysical Research: Solid Earth} , \emph{124}
(4), 3918--3932. https://doi.org/10.1029/2018JB017226

De Barros, L., Cappa, F., Deschamps, A., \& Dublanchet, P. (2020).
Imbricated aseismic slip and fluid diffusion drive a seismic swarm in
the Corinth Gulf, Greece. \emph{Geophysical Research Letters} ,
47(9).\emph{https://doi.org/10.1029/2020GL087142}

Dodge, D. A., Beroza, G. C., \& Ellsworth, W. L. (1996). Detailed
observations of California foreshock sequences: Implications for the
earthquake initiation process. \emph{Journal of Geophysical Research:
Solid Earth} , \emph{101} (B10),
22371--22392.\emph{https://doi.org/10.1029/96JB02269}

Ebel, J. E., \& Chambers, D. W. (2016). Using the locations of M \selectlanguage{english}[?] 4
earthquakes to delineate the extents of the ruptures of past major
earthquakes. \emph{Geophysical Journal International} , \emph{207} (2),
862--875.\emph{https://doi.org/10.1093/gji/ggw312}

Ellsworth, W. L. (2013). Injection-induced earthquakes. \emph{Science}
,\emph{341} (6142),
1225942.\emph{https://doi.org/10.1126/science.1225942}

Fischer, T., \& Horalek, J. (2003). Space-time distribution of
earthquake swarms in the principal focal zone of the NW Bohemia/Vogtland
seismoactive region: Period 1985-2001. \emph{Journal of Geodynamics}
,\emph{35} (1--2), 125--144.
https://doi.org/10.1016/S0264-3707(02)00058-3

Goebel, T. H. W., Weingarten, M., Chen, X., Haffener, J., \& Brodsky, E.
E. (2017). The 2016 Mw5.1 Fairview, Oklahoma earthquakes: Evidence for
long-range poroelastic triggering at \textgreater{}40 km from fluid
disposal wells. \emph{Earth and Planetary Science Letters} , \emph{472}
, 50--61. https://doi.org/10.1016/j.epsl.2017.05.011

Guglielmi, Y., Cappa, F., Avouac, J. P., Henry, P., \& Elsworth, D.
(2015). INDUCED SEISMICITY. Seismicity triggered by fluid
injection-induced aseismic slip. \emph{Science} , \emph{348} (6240),
1224--1226.\emph{https://doi.org/10.1126/science.aab0476}

Gutenberg, B., \& Richter, C. F. (1944). Frequency of earthquakes in
California. \emph{Bulletin of the Seismological Society of America}
,\emph{34} (4), 185--188.

Hainzl, S., \& Ogata, Y. (2005). Detecting fluid signals in seismicity
data through statistical earthquake modeling. \emph{Journal of
Geophysical Research} , \emph{110} (B5),
1--10.\emph{https://doi.org/10.1029/2004JB003247}

Hasegawa, A. (2017). Role of H\textsubscript{2}O in generating
subduction zone earthquakes. \emph{Monographs on Environment, Earth and
Planets} , \emph{5} (1),
1--34.\emph{https://doi.org/10.5047/meep.2017.00501.0001}

Hasegawa, A., Nakajima, J., Umino, N., \& Miura, S. (2005). Deep
structure of the northeastern Japan arc and its implications for crustal
deformation and shallow seismic activity. \emph{Tectonophysics}
,\emph{403} (1--4),
59--75.\emph{https://doi.org/10.1016/j.tecto.2005.03.018}

Helmstetter, A., \& Sornette, D. (2002). Diffusion of epicenters of
earthquake aftershocks, Omori's law, and generalized continuous-time
random walk models. \emph{Physical Review E - Statistical Physics,
Plasmas, Fluids, and Related Interdisciplinary Topics} , \emph{66} (6),
24. https://doi.org/10.1103/PhysRevE.66.061104

Horiuchi, S. S., \& Iwamori, H. (2016). A consistent model for fluid
distribution, viscosity distribution, and flow-thermal structure in
subduction zone. \emph{Journal of Geophysical Research: Solid Earth}
,\emph{121} (5), 3238--3260. https://doi.org/10.1002/2015JB012384

Hubbert, M. K., \& Rubey, W. W. (1959). \emph{Role of fluid overpressure
in the mechanics of overthrust faulting} \emph{{[}Geological Society of
America bulletin{]}} , \emph{70} (pp.
167--206).\emph{https://doi.org/10.1130/0016-7606(1959)70}

Imanishi, K., \& Ellsworth, W. L. (2006). Source scaling relationships
of microearthquakes at Parkfield, CA, determined using the SAFOD pilot
hole seismic array \emph{{[}Geophysical monograph series{]}}
.\emph{Geophysical Monograph Series} ,
81--90.\emph{https://doi.org/10.1029/170GM10}

Italiano, F., Martinelli, G., \& Nuccio, P. M. (2001). Anomalies of
mantle-derived helium during the 1997-1998 seismic swarm of
Umbria-Marche, Italy. \emph{Geophysical Research Letters} , \emph{28}
(5), 839--842. https://doi.org/10.1029/2000GL012059

Iwamori, H. (1998). Transportation of H2O and melting in subduction
zones. \emph{Earth and Planetary Science Letters} , \emph{160} (1--2),
65--80. https://doi.org/10.1016/S0012-821X(98)00080-6

Kaneko, Y., \& Shearer, P. M. (2014). Seismic source spectra and
estimated stress drop derived from cohesive-zone models of circular
subshear rupture. \emph{Geophysical Journal International} ,\emph{197}
(2), 1002--1015. https://doi.org/10.1093/gji/ggu030

Kato, A., Obara, K., Igarashi, T., Tsuruoka, H., Nakagawa, S., \&
Hirata, N. (2012). Propagation of Slow Slip Leading Up to the 2011 Mw
9.0 Tohoku-Oki Earthquake. \emph{Science} , \emph{335} (6069), 705--708.
https://doi.org/10.1126/science.1215141

Kato, A., \& Nakagawa, S. (2014). Multiple slow-slip events during a
foreshock sequence of the 2014 Iquique, Chile Mw 8.1
earthquake.\emph{Geophysical Research Letters} , \emph{41} (15),
5420--5427.

Kato, A., Fukuda, J., Nakagawa, S., \& Obara, K. (2016). Foreshock
migration preceding the 2016 MW 7.0 Kumamoto earthquake,
Japan.\emph{Geophysical Research Letters} , \emph{43} (17),
8945--8953.\emph{https://doi.org/10.1002/2016GL070079}

Kumazawa, T., \& Ogata, Y. (2013). Quantitative description of induced
seismic activity before and after the 2011 Tohoku-Oki earthquake by
nonstationary ETAS models. \emph{Journal of Geophysical Research: Solid
Earth} , \emph{118} (12),
6165--6182.\emph{https://doi.org/10.1002/2013JB010259}

Kumazawa, T., Ogata, Y., \& Tsuruoka, H. (2019). Characteristics of
seismic activity before and after the 2018 M6. 7 Hokkaido Eastern Iburi
earthquake. \emph{Earth, Planets and Space} , \emph{71} (1), 1--17.

Lay, T., \& Kanamori, H. (1981). An asperity model of large earthquake
sequences. In. \emph{Maurice Ewing Series} . Maurice Ewing Series.
American Geophysical Union,
(579--592).\emph{https://doi.org/10.1029/ME004p0579}

Llenos, A. L., \& Michael, A. J. (2013). Modeling earthquake rate
changes in Oklahoma and Arkansas: Possible Signatures of induced
seismicity. \emph{Bulletin of the Seismological Society of America}
,\emph{103} (5), 2850--2861.\emph{https://doi.org/10.1785/0120130017}

Lohman, R. B., \& McGuire, J. J. (2007). Earthquake swarms driven by
aseismic creep in the Salton Trough, California. \emph{Journal of
Geophysical Research: Solid Earth} ,
112(B4).\emph{https://doi.org/10.1029/2006JB004596}

Madariaga, B. Y. R. (1976). Dynamics of an expanding circular
fault.\emph{Bulletin of the Seismological Society of America} ,
\emph{66} ,
639--666.\emph{https://doi.org/10.1111/j.1461-0248.2009.01352.x}

Mendoza, C., \& Hartzell, S. H. (1988). Aftershock patterns and main
shock faulting. \emph{Bulletin of the Seismological Society of America}
,\emph{78} (4), 1438--1449.

Mogi, K. (1969). Some features of recent seismic activity in and near
Japan (2): Activity before and after Great Earthquakes. \emph{Bulletin
of the Earthquake Research Institute} , \emph{47} , 395--417.

Mogi, K. (1989). The mechanism of the occurrence of the Matsushiro
earthquake swarm in central Japan and its relation to the 1964 Niigata
earthquake. Tectonophysics, 159(1--2), 109--119.

Nanjo, K. Z., Miyaoka, K., Tamaribuchi, K., Kobayashi, A., \& Yoshida,
A. (2018). Related spatio-temporal changes in hypocenters and the b
value in the 2017 Kagoshima Bay swarm activity indicating a rise of hot
fluids. \emph{Tectonophysics} , \emph{749} ,
35--45.\emph{https://doi.org/10.1016/j.tecto.2018.10.023}

Nur, A., \& Booker, J. R. (1972). Aftershocks caused by pore fluid
flow?\emph{Science} , \emph{175} (4024),
885--887.\emph{https://doi.org/10.1126/science.175.4024.885}

Ogata, Y. (1988). Statistical models for earthquake occurrences and
residual analysis for point processes. \emph{Journal of the American
Statistical Association} , \emph{83} (401),
9--27.\emph{https://doi.org/10.1080/01621459.1988.10478560}

Ogata, Y. (1992). Detection of precursory relative quiescence before
great earthquakes through a statistical model. \emph{Journal of
Geophysical Research} , \emph{97} (B13),
19845--19871.\emph{https://doi.org/10.1029/92JB00708}

Ogata, Y. (2006). \emph{Statistical analysis of seismicity: Updated
version} (p. SASeis2006). Institute of Statistical Mathematics.

Okada, T., Matsuzawa, T., Umino, N., Yoshida, K., Hasegawa, A.,
Takahashi, H. et al. (2016). Hypocenter migration and crustal seismic
velocity distribution observed for the inland earthquake swarms induced
by the 2011 Tohoku-Oki earthquake in NE Japan: Implications for crustal
fluid distribution and crustal permeability. In \emph{Crustal
Permeability} ,
(307--323).\emph{https://doi.org/10.1002/9781119166573.ch24}

Parotidis, M., Rothert, E., \& Shapiro, S. A. (2003). Pore-pressure
diffusion: A possible triggering mechanism for the earthquake swarms
2000 in Vogtland/NW-Bohemia, central Europe. Geophysical Research
Letters, 30(20), 10--13. https://doi.org/10.1029/2003GL018110

Prieto, G. A., Parker, R. L., \& Vernon III, F. L. (2009). A Fortran 90
library for multitaper spectrum analysis. \emph{Computers and
Geosciences} , \emph{35} (8),
1701--1710.\emph{https://doi.org/10.1016/j.cageo.2008.06.007}

Rice, J. R. (1992). Chapter 20 Fault stress states, pore pressure
distributions, and the weakness of the San Andreas fault.
In.\emph{International Geophysics} . Elsevier, \emph{51}
.\emph{https://doi.org/10.1016/S0074-6142(08)62835-1}

Roland, E., \& McGuire, J. J. (2009). Earthquake swarms on transform
faults. \emph{Geophysical Journal International} , \emph{178} (3),
1677--1690.\emph{https://doi.org/10.1111/j.1365-246X.2009.04214.x}

Ross, Z. E., Kanamori, H., \& Hauksson, E. (2017b). Anomalously large
complete stress drop during the 2016 Mw5.2 Borrego Springs earthquake
inferred by waveform modeling and near-source aftershock
deficit.\emph{Geophysical Research Letters} , \emph{44} (12),
5994--6001.\emph{https://doi.org/10.1002/2017GL073338}

Ross, Z. E., Rollins, C., Cochran, E. S., Hauksson, E., Avouac, J. P.,
\& Ben-Zion, Y. (2017a). Aftershocks driven by afterslip and fluid
pressure sweeping through a fault-fracture mesh. \emph{Geophysical
Research Letters} , \emph{44} (16),
8260--8267.\emph{https://doi.org/10.1002/2017GL074634}

Ross, Z. E., Kanamori, H., Hauksson, E., \& Aso, N. (2018). Dissipative
intraplate faulting during the 2016 Mw6.2 Tottori, Japan
earthquake.\emph{Journal of Geophysical Research: Solid Earth} ,
\emph{123} (2), 1631--1642.\emph{https://doi.org/10.1002/2017JB015077}

Rossi, G., Zuliani, D., \& Fabris, P. (2016). Long-term GNSS
measurements from the northern Adria microplate reveal fault-induced
fluid mobilization. \emph{Tectonophysics} , \emph{690} ,
142--159.\emph{https://doi.org/10.1016/j.tecto.2016.04.031}

Rossi, G., Fabris, P., \& Zuliani, D. (2018). Overpressure and Fluid
Diffusion Causing Non-hydrological Transient GNSS
Displacements.\emph{Pure and Applied Geophysics} , \emph{175} (5),
1869--1888. https://doi.org/10.1007/s00024-017-1712-x

Ruhl, C. J., Abercrombie, R. E., Smith, K. D., \& Zaliapin, I. (2016).
Complex spatiotemporal evolution of the 2008 Mw 4.9 Mogul earthquake
swarm (Reno, Nevada): Interplay of fluid and faulting. \emph{Journal of
Geophysical Research: Solid Earth} , \emph{121} (11), 8196--8216.

Saiga, A., Matsumoto, S., Uehira, K., Matsushima, T., \& Shimizu, H.
(2010). Velocity structure in the crust beneath the Kyushu
area.\emph{Earth, Planets and Space} , \emph{62} (5),
449--462.\emph{https://doi.org/10.5047/eps.2010.02.003}

Sato, T., \& Hirasawa, T. (1973). Body wave spectra from propagating
shear cracks. \emph{Journal of Physics of the Earth} , \emph{21} (4),
415--431.\emph{https://doi.org/10.4294/jpe1952.21.415}

Scholz, C. H. (1998). Earthquakes and friction laws. \emph{Nature}
,\emph{391} (6662), 37--42.\emph{https://doi.org/10.1038/34097}

Segall, P. (1989). Earthquakes triggered by fluid
extraction.\emph{Geology} , \emph{17} (10), 942--946.
https://doi.org/10.1130/0091-7613(1989)017\textless{}0942:ETBFE\textgreater{}2.3.CO;2

Shapiro, S. A., Huenges, E., \& Borm, G. (1997). Estimating the crust
permeability from fluid-injection-induced seismic emission at the KTB
site. \emph{Geophysical Journal International} , \emph{131} (2),
F15--F18.\emph{https://doi.org/10.1111/j.1365-246X.1997.tb01215.x}

Shelly, D. R., Taira, T., Prejean, S. G., Hill, D. P., \& Dreger, D. S.
(2015). Fluid-faulting interactions: Fracture-mesh and fault-valve
behavior in the February 2014 Mammoth Mountain, California, earthquake
swarm. Geophysical Research Letters, 42(14), 5803--5812.
https://doi.org/10.1002/2015GL064325

Shelly, D. R., Ellsworth, W. L., \& Hill, D. P. (2016). Fluid-faulting
evolution in high definition: Connecting fault structure and
frequency-magnitude variations during the 2014 Long Valley Caldera,
California, earthquake swarm. \emph{Journal of Geophysical Research:
Solid Earth} , \emph{121} (3),
1776--1795.\emph{https://doi.org/10.1002/2015JB012719}

Sibson, R. H. (1992). Implications of fault-valve behaviour for rupture
nucleation and recurrence. \emph{Tectonophysics} , \emph{211} (1--4),
283--293.\emph{https://doi.org/10.1016/0040-1951(92)90065-E}

Sibson, R. H. (2020). Preparation zones for large crustal earthquakes
consequent on fault-valve action. \emph{Earth, Planets and Space}
,\emph{72} (1), 1--20.

Talwani, P., \& Acree, S. (1985). Pore pressure diffusion and the
mechanism of reservoir-induced seismicity. In. \emph{Pure and Applied
Geophysics PAGEOPH} . Berlin: Springer, 122(6),
(947--965).\emph{https://doi.org/10.1007/BF00876395}

Talwani, P., Chen, L., \& Gahalaut, K. (2007). Seismogenic permeability,
ks. \emph{Journal of Geophysical Research} , \emph{112} (B7),
1--18.\emph{https://doi.org/10.1029/2006JB004665}

Yoshida, K., Saito, T., Emoto, K., Urata, Y., \& Sato, D. (2019b).
Rupture directivity, stress drop, and hypocenter migration of small- and
moderate-sized earthquakes in the Yamagata--Fukushima border swarm
triggered by upward pore-pressure migration after the 2011 Tohoku-Oki
earthquake. \emph{Tectonophysics} , \emph{769} . PubMed:\emph{228184}

Terakawa, T., Hashimoto, C., \& Matsu'ura, M. (2013). Changes in seismic
activity following the 2011 Tohoku-oki earthquake: Effects of pore fluid
pressure. \emph{Earth and Planetary Science Letters} , \emph{365} ,
17--24.\emph{https://doi.org/10.1016/j.epsl.2013.01.017}

Ueno, H., Hatakeyama, S., Aketagawa, T., Funasaki, J., \& Hamada, N.
(2002). Improvement of hypocenter determination procedures in the Japan
Meteorological Agency. \emph{Quarterly Journal of Seismology} ,\emph{65}
, 123--134.

Utsu, T. (1961). A statistical study on the occurrence of
aftershocks.\emph{Geophysical Magazine} , \emph{30} , 521--605.

Utsu, T., Ogata, Y., S, R., \& Matsu'ura. (1995). The centenary of the
Omori formula for a decay law of aftershock activity. \emph{Journal of
Physics of the Earth} , \emph{43} (1),
1--33.\emph{https://doi.org/10.4294/jpe1952.43.1}

Vidale, J. E., \& Shearer, P. M. (2006). A survey of 71 earthquake
bursts across southern California: Exploring the role of pore fluid
pressure fluctuations and aseismic slip as drivers. \emph{Journal of
Geophysical Research: Solid Earth} , \emph{111}
(B5),\emph{https://doi.org/10.1029/2005JB004034}

Waite, G. P., \& Smith, R. B. (2002). Seismic evidence for fluid
migration accompanying subsidence of the Yellowstone
caldera.\emph{Journal of Geophysical Research: Solid Earth} , \emph{107}
(B9), ESE 1--ESE 1.\emph{https://doi.org/10.1029/2001JB000586}

Wada, I., \& Behn, M. D. (2015). Focusing of upward fluid migration
beneath volcanic arcs: Effect of mineral grain size variation in the
mantle wedge. \emph{Geochemistry, Geophysics, Geosystems} ,\emph{16}
(11), 3905--3923. https://doi.org/10.1002/2015GC005950

Waldhauser, F. (2002). Fault structure and mechanics of the Hayward
Fault, California, from double-difference earthquake
locations.\emph{Journal of Geophysical Research} , \emph{107} (B3),
2054.\emph{https://doi.org/10.1029/2000JB000084}

Waldhauser, F., \& Ellsworth, W. L. (2000). A double-difference
earthquake location algorithm: Method and application to the Northern
Hayward fault, California. \emph{Bulletin of the Seismological Society
of America} , \emph{90} (6),
1353--1368.\emph{https://doi.org/10.1785/0120000006}

Wessel, P., \& Smith, W. H. F. (1998). New, improved version of generic
mapping tools released. \emph{Eos, Transactions American Geophysical
Union} , \emph{79} (47),
579--579.\emph{https://doi.org/10.1029/98EO00426}

Wetzler, N., Lay, T., Brodsky, E. E., \& Kanamori, H. (2018). Systematic
deficiency of aftershocks in areas of high coseismic slip for large
subduction zone earthquakes. \emph{Science Advances} , \emph{4} (2),
eaao3225.\emph{https://doi.org/10.1126/sciadv.aao3225}

Woessner, J., Schorlemmer, D., Wiemer, S., \& Mai, P. M. (2006). Spatial
correlation of aftershock locations and on-fault main shock
properties.\emph{Journal of Geophysical Research} , \emph{111} (B8),
(B8).\emph{https://doi.org/10.1029/2005JB003961}

Yabe, S., \& Ide, S. (2018). Variations in precursory slip behavior
resulting from frictional heterogeneity. \emph{Progress in Earth and
Planetary Science} , \emph{5} (1),
43.\emph{https://doi.org/10.1186/s40645-018-0201-x}

Yoshida, K., \& Hasegawa, A. (2018a). Sendai-Okura earthquake swarm
induced by the 2011 Tohoku-Oki earthquake in the stress shadow of NE
Japan: Detailed fault structure and hypocenter
migration.\emph{Tectonophysics} , \emph{733} ,
132--147.\emph{https://doi.org/10.1016/j.tecto.2017.12.031}

Yoshida, K., \& Hasegawa, A. (2018b). Hypocenter migration and
seismicity pattern change in the Yamagata--Fukushima Border, NE Japan,
caused by fluid movement and pore pressure variation. \emph{Journal of
Geophysical Research: Solid Earth} , \emph{123} (6),
5000--5017.\emph{https://doi.org/10.1029/2018JB015468}

Yoshida, K., Hasegawa, A., \& Yoshida, T. (2016a). Temporal variation of
frictional strength in an earthquake swarm in NE Japan caused by fluid
migration. \emph{Journal of Geophysical Research: Solid Earth}
,\emph{121} (8), 5953--5965.\emph{https://doi.org/10.1002/2016JB013022}

Yoshida, K., Hasegawa, A., \& Okada, T. (2016b). Heterogeneous stress
field in the source area of the 2003 M6.4 Northern Miyagi Prefecture, NE
Japan, earthquake. \emph{Geophysical Journal International} ,\emph{206}
(1), 408--419.\emph{https://doi.org/10.1093/gji/ggw160}

Yoshida, K., Hasegawa, A., Yoshida, T., \& Matsuzawa, T. (2019a).
Heterogeneities in stress and strength in tohoku and its relationship
with earthquake sequences triggered by the 2011 M9 Tohoku-Oki
earthquake. \emph{Pure and Applied Geophysics} , \emph{176} (3),
1335--1355.\emph{https://doi.org/10.1007/s00024-018-2073-9}

Yoshida, K., Uchida, N., Hiarahara, S., Nakayama, T., Matsuzawa, T.,
Okada, T., et al. (2020a). 2019 M6. 7 Yamagata-Oki earthquake in the
stress shadow of 2011 Tohoku-Oki earthquake: Was it caused by the
reduction in fault strength? Tectonophysics, 793, 228609.

Yoshida, K., Taira, T., Matsumoto, Y., Saito, T., Emoto, K., \&
Matsuzawa, T. (2020b). Stress release process along an intraplate fault
analogous to the plate boundary: A case study of the 2017 M5. 2
Akita-Daisen earthquake, NE Japan. \emph{Journal of Geophysical Research
-- Solid Earth} , e2020JB019527.

Yukutake, Y., Ito, H., Honda, R., Harada, M., Tanada, T., \& Yoshida, A.
(2011). Fluid-induced swarm earthquake sequence revealed by precisely
determined hypocenters and focal mechanisms in the 2009 activity at
Hakone volcano, Japan\emph{. Journal of Geophysical Research}
,\emph{116} (B4).\emph{https://doi.org/10.1029/2010JB008036}

Zhao, D., Yanada, T., Hasegawa, A., Umino, N., \& Wei, W. (2012).
Imaging the subducting slabs and mantle upwelling under the Japan
Islands. \emph{Geophysical Journal International} , \emph{190} (2),
816--828.

\selectlanguage{english}
\FloatBarrier
\end{document}

