Electrical conductive fluid-rich zones and their influence on the earthquake initiation, growth, and arrest processes: observations from the 2016 Kumamoto earthquake sequence, Kyushu Island, Japan

Abstract

Crustal earthquake ruptures tend to initiate near fluid-rich zones. However, it is relatively unknown whether fluid-rich zones can further promote or arrest these ruptures. We image the electrical resistivity structure around the focal area of the 2016 Kumamoto earthquake sequence by using 200 sites broadband magnetotelluric data, and discuss its quantitative relationship to earthquake initiation, growth, and arrest processes. The ruptures that initiated along the outer edge of the low-resistivity fluid-rich zones (< 30 Ωm) tended to become large earthquakes, whereas those that initiated either distal to or within the fluid-rich zones did not. The ruptures were arrested by high-temperature (> 400 °C) fluid-rich zones, whereas shallower low-temperature (200–400 °C) fluid-rich zones either promoted or arrested the ruptures. These results suggest that the distribution of mid-crustal fluids contributes to the initiation, growth, and arrest of crustal earthquakes. The pre-failure pressure/temperature gradient (spatial difference) of the pore fluids may contribute to the rupture initiation, propagation, and arrest.

Introduction

Electrically conductive fluid-rich zones have been detected in the mid-crust of several island arc settings using broadband magnetotelluric (MT) observations, with numerous studies suggesting that the spatial relationship between these 10- to 30-km-wide low-resistivity zones and vigorous crustal seismicity highlights the importance of crustal fluid migration and accumulation on the onset of earthquake rupture (e.g., Ogawa et al. 2001; Tank et al. 2005; Wannamaker et al. 2009; Yoshimura et al. 2009; Becken et al. 2011; Ichihara et al. 2014; Ogawa et al. 2014; Cai et al. 2017; Bedrosian et al. 2018). However, it is relatively unknown whether localized fluids also play a role in the growth and/or arrest of earthquake rupture, largely because the final earthquake magnitudes and their locations have not been analyzed to determine whether they are related to the subsurface fluids, or if the observed earthquake/magnitude distribution is due to a purely random process. It is therefore useful to investigate both the spatial relationship between low-resistivity zones and hypocenter locations, and the spatial relationship between low-resistivity zones, and the earthquake magnitudes and their slip distribution to elucidate the potential importance of crustal fluids in earthquake processes. The 2016 Kumamoto earthquake sequence provides the rare opportunity for such an investigation because numerous aftershocks that span broad magnitude range earthquakes, have been detected during the 14 April 2016–31 December 2019 period. Furthermore, the seismicity before the 2016 Kumamoto earthquake sequence has been captured by a densely deployed seismometer network (Matsumoto et al. 2015, 2018). In this study, we precisely relocate the detected earthquakes during the January 1993–June 2019 via the double-difference method and a detailed three-dimensional (3-D) seismic velocity structure (Shito et al. 2017). We choose the earthquakes that were detected via manual picking of the P-wave first arrivals, with 41,727 events (25,785 and 15,942 before and after the Kumamoto earthquake sequence, respectively) used in the analysis.

The Mw 7.0 mainshock, which occurred on 16 April 2016, was preceded by an Mw 6.2 foreshock on 14 April 2016. The ruptures of these two earthquakes initiated near the junction of the Hinagu and Futagawa faults at approximately 13 km below sea level (b.s.l., see Fig. 1). The observed surface displacements and surface expressions of the active faults, and the aftershock distribution, have indicated that the mainshock rupture propagated mainly ENE along the Futagawa Fault and terminated around Aso Volcano (Fujiwara et al. 2016; Shirahama et al. 2016) (Fig. 1). The mainshock rupture also propagated along the westward-dipping plane of the Hinagu Fault, whereas the foreshock rupture propagated along a vertical plane beneath the fault.

Fig. 1

Map of the study region showing (a) the broadband magnetotelluric (MT; 153 open squares) and telluric-only sites (47 open circles) used in this study are shown. The blue star and line, the red star, and two dashed boxes indicate the epicenters and fault planes of the Mw 6.2 foreshock and Mw 7.0 mainshock of the 2016 Kumamoto earthquake sequence, respectively. Black open circles are the epicenters of the crustal earthquakes (M > 2, < 30 km depth) that occurred after the Mw 6.2 foreshock of the 2016 Kumamoto earthquake sequence (14 April 2014–30 June 2019). White open circles are the epicenters of the earthquakes before the Mw 6.2 foreshock of the 2016 Kumamoto earthquake sequence (1 January 1993–14 April 2014). Large and small triangles show the current active and extinct Quaternary volcanoes, respectively. Faults are represented by green lines. HF: Hinagu Fault, FF: Futagawa Fault, MF: Midorikawa Fault, KP: Kumamoto Plain. Note that there are an additional 14 sites used in this study that lie outside of the presented study region (see Additional file 1: Fig. S1)

Magnetotelluric (MT) observations and analysis

We conducted broadband MT surveys at 90 sites around the focal region of the 2016 Kumamoto earthquake sequence in September–October 2017, March 2018, and October–December 2018. Compiled with the broadband MT data that estimated the one-dimensional (1-D) resistivity structure of this region (Asaue et al. 2012; Aizawa et al. 2017), and the two-dimensional (2-D) and 3-D resistivity structures of Aso Volcano (Takakura et al. 2000; Hata et al. 2018), we used 200 sites (magnetic and telluric; 153 sites, telluric-only; 47 sites), which were distributed at 2- to 5-km intervals around the focal area of the 2016 Kumamoto earthquake sequence (Fig. 1).

The MT response functions were calculated from the obtained electric and geomagnetic time series data using a robust estimation code (Chave and Thomson 2004). We obtained the MT response functions across a broad period range (0.005–3000 s) using remote-reference processing (Gamble et al. 1979) by using the magnetic field data recorded at least 100 km away from the sites. We used the magnetic data from nearby sites to calculate the MT response functions at the sites where only the telluric field variations were recorded.

We inverted the MT response functions using a 3-D code with an algorithm that transforms the model renewal equation from the model space to the data space (Siripunvaraporn and Egbert 2009). The full impedance tensor (four complex components: Zxx, Zxy, Zyx, and Zyy) and geomagnetic transfer functions (two complex components: Tx and Ty) across a broad period range (0.0125–3278 s period range) were used as inputs for the 3-D inversion. The different locations of electric and magnetic fields at certain observation sites were taken into account in the inversion. The horizontal mesh size was set to 1500 m in the region surrounding the focal area of the 2016 Kumamoto earthquake sequence, and logarithmically increased with increasing distance from the target region. The minimum vertical mesh size was set at 100 m from the highest point (1500 m altitude) to 0 km b.s.l., and the air was approximated by 10⁸-Ωm blocks. We used a 50, 100, 100, 200, 200 m, … vertical mesh size distribution to represent the ocean. The total mesh consisted of 78 × 78 × 59 units in the x, y, and z directions, respectively, which included seven vertical air layers as the default setting (Additional file 1: Fig. S1). The model also took topography and bathymetry into account. The calculation area was 450 × 450 × 433 km, and the seawater above the bathymetry was represented by 0.33-Ωm blocks. The initial resistivity of the land was set at 100 Ωm. We removed the noisy outlier data via visual inspection. We then assumed 10% of |SQRT(Zxy × Zyx)| as the error floor of the impedance tensor and 10% of SQRT(Tx² + Ty²) as the error floor geomagnetic transfer function. We reduced the root-mean-square (RMS) misfit of the model from 4.08 to 1.57 when we used a smoothing parameter (τ = 10 with δx = δy = δz = 0.1) in the model covariance matrix. Comparisons of the observed and modeled data fits are provided in the supplementary material (Additional file 1: Figs. S2–S5).

Another hyper-parameter λ that balance the trade-off between RMS misfit and the model smoothness is important for the structure, because inappropriate λ sometimes produce too rough or to smooth structure. We checked the consistency of the inverted resistivity structure using well-log data to avoid unrealistic resistivity values due to the inversion overfitting for the MT data. We used the well-log resistivity data measured at three locations near the Futagawa Fault (NEDO 1995; NRA 2018). The drilling depths spanned 750–1000 m and 0–1500 m depth from the surface. There is good consistency between the inverted resistivity structure beneath the drilling sites and well-log data (Additional file 1: Fig. S6). In the final iteration, we selected the structure that produced the minimum RMS between the well-log data and MT inversion. The corresponding hyper-parameter λ is 1.0.

Three-dimensional resistivity structure and interpretation

Horizontal slices of the final 3-D resistivity structure (Fig. 2; Additional file 1: S7) indicate that the focal area of the 2016 Kumamoto earthquake sequence is characterized by a mixture of low- and high-resistivity zones that possess 10-km-scale spatial lengths. The two dominant low-resistivity zones are located on the western side of the Futagawa–Hinagu fault zone (C1) and on the northeastern side of the Futagawa Fault (C2). C1 occupies a large volume at 18 km depth, and gradually branches into two moderate conductivity zones at shallower depths that are located along the southern (C1S) and northern (C1N) sections of the mainshock and foreshock hypocenters, respectively (Fig. 2; Additional file 1: S7). The presence of the four low-resistivity zones (C1, C1N, C1S, and C2) is verified via a sensitivity test (Additional file 1: Fig. S8).

Fig. 2

Horizontal slices of the 3-D resistivity structure at various depths. Hypocenters within 2 km of each slice are shown by circles (white and black circles show the hypocenters before and after the Mw 6.2 foreshock of the 2016 Kumamoto earthquake sequence, respectively). Large and small stars mark the mainshock and foreshock hypocenters, respectively. The two dashed lines across the Futagawa and Hinagu faults indicate the locations of the vertical slices shown in Fig. 3d–e. White solid line denotes the caldera rim of Aso Volcano. Other notations are the same as in Fig. 1

We interpret C1 as a zone of high-temperature magmatic fluids, even though there are no active volcanoes or geothermal areas around C1, because two Quaternary volcanoes, Kinpo and Akai (most recent eruption occurred at approximately 0.2 Ma), are located to the east and the west of C1, respectively. Approximately 2 km³ of lava was extruded from Akai Volcano around 200 ka, with widespread flow deposits present beneath Kumamoto City (Watanabe et al. 1979). The very low seismicity inside C1 over the past 25 years (Fig. 2) also supports the interpretation that C1 corresponds to a high-temperature (> 400 °C) ductile zone (Fournier 1999). However, many earthquakes occurred inside C1S and C1N; these two branches are interpreted as low-temperature zones (< 400 °C) with developed fracture networks that are filled with pore water. These fluid pathways allow the upward migration of deep magmatic volatiles. This interpretation is supported by observations of magmatic helium gas in the hot springs around C1S and C1N (Horiguchi and Matsuda 2013), and a radon-222 concentration anomaly under non-equilibrium conditions along Hinagu Fault in C1S (Koike et al. 2014). C2 is located to the north of the central cone of Aso Volcano, which is consistent with previous results (Hata et al. 2018; Matsushima et al. 2020), and is interpreted as a zone of high-temperature magmatic fluids due to the shallow extension of C2 to the active crater of Aso Volcano, which constantly emits high-temperature (> 400 °C) volcanic gases.

Figure 3b, e shows the slice of the 3-D resistivity structure along (two dashed red boxes in Fig. 2) and across (B–B′ in Fig. 2) the fault planes of the mainshock ruptures. A close-up of the C1N resistivity structure and its relationship to the located hypocenters show that there are clear upper bounds to the regional seismicity in the shallow sections of C1S and C1N, with C1N possessing a lower resistivity above these upper bounds (Fig. 3b, e). We interpret the low-resistivity zones above the upper bounds of seismicity as clay-rich, electrically conductive zones that favor a low-temperature environment (< 200 °C). Drilling along the eastern section of the Futagawa Fault showed that the low-resistivity zone (~ 10 Ωm) corresponded to a smectite layer (NEDO 1995) (Fig. 3). This clay layer may act as a cap that prohibits the upward migration of deep hydrothermal fluids, with the generation of a high-pore-pressure zone and enhanced seismicity beneath it (Fig. 3b, e). The pore-fluid pressure difference constrained by the permeability structure is essentially the same as those in widely accepted resistivity models for geothermal zones and volcanoes (e.g., Fournier 1999; Tsukamoto et al. 2018), and has also been proposed for active faults (Sibson 2007).

Fig. 3

Dipping and vertical slices of the electrical resistivity structure and seismicity along and across the faults. The dip angle of each slice is shown at the top of each figure. The hypocenters within 2 km of each slice are plotted on the resistivity profiles. a Near-vertical (dip angle = 89°) slice along Hinagu Fault, which slipped during the Mw 6.2 foreshock. Blue star and contours indicate the foreshock hypocenter and spatial slip distribution on the fault plane (0.3-m contour interval). b Northward and westward-dipping slices along the mainshock fault planes (Figs. 1, 2). The mainshock rupture propagated along the Futagawa and Hinagu faults; these two fault planes merged at the location of the vertical dashed line. Red star and contours indicate the mainshock hypocenter and spatial slip distribution on the fault plane (1.0-m interval). Dashed white lines show the inferred 200 °C and 400 °C isotherms, which correspond to the base of the altered clay (smectite) layer and brittle–ductile transition zone, respectively. c Vertical slice along the Futagawa Fault. d Vertical slice across the foreshock and mainshock hypocenters (A–A′ profile in Fig. 2). e Vertical slice across the midpoint of Futagawa Fault (B–B′ profile in Fig. 2), where the largest slip due the mainshock is estimated

Relationship between resistivity structure and earthquake rupture process

Earthquake rupture arresting process

We investigated the resistivity structure from the viewpoint of the earthquake arresting process to search for a spatial relationship between the final slip distribution and resistivity structure. The major rupture of the mainshock of the 2016 Kumamoto earthquake sequence occurred in the zone between C1 and C2 along the Futagawa Fault (Fig. 2). Figure 3a, b shows the slice of the 3-D resistivity structure along the fault planes of the foreshock and mainshock ruptures (Figs. 2, 3). The slip distribution along the fault planes was contoured from strong-motion inversion results (Asano and Iwata 2016). The mainshock rupture initiated near the upper edge of C1 (Fig. 3b–d), and mainly propagated northeastward from C1. According to the temporal slip progression model of Asano and Iwata (2016), the slip rate was not high during the first 4 s after the onset of the earthquake, with rupture occurring in the relatively high-resistivity zone; the slip rate then gradually accelerated around C1N between 4 and 10 s after the onset. The existence of C1N along a dipping plane of the Futagawa Fault (Fig. 3b), and not along a vertical plane beneath the Futagawa Fault (Fig. 3c), indicates that the rupture propagation was guided along the moderately conductive C1N (dipping plane) rather than the resistive zone (vertical plane beneath the Futagawa Fault). The rupture terminated near Aso Volcano between 10 and 20 s after the onset. The slip distribution clearly shows that the rupture arrested along the western edge of C2, beneath Aso Volcano. It is reasonable that this brittle rupture cannot propagate into C2, since C2 is interpreted as a high-temperature (> 400 °C) ductile zone. Here we note that the minor rupture of the mainshock along the Hinagu Fault terminated along the northern edge of another low-resistivity zone, C1S (Figs. 2, 3b). Furthermore, the foreshock rupture also occurred in the zone between C1S and C1N (Figs. 2, 3a), both of which are interpreted as fluid-rich brittle zones (< 400 °C). It appears that C1N acted to arrest the foreshock rupture, whereas it acted to guide the mainshock rupture. These results mean that the high-temperature (> 400 °C) magmatic fluid zone (C2) arrested the rupture of the crustal earthquakes due to its ductility, whereas the shallower low-temperature (200–400 °C) fluid-rich zone (C1N) played a key, but variable, role in crustal earthquake rupture, both allowing the Mw 7.0 mainshock rupture to propagate and arresting the Mw 6.2 foreshock rupture.

Earthquake rupture growth process

We also quantitatively investigated the resistivity structures from the viewpoint of the earthquake rupture growth process to search for a spatial relationship between the electrical resistivity structures, and the earthquake rupture initiation locations (hypocenters) and their final magnitudes. We defined the < 30-Ωm structures as low-resistivity zones; 30 Ωm corresponds to 0.5 volume% good-connectivity pore networks (Hashin and Shtrikman 1962) that are filled with 0.1-Ωm brine, which are typical structures in the mid-to-lower crust (Nesbitt 1993; Sakuma and Ichiki 2016). We then calculated the minimum distance from the hypocenters to the nearest low-resistivity zones using the hypocenters within 40 km of the mainshock hypocenter. We defined the hypocenters that were surrounded by seven low-resistivity blocks (middle, upper, lower, north, south, east, and west) as earthquakes that initiated within a low-resistivity zone. The other hypocenters were defined as the earthquakes that initiated outside of a low-resistivity zone. We used the HYPOMH program (Hirata and Matsuura 1987) to calculate the earthquake magnitudes.

Figure 4 shows the relationship between earthquake magnitude and distance to a low-resistivity zone. Because the hypocenter is the point where an earthquake rupture initiates, the notable observation in Fig. 4 is that all of the large (M > 5.3) earthquakes occurred within 5 km of the low-resistivity zones; however, their ruptures did not initiate within these zones. The average observation site spacing (2–5 km) and spatial smoothness constraints in the inversion indicate that the ruptures of the large earthquakes must initiate at the outer edge of the low-resistivity zones. In contrast, the earthquakes that initiate either distal to or within the low-resistivity zones do not grow into large earthquakes. Such a relationship is unexpected if the rupture growth is a random non-deterministic process. The averaged distances between the hypocenters and low-resistivity zones within different two-magnitude ranges also show that the large earthquakes tend to occur along the outer edge of the low-resistivity zones (Fig. 4). This tendency did not change before or after the 2016 Kumamoto earthquake sequence. The results suggest that the large earthquakes tended to selectively initiate ruptures along the outer edge of the low-resistivity zones, whereas the smaller earthquake ruptures initiated everywhere.

Fig. 4

Relationship between earthquake magnitude and distance to a low-resistivity zone (LRZ). Crosses and open circles represent earthquakes before and after the 2016 Kumamoto earthquake sequence, respectively (red: within LRZ, black: outside of LRZ). Here, we defined the low-resistivity zones as < 30 Ωm structures. Blue and red stars indicate the Mw 6.2 foreshock and Mw 7.0 mainshock, respectively. We used only the hypocenters within a 40-km radius of the mainshock. Yellow circles indicate the average distance values for three magnitude ranges (M1–M3, M3–M5, and M5–M7), which were calculated for the hypocenters that were located outside of LRZs. Histograms of the earthquake magnitudes (upper) and earthquake distances to the nearest LRZ (right) are also shown. The red area in the histogram indicates the hypocenters within the LRZ

Discussion

Figure 5a shows the schematic relationship between the low-resistivity zones and earthquake magnitude. A localized stress accumulation around the mechanically weak low-resistivity zones, and/or fluid supply from the low-resistivity zone have been proposed as potential crustal earthquake initiation mechanisms around low-resistivity zones (e.g., Ogawa et al. 2001; Ichihara et al. 2014; Aizawa et al. 2017; Cai et al. 2017). However, these proposed mechanisms might not account for the relationship between the low-resistivity zones and the final earthquake magnitudes. The large earthquake did not occur distal to the low-resistivity zones, even though the stress field was estimated high there (Matsumoto et al. 2016). Therefore, we argue that the pore-fluid pressure—which is considered to be high inside and near low-resistivity zones, such as deep magmatic fluid zones (Fournier 1999; Lee et al. 2020) or a fracture zone that transports magmatic volatiles (Aizawa et al. 2016; Lee et al. 2020)—plays an important role in the evolution of crustal earthquake rupture. We hypothesize that the pre-failure pressure/temperature (PT) gradient (spatial difference) of the pore fluids contributes to the propagation and arrest of earthquake rupture. This proposed mechanism is based on Griffith’s criterion, which states that real rocks must contain flaws, such as cracks and voids. Note that the observed resistivity values (1–20,000 Ωm range; Figs. 2, 3) suggest that fluid-filled cracks can exist anywhere in the crust because dry granite and gabbro possess electrical resistivities of > 100,000 Ωm at < 400 °C (Kariya and Shankland 1983; Fuji-ta et al. 2004). When the rupture initiates at one of the cracks, it can propagate along other cracks, and these rupturing cracks can coalesce. Such crack coalescence can enhance both the rupture magnitude and propagation speed, culminating in macroscopic rupture (Ashby and Hallam 1986; Kame and Yamashita 1997). Tensile cracking may be generated when the rupture (slip) initiates at one of the cracks (Fig. 5b) (Ashby and Hallam 1986; Liu et al. 2017), which is consistent with the novel view of crustal earthquakes, whereby rupture initiation not only induces slip, but also accompanies the opening of tensile cracks due to the presence of localized high-pressure fluids (Hayashida et al. 2020). If the rupture initiates at a crack near the edge of a low-resistivity zone and crack coalescence occurs outward from the low-resistivity zones, then the pre-failure PT gradient in the pore fluids may promote pore-fluid migration into the crack, which may enhance crack propagation and widening at the tip of the coalesced crack (Fig. 5b). This process likely occurs successively at various spatial scales, and subsequently tends to advance the rupture front; it also yields a high probability of resulting in a large earthquake. Conversely, if the rupture initiated either within or distal to the low-resistivity zone, then the small PT gradient in the pore fluids may be less likely to promote rupture growth, resulting in a small earthquake. The relative pore-pressure difference would be maintained based on the phase diagram of water (Chen et al. 2017), even if water vaporization occurs along the edge of the newly opened crack.

Fig. 5

Schematic model for the relationship between the electrical resistivity structure and local seismicity. a Interpretation of the resistivity structure in terms of fluids and isotherms. Red areas indicate relatively low-resistivity zones. Red and blue arrows represent hot- and cold-fluid movements, respectively. Large and small crosses show the hypocenters of large and small earthquakes, respectively; large earthquakes tend to initiate along the outer edge of the low-resistivity zones. b Schematic model for rupture initiation and its propagation outward from the outer edge of a low-resistivity zone. From left to right, the three figures represent the temporal evolution of rupture propagation. The model is based on the hypothesis that the pre-failure pressure/temperature (PT) gradient (spatial difference) of the pore fluids contributes to crack propagation and coalescence. One possible mechanism for crack coalescence is shown here. The solid arrows indicate the slip directions and the opening of the tensile crack. The spatial difference in the pre-failure pore-fluid pressure tends to both promote fluid movement (red arrow) and enhance outward rupture propagation (open arrow) from the low-resistivity zones. Note that the direction of rupture growth is constrained by both the orientation of the crack that initially slipped and the regional stress field

The PT gradient may also contribute to the arrest of the rupture. The presence of high pore-fluid pressure inherently acts to slow rupture growth due to dilatant hardening (French and Zhu 2017). Furthermore, rupture propagation from a zone of low PT conditions to a zone of high PT conditions is considered to be reduced by the opposite mechanism, as shown in Fig. 5b. This means that the cracks to the high PT fluids act as barriers to rupture propagation. However, the high PT condition may act to advance the rupture front away from the low-resistivity zones once the rupture is well within these zones. C1N arrested the rupture of the Mw 6.2 foreshock of the 2016 Kumamoto earthquake sequence, but this zone was highly damaged by the foreshock and subsequent aftershocks, which reduced the pore-fluid pressure around C1N and potentially placed it in a vulnerable state for the next large rupture propagation. Furthermore, the spatial scale and slip amount at the mainshock rupture front (10 km and 2 m, respectively) were far larger than those at the foreshock rupture front (5 km and 0.5 m) when the rupture reached C1N; therefore, the mainshock rupture might be able to propagate deep into C1N, resulting in the promotion of the mainshock.

We found changes in the seismicity before and after the 2016 Kumamoto earthquake sequence. Large (M > 5) earthquakes only occurred at the top of C1 and C2 before the 2016 Kumamoto earthquake sequence, whereas they also occurred along the outer edges of the shallow sections of C1S and C1N after the earthquake sequence. The observed changes in the hypocenter locations of the large earthquakes may indicate that pore-pressure diffusion occurred via the mechanism proposed in Fig. 5b. Another change in the seismicity before and after the 2016 Kumamoto earthquake sequence was observed at the eastern edge of C1. No earthquakes occurred at > 15 km depth before the 2016 Kumamoto earthquake sequence, whereas the Mw 6.2 foreshock marked a shift to the occurrence of deeper earthquakes (Mitsuoka et al. 2020; Shito et al. 2020) (18-km-depth slice in Fig. 2). The deep aftershocks are actually located along the eastern edge of C1 (Fig. 2), and magmatic fluids likely moved horizontally along the rupture formed by these deep aftershocks. This may indicate that the lateral expansion of C1 in the lower crust was stimulated by the rupture of a large crustal earthquake. Deep low-resistivity zones, similar to C1, have been found beneath active faults (e.g., Ogawa et al. 2001; Yoshimura et al. 2009; Becken et al. 2011; Cai et al. 2017). These results suggest that there may be a positive feedback, whereby deep fluid-rich zones generate large crustal earthquakes and induce the evolution of high-temperature fluid-rich zones in the mid-to-lower crust.

The low-resistivity zones and their observed relationship to local seismicity suggest that the mid-crustal fluid distribution controls the initiation, growth, and arrest of crustal earthquake rupture. These results demonstrate that the 3-D imaging of electrical low-resistivity zones using spatially dense MT observations provides valuable information for assessing the locations and maximum magnitudes of future earthquakes, particularly since some large (M 6–7) intraplate earthquakes have occurred in zones with no significant surface fault traces (Semmane et al. 2005). It has recently been suggested that the strong frictional coupling (plate locking) of megathrust interplate earthquakes occurs in relatively high-resistivity zones, which are interpreted to be fluid- and sediment-sparse zones (Wannamaker et al. 2014; Heise et al. 2017). Therefore, evaluations of the potential for large earthquakes and their locations in a given region may be refined when we take into account the locations of low-resistivity fluid-rich zones, regardless of earthquake type.

Conclusions

The relationship between resistivity structure and earthquake rupture process suggested that the distribution of mid-crustal fluids contributes to the initiation, growth, and arrest of crustal earthquakes. In addition to the pre-failure PT gradient, heterogeneities in the stress field, fracture strength, and crack density may also contribute to the initiation and arrest of earthquakes (e.g., Umeda et al. 1996; Kame and Yamashita 1999; Matsumoto et al. 2018). Furthermore, all of the parameters depend upon each other, and can undergo temporal changes during rupture growth. Modeling approaches that explore the evolution of earthquake rupture over a range of spatial and time scales are necessary to elucidate the key factors that control earthquake rupture growth and arrest. Constraints on the locations of fluid-rich crustal zones may provide essential clues to address this problem.

Availability of the data

The broadband MT data from Aso Volcano are available at https://www.gsj.jp/researches/openfile/openfile2018/openfile0654.html, although the data ownership belongs to the AIST in any case. The other MT data are available from the author upon request. The digital data for the 3-D electrical resistivity model is available on the SEVO website (http://www.sevo.kyushu-u.ac.jp/open_file/RHO_Kumamoto_Aizawa_et_al_2020.txt).