Inelastic strain in the hypocentral region of the 2000 Western Tottori earthquake (M 7.3) inferred from aftershock seismic moment tensors

Abstract

Inelastic deformation due to seismic activity is an important signal that reflects fault evolution. In particular, aftershock sequences indicate the evolution of damage in a medium that has experienced a large earthquake. Herein, we discuss the inelastic strain rate surrounding the fault that produced the M 7.3 Western Tottori earthquake in 2000 using long-term aftershock analysis. To obtain high-resolution focal mechanisms 18 years after the earthquake occurrence, we conducted dense seismic observations in the focal area. The inelastic strain rate estimated from the aftershock seismic moment tensor data showed spatial variations within a range of 10⁻⁷–10⁻¹¹ per year, 18 years after the main shock. By comparing the inelastic strain rates from immediately after the earthquake and 18 years later, we detected the increase in the spatial variations in the inelastic strain rate; the variations are as small as 10² (= 10⁻⁵/10⁻⁷) for the early stage but as large as 10⁴ (= 10⁻⁷/10⁻¹¹) for the later period. In addition, the decay of the rate during these two periods varied spatially from spatial bin to bin. Certain bins in the northern segment of the earthquake fault, southern edge of the fault, and surrounding the location of the preceding swarm activity to the M 7.3 event showed slower decay rates than the inverse of the lapse time since the occurrence of the M 7.3 earthquake. We modeled this decay rate change as the relaxation response of a power-law fluid to an elastic strain input from the large earthquake. Most parts of the fault can be explained by this model. However, the areas with low decay rates suggest the presence of a dragging mechanism, such as aseismic slip, at or around these locations.

Main text

Introduction

Crustal strength is the key to understanding the development of the crust, including large earthquakes, because it controls the medium behavior under tectonic stress loading in the crust. Stress and strain fields in the crust are parameters that can be detected via geophysical observation. Recent seismic network and global navigation satellite system (GNSS) observations have provided high-quality data that can be used to investigate stress and strain rates in the crust. GNSS data provide the total strain rate at the ground surface, which consists of elastic and inelastic strain. Earthquake activity also contributes to inelastic deformation in a seismogenic zone. Kostrov (1974) proposed a method to determine the average strain rate in a region by summing the moment tensors released in the affected volume. Numerous studies have applied this method to analyze seismic activity data during investigations of crustal deformation and tectonic evolution (Wesnousky et al. 1982; England and Molnar 2005). Earthquakes are an inelastic phenomenon in a medium, where a large number of earthquakes in a volume (ΔV) create inelastic strain in that volume. In this study, we discuss the inelastic strain produced by earthquakes and estimate the inelastic strain field in the hypocentral area of the M 7.3 Western Tottori Earthquake (W-Tottori EQ) that occurred in 2000.

The earthquake occurred in the center of the San-in district in southwestern Japan on October 6, 2000, whose aftershock activity continues to this day. According to a total strain rate study based on GNSS data, Nishimura and Takada (2017) found a strain concentration zone in the north of the W-Tottori EQ’s hypocentral area. They identified this zone as part of the Northern Chugoku shear zone (Gutscher and Lallemand 1999). The earthquake occurred due to left-lateral strike-slip faulting and consisted of two major and six minor faults, inferred from precisely determined aftershock distributions (Yukutake and Iio 2017). In this area, seismic activity that preceded the main shock was observed in 1989, 1990, and 1997 (Shibutani et al. 2002). The large co-seismic slip section of the fault appears to be distributed around the preceding seismic activity, as estimated from strong motion records (Iwata and Sekiguchi 2002). This suggests that the strength of this area is lower than other areas. Yukutake et al. (2007) estimated, in detail, the deviatoric stress field derived from focal mechanism solutions, showing that the stress field is a strike-slip regime with the principal compression direction-oriented WNW–ESE. The fault also showed remarkable spatial heterogeneity around the large co-seismic slip area, which was attributed to stress changes caused by the co-seismic slip. These stress and strength heterogeneities detected by previous studies can be interpreted as time slices of the fault system maturity. Detailed analyses of the aftershock activity from an inelastic strain evolution perspective may also provide important information on fault development. Therefore, we investigate how aftershock activity in the Tottori area damaged the crust with a well-known detailed stress field and tectonic structure. However, we note that the maturing process discussed in this study is not always applicable on a geologic time scale. Based on geologic surveys performed in this area, numerous small faults show slip directions consistent with the W-Tottori EQ (e.g., Kobayashi et al. 2003). However, we have not observed large mature faults in the aftershock area, which suggests that the fault system in this area may be at an earlier developmental stage. Therefore, this study is able to document a mechanism for younger fault system development.

The previous geodetic studies described above reported that there was a concentration of strain in the San-in area, which also includes the main shock hypocenter location. These studies investigated the total strain at the surface using GNSS observations and modeled the velocity field of a shear zone running in the E–W direction. The detected total strain rates were on the order of 10⁻⁹ to 10⁻⁷ per year, with scale lengths greater than 20 km due to station separation. The inelastic strain estimated in this study corresponds to the hypocentral zone at depths ranging from 0 to ~ 20 km. A comparison of these sets of observations is important when considering crustal deformation after the earthquake. However, inelastic strain due to seismic activity creates elastic strain at the surface. This elastic strain may be small as no geodetic study has detected significant deformation in the hypocentral area after the main shock. In addition, the length scale reported in Nishimura and Takada (2017) is larger than that used in our target area (see below). Therefore, in this study, we only discuss inelastic behavior in the hypocentral area.

Data and observations

In this study, we attempted to estimate the inelastic deformation due to aftershocks because aftershock activity is an inelastic response to step-like stress loading in the Yukutake and Iio 2017 surrounding the fault caused by the main shock. To evaluate this response, we required temporal changes in the inelastic strain rate. As described above, seismic moment tensors are essential data to estimate the rate after the occurrence of the main shock. For the W-Tottori EQ, dense seismic observations were performed immediately after the main shock (Yukutake and Iio 2017). However, we also required seismic moment tensor data for earthquakes during a second period to estimate the rate change. A high detectability of aftershocks and focal mechanisms was required for this study due to the decline in aftershock activity based on the Modified Omori’s law (Utsu 1957). We conducted a dense seismic observation (referred to as a “0.1 Manten seismic observation” in the Crustal Dynamics research project supported by the Ministry of Education, Culture, Sports, Science, and Technology [MEXT], Japan) in the hypocentral area of the W-Tottori EQ using the station distribution shown in Fig. 1. The word “Manten” refers to 10,000 points in Japanese. Thus, “0.1 Manten” indicates that 1,000 stations were used. The 1,000 stations were deployed throughout the study area for an approximately 1-year observation period from March 2017 to April 2018. A seismometer (natural frequency of 2 or 4.5 Hz) and recording system (composed of a memory card and 48 type-D alkaline batteries) were installed at each station. Specifics of the observation method are described in detail in Matsumoto et al. (2020). Vertical component seismometers were used for this deployment. Data collected during this deployment were merged with three-component seismograms recorded by Hi-net, which is operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). Automatic picking and focal mechanism determination processes were applied to the continuous seismic records. The hypocenters of 7,315 events, which were determined with root mean square (RMS) residuals for P-wave arrival times smaller than 0.15 s in the hypocenter locating inversion, were used in this study. A smoothed depth-dependent velocity structure was adopted in this study, as in previous studies (Yukutake and Iio 2017), which was used for hypocenter determination. We only used events that were in the area covered by the 0.1 Manten stations. Focal mechanisms for the pick data were also determined using the HASH algorithm described in Hardebeck and Shearer (2002). Figure 2 shows the epicenter map and a magnitude–frequency diagram for the events determined in this study. The diagram shows that the completeness of the detection processing was at an approximate magnitude of − 0.5 for 3641 events within the network. High-quality focal mechanism solutions, which comprise ranks A and B under the criteria defined by HASH, were also determined. We successfully determined the focal mechanisms of most small earthquakes, i.e., more than 80% of the mechanisms were determined even in the magnitude range from 0.5 to 1.0. We used focal mechanism data from 2315 rank A and B events in the analyses.

Fig. 1

Left) Target area and moment tensor solutions since the 2000 Western Tottori Earthquake, until Apr 2018, recorded by F-net, NIED. The beach balls show the moment tensor solutions plotted on the lower hemisphere. Crosses show Hi-net and NIED stations used in this study. Triangle indicates an active volcano. The red star indicates the epicenter of the main shock. Dotted lines indicate prefectural boundaries in the region. Right) Station distribution of the “0.1 Manten” dense seismic observations in blue circles. The location of the right map corresponds to the dashed rectangle in the left figure. Other symbols are the same as in the left figure

Fig. 2

Left) Map of the epicenters of the events detected during the observation period. The white dotted line indicates the area covered by the network. Right) Histogram of the magnitude of the detected earthquakes. Light blue bars show the number of earthquakes detected. Detection completeness is an approximate magnitude of − 0.5. Dark blue bars indicate the number of focal mechanisms with enough accuracies

To evaluate the accuracy and reliability of the automatic hypocenter and focal mechanism determinations, we compared the results with manually picked data. We then calculated the differences in the hypocenter locations and focal mechanism solutions for the results determined using manual and automatic picking. Fifty events with various magnitudes were selected from the hypocenter data, followed by a comparison. Figure 3 shows the differences in the locations and tensors calculated from focal mechanisms with magnitudes of − 0.5 to 3.2. We only used the auto-picked hypocenter results that had an RMS residual of less than 0.15 s and focal mechanisms with A and B HASH ranks. The colors plotted in Fig. 3 are classified by the focal mechanism HASH rank using the manually picked data (reference data). For reference focal mechanisms with A and B ranks, the horizontal differences in the hypocenters are less than 0.4 km for magnitudes greater than 0. The vertical differences in location are less than 1 km. The reliability of the focal mechanism was evaluated by calculating the tensor dot product (DTensor, dT) of the automatic and manual solutions. The tensor dot product between two order tensors, i.e., m and m’, can evaluate the similarity of the two tensors (Terakawa 2017), which is expressed as follows:

$${\text{dT}}\left( {m,m{\prime}} \right) = \sum\nolimits_{i = 1}^{3} {\sum\nolimits_{j = 1}^{3} {m_{ij} } } m{\prime}_{ij} /\sqrt {\sum\nolimits_{i = 1}^{3} {\sum\nolimits_{j = 1}^{3} {m_{ij}^{2} } } } \sqrt {\sum\nolimits_{i = 1}^{3} {\sum\nolimits_{j = 1}^{3} {m{\prime}_{ij}^{2} } } }$$

(1)

Fig. 3

Results of the reliability test for the auto-picking processing for various magnitude ranges. Manually picked phases and polarity data were used as a reference. Top) Tensor product showing the discrepancies between tensors by the auto-picking and reference focal mechanisms. Middle) Location differences in horizontal and Bottom) depth directions. The symbol colors indicate the quality of the focal mechanism data for the reference event using the HASH rank

The products of the focal mechanisms for the reference events with A and B ranks are above 0.95. This implies that the focal mechanisms are highly similar to the reference mechanism. Therefore, in this study, we adopted a focal mechanism solution for an event with a magnitude greater than 0.5 using the automatic picking procedure.

We only used the focal mechanisms of events that occurred within the “0.1 Manten” network-covered area. In addition to the data described above, for a period from October 15, 2000, to November 30, 2000, we used focal mechanisms that were estimated using the HASH algorithm for manually picked polarity data (Yukutake and Iio 2017). Consequently, we analyzed these datasets and their corresponding focal mechanisms immediately after and 18 years after the main shock of the W-Tottori EQ. Figure 4 shows a map view and vertical cross-section of the P-axis distribution for the events immediately after and 18 years after the occurrence of the main shock. This cross-section is 5 km wide, taken from the earthquakes distributed along the main shock fault (striking N30 °W), and includes the hypocenter. Most events occurred around the fault. The general tendency of the P-axis directions obtained by the “0.1 Manten” observations is similar to that of the events that occurred immediately after the main shock. However, the aftershock activity differs along the fault sections.

Fig. 4

Map view and vertical cross-section of P-axis distribution along the main shock fault. P-axes of the events located in the dashed rectangle on the map are plotted on the cross-section; left and right figures show data immediately after and 18 years. after the main shock, respectively. The star indicates the hypocenter of the main shock. The area where co-seismic slip is larger than 2.5 m (Iwata et al. 2002) is indicated by a pink dashed line

Inelastic strain due to seismic activity

Matsumoto et al. (2016) estimated the inelastic strain rate around the hypocentral region of the 2016 Kumamoto earthquake (M 7.2), suggesting that the inelastic strain was related to stress loading on a fault related to the large earthquake. In this study, we adopted this approach, following their method to estimate the inelastic strain and thus detect temporal changes in inelastic behavior around the main shock of the W-Tottori EQ.

A region with many faults displaced by earthquakes under triaxial stress can be treated as an inelastic deformed body from a macroscopic perspective. Matsumoto (2016) discussed the relationship between the deviatoric stress field and seismic moment tensors through inelastic deformation. Considering a body with numerous earthquakes in a ΔV deformed by earthquake faulting, we can express the inelastic deformation using moment tensors distributed throughout the volume (Aki and Richards 1980; Ichihara et al. 2016; Kusakabe et al. 2016). The relationship between inelastic strain, $\varepsilon_{ij}^{p} ,$ in the volume and the moment density tensor can be expressed using the compliance, C_ijkl, of the volume (Noda and Matsu’ura 2010) as follows:

$$\varepsilon_{ij}^{p} = \varvec{C}_{ijkl} m_{kl} ,$$

(2)

where m_kl is the moment density tensor. The term $\varepsilon_{ij}^{p}$ corresponds to the increment of inelastic strain in the medium after the events occurred in ΔV (the increment of strain). This inelastic strain is identical to the plastic strain in a medium because earthquake faulting is macroscopically a plastic phenomenon, as described in Kostrov (1974) and Matsumoto (2016). Assuming a uniform distribution of earthquake point sources in the volume, the average moment density tensor in the volume can be written as follows:

$$m_{ij} = \sum\nolimits_{k=1}^{K} {M_{ij}^{k} /\Delta V} ,$$

(3)

where the seismic moment tensor of the kth earthquake is $M_{ij}^{k}$ and K is the number of earthquakes. Then, we can obtain the following:

$$\varepsilon_{ij}^{p} = \varvec{C}_{ijpq} \sum\nolimits_{k=1}^{K} {M_{pq}^{k} /\Delta V}$$

(4)

This formula was originally proposed by Kostrov (1974) to evaluate strain in a volume. Equation (4) indicates that the average strain increment can be estimated by summing the seismic moment tensors in a volume for the duration of the target period. In this study, we attempted to obtain the inelastic strain rate calculated from the summation of the seismic moment tensors divided by the period of the dataset. The period length was 47 days for the dataset collected immediately after the main shock and 395 days for the data collected 18 years after the main shock. The procedure for obtaining the inelastic strain rate and its change in a spatial bin from the real data is as follows:

• 1) Determine the focal mechanisms of events that occurred in the bin from the polarities of the first P-wave arrival onsets.

• 2) Obtain the moment tensors from the focal mechanisms and magnitudes of the events.

• 3) Sum the moment tensors of the events in the bin.

• 4) Calculate the inelastic strain from the total moment tensor, compliance of the medium, and volume of the bin (Eq. (4)).

• 5) Calculate the inelastic strain rate by dividing the strain by the duration of the observation period.

• 6) Compare the rate for 18 years after the main shock with the rate for immediately after the main shock.

Spatial and temporal inelastic strain distributions

We generated a moment tensor dataset from the focal mechanisms and earthquake magnitudes (Ma). The tensor shape was obtained from the estimated focal mechanism. The scalar moment (M0) was estimated using an empirical relationship used in Matsumoto et al. (2016): (log(M0) = 1.15 Ma + 10.548). The strain rate tensors were estimated in 5-km-deep spatial blocks set at horizontal intervals of 3 km. In addition, the same procedure was applied to a block distribution shifted by a half block to the north and east from the origin to smooth the results. The sum of the seismic moment tensors in each bin volume, ΔV (3 km × 3 km × 5 km), was calculated for any bin that obtained more than ten focal mechanisms. We assumed that the compliance of the homogeneous medium was the inverse of the rigidity (shear modulus) of 33 GPa, i.e., ordinary crustal rocks. Figure 5 shows the estimated inelastic strain rate for the two observation periods. The strain tensor shape for 18 years after the main shock exhibits a spatial pattern that is similar to the pattern present immediately after the main shock.

Fig. 5

Map showing the inelastic strain rate tensors estimated from the two periods of data in spatial blocks. The tensor is plotted on the lower hemisphere in a similar way to the seismic moment tensor. The upper and lower maps show the results from Oct 15, 2000 to Nov 30, 2000 and from Mar 2017 to Apr 2018, respectively. The figures on the left, middle, and right correspond to the depth ranges shown at the top of the figure. The symbol sizes and colors reflect the magnitude of the inelastic strain rate, as shown in the legend at the lower left of each map. Gray dots are epicenters within the depth range for each period

The “0.1 Manten” observations were endowed with the high detectability of small earthquakes, as described above. Therefore, the inelastic rate change from several tens of days to 18 years after the main shock can be estimated in the aftershock area. We extracted the spatial blocks where the inelastic strain had been obtained for both periods. All events we selected here had a distance from the W-Tottori EQ fault plane over 1.5 km. Here, we calculated the rate from the tensor magnitude defined as $\sqrt {\frac{{\dot{\varepsilon_1 }^{2} + \dot{\varepsilon_2 }^{2} + \dot{\varepsilon_3 }^{2} }}{2}}$, where $\dot{\varepsilon_1 }, \dot{\varepsilon_2 },$ and $\dot{\varepsilon_3 },$ are the maximum, moderate, and minimum principal strain rates, respectively. Figure 6 shows the range of inelastic strain rates for the time lapse after the main shock in the derived block solutions (see also Additional file 1: Fig. S2b). The rates in respective blocks show spatial variations within a range of 10⁻⁵–10⁻⁷ per year for data from 2000 and within a range of 10⁻⁷–10 ⁻¹¹ per year for 2017–2018 data. The range of the strain rates broadened to the order of 10⁴ over 18 years, implying that the change in the strain rate varies by section around the earthquake fault. To model the temporal change in the strain rate, we estimated the average value within the entire target area. The F-net centroid moment tensor (CMT) catalog from NIED provided the moment tensor solutions for events with relatively larger magnitudes than those determined by this study in the target area. Therefore, obtaining spatial variations from the F-net CMT data was difficult. The average strain rate was estimated using the tensor magnitude of the CMT data divided by the lapsed time and volume size of the entire region (30 km along the earthquake fault × 10 km wide × 10 km deep). We adopted a lapsed time from the event origin time to the preceding event that obtained the CMT solution. Figure 6 shows the average rates, which declined with the elapsed time from the main shock. A power-law function fitted to the average inelastic strain rate ($\dot{\varepsilon_i } \propto t^{ - 1.2 \pm 0.1}$, where t is the lapse time and the power-law exponent is − 1.2, with a standard error of 0.1) can explain the decline, as shown in Fig. 6. We calculated the decay rate at each spatial bin from the inelastic strain rates in the two periods (i.e., immediately and 18 years after the main shock). Here, we defined the “P value” in the power-law relation, $\dot{\varepsilon_i } \propto t^{ - P}$.

Fig. 6

Inelastic strain rate versus time elapsed since the main shock occurred. Red and green symbols indicate the estimated inelastic strain rates for the period immediately after and 18 years after the main shock, respectively. These symbols are plotted in the middle of the time-span used for calculating inelastic strain. The earlier time lapse window is indicated by a red arrow. The latter window corresponds to the width of the symbol. The blue circle indicates the average strain rate for the entire target area, estimated from the F-net moment tensor solution. The solid line is a regression line for the average rate change data

Figure 7 shows the variations in the decay rate along the earthquake fault, where we plot the P value along the fault. The P value was taken as the average of the three blocks closest to the fault in the perpendicular direction.

Fig. 7

a Vertical cross-section of the P value distribution along the fault plane. The shaded area shows a block where the number of moment tensors used in the analysis is less than 10 in either period (2000 or 2017). The pink dashed line indicates the contour line of co-seismic slip larger than 2.5 m, from Iwata et al. (2002). The star indicates the location of the W-Tottori EQ hypocenter. The activity that preceded the W-Tottori EQ from Shibutani et al. (2002) is plotted as open circles. b The 90% confidence interval of the inelastic strain rate distribution for the period in 2000 (left) and in 2017–2018 (right). Other symbols are the same as in Fig. 7a

To smooth the results, we applied a “surface” function in the GMT software package (Smith and Wessel 1990) to the results. Reliable values can be observed in unshaded areas and the blue “DTensor” color. We only discuss results that satisfy the following conditions: (1) the number of earthquakes used for the estimation is greater than 10 for both periods and (2) the 90% confidence interval of the inelastic tensor estimation is greater than 0.6 as the tensor dot product (DTensor). The confidence intervals were evaluated for each spatial bin in the range of 90% in the descending tensor product order distribution between the optimal value and trial solutions via bootstrap resampling of the moment tensor data. Small decay areas of the inelastic strain rate are observed along the central, northern, and southern edges of the fault plane. The low P value at the fault edge can be attributed to a complex stress field at the edge of the co-seismic slip. However, low P values are also observed approximately 5 and 10 km north of the hypocenter. The aftershock activity was scattered along the northern segment of the fault, as shown in Fig. 4, suggesting that the rupture strength was low in this section. Shibutani et al. (2002) showed the swarm activity that preceded (circles in Fig. 7) the W-Tottori EQ located around the hypocenter of the main shock. This swarm is in a low P value area, which can also imply a low strength in the medium. However, the state of stress in the northern segment of the fault is higher than in the south, such that seismic activity continued for a long period. Therefore, we attempted to model this activity.

Discussion

Numerous studies have modeled the aftershock sequence as a behavior of continuum damage material. For example, Ben-Zion and Lyakhovsky (2006) consider a model for an aftershock sequence using damage rheology. Nanjo and Turcotte (2005) show that an aftershock sequence can be explained by a medium with power-law creep properties. Nanjo (2007) proposes a model introducing yield stress, which can also explain aftershocks as a response to a fiber-bundle material. These models consider aftershock activity as a phenomenon caused by elastic strain loaded by a large earthquake in the surrounding media, which is then replaced by inelastic strain due to aftershock failure. Previous studies have typically attempted to explain the decay in the number of aftershocks using a modified form of well-known Omori’s law based on rock rheology. Inelastic strain estimated in this study from the seismic moment tensors can directly connect to a physical model because transformation from a physical parameter (such as strain) to the number of earthquakes requires strong assumptions. Based on previous studies, the power-law decay in the inelastic strain rate (i.e., $\dot{\varepsilon_i } \propto t^{ - p}$), which was observed in this study, can be modeled by an inelastic medium. Therefore, following Nanjo (2007), we considered that aftershock activity was a response of the inelastic medium to elastic strain loaded by an earthquake. We also assumed that inelastic strain due to aftershock activity followed the constitutive relation for a power-law fluid. Inelastic deformation of Earth’s interior is often modeled by a rock creep mechanism, such as dislocation, diffusion, or pressure solution creep. The constitutive relation for dislocation creep is similar to that of a power-law fluid (Rutter 2010). However, this can only be observed under high-temperature mantle conditions. The pressure solution creep model (Rutter 1976; Shimizu 1995) is applicable for inelastic deformation in the crust. This model is not able to explain the long-term tail observed in the aftershock sequence. Diffusion creep, which is a special case of a power-law fluid, also shows a rapid decline with lapsed time. Thus, the power-law fluid model is appropriate as one of the physical models that can explain inelastic phenomena in the Earth. Therefore, we adopted the relation for a power-law fluid in this study.

Relaxation of the medium due to inelastic strain can be expressed under the conditions that (1) elastic strain is suddenly applied to a volume at t = 0 by a main shock, (2) no tectonic stress above the yield stress is loaded onto the volume, and (3) total strain is the sum of the elastic and inelastic strains, which is constant for t > 0. A relationship between the inelastic strain rate and time then becomes identical to the observed relation (see detailed expression in the Additional files 1):

$$\dot{\varepsilon }_{i} \propto t^{ - p} ,$$

(5)

where $\dot{\varepsilon }_{i}$ is the inelastic strain rate and t is the lapsed time after the main shock. The P value in Eq. (5) relates to the stress–strain relationship in the medium:

$$\dot{\varepsilon }_{t} = \dot{\varepsilon }_{e} + \frac{1}{\tau }\left( {\frac{{\sigma - \sigma_{y} }}{E}} \right)^{n} ,$$

(6)

$$P = \frac{n}{n - 1},n > 1,$$

(7)

where $\dot{\varepsilon }_{t}$ and $\dot{\varepsilon }_{e}$ are the total and elastic strain rate, respectively, $\sigma$ and $\sigma_{y}$ are the stress and yield stress, respectively, E is Young’s modulus, and τ is the characteristic time of the medium (see Additional file 1). Equations (6) and (7) are equivalent to a case in which there is a sudden decrease in the yield stress. For example, the strength of the fault can be reduced due to fluid injection in the focal area, as suggested by Sibson (1990). In this case, $\sigma - \sigma_{y}$ is also positive. The n term in Eq. (6), which is a parameter for a power-law fluid, expresses characteristics of the fluid behavior. For example, when n = 1 and n > 1, Eq. (6) becomes the equation for a medium with diffusion and dislocation creep properties, respectively. In general, τ is related to the medium properties and conditions, such as temperature. However, discussing these parameters is out of the scope of this study because we only consider the temporal decay of inelastic strain. Therefore, we will not discuss the characteristic time or other medium properties further in this paper.

The decay of the number of aftershocks with elapsed time from an earthquake can be expressed by the modified Omori law (i.e., n (t)= K(t + c)^−p, where n (t) is the number of earthquakes at the elapsed time, t, and K, c, and p are constants). Nanjo (2007) shows that p is related to n based on power-law fluid modeling, i.e., p = n/(n−1). Therefore, by comparing Eq. (7) with this relationship, we can obtain the relationship among the P, p, and n values as follows:

$$p = P = \frac{n}{n - 1},n > 1.$$

(8)

Nanjo et al. (2007) investigate the p value for numerous inland earthquakes, obtaining p values from 0.08 to 1.25. They also evaluated the n value in Eq. (6), which ranged from 5 to 13 for the aftershock sequences described in Nanjo et al. (2007). They suggest that this range is also possible for the behavior of a medium with a power-law constitutive relationship between stress and strain rate in the crust, despite being higher than the case for a medium with dislocation and diffusion creep properties.

In this study, we obtained similar ranges for their estimations. Most fault sections appear to behave as a power-law fluid with n values between 2 and 10 (see description in Additional file 1). In contrast, we observe that the P values of certain fault sections were small based on the standard error estimation (see Additional file 1). As the n value becomes infinitely large for P = 1, P > 1 can be explained by elastic strain step response for the power-law fluid medium (black to dark blue region in Fig. 7a). However, we must consider other factors that can cause the slow decay rate observed along the central and northern sections of the fault.

For the slow decay of the inelastic strain rate, we considered the stress loading on the medium. The slow decay of the inelastic strain rate (i.e., P value smaller than 1.0) was observed in limited areas around the fault, suggesting that it was not caused by regional tectonic stress loading. We observed slow decay in the northern and central sections of the fault. The inelastic strain rate in these sections remained high, despite being 18 years after the main shock. One possible cause of this is that faults present in, around, or beneath the earthquake experienced aseismic slip after the main shock of the W-Tottori EQ. In this case, the decay may appear slow because the rate was estimated using inelastic strain rates at only two periods (see Additional file 1). In addition, the post-seismic slip may occur in small areas that are insufficient to affect a wide area of the fault, which may explain the low P value observed for limited sections of the fault. Although this possibility requires low strength for fault slip in the post-seismic slip area, as an example, swarm activity (Shibutani et al. 2002) that occurred in 1989, 1990, and 1997 may have been a response to a small fault slip. The area of the preceding activity was in the northwestern region of a large slip area during the W-Tottori EQ and corresponded to the low P area. This suggests that an easily slipping fault is present in this area. In this case, we consider that the activity is a response to either the larger earthquake or localized aseismic slip. For further studies, investigating spatiotemporal variations in the total inelastic strain for a longer duration, including the co-seismic slip and preceding swarm, may enable us to model the physical properties of the fault, despite the requirement of long-term accurate moment tensor data.

Conclusions

In this study, we investigated the inelastic strain rate surrounding the fault associated with the W-Tottori EQ by performing dense seismic observations. We obtained high-quality focal mechanisms in the aftershock area of the 2000 Western Tottori earthquake (M7.3). Our conclusions are as follows:

• By introducing a transform seismic moment tensor to the inelastic strain, we can model the aftershock deformation process as a power-law medium.

• The inelastic strain rate estimated from the aftershock seismic moment tensor data shows spatial variations ranging from 10⁻⁷ to 10⁻¹¹ per year 18 years after the main shock.

• Comparison of the above inelastic strain rate with the rate immediately after the main shock shows that the $\dot{\varepsilon }_{i}$ decay rate is low at the northern end of the co-seismic large slip area, where the swarm that preceded the main shock occurred.

• Most fault sections behave as power-law fluids with n = 2–10. The small rate observed in the north can be attributed to behavior caused by loading from small faults.

These results show that aftershock activity and fault slip caused the medium to deform and release elastic strain around a large earthquake fault. The deformation was spatially heterogeneous, which may be related to the fault development process. To capture this development process, we must combine geodetic observations and discuss ancient fault distributions detected by geologic surveys in future studies.