Electrical Conductive Fluidized Zones and Their In uence on the Earthquake Nucleation, Growth, and Arrest Processes of the 2016 Kumamoto Earthquake Sequence, Kyushu Island, Japan


 Crustal earthquake ruptures tend to nucleate near fluidized zones. However, it is relatively unknown whether fluidized 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 broad-band magnetotelluric data, and discuss its quantitative relationship to earthquake nucleation, growth, and arrest processes. The result shows that the earthquake hypocenters are all located within 10 km from low-resistivity fluidized zones < 30Ωm. The ruptures that nucleated along the outer edge of the low-resistivity fluidized zones tended to become large earthquakes, whereas those that initiated either distal to or within the fluidized zones did not. The ruptures were arrested by high-temperature (>400°C) fluidized zones, whereas shallower low-temperature (200°C–400°C) fluidized zones either promoted or arrested the ruptures. These results suggest that the distribution of mid-crustal fluids contributes to the nucleation, growth, and arrest of crustal earthquakes.


Introduction
Electrically conductive uidized 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-30-km-wide low-resistivity zones and vigorous crustal seismicity highlights the importance of crustal uid migration and accumulation on the onset of earthquake rupture relatively unknown whether localized uids also play a role in the growth and/or arrest of earthquake rupture, largely because the nal earthquake magnitudes and their locations have not been analyzed to determine whether they are related to the subsurface uids, 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 uids 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;Matsumoto et al., 2018). In this study, we precisely relocate the detected earthquakes during the January 1993-June 2019 via the doubledifference 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 rst arrivals, with 41,727 events (25, 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 eld 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 eld 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 to 3278s period range) were used as inputs for the 3-D inversion. The different locations of electric and magnetic elds 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 8 -Ωm blocks. The total mesh consisted of 78 × 78 × 56 units in the x, y, and z directions, respectively, which included seven vertical air layers as the default setting (Fig. S1). The model also took topography and bathymetry into account. The calculation area was 450 × 450 × 275 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 for the impedance tensor and 10% of SQRT(Tx 2 + Ty 2 ) as the error for the geomagnetic transfer functions. We reduced the root-mean-square (RMS) mis t of the model from 4.08 to 1.56 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 ts are provided in the supplementary material (Figs S2-S4).
Another hyper parameter λ that balance the trade-off between RMS mis t 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 over tting 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 (Fig. S5). The minimum RMS between the well-log data and MT inversion is obtained from the nal model results. The corresponding hyper-parameter λ is 0.33.

Three-dimensional Resistivity Structure And Interpretation
Horizontal slices of the nal 3-D resistivity structure ( Fig. 2) 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. The presence of the four low-resistivity zones (C1, C1N, C1S, and C2) are veri ed via a sensitivity test (Fig. S6).
We interpret C1 as a zone of high-temperature magmatic uids, 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 3 of lava was extruded from Akai Volcano around 200 ka, with widespread ow 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 lled with pore water. These uid 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  Matsushima et al., 2020), and is interpreted as a zone of high-temperature magmatic uids due to the shallow extension of C2 to the active crater of Aso Volcano, which constantly emits high-temperature (> 400 °C) volcanic gases.
The largest surface displacements along the Futagawa Fault (> 2 m) were observed above C1N (Shirahama et al., 2016); therefore, we interpret C1N as the damage zone formed by repeated faulting and fracturing over time. Figures 3b and d show the slice of the 3-D resistivity structure along and across the fault planes of the mainshock ruptures. A close-up of the C1N resistivity structure and its relationship to the located hypocenters shows that there are clear upper bounds to the regional seismicity in the shallow sections of C1S and C1N, with C1N possessing a higher conductivity above these upper bounds ( Fig. 3b  and d). 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 uids, with the generation of a high-pore-pressure zone and enhanced seismicity beneath it (Figs. 3b and d). These triple-layered structures are constrained by the 200 °C and 400 °C isotherms, and the corresponding lithologies (C1 − C1N − C1Ntop, and C1 − C1S − C1Stop) are essentially the same as those in widely accepted resistivity models for geothermal zones and volcanoes (Fournier, 1999;Tsukamoto et al., 2018), and have also been proposed for active faults (Sibson, 2007 We investigated the resistivity structure from the viewpoint of the earthquake arresting process to search for a spatial relationship between the nal 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). Figures 3a and b show the slice of the 3-D resistivity structure along the fault planes of the foreshock and mainshock ruptures (Figs 2 and 3). The slip distribution along the fault planes was contoured from strong-motion inversion results (Asano and Iwata, 2016). The mainshock rupture nucleated near the upper edge of C1 ( Fig. 3a and c), and mainly propagated northeastward from C1. The temporal slip progression indicated that the slip rate was not high during the rst 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 fault (Fig. S7), 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 and 3b). Furthermore, the foreshock rupture also occurred in the zone between C1S and C1N (Figs 2 and 3a), both of which are interpreted as uidized 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 uid zone (C2) arrested the rupture of the crustal earthquakes due to its ductility, whereas the shallower uidized zone (C1N) played a key, but variable, role in crustal earthquake rupture, 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 nucleation locations (hypocenters) and their nal magnitudes. We de ned 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 lled with 0.1-Ωm brine, which are typical structures in the mid-tolower crust (Nesbitt, 1993). 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 de ned the hypocenters that were surrounded by seven low-resistivity blocks (middle, upper, lower, north, south, east, and west) as earthquakes that nucleated within a low-resistivity zone. The other hypocenters were de ned as the earthquakes that initiated outside of a low-resistivity zone. We de ned the hypocenters located within 5 km of a low-resistivity zone as the earthquakes that nucleated along the outer edges of a low-resistivity zone. We used the HYPOMH program (Hirata and Matsuura, 1987) to calculate the earthquake magnitudes. Fig. 4 shows the relationship between earthquake magnitude and distance to a low-resistivity zone. All of the hypocenters are located within 10 km of the low-resistivity zones due to the prevalence of lowresistivity zones in the study area. The notable feature is that the large (M > 5.3) earthquakes all nucleated within 5 km of the low-resistivity zones, but their ruptures did not nucleate 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 nucleate at the outer edge of the low-resistivity zones. In contrast, the earthquakes that nucleate 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 nondeterministic 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 nucleate ruptures along the outer edge of the low-resistivity zones, whereas the smaller earthquakes nucleated everywhere. 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 uid supply from the low-resistivity zone have been proposed as potential crustal earthquake nucleation 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 lowresistivity zones and the nal earthquake magnitudes. The large earthquake did not occur distal to the low-resistivity zones, even though the stress eld was estimated high there (Matsumoto et al., 2016). Therefore, we argue that the pore-uid pressure-which is considered to be high inside and near lowresistivity zones, such as deep magmatic uid 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 uids contributes to the propagation and arrest of earthquake rupture. This proposed mechanism is based on Gri th's criterion, which states that real rocks must contain aws, such as cracks and voids. Note that the observed resistivity values (1-20,000 Ωm range; Figs. 2 and 3) suggest that uid-lled 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 nucleates 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) nucleates 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 nucleation not only induces slip but also accompanies the opening of tensile cracks due to the presence of localized high-pressure uids (Hayashida et al., 2020). If the rupture nucleates at a crack near the edge of a low-resistivity zone and crack coalescence occurs outward from the lowresistivity zones, then the pre-failure PT gradient in the pore uids may promote pore-uid 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 nucleated either within or distal to the low-resistivity zone, then the small PT gradient in the pore uids 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 , even if water vaporization occurs along the edge of the newly opened crack.

Discussion
The PT gradient may also contribute to the arrest of the rupture. The presence of high pore-uid 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 with high PT uids 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 uid 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, resulted in the promotion of the mainshock.
We found changes in the seismicity before and after the 2016 Kumamoto earthquake sequence. The number of earthquakes that nucleated within the low-resistivity zones before the earthquake sequence was approximately 16% of the total number of the earthquakes; this number decreased to 8% after the earthquake sequence. This observation may indicate that pore pressure in the low-resistivity zones decreased during the 2016 Kumamoto earthquake sequence due to pore-uid pressure diffusion via the mechanism proposed in Fig. 5b. 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 are also consistent with the speculation that the high poreuid pressure diffused from the low-resistivity zones to the surrounding high-resistivity zones.
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 (18 km depth slice in Fig. 2) marked a shift to the nucleation of deeper earthquakes. The deep aftershocks are actually located along the eastern edge of C1 (Fig. 2), and magmatic uids 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 uidized zones generate large crustal earthquakes and induce the evolution of high-temperature uidized zones in the mid-to-lower crust.
The low-resistivity zones and their observed relationship to local seismicity suggest that the mid-crustal uid distribution controls the nucleation, 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 signi cant 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 highresistivity zones, which are interpreted to be uid-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 re ned when we take into account the locations of low-resistivity uidized zones, regardless of earthquake type.

Conclusions
The relationship between resistivity structure and earthquake rupture process suggested that the distribution of mid-crustal uids contributes to the nucleation, growth, and arrest of crustal earthquakes. In addition to the pre-failure PT gradient, heterogeneities in the stress eld, fracture strength, and crack density may also contribute to the nucleation and arrest of earthquakes (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 uidized crustal zones may provide essential clues to address this problem. lie outside of the presented study region (see Fig. S1).

Figure 2
Horizontal slices of the 3-D resistivity structure at various depths. Hypocenters within 2 km of each slice are shown by circles (black and white circles show the hypocenters before and after the Mw 6.   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 de ned 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 average 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.

Figure 5
Schematic model for the relationship between the electrical resistivity structure and local seismicity. (a) Interpretation of the resistivity structure in terms of uids and isotherms. Red areas indicate relatively lowresistivity zones. Red and blue arrows represent hot-and cold-uid movements, respectively. Large and small crosses show the hypocenters of large and small earthquakes, respectively; large earthquakes tend to nucleate along the outer edge of the low-resistivity zones. (b) Schematic model for rupture nucleation and its propagation along the outer edge of a low-resistivity zone. The model is based on crack coalescence, with the gure showing one possible mechanism for crack coalescence. The spatial difference in pre-failure pore-uid pressure tends to both promote uid movement (red arrow) and enhance outward rupture propagation 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 eld.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.