Geophys. J. Int. (2020) 221, 289–304 doi: 10.

1093/gji/ggaa008
Advance Access publication 2020 January 09
GJI Geodynamics and tectonics

Inelastic strain rate and stress fields in and around an aseismic zone
of Kyushu Island, Japan, inferred from seismic and GNSS data

Y. Yuasa,1 S. Matsumoto,2 S. Nakao,3 T. Matsushima2 and T. Ohkura4

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
1 Graduate School of Science, Kyushu University, Fukuoka 819-0395, Japan. E-mail; yuasa@sevo.kyushu-u.ac.jp
2 Institute
of Seismology and Volcanology, Faculty of Science, Kyushu University, Shimabara 855-0843, Japan
3 Graduate School of Science and Engineering, Kagoshima 890-0065, Japan
4 Aso Volcanological Laboratory, Kyoto University, Kumamoto 869–1404, Japan

Accepted 2020 January 1. Received 2019 December 30; in original form 2019 January 10

SUMMARY
Understanding earthquake processes and crustal deformation requires knowledge of the stress
concentration process in the crust. With the enhancement of observation networks, it has
become possible to consider in detail the relationships between localized deformation and
seismic activity in island arcs and the process of stress concentration. According to previous
studies, inelastic deformation in localized weak zones in the crust is considered to play an
important role in the stress concentration process. Kyushu, located in southwest Japan, has
a 20–30 km band-like active seismic activity and an enclosed aseismic zone. In particular, a
part of the seismic active region called the Beppu-Simahara Graben, which is dominated by
north–south extensional deformation, is characterized by high seismic activity and a remark-
able aseismic zone. We identified the relationship between inelastic deformation and stress
concentration processes in this area by using analyses of geodetic and seismic data. The results
inverted from both the strain rate field obtained by the geodetic observations and the devia-
toric stress field estimated from focal mechanism data reveal a large inelastic deformation zone
(∼ 10−7 yr−1 ) beneath the area of active seismicity. From comparison with previous works, the
inelastic deformation zone in the lower crust may correspond to an area with high temperature
and/or fluid. This may suggest that inelastic deformation is in progress in the area where the
strength of lower crustal rocks has reduced due to the presence of geothermics and/or fluids.
Furthermore, we confirmed that this inelastic deformation causes stress concentrations of up
to 10 kPa yr−1 in the upper crust. These results show that stress concentration occurs locally
in the upper crust, above the inelastic deformation zone in the weakened lower crust, owing to
the presence of geothermal and/or fluid; this stress concentration induces seismic activity and
crustal deformation.
Key words: Elasticity and anelasticity; Japan; Seismicity and tectonics; Dynamics: seismo-
tectonics; Intra-plate processes.

of the lower crust caused a stress concentration on the intraplate

1 I N T RO D U C T I O N
earthquake fault (Stuart et al. 1997; Kenner & Segall 2000; Fang
Inelastic deformation, including earthquakes, releases stress accu- & Hou 2019). Moreover, for earthquakes occurring in an island
mulated in the crust, which causes a stress change in the surround- arc, similar models have been suggested. From analysis using finite
ing medium. Therefore, it is necessary to understand the process element modelling, it was suggested that the viscous heterogeneity
of stress concentration and relaxation to understand the processes structure of the crust (Hyodo & Hirahara 2003) and/or the mantle
of earthquake occurrence and crustal deformation. In particular, (Yamasaki & Seno 2005) causes stress concentrations in the upper
for intraplate earthquakes, the stress concentration on faults due to crust, 200–300 km away from plate boundaries. Iio et al. (2002)
crustal and/or mantle heterogeneity and inelastic deformation has inferred that the lower crust and the uppermost mantle are weak-
been discussed. ened by water from a slab. Geodetic observations using the Global
For example, for earthquakes occurring in the continental crust, Navigation Satellite System (GNSS) have revealed a strain con-
it was demonstrated that aseismic slip on the detachment fault be- centration zone in Japan (Sagiya et al. 2000; Miura et al. 2004).
neath the earthquake fault or ductile deformation in the weak zone Meneses-Gutierrez et al. (2018) showed that strain concentration


C The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. 289
290 Y. Yuasa et al.

owed to inelastic deformation in the crust and a low elastic constant We estimate the daily positions of each station using the Bernese
layer in the upper crust. Hasegawa et al. (2005) further argued that GNSS Software Version 5.2 (Dach et al. 2015). We use precise
concentrated deformation is caused by inelastic deformation in the ephemerides and rotation parameters from the Center for Orbit De-
crust, for which the strength had been reduced by a fluid. termination in Europe. The coordinates of the GNSS sites are esti-
Such localized seismic activity and deformation, implying stress mated with respect to ITRF2014 (Altamimi et al. 2016). Moreover,
concentration, are also observed in Kyushu Island, southwest Japan. to improve the accuracies of the station coordinates, we estimate
In the Beppu-Shimabara Graben (BSG), high seismic activity the tropospheric delay every 1 hr and their horizontal gradients
(Fig. 1) and NS extension deformation are observed over an area every 24 hr. The observed horizontal velocity vectors relative to
with a width of 20–30 km (Tada 1984, 1985). In this region, Savage a reference point are shown in Fig. 3. We adopted GEONET site

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
et al. (2016) estimated a low stress ratio (tension-dominated stress), 950456 (Kamitsushima, 34.656◦ N, 129.482◦ E) as the reference sta-
a low compressional strain rate and strong seismic anisotropy using tion. The error ellipse is very small because the vector uncertainties
seismic and GNSS data. Moreover, the 2016 Kumamoto earthquake are less than 0.25 mm yr−1 at most. The vectors at the observation
(Mw = 7.3; Asano & Iwata 2016; Shimizu et al. 2016) occurred points located north of 32.5◦ N are essentially westward, whereas
in this area. Along the Yatsushiro sea (YS) and Amakusanada sea those south of this parallel exhibit counter-clockwise rotation with
(AS), many earthquakes are distributed in a band of 20–30 km southing.
width. However, there is a remarkable deficit in the seismic activ- We calculate the total strain rate at each grid from the displace-
ity around the Amakusa Islands (AI; Fig. 1). In addition, Wallace ment rate vectors, using a method developed by Shen et al. (1996)
et al. (2009) showed a left lateral shear zone in southern Kyushu and Sagiya et al. (2000). In this method, the relationship between
(Fig. 1). Previous studies estimated that the tectonic stress field an observed displacement rate (U, V) at a GNSS site X = (X, Y ),
of most of Japan is typified by reverse-fault and strike-slip types horizontal displacement rate (u, v), total strain rate (ε̇x x , ε̇x y , ε̇ yy )
(Townend & Zoback 2006; Terakawa & Matsu’ura 2010), while and rotation rate ω at a given point x i = (x i , yi ) is written as
the BSG has a normal-fault stress field, which is one of the char- follows:
acteristics of this region. Furthermore, Matsumoto et al. (2015) ⎛ ⎞
estimated that central Kyushu has a tension-dominated stress field. u
⎜v ⎟
In this study, we consider the possibility that localized seismic     ⎜ ⎟  i
1 0 xi yi 0 yi ⎜ε̇x x ⎟
activity and crustal deformation is attributed to inelastic deforma- U
= ⎜ ⎟ + eix , (1)
tion in the lower crust. In addition, we consider the stress con- V 0 1 0 xi yi −xi ⎜ ⎟
⎜ε̇x y ⎟ ey
⎝ε̇ yy ⎠
ditions in and around the aseismic zone using data on inelastic
deformation. ω

where xi and yi are X − x i and Y − yi , respectively. (eix , eiy ) are
observation errors, which are weighted depending on the distance
2 D ATA A N D M E T H O D O L O G Y
between the GNSS site X and the given point x i . It is written as
We consider the seismic activity and crustal deformation in this re- follows:
gion using a mechanical model as shown in Fig. 2. We assume that  
Ri2
the upper crust and lower crust are elastic and elastic-plastic bodies, eix, y = σx,y
i
exp , (2)
respectively. In this model, constant stress loading from outside the 2D 2
region causes spatial long wavelength deformation, while localized
where σx,yi
is the observation error of the displacement rate com-
inelastic deformation in the lower crust causes short wavelength
ponent, Ri = |X − x i | and D is the parameter which controls the
deformation. The sum of both is obtained from geodetic observa-
spatial variation of strain rate. D is termed the distance decay con-
tions at the ground surface. We therefore attempt to extract short
stant (DDC). We can estimate the strain rate by solving eq. (1). In
wavelength components from observed crustal deformation data to
this study, the spatial grid interval is 0.02◦ and we adopt a DDC of
estimate inelastic deformation in the lower crust using the method
25 km, which is consistent with the GNSS station distribution. We
proposed by Meneses-Gutierrez & Sagiya (2016). We then estimate
use the stations within 2D from the point of strain rate calculation.
the inelastic deformation in the lower crust using the obtained short
If DDC is small, strain rates reflecting local effects can be obtained,
wavelength component. In this study, we calculate the total strain
but the calculable points are reduced.
rate at the ground surface from observed displacement velocity
The total strain rate fields are shown in Fig. 4(a). In most ar-
data provided by a GNSS network. In the following subsections, we
eas, the EW and NS components of the total strain rate dominate
present the details of this method and the data used to validate this
the compression and extension, respectively. In addition, the total
approach.
strain rate field is perturbed within the order of 10−7 yr−1 . For exam-
ple, the total strain rate in the NS component around the Shimabara
Peninsula, located north of Amakusa, exhibits compression and that
2.1 GNSS data and wavelength decomposition
in the EW component along the YS exhibits extension, which are
opposite to the general trend. These spatial characteristics can be
2.1.1 Calculation of the strain rate field from GNSS data
attributed to the following factors: (1) the elastic strain rate is caused
Time-series data acquired by the GNSS Earth Observation Net- by mechanical interaction from outside of the model region, (2) the
work system (GEONET) of the Geospatial Information Authority heterogeneous properties of the crust and (3) inelastic deformation
of Japan and temporal stations installed by Kyushu, Kyoto, and in the lower crust. The first factor creates a nearly constant strain
Kagoshima Universities are used to calculate the displacement rate rate field over the target region, while the latter two create pertur-
from 2004 January 1 to 2014 December 31. Note that we use data bations of the strain rate field. The heterogeneity of the observed
before the 2016 Kumamoto earthquake and not those data contain- strain rate field is extremely large compared with the expected one
ing large stress and strain field changes induced by the earthquake. from the elastic properties (Saiga et al. 2010). The seismic velocity
 Inelastic strain rate and stress fields of Kyushu Island, Japan 291

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure 1. Tectonic setting of the study area and epicentre distribution. (a) The black solid lines in the map indicate the active fault traces from the active
fault database of the National Institute of Advanced Industrial Science and Technology. The red rectangles are source faults of the Kumamoto earthquake
(Asano & Iwata 2016). They are called Futagawa fault (FF) and Hinagu fault (HF). The red triangles are active volcanoes. TB—Tachibana Bay; KM—Kyushu
Mountains; AI—Amakusa Islands; AS—Amakusanada Sea: YS—Yatsushiro Sea. (b) The events plotted in the map using dots occurred from 1995 to 2014
and had magnitudes exceeding 1.0 and depths of less than 20 km. The blue rectangle is the study area. The red ellipse shows the Beppu-Shimabara Graben
(BSG). The pink ellipse shows the left lateral shear zone pointed out by Wallace et al. (2009). The star shows the location of the Kumamoto earthquake (2016,
Mw = 7.3). The Philippine Sea slab (PHS) subducts with an NW direction.

Figure 2. Schematic illustration of crustal deformation in the target region.

perturbation was at most ∼5 per cent, which is insufficient to gener- 2.1.2 Extraction of the short-wavelength component
ate a drastic change in the strain field (i.e. assuming that the density
We attempt to extract the short-wavelength component of the ob-
is constant, the modulus of rigidity changes is expected to be ∼10
served strain rate field based on a method developed by Meneses-
per cent). Thus, it is plausible that local inelastic deformation of the
Gutierrez & Sagiya (2016) and Menese-Gutierrez et al. (2018). For
lower crust generates this spatial perturbation of the strain rate field.
the case considered here, we apply a moving average filter to the
Therefore, in this study, we assume that the components of the strain
observed total strain rate field in the study area to eliminate the
rate with long and short wavelengths can be attributed to elastic de-
long-wavelength component. The removed long-wavelength com-
formation due to external (i.e. Philippine Sea slab, PHS, subduction)
ponents are shown in Fig. 4(b). We use filtered data with a smoothing
and inelastic deformation in the lower crust, respectively. Meneses-
radius of 60 km. This radius was selected to extract the optimum
Gutierrez & Sagiya (2016) validate this assumption and the idea
characteristic of the strain rate field in the short-scale length range
of wavelength decomposition through a comparison between inter-
(see Fig. A1 for the radius dependency). Contraction and exten-
seismic and post-seismic deformation. However, we cannot validate
sion of 40–50 km in width are obtained with the short-wavelength
those parameters easily in this study because a large earthquake
component (Fig. 4c). We found Contraction around the Shimabara
that could be used for verification did not occur around the target
Peninsula, and may correspond to a magma chamber after eruption
region.
292 Y. Yuasa et al.
p
Therefore, we can estimate ε̇ jk by solving the observation eq. (5)
using a least-squares method, taking the short-wavelength compo-
nent of the strain rate at the surface as the input data.
In our analysis of Kyushu Island, we arrange the boxes as shown in
Fig. 5. The size of the box is 0.28◦ in the latitudinal direction, 0.30◦
in the longitudinal direction and 10 km in the vertical direction. The
depth of the box is 15–25 km corresponding to the lower crust. In
addition, referring to previous studies (Abe et al. 2010; Kohno et al.
2008), we arranged boxes at the locations of the Unzen and Aso vol-

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
cano magma storage regions, respectively. Where a volcano’s box
and the normal box of the lower crust overlap, the volcano box is
prioritized. Thus, we arrange 30 boxes in the target region. Gener-
ally, the inelastic strain rate tensor has six independent parameters
under assumption of an isotropic medium. Therefore, there are six
unknowns in each box, and the total number of unknowns is too
large (6 × 30 boxes = 180) for the least-squares solver to stabilize.
Reducing the number of parameters is necessary for solving this
problem. To reduce the number of parameters, there is a method
of solving by giving physically appropriate constraints to the solu-
tion. For example, in the estimation of slip distribution on the fault
Figure 3. Displacement rate vector distribution from 2004 to 2014.
The reference point is GEONET site 950 456 (Kamitsushima 34.656◦ N, plane, a constraint, such as smoothness of the slip distribution, is
129.482◦ E). The blue and red vectors are at the GEONET and temporal used. However, we cannot use this constraint because the inelastic
station, respectively. deformation zone may be localized in the lower crust. To reduce the
number of parameters, we constrain the inelastic strain rate tensor
of the Unzen Volcano, located in the centre of the Shimabara Penin- using a deviatoric stress field because, according to a flow rule,
sula (Kohno et al. 2008). This is the major chamber located at a plastic strain rate is proportional to deviatoric stress (Jaeger et al.
depth of 15 km beneath Tachibana Bay to the west of the Shimabara 2007). In particular, we apply the Prandtl–Reuss formula, which is
Peninsula. the constitutive equation for isotropic media:
p
εi j = si j dλ, (6)
where dλ is a positive scalar. We cannot evaluate dλ accurately
2.2 Inversion method and model concept
because it relates to the material’s plastic properties. However, the
Barbot et al. (2017) developed a representation of displacement relative relationship between each component of the tensor does
and its spatial derivatives in an elastic half space due to uniformly not depend on this value. Therefore, we can know the relative mag-
distributed inelastic deformation over a cuboid. We cannot find a nitude of each component of the inelastic strain rate tensor from
localized distinct weak structure as fault in the lower crust in the the deviatoric stress tensor. In other words, since the direction of
target area; therefore, we do not use the Okada (1992) model but use the inelastic strain rate can be constrained by the deviatoric stress
the model developed by Barbot et al. (2017) instead. As reported by tensor, the unknown parameter is only its magnitude. Thus, we can
Barbot et al. (2017), the total strain change εi at the ground surface transform eq. (5) into the following form, when we can constrain
caused by distributed inelastic deformation can be calculated from the form of the tensor of the box:
the Green’s function G i jk for the cuboid source developed by Barbot ⎛ ⎞ ⎛ ⎞
ε̇1 aG 111 + bG 112 + cG 113 + dG 122 + eG 123 + f G 133
et al. (2017), such that ⎝ε̇2 ⎠ = ⎝aG 211 + bG 212 + cG 213 + dG 222 + eG 223 + f G 233 ⎠
p
εi = G i jk ε jk ( i = 1, 2, 3 ) , (3) ε̇3 aG 311 + bG 312 + cG 313 + dG 322 + eG 323 + f G 333
p × ε̇0 , (7)
where ε jk is the inelastic strain tensor in the cuboid. Note that strain
changes are tensors, which are represented as vectors to simplify where coefficients a–f are the ratios of each component of the tensor
the notation of eq. (3). The component of ε is (εx x , εx y, ε yy ). to the size of the tensor value ε̇0 = 12
p p
ε̇i j ε̇i j , respectively.
p
Furthermore, in this case, if an inelastic strain change δεkl occurs
In contrast, since we cannot estimate the deviatoric stress in the
in the crust for a time interval T , eq. (3) can be converted as
aseismic zone, owing to the absence of earthquakes, we cannot use
follows:
this constraint. However, using the constraint of no volume change
p p p
p
ε̇i = G i jk ε̇ jk , (4) ( ε̇33 = −ε̇11 − ε̇22 ), eq. (5) can be modified as follows, and the
p
unknown parameters can be reduced from 6 to 5:
where ε̇i and ε̇kl
are the total strain rate εi /T and the inelastic ⎛ p⎞
⎞ ε̇11
p
strain rate δεkl /T , respectively. Eq. (4) written in matrix form is ⎛ ⎞ ⎛
G 11 − G 16 G 12 G 13 G 14 − G 16 G 15 ⎜
p ⎟
as follows: ε̇1 ⎜ε̇12 ⎟
⎛ p⎞ ⎝ε̇2 ⎠ = ⎝G 21 − G 26 G 22 G 23 G 24 − G 26 G 25 ⎠ ⎜ε̇13 p ⎟
ε̇11 ⎜ p ⎟. (8)
ε̇3 G 31 − G 36 G 31 G 33 G 34 − G 36 G 35 ⎝ε̇22 ⎠
⎛ ⎞ ⎛ ⎞ ⎜ε̇12 p ⎟
G 11 G 12 G 13 G 14 G 15 G 16 ⎜ p ⎟
⎜ p
ε̇1 ⎟ ε̇23
⎝ε̇2 ⎠ = ⎝G 21 G 22 G 23 G 24 G 25 G 26 ⎠ ⎜ε̇13 ⎟
⎜ε̇ p ⎟ . (5)
Finally, from eqs (5), (7) and (8), the total strain rate on the ground
ε̇3 G 31 G 31 G 33 G 34 G 35 G 36 ⎜ 22 ⎟
⎝ε̇ p ⎠
23
surface caused by the inelastic deformation in the aseismic zone,
p
ε̇33 the magma storage region and the constrained box can be written
 Inelastic strain rate and stress fields of Kyushu Island, Japan 293

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure 4. Strain rate fields calculated at ground surface for components for the EW, shear and NS directions. Blue and red colours correspond to contraction
and extraction, respectively. The rightmost column shows the direction of the principal axis of each strain rate field. (a) Strain rate field calculated using the
method developed by Shen et al. (1996) and Sagiya et al. (2000). (b) Long-wavelength and (c) short-wavelength components of the calculated strain rate.

in matrix form as follows: where d is the short-wavelength component of the total strain rate
⎡ sg ⎤ at each grid.
 sg vol con  m
d = G G G ⎣ mvol ⎦ (9)
mcon 2.3 Estimation of deviatoric stress field and use for
with constraints
⎛p ⎞ To perform the constraint with deviatoric stress, we estimate the
⎛ ⎞ ε̇11
G 15 ⎜
p ⎟ deviatoric stress field in the target region. Fig. 6(a) is the dataset used
G 11 − G 16 G 12 G 13 G 14 − G 16 ⎜ε̇12 ⎟
G sg msg = ⎝G 21 − G 26 G 22 G 23 G 24 − G 26 ⎠ ⎜
G 25 ⎜ε̇13 ⎟ for estimating the deviatoric stress field; this dataset was obtained by
p
⎟
G 35 ⎝ε̇22 ⎠ Matsumoto et al. (2015) and represents earthquakes that occurred
p
G 31 − G 36 G 31 G 33 G 34 − G 36
p
ε̇23 at depths of 5–15 km with magnitudes of 1.0 or more from 2004
to 2014. Most of the earthquakes in this region are shallower than
⎛ p ⎞
15 km, and, thus, we cannot directly estimate the deviatoric stress
ε̇11 of the lower crust. However, we can estimate the deviatoric stress
⎛ p ⎟
⎞ ⎜ε̇12
G 11 G 12 G 13 G 14 G 15 G 16 ⎜ ⎜ε̇ p ⎟
⎟ in the lower crust using the continuity of stress. Fig. 6(b) shows the
G vol mvol = ⎝G 21 G 22 G 23 G 24 G 25 G 26 ⎠ ⎜ 13 ⎟
⎜ε̇ ⎟
p
estimated deviatoric stress field in the upper crust using the method
G 31 G 31 G 33 G 34 G 35 G 36 ⎜ 22 ⎟
⎝ε̇ p ⎠
of Matsumoto (2016). We performed this estimation only for boxes
23
p
with eight or more events. Owing to this limitation, we cannot
ε̇33 estimate the deviatoric stress of the southwestern box. Therefore,
we excluded this box from the estimation of inelastic strain rates.
⎛ ⎞
aG 111 + bG 112 + cG 113 + dG 122 + eG 123 + f G 133 The deviatoric stress field in this area was found to be a strike-
G con
m con ⎝
= aG 211 + bG 212 + cG 213 + dG 222 + eG 223 + f G 233 ⎠ slip or normal slip with a tensile axis in the northwest–southeast
aG 311 + bG 312 + cG 313 + dG 322 + eG 323 + f G 333 or south-south direction. Moreover, the compression axis was not
× ε̇0 , found necessarily parallel to the orientation of the PHS subduction.
294 Y. Yuasa et al.

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure 5. Horizontal position and depth of the boxes. The yellow rectangle shows the box beneath the aseismic zone around the Amakusa Islands. The red
boxes correspond to the magma storage regions of Unzen and Aso, respectively. Boxes with grey shades are excluded because they are covered with red boxes.
The blue inverted triangle indicates the location of the Global Navigation Satellite System (GNSS) sites.

This is because the subduction faults, such as the Nankai Trough, standard error as a weight. Thus, the observation equation to be
are weak and the stress change due to locking of the subduction zone solved is
is lower than the background absolute stress (Wang 2000). From the
continuity of stress, we then assume that the deviatoric stress in the W 1/2 d = W 1/2 Gm, (11)
lower crust is similar to the estimated deviatoric stress in the upper where W is a matrix containing 1/σi2 in diagonal elements.
crust. From this assumption, we can estimate the inelastic strain rate
tensor in the lower crust using the flow rule.
In Fig. 5, boxes with five independent components of the inelastic 3 I N E L A S T I C D E F O R M AT I O N I N A N D
strain rate are coloured yellow, while red boxes show the Unzen and A RO U N D T H E A S E I S M I C Z O N E
Aso volcano magma storage regions. Therefore, the inelastic strain
rate of the red boxes can have an isotropic component (i.e. corre- We estimate that inelastic deformation occurs from the Kumamoto
sponding to volumetric change) because of the volcanic activity, area to the AS (region A in Fig. 8) and from the AS to Hokusatsu
yielding six independent components. Adding the parameters to the (region B in Fig. 8). The value is approximately 10−7 yr−1 , which is
remaining 26 boxes with the constraint, we obtain 43 total parame- several tens of times larger than that released by only earthquakes in
ters ( = 5 × 1 boxes + 6 × 2 box + 26 boxes). The grid number at the upper crust (Matsumoto et al. 2016). It is estimated by Mazzotti
the surface, which corresponds to the points at which the strain rate et al. (2000) that the residual strain rate, excluding the strain change
was calculated, was 6077; hence, d = 6077 × 3 = 18 231. Note due to plate subduction, is high around region A. They pointed out
that Ṁ0 was positive in this case. Finally, the observation equation that this area is an internal deformation concentration zone. The re-
to be solved is as follows: sults in this study are consistent with their point of view. In contrast,
no inelastic deformation was found around the Kyushu Mountains
d = Gm (10) (KM). In this study, normal fault type inelastic deformation with
a tensile axis in the NS direction is detected beneath the aseismic
with
zone (yellow box in Fig. 5), and the magnitude of the value is the
T
d = ε̇11 ε̇21 ε̇31 · · · ε̇26077 ε̇36077 same as for the surrounding area. This suggests that the area be-
neath the aseismic zone is deformed as well as the surrounding area.
The volume change rates of Unzen and Aso are −3.73 × 106 and
 
G = G sg G vol(Unzen) G vol(Aso) G 1con · · · G 26con −1.06 × 106 m3 yr−1 , respectively (Table 1). This is consistent with
the results of previous studies (Kohno et al. 2008; Mochizuki &
Mitsui 2016).
 T Here, a bootstrap test is performed by randomly resampling the
m = msg mvol(Unzen) mvol(Aso) m1con · · · m26con
dataset of the seismic moment tensor to verify the effect of the ac-
We can calculate the standard error σi of the total strain rate at curacy of the deviatoric stress tensor estimated. The estimated de-
the ground surface (Fig. 7), using the inverse of the square of the viatoric stress field of the 95 per cent confidence interval is shown
 Inelastic strain rate and stress fields of Kyushu Island, Japan 295

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure 6. (a) P- and T-axes of earthquakes used for our analysis. The blue and green lines in the left-hand side figure show the P-axis plunging ≤ 45◦ and
> 45◦ from the horizontal, respectively. The blue and red lines in the right-hand side figure show the T-axis plunging ≤ 45◦ and > 45◦ from the horizontal,
respectively. (b) Estimated deviatoric stress field in the upper crust.

Figure 7. Standard error distribution of estimated strain rate. The error is displayed in the colour scale shown below the maps. There are some places in which
the standard error exceeds 1.0 × 10−7 yr−1 , but it is generally less than 3.0 × 10−8 yr−1 .

in Fig. 9(a). The T- and B-axis confidence intervals appear to be is found to be on the order of 10−7 /yr . Furthermore, no inelastic
extending. We suggest that this is caused by the stress ratio (Fig. 9a) deformation is observed around the KM.
because the minimum compressive stress is superior, and the mag- To evaluate the reliability of model, we examine the influence of
nitudes of the intermediate stress and maximum compressive stress the spatial arrangement of boxes on the results (Fig. 10). We perform
are almost the same. Note that the deviatoric stress used to calculate this test by estimating the inelastic strain rate using a model, in which
the stress ratio is estimated by the method of Matsumoto (2016). the box arrangement is shifted by ± 0.03◦ in the east–west direction
The result of the inversion performed within this confidence inter- and 0.025◦ in the north–south direction (i.e. shifts of ∼10 per cent
val is shown in Fig. 9(b). The upper limit of the inelastic strain rate of the box size). Since the shift of 0.03◦ to the west causes a lack
296 Y. Yuasa et al.

4 DISCUSSION
According to Saiga et al. (2010), it is estimated that there is a low
velocity area around the AI at a depth of 20 km. Furthermore, this
distribution is consistent with relatively high attenuation zones in
the lower crust (Wang et al. 2017). These suggest that deep parts
of the BSG have a relatively high temperature or fluid. Tada (1984,
1985) showed that the BSG is an extension of the Okinawa Trough.
It is possible that high temperature fluid may be supplied from the

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
deep part by crustal thinning. In Hokusatsu, a low resistivity area
is estimated (Hata et al. 2017). There is a gold deposit in this area,
and Hata et al. (2017) showed that the low resistivity reflects a
hydrothermal or hydrothermal alteration area associated with the
deposit. Although the causes of this low velocity and low resistiv-
ity could not be identified from these previous studies alone, these
structures suggest the presence of water and wet lower crustal rocks
in these region. In addition, the higher geothermal gradient (30–40
◦
C km-1 ) in these regions than that in other areas (Tanaka et al. 2004;
Pollitz et al. 2017) suggests that the temperature reaches >600 ◦ C
at 20 km depth. It is supposed that feldspars are the most abundant
mineral constituents of the lower crustal rocks, and an experimental
study (Rybacki & Dresen 2004) showed that wet feldspars, under a
Figure 8. Estimated inelastic strain rate in the lower crust. The beach balls
similar environment in the lower crust to it in this study (>600 ◦ C,
of Unzen and Aso show the deviatoric component of the estimated inelastic strain rate = 10−14 s−1 , shear stress = 10 MPa), could plastically de-
strain rate tensor (see Table 1 for volume change rate). A and B show form. On the contrary, around the KM, where inelastic deformation
inelastic deformation zones corresponding to the lower part of the BSG and was not detected, there is no observation suggesting the presence
high shear zone, respectively. of high temperature fluids. Therefore, the estimated high inelastic
strain rate in these regions may be attributed to the decrease in rock
strength due to the presence of geothermics and water. This is con-
sistent with the lower crust in these areas estimated to be relatively
low viscosity (Moore et al. 2017; Pollitz et al. 2017).
It is important to examine the stress change due to inelastic
deformation because inelastic deformation relaxes the stress in the
of moment tensor data for the west end of the target region, we can
medium and this process loads stress on the surrounding medium.
estimate only one deviatoric stress in the western box of the target
To investigate how inelastic deformation affects the medium in
area. Thus, we cannot not perform the test in this case. As similar
the other parts, we calculate the tensor product dT between the
results were obtained for all models, we conclude that the result does
background stress tensor σi j (Fig. 12a) and the stress change tensor
not significantly depend on this degree of model selection. Thus, we
σi j caused by inelastic deformation in the lower crust. Fig. 12(a)
adopt the model of Fig. 5 as the optimal model based on the GNSS
is similar to Fig. 6(b), but the background stress field (Fig. 12a)
site distribution, the spatial scale of the observed strain rate field and
has been estimated including the shaded grey area, which has been
the epicentre distribution. Note that the parameter settings for both
removed in Fig. 5. The tensor product is given by the following
procedures for removing the long-wavelength component of the
equation:
strain rate and for the DDC affect the inverted results. We compare
the long-wavelength component removed in our analysis with the
3 3
strain rate distribution calculated from the slip deficit distribution i =1 j =1 σi j σi j
dT = . (12)
on the plate boundary (Noda et al. 2018). The similarity of both 3 3 3 3
i =1 j =1 σi2j i =1 j =1 σi2j
rate fields reveals that our procedure correctly removes the elastic
strain rate from the PHS subduction.
In the inversion process, we calculated the strain rate at the ground A negative value for dT indicates stress relaxation, whereas a pos-
surface using the modelled inelastic strain rate. The residual be- itive value indicates stress concentration. Figs 12(b) and (c) show
tween the observed and calculated strain rates is shown in Fig. 11. the stress change and dT at 10 km and 20 km depth, respectively.
The residual at the edge of the target area is large, but the residuals The dT at 20 km depth is negative in almost all regions because
in most areas are less than 3.0 × 10−8 yr−1 . The standard error in stress relaxation due to the inelastic deformation occurs. In con-
the strain rate field estimation performed using the method of Shen trast, the dT at 10 km depth spatially varies. Fig. 12(b) indicates
et al. (1996) is 1.0 − 3.0 × 10−8 yr−1 (Fig. 7). This result is consid- that stress concentration is progressing in region A and region B,
ered to be comparatively well reproduced by our model. In contrast, while stress relaxation is in progress around the KM. As mentioned
areas with large residuals are found around the centre of the target above, we constrain the direction of the inelastic strain rate in the
area (e.g. Kumamoto and Shimabara). These large residuals may be lower crust by the shape of the deviatoric stress tensor in the upper
caused by the elastic heterogeneity. In fact, a low seismic velocity crust located directly above. Therefore, the large inelastic strain rate
anomaly can be seen in these areas (Saiga et al. 2010). In contrast, in the lower crust results in a stress concentration in the crust just
there are possible errors due to the removal of long wavelength pat- above. In this case, the dT value will be positive. However, the dT at
tern too. However, we do not have a technique to evaluate this error a block could be negative if the inelastic strain rate of the surround-
quantitatively in this study. ing lower crust is relatively stronger than one under the block. In
 Inelastic strain rate and stress fields of Kyushu Island, Japan 297

Table 1. Volume change rate of the major magma storage regions of the Unzen and Aso volcanoes.
95% confidence interval Source depth
Volume change rate (m3 yr−1 ) (m3 yr−1 ) (km)
Unzen −3.73 × 106 − 3.30 × 106 to −4.20 × 106 15
Aso −1.06 × 106 − 0.986 × 106 to −1.13 × 106 12.5

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure 9. Bootstrap test results. (a) 95 per cent confidence interval of the deviatoric stress field for restraint obtained by bootstrap. Left: direction of the
principal axis. Red, green and blue dots are the T-axis, B-axis and P-axis, respectively. The grey shaded beach ball shows the optimum solution. Right: stress
ratio R = (S2 − S3 )/(S1 − S3 ), where S1 , S2 and S3 are the maximum compressive stress, intermediate stress and minimum compressive stress, respectively.
Solid lines are the optimum solution, and broken lines are the upper and lower limits of the confidence interval. (b) Inversion results using deviatoric stress in
the confidence interval obtained by bootstrap. Left: 95 per cent confidence interval of the magnitude of inelastic strain rate. Right: inelastic strain rate tensor
within the 95 per cent confidence interval beneath the aseismic zone.

addition, comparing the dT at 10 km depth with the epicentre dis- sense to the stress change (negative dT area), which could be at-
tribution (Fig 12d), we find that earthquakes occurred in the stress tributed that the stress field in the area was formed by something
concentration area. This suggests that the inelastic deformation may other than inelastic deformation in the lower crust. However, this
contribute to the earthquake occurrences in the upper crust. In con- study cannot determine the cause of the stress field in the negative
trast, the area where the background stress field reveals opposite dT area.
298 Y. Yuasa et al.

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure 10. Top: short-wavelength component of observed strain rate. Centre: calculation results for the best-fit model. Bottom: residual between observation
and the model.

We observe a stress relaxation area where seismicity is high from has spread more extensively than the HF and that a maximum slip
the YS to Kumamoto area, including the epicentre of the 2016 Ku- of 5 m or more has been estimated in the FF. Moreover, to evalu-
mamoto earthquake. This result implies that the stress increment ate the amount of stress change quantitatively, the Mises equivalent
due to inelastic deformation in the lower crust did not contribute stress, which is a scalar quantity, is calculated (Fig. 12e). The stress
to initiate the rupture of the Kumamoto earthquake. This may indi- is loading at 2 kPa yr−1 or higher over the entire target area, while
cate that the state of stress at the Kumamoto earthquake was high the stress loading is ∼5 KPa yr−1 in the aseismic zone. Assuming
enough before the occurrence of the earthquake, and the initiation that the stress drop of inland earthquakes is ∼10 MPa, we suggest
of the 2016 Kumamoto earthquake may be attributed to a reduction that an earthquake is expected to occur at least once every several
in strength or stress perturbation caused by something other than thousand years. On the other hand, the current low seismic activity
inelastic deformation in the lower crust. In addition, the dT is posi- is considered to be caused by the aseismic zone having a higher
tive around the Futagawa fault (FF) and negative around the Hinagu strength than the surrounding area. A high velocity anomaly is seen
fault (HF). This could be related to the fact that the failure on the FF at 5–15 km depth in the AI, and the Bouguer gravity anomaly (with
 Inelastic strain rate and stress fields of Kyushu Island, Japan 299

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure 11. Estimated inelastic strain rate field for each model. (a) Result using the model in Fig. 5 shifted 0.03◦ to the east. (b and c) Results of using the
model in Fig. 5 shifted ± 0.025◦ to the N or S.

Figure 12. The background stress field, stress change at 10 and 20 km depths due to inelastic deformation in the lower crust. (a) The background stress field in
the upper crust estimated using the earthquakes in Fig. 6(a). (b) The tensor product (dT) between background stress and stress change at a depth of 10 km due
to inelastic deformation. Beach balls are the stress change tensor of each grid at the depth. Red and blue colours correspond to stress concentration zones and
stress relaxation zones, respectively. The grey areas are where the background stress field could not be estimated. (c) dT between background stress and stress
change at a depth of 10 km due to inelastic deformation. (d) Relationship between earthquakes (green dot) and the value of dT at 10 km depth. The yellow star
and black rectangles indicate the failure initiation point and the faults of the 2016 Kumamoto earthquake, respectively. (e) Mises stress change rate at 10 km
depth. The yellow star and black rectangles indicate the failure initiation point and the faults of the 2016 Kumamoto earthquake, respectively.

an assumed density of 2670 kg m−3 ) is also high, suggesting that and inelastic deformation is related to the absolute stress value, it
high-strength plutonic rocks are probably present. is important to analyse crustal deformation data in the co-seismic
The present dataset can only discuss deformation in the steady and post-seismic period, rather than performing an analysis using
state. However, as elastic deformation is caused by stress change the data in the interseismic period only.
300 Y. Yuasa et al.

5 C O N C LU S I O N S Barbot, S., Moore, J.D.P. & Lambert, V., 2017. Displacement and stress
associated with distributed anelastic deformation in a half space, Bull.
We investigate the inelastic strain rate in and around an aseismic seism. Soc. Am., 107(2), 821−855.
zone in the region of Kyushu Island, Japan, based on geodetic and Dach, R., Lutz, S., Walser, P. & Fridez., P., 2015. Bernese GNSS Software
seismic moment tensor data. The strain rate observed from the Version 5.2, User Manual, Astronomical Institute, University of Bern,
geodetic observation is the sum of the elastic and inelastic defor- Bern Open Publishing, 852.
mation in the crust. We assume that the short-wavelength component Fang, P. & Hou, G., 2019. Channel flow and fault segmentation with impli-
of the strain rate at the ground surface is induced by localized in- cations for the generation of earthquakes in the Longmenshan fault zone,
elastic deformation and that the long-wavelength component is the eastern Tibetan Plateau, J. Asian Earth Sci., 177, 107–116.

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
elastic strain rate caused by mechanical interaction from outside Hasegawa, A., Nakajima, J., Umino, N. & Miura, S., 2005. Deep structure
of the northeastern Japan arc and its implications for crustal deformation
of the model region (e.g. the slip deficit at the plate interface). In
and shallow seismic activity, Tectonophysics, 403, 59–75.
addition, the stress continuity and flow law are used to constrain Hata, M. et al., 2017. 3-D electrical resistivity structure based on geomag-
the relative magnitudes of the components of the inelastic strain netic transfer functions exploring the features of arc magmatism beneath
rate tensor in the lower crust from the deviatoric stress in the upper Kyushu, Southwest Japan Arc, J. geophys. Res., 122, 172–190.
crust. Under this assumption and the constraint, we implement an Hyodo, M. & Hirahara, K., 2003. A viscoelastic model of interseismic strain
inversion method to determine the inelastic strain rate in the lower concentration in Niigata-Kobe Tectonic Zone of central Japan, Earth
crust. As a result of analysis, an inelastic deformation zone is es- Planets Space, 55, 667–675.
timated in an area where high temperature and/or the presence of Iio, Y., Sagiya, T., Kobayashi, Y. & Shinozaki, I., 2002. Water-weakened
water is suggested. We find that the stress concentration is progress- lower crust and its role in the concentrated deformation in the Japanese
ing above this inelastic deformation zone. In addition, seismicity is Islands, Earth planet. Sci. Lett., 203, 245–253.
Jaeger, J.C., Cook, N.G.W. & Zimmerman., R.W., 2007. Fundamentals of
active in this area. The estimated inelastic strain rate tensor beneath
Rock Mechanics, 4th edn, Blackwell.
the aseismic zone exhibits the same tendency as the surrounding Kenner, S.J. & Segall., P., 2000. A mechanical model for intraplate earth-
tensors. Also, as a result of calculation, the stress change rate in quakes: Application to the New Madrid seismic zone, Science, 289, 2329–
the aseismic zone due to the inelastic deformation is 5 kPa yr−1 . 2332.
The magnitude of this value is similar to the surrounding upper Kohno, Y., Matsushima, T. & Shimizu, H., 2008. Pressure sources beneath
crust; however, we conclude that earthquakes are not occurring be- Unzen Volcano inferred from leveling and GPS data, J. Volc. Geotherm.
cause the strength is high. Both the velocity structure and gravity Res., 175(1–2), 100−109.
anomalies in the aseismic zone have positive anomalies, which may Matsumoto, S., 2016. Method for estimating the stress field from seismic
suggest the existence of high-strength plutonic rocks. moment tensor data based on the flow rule in plasticity theory, Geophys.
Our results reveal the possibility of an aseismic zone around the Res. Lett., 43, 8928–8935.
Matsumoto, S. et al., 2015. Spatial heterogeneities in tectonic stress in
AI. There is a method to detect inelastic strain change using only
Kyushu, Japan and their relation to a major shear zone, Earth Planets
interseismic period data proposed by Noda & Matsu’ura (2010). Space, 67, 172, doi:10.1186/s40623-015-0342-8.
However, for more accurate separation of both, because inelastic Matsumoto, S., Nishimura, T. & Ohkura, T., 2016. Inelastic strain rate in the
deformation is related to the absolute value of stress and elastic de- seismogenic layer of Kyushu Island, Japan, Earth Planets Space, 68, 207,
formation is caused by stress change, the use of not only interseismic doi:10.1186/s40623-016-0584-0.
data but also co-seismic and post-seismic data is required. Mazzotti, S., Le Pichon, X., Henry, P. & Miyazaki, S.-I., 2000. Full inter-
seismic locking of the Nankai and Japan-west Kurile subduction zones:
an analysis of uniform elastic strain accumulation in Japan constrained
by permanent GPS, J. geophys. Res., 105, 13 159–13 178.
AC K N OW L E D G E M E N T S ,Meneses-Gutierrez , A. & Sagiya, T., 2016. Persistent inelastic deformation
in central Japan revealed by GPS observation before and after the Tohoku-
We thank Prof. Sagiya, the anonymous reviewer and the editor for oki earthquake, Earth planet. Sci. Lett, 450, 366–371.
insightful comments that significantly improved the manuscript. We Meneses-Gutierrez, A., Sagiya, T. & Sekine, S., 2018. Crustal deformation
also thank Dr Moore for many useful comments. This study was process in the Mid-Niigata region of the Niigata-Kobe Tectonic Zone as
supported by the Ministry of Education, Culture, Sports, Science observed by dense GPS network before, during, and after the Tohoku-oki
and Technology (MEXT) of Japan, under its Earthquake and Vol- earthquake, J. geophys. Res., 123, 6072–6085.
Miura, S., Sato., T., Hasegawa, A., Suwa, Y., Tachibana, K. & Yui, S.,
cano Hazards Observation and Research Program.
2004. Strain concentration zone along the volcanic front derived by GPS
YY analysed the data and designed the research with help from observations in NE Japan arc, Earth Planets Space, 56(12), 1347–1355.
SM. SN provided displacement rate data. TM and TO acquired the Mochizuki, K. & Mitsui, Y., 2016. Crustal deformation model of the Beppu-
data. All authors have read and approved the final manuscript. Shimabara graben area, central Kyushu, Japan, based on inversion of
three-component, Earth Planets Space, 68(1), 177, doi:10.1186/s40623-
016-0550-x.
Moore, J.D.P. et al., 2017. Imaging the distribution of transient viscosity
REFERENCES after the 2016 Mw 7.1 Kumamoto earthquake, Science, 356, 163–167.
Abe, Y., Ohkura, T., Shibutani, T., Hirahara, K. & Kato, M., 2010. Crustal Noda, A. & Matsu’ura, M., 2010. Physics-based GPS data inversion to
structure beneath Aso Caldera, Southwest Japan, as derived from receiver estimate 3-D elastic and inelastic strain fields, Geophys. J. Int., 182, 513–
function analysis. J. Volc. Geotherm. Res., 195(1), 1–12. 530.
Altamimi, Z., Rebischung, P., Metvier, L. & Collilieux, X., 2016. ITRF2014: Noda, A., Saito, T. & Fukuyama, E., 2018. Slip-deficit rate distribution
a new release of the International Terrestrial Reference Frame modeling along the Nankai trough, Southwest Japan, with elastic lithosphere and
nonlinear station motions, J. geophys. Res., 121, 6109–6131. viscoelastic asthenosphere, J. geophys. Res., 123, 8125–8142.
Asano, K. & Iwata, T., 2016. Source rupture processes of the foreshock and Okada, Y., 1992. Internal deformation due to shear and tensile faults in a
mainshock in the 2016 Kumamoto earthquake sequence estimated from half-space, Bull. seism. Soc. Am., 82, 1018–1040.
the kinematic waveform inversion of strong motion data, Earth Planets Pollitz, F.F., Kobayashi, T., Yarai, H., Shibazaki, B. & Matsumoto, T., 2017,
Space, 68, 147 , doi:10.1186/s40623-016-0519-9. Viscoelastic lower crust and mantle relaxation following the 14–16 April
 Inelastic strain rate and stress fields of Kyushu Island, Japan 301

2016 Kumamoto, Japan, earthquake sequence, Geophys. Res. Lett., 44, Wang, Z., Zhao, D., Liu, X. & Li, X., 2017. Seismic attenuation tomography
8795–8803. of thesource zone of the 2016 Kumamoto earthquake (M7.3), J. geophys.
Rybacki, E. & Dresen, G., 2004. Deformation mechanism maps for feldspar Res., 122, 2988–3007.
rocks, Tectonophysics, 382, 173–187.
Saiga, A., Matsumoto, S., Uehira, K., Matsushima, T. & Shimizu, H., 2010.
Velocity structure in the crust beneath the Kyushu area, Earth Planets
A P P E N D I X A : F I LT E R R A D I U S
Space, 62, 449–462.
Sagiya, T., Miyazaki, S. & Tada, T., 2000. Conyinuous GPS array and
DEPENDENCY
present-day crustal deformation of Japan, Pure appl. Geophys., 157, We confirmed the long-wavelength component of the strain rate

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
2303–2322 obtained for each different filter radius. The tested radii were 50,
Savage, M. et al., 2016. Stress, strain rate and anisotropy in Kyushu, Japan,
60, 70, 80 and 90 km. The results obtained for 50 and 90 km are
Earth planet. Sci. Lett., 439, 129–142.
shown in Fig. A1. We can confirm that there is almost no difference
Shen, Z.K., Jackson, D.D. & Ge, B.X., 1996. Crustal deformation across and
beyond the Los Angeles basin from geodetic measurements, J. geophys. between the results obtained for both radii. Therefore, it is expected
Res., 101, 27 957–27 980. that the filter radius selection does not significantly affect the result.
Shimizu, H. et al., 2016. Urgent joint seismic observation of the 2016
Kumamoto earthquake seismic activities and their background, in JPGU
2016 Ann. Meeting, MIS34-02, Chiba, Japan. A P P E N D I X B : M O M E N T T E N S O R D ATA
Stuart, W.D., Hildenbrand, T.G. & Simpson, R.W., 1997. Stressing of the SET
New Madrid Seismic Zone by a lower crust detachment fault, J. geophys.
Res., 102, 27 623–27 633 We confirmed the magnitude–frequency of the moment tensor cata-
Tada, T., 1984. Spreading of the Okinawa trough and its relation to the crustal logue used for estimating the deviatoric stress field. In the method of
deformation in Kyushu, J. Seismol. Soc. Jpn., 37, 407–415 (in Japanese Matsumoto (2016), since the deviatoric stress is estimated by sum-
with English abstract). ming the moment tensors, the estimated stress resembles the mo-
Tada, T., 1985. Spreading of the Okinawa trough and its relation to the ment tensor of a large earthquake if the number of microearthquakes
crustal deformation in Kyushu (2), J. Seismol. Soc. Japan, 38, 1–12 (in
is small. Fig. B1 shows the magnitude–frequency of each box
Japanese with English abstract).
Tanaka, A., Yamano, M., Yano, Y. & Sasada, M., 2004. Geothermal gradient
(Fig. 5). We could use many microearthquakes for estimating the
and heat flow data in and around Japan (I): appraisal of heat flow from deviatoric stress. Furthermore, Fig. B2 shows the estimated devi-
geothermal gradient data, Earth Planets Space, 56, 1193–1196. atoric stress using earthquakes less than M 2.5. Similar results to
Terakawa, T. & Matsu’ura, M., 2010. The 3-D tectonic stress fields in Fig. 12(a) were obtained. Therefore, we can conclude that the de-
and around Japan inverted from centroid moment tensor data of seismic viatoric stress estimated in this study does not depend on a specific
events, Tectonics, 29, TC6008. earthquake alone.
Townend, J. & Zoback, M.D., 2006. Stress, strain, and mountain building in In contrast, especially in the sea area, there are some earthquakes
central Japan, J. geophys Res., 111, B03411. for which we could not determine the mechanism. Fig. B3 shows
Yamasaki, T. & Seno, T., 2005. High strain rate zone in central Honshu a comparison of the cumulative moment tensor of all the detected
resulting from the viscosity heterogeneities in the crust and mantle, Earth
earthquakes and the earthquake for which we could estimate the
planet. Sci. Lett., 232, 13–27
Wallace, L.M., Ellis, S., Miyao, K., Miura, S., Beavan, J. & Goto, J., 2009.
mechanism. It is expected from Fig. B3 that earthquakes for which
Enigmatic, highly active left-lateral shear zone in southwest Japan ex- we could not determine a mechanism exist regardless of magni-
plained by aseismic ridge collision, Geology, 37(2), 143–146. tude. However, since the same result is obtained by the estimation
Wang, K., 2000. Stress–strain ‘paradox’, plate coupling, and forearc seis- excluding earthquakes larger than M 2.5, it is expected that the in-
micity at the Cascadia and Nankai subduction zones, Tectonophysics, 319, fluence of the lack of earthquakes for which we couldnot estimate
321–338. the mechanism is small.
302 Y. Yuasa et al.

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020

Figure A1. Short-wavelength (bottom) and long-wavelength (top) components of strain field for moving average filter radii of 50 and 90 km. For analysis,
60 km was used.
 Inelastic strain rate and stress fields of Kyushu Island, Japan 303

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure B1. Magnitude–frequency for each box of the mechanism solution data set used to estimate the deviatoric stress field. The place of each figure
corresponds to the place of each box in Fig. 5.
304 Y. Yuasa et al.

Downloaded from https://academic.oup.com/gji/article-abstract/221/1/289/5698808 by Kyushu University Central Library (S) user on 07 July 2020
Figure B3. Cumulative moment tensor from 2004 January 1 to 2014 March
Figure B2. The estimated deviatoric stress field using only earthquakes 31. The green and red lines are the cumulative moment tensors of all the
less than M 2.5. We could not estimate the deviatoric stress southwest of the detected earthquakes and the earthquake for which we could estimate the
target area, owing to the absence of earthquakes. mechanism, respectively.