Three‐Dimensional Electrical Resistivity Structure Beneath a Strain Concentration Area in the Back‐Arc Side of the Northeastern Japan Arc

Intraplate earthquakes occur more frequently in the Japanese islands than in other regions. Large intraplate earthquakes in the island arc preferentially occur in strain concentration zones detected by geological and geodetic studies. Crustal heterogeneity plays a crucial role in generating large intraplate earthquakes and strain concentrations. Thus, we elucidated the crustal heterogeneities beneath a strain concentration area on the back‐arc side of the northeastern Japan Arc based on electrical resistivity, which is sensitive to weak zones in the crust. By deploying wideband magnetotelluric surveys, we revealed the three‐dimensional electrical resistivity structure in the crust, suggesting the coexistence of two different types of strain‐concentration mechanisms in the strain‐concentration area. The shallow conductive layers and lower‐crustal conductors appear to act as low‐elastic‐modulus and low‐viscosity areas, respectively, and are responsible for the strain concentration. We found shallow and lower‐crustal conductors in the strain concentration zone revealed by geological studies. Those conductive areas are considered to act as mechanically weak zones and cause stress loading on the brittle pars of pre‐existing faults, resulting in large intraplate earthquakes. The resultant earthquakes presumably contribute to strain accumulation on a geological timescale. In addition, a spatial correlation between the epicenters of large intraplate earthquakes and edges of lower‐crustal conductors implies the contribution of the fluids in the lower crust to the generation of large earthquakes. We also identified vertical conductors ranging from the lower crust to Quaternary volcanoes, which may indicate trans‐crustal magmatic systems under these volcanoes.


Introduction
Shallow intraplate earthquakes are the most destructive type of natural disaster.They can cause devastating damage because of their shallow hypocenters and proximity to inhabited areas.Intraplate earthquakes in the Japanese islands are more frequent than in other regions (e.g., Kato, 2014).Beneath northeastern Japan arc, the Pacific Plate subducts under the North American (or Okhotsk) Plate at a convergence rate of approximately 8 cm/ yr in a nearly WNW direction (e.g., Loveless & Meade, 2010).Many large intraplate earthquakes (M > 7) have occurred in northeastern Japan because of the stress field caused by plate convergence (Figure 1a).
The distribution of the epicenters of historical large earthquakes (M > 7) is not uniform in northeastern Japan, but is preferentially located in the strain concentration zones detected by geological studies (Okamura, 2002;Okamura et al., 2007) (yellow areas in Figure 1).Okamura (2002) and Okamura et al. (2007) defined the geological strain concentration zone, hereafter denoted as SCZGS (Strain Concentration Zone detected by Geological Studies), based on the distribution of faults and folds.The SCZGS is a strain concentration zone on a geological time scale (in the past approximately 3 million years), and inelastic deformation (such as displacement by large intraplate earthquakes) contributes to its strain concentration.In a major part of the SCZGS, geodetic measurements using the global positioning system (GPS) have also detected the strain concentration on the time scale of 10 years (Miura et al., 2004;Sagiya, 2002Sagiya, , 2004)).Miura et al. (2004) revealed a strain concentration zone along the Ou Backbone Range (OBR), which is hereafter denoted as SCZOBR.Sagiya (2002) suggested another strain concentration zone along the eastern margin of the Sea of Japan, denoted as SCZEJS.The southernmost part of the SCZEJS is shared with the Niigata-Kobe tectonic zone (NKTZ) (Sagiya, 2004), which is another strain concentration zone continuing to southwest Japan.
Seismic tomography images show low-velocity and low-Q (high attenuation) anomalies under the OBR (Hasegawa & Nakajima, 2004;Hasegawa et al., 2005;Nakajima et al., 2013).Based on those findings, Hasegawa et al. (2009) attributed the high strain rate of the SCZOBR to the locally weak crust due to the presence of the slabderived fluid.Muto (2011) also concluded that the lower crust under the OBR has low viscosity by combining a melt fraction inferred based on rock velocity measurements of xenoliths (Nishimoto et al., 2008) and the flow law for the partially molten plagioclase (Dimanov et al., 1998).However, the strain-concentration mechanisms of SCZEJS and SCZGS, which cover wider areas than SCZOBR, have rarely been investigated.Based on dense GPS measurement data, Ohzono et al. (2012) suggested that the strain concentration in the northern-most part of the NKTZ is due to the low elastic modulus of the thick sedimentation layer.Because the northern-most part of the NKTZ is shared with the southern-most part of the SCZEJS (Figure 1a), their interpretation is applicable to the southern-most part of the SCZEJS.However, their model, suggesting different strain concentration mechanisms between SCZOBR and SCZEJS, was not examined in subsequent studies.
Crustal heterogeneity is believed to play a crucial role in the generation of intraplate earthquakes and strain concentrations (Hasegawa et al., 2009;Iio et al., 2004;Kato, 2014;Ogawa et al., 2001).Thus, to identify the mechanism of strain concentration and the generation process of large intraplate earthquakes, it is necessary to reveal the subsurface heterogeneity from the near-surface to the deep crust.One of the most effective approaches for imaging subsurface heterogeneity is probing the subsurface electrical resistivity structure using electromagnetic (EM) induction surveys.Electrical resistivity is sensitive to temperature and the abundance of liquid phases (such as melts and aqueous fluids) in rocks.It has been recognized that dehydration on subducting plates releases large amounts of fluids in subduction zones (e.g., Iwamori, 1998;van Keken et al., 2011).EM induction surveys have been conducted to estimate the fluid distribution in a number of subduction zones, including Cascadia (e.g., Bedrosian et al., 2018;Wannamaker et al., 2014), New Zealand (e.g., Heise et al., 2024;Wannamaker et al., 2009), Costa Rica (e.g., Worzewski et al., 2011), Central Andes (e.g., Araya Vargas et al., 2019), Central Chile (e.g., Cordell et al., 2019) and also Japan (as mentioned below).Those studies have found mid-and lower-crustal fluid-rich areas that are associated with earthquakes and volcanism.
The EM induction survey allows us to clarify the distribution of sedimentary layers (Ichihara et al., 2013;Uyeshima et al., 2005) and the presence of fluid-rich areas in the mid-and lower crust (Ogawa & Honkura, 2004;Usui et al., 2021).Both are expected to act as mechanically weak regions.Because the electrical resistivity can be interpreted as a proxy for mechanical strength in the lithosphere (e.g., Liu & Hasterok, 2016;Pommier et al., 2013), it has been used to constrain lithospheric rheological properties in various areas worldwide (e.g., Comeau et al., 2016;Dong et al., 2020;Liu & Hasterok, 2016;Murphy et al., 2019;Sheng et al., 2022).For example, Liu and Hasterok (2016) converted the resistivity structure to the viscosity structure by an empirical relationship and improved the prediction of surface topography and lithosphere deformation in the Western United States.In a similar way, in the central part of the Tibetan Plateau, Dong et al. (2020) revealed contrasting lateral viscosity differences in the lower crust that seem responsible for complex surface deformation patterns.
To investigate the mechanisms of strain concentration and the generation of intraplate earthquakes, the present work estimates the 3-D crustal resistivity structure beneath a strain concentration area on the back-arc side of northeastern Japan (Figure 1c), where the SCZGS covers a wide area, and its eastern and western sides are shared with SCZOBR and SCZEJS, respectively.A number of Quaternary volcanoes (e.g., Chokaisan, Hijiori, and Gassan) are situated in the study area (Figure 1c).Because volcanic activities involve subsurface fluids, the presence of these volcanoes implies the important role of fluids in the tectonics of the study area.In parts of the study area, two-dimensional (2-D) resistivity structures along single survey lines were estimated (Asamori et al., 2011;Ichihara et al., 2011).However, the crustal resistivity structure in the area is not strictly 2-D, as illustrated in subsequent sections.Therefore, their 2-D inversion could yield false structures.Elucidating the regional-scale 3-D crustal structure is essential for investigating the relationship between the resistivity structure, strain concentration, and the occurrence of earthquakes.
In the following sections, we provide information on our EM induction survey and obtained response functions.Next, we show the 3-D resistivity structure estimated by inversion.Based on the 3-D electrical resistivity structure, we discuss the contribution of crustal heterogeneity to the strain concentration and the generation of large intraplate earthquakes.In addition, we discuss crustal conductors that are possibly associated with volcanic activity.

Magnetotelluric Measurement and Data Analysis
We conducted wideband magnetotelluric (MT) measurements in the study area (Figure 2).The wideband MT survey is the most appropriate method for estimating the near-surface structure to that in the deep crust because it can simultaneously measure wideband electromagnetic signals (e.g., from several hundred Hz to over 10,000 s).We deployed 91 stations in the study area (Figure 2) from 2006 to 2011.Among the 91 stations, 55 were measured using the MTU systems (Phoenix Geophysics).The sampling frequencies of the MTU systems were 2,400, 150, and 15 Hz.We measured the 2,400 and 150 Hz data intermittently, whereas the 15 Hz data were seamlessly measured.The remaining 36 stations were measured using the ADU systems (Metronix).The sampling frequencies of the ADU systems were 1,024 and 32 Hz.We measured the 1,024 Hz data only at midnight, whereas the 32 Hz data were measured seamlessly.The observation stations along line FF' (Figure 2) were deployed by the Japan Atomic Energy Agency (Asamori et al., 2011).The other stations were deployed in the "Multidisciplinary Research Project for High Strain Rate Zones" (2008)(2009)(2010)(2011)(2012).The measured data, except for those of lines CC' and FF', were analyzed for the first time in this study.The observation periods differed among stations, ranging from several days to a month.At the 12 stations denoted by inverted triangles in Figure 2, we deployed the equipment only to measure the electric field, or the magnetic field was not properly measured.For these stations, the magnetic field of a neighboring station was used to estimate the MT response functions.Such inter-station technique has been used in numerous MT surveys (e.g., Comeau et al., 2020;Käufl et al., 2020;Muñoz & Ritter, 2013).Different locations of the telluric and magnetic fields were considered for subsequent MT inversions.Detailed information concerning each observation station is listed in Supporting Information S2 (Tables S1-S3).The Sawauchi (Nittetsu Mining Consultants Co., Ltd.) and Esashi stations (Geospatial Information Authority of Japan) (Figure 1a) were used as the reference stations for the remote reference method (Gamble et al., 1979).MTU systems were used at the remote reference stations.
Using the observed EM field time series data, we estimated the impedance tensor, vertical magnetic transfer function (VMTF), and interstation horizontal magnetic transfer function (HMTF).We used the two-stage bounded influence method (Chave & Thomson, 2004) to estimate the response functions and the subset deletion jackknife method (Eisel & Egbert, 2001) to estimate the error of the response functions.Prior to data analysis, we downsampled the 15 Hz data of the MTU stations and the 32 Hz data of the ADU stations to 1 Hz.The primary reason for the down-sampling was to use the Sawauchi or Esashi stations (Figure 1a) as the reference stations of the remote reference method (Gamble et al., 1979) for long-period responses (≥30 s) of the ADU stations.To estimate the short-period response functions (<30 s) of the ADU stations, the magnetic field of a neighboring station was used in the remote reference method (Gamble et al., 1979).For the MTU stations, we used either the Sawauchi or Esashi stations as remote reference stations for all periods.Information regarding the remote reference stations for the respective stations is summarized in Supporting Information S2 (Tables S1-S3).Tables S1-S3 in Supporting Information S2 present the reference stations for the HMTF.We selected the cleanest station among those that shared an observation period as the reference station.Because the sampling frequencies differed between the MTU and ADU stations, except for the 1 Hz data, the periods in which the response functions were estimated were not the same.The periods for the MTU stations were 12 periods from 0.03048 (32.81 Hz) to 2731 s, whereas those for the ADU stations were 12 periods from 0.03125 (32.00 Hz) to 2731 s.The five longest periods, 42.67-2731 s, were common between the MTU and ADU stations because we used the down-sampled  8).The observation stations along line FF' were deployed by Japan Atomic Energy Agency (Asamori et al., 2011), whereas the other stations were deployed by the "Multidisciplinary Research Project for High Strain Rate Zones" (2008)(2009)(2010)(2011).Active faults (Nakata & Imaizumi, 2002)  Figure 3 shows the phase tensor (PT) ellipses (Caldwell et al., 2004) and induction vectors (Parkinson convention).Because the periods at which the response functions were estimated were not the same between the MTU and the ADU stations, except for the long-period responses (≥42.67 s), periods slightly differ among the  et al., 2004), respectively.Black arrows in (c) denote the real part of the Parkinson vectors.Because the periods at which the response functions were estimated were not the same between the MTU and the ADU stations, except for the long-period responses (≥42.67 s), periods slightly differ among the two measurement systems in the two left columns.Red triangles indicate Quaternary volcanos.Shaded areas in the leftmost columns are the schematic illustration of sedimentary plains and basins (Geological Survey of Japan, 2022).The response functions that appeared to be heavily affected by cultural noises were not plotted.two measurement systems.At periods of 1.07 (MTU) and 1.00 s (ADU), the oblateness of the PT ellipses are relatively small, and both Φ min and Φ max are higher than 45°at the northwest area near the coast, where the sedimentary Shonai plain is located.During these periods, the induction vectors pointed toward the sedimentary plains and basins (Figure 3c).At periods longer than 10 s, particularly at 128 s, the PT ellipses were flattened in the northeastern part of the study area.In the period range, Φ min (<45°) is significantly lower than Φ max (>45°).Because the major axes of the PT ellipses, corresponding to Φ max , are aligned perpendicular to the island arc and coastline, in the northeastern part, the phase between the along-arc electric field and across-arc magnetic field is relatively large.This feature was confirmed in the sounding curves of the phases (SKT150 and SKT160 in Figure S11 in Supporting Information S1).The Sea of Japan and Syonai Plain contribute to the flattened PT ellipses because the strong spatial difference in the surface electrical resistivity can cause splitting of the impedance tensor (Fischer & Weaver, 1986).However, although the coastal effect makes the induction arrow point toward the sea, the induction arrows in the eastern half of the survey area are aligned in the along-arc direction at periods longer than 100 s.Thus, it appears that the inland sedimentary basins and possibly deep crustal resistivity anomalies also contribute to the flattened PT ellipses.

MT Inversion
Using the estimated response functions described above, we estimated the resistivity structure beneath the survey area using 3-D MT inversion.In the modeling, we used the hexahedral mesh version of the FEMTIC code (Usui, 2015(Usui, , 2021;;Usui et al., 2017), in which the finite-element method was used to compute the EM field and response functions.A more detailed description of the scheme of the FEMTIC code is provided in Appendix A. Figure 4 shows the computational mesh used in this study.We set the positive x direction to N20°E, which is nearly parallel to the coastline and the along-arc direction.Using the deformed hexahedral elements, we incorporated the topography and land-sea distribution into the computational mesh.As shown in Figure 4a, we locally refined the elements around each observation station to achieve a higher resolution of the topography and resistivity distribution in the vicinity of the observation points.The minimum horizontal edge length was approximately 300 m.We used the 10-m grid elevation of the Geospatial Information Authority of Japan and the JODC-Expert Grid data for Geography-500 m (J-EGG500) (Japan Oceanographic Data Center) as altitude and bathymetry data, respectively, around the survey area.For the other areas, we used ETOPO1 (Amante & Eakins, 2009).Sea resistivity was kept at 0.3 Ωm in the inversion.The total number of elements used in the mesh is 650,220, and the number of parameter cells is 464,643.The parameter cell denotes a group of elements that share the same electrical resistivity (Usui, 2015).The initial model was the 10-Ωm half-space except for the sea layer.Even when we used the 100-Ωm half-space as the initial model, no significant difference in the resultant subsurface structure was noticeable (Figures S2 and S3 in Supporting Information S1).
For the inversion, we used the impedance tensor, VMTF, and HMTF.Campanyà et al. (2016) confirmed the effectiveness of the combined inversion of the impedance tensor, VMTF, and HMTF in constraining the resistivity structure.Prior to inversion, the apparent outliers in the sounding curves were removed from the data set.Because the VMTF and HMTF are less sensitive to the absolute values of the subsurface resistivity than the impedance tensor, we first inverted the impedance tensor only.We then started the combined inversion of all the response functions from the resulting model of the former inversion.In the former impedance tensor inversion, we first estimated the resistivity structure only as model parameters until convergence and then restarted the inversion by adding the galvanic distortion parameters to the unknowns to prevent overestimation of the galvanic distortion parameters.The FEMTIC code seeks a solution that minimizes the objective function consisting of the data-misfit term, model-roughness term, and the term of the norm of galvanic distortion (Usui, 2015;Usui et al., 2017; see Appendix A).In the modeling, as in Usui et al. (2017) and Yoshimura et al. (2018), we selected the preferred model using the L-curve criteria from the results obtained by several different α values (1.0, 1.8, 3.0, 5.6, 10, 18, 30, and 100 in the present work) while keeping trade-off parameter β as 0.1.Figure 5 shows the L-curve drawn from the data misfit and model roughness for respective α values.From the figure, the points corresponding to α = 5.6 and α = 10 appear to be the "knee" of the curve, that is, the preferable models.The ratio of the variances, which were calculated from the normalized differences between the calculated and observed response functions, between these two models was 1.62.The variance ratio was regarded as an F-test statistic (Ichiki et al., 2021;Toh et al., 2006).Because the variance ratio is greater than the critical value of the 5% significance level (1.03), the Ftest suggests that the data misfit for α = 5.6 is significantly smaller than the data misfit for α = 10.However, the resultant electrical resistivity structures for α = 5.6 and α = 10 exhibit similar features (Figures S4 and S5 in Supporting Information S1).Thus, we selected the model for α = 5.6 as the preferred model because it explains the observed response functions significantly better.

Resultant Electrical Resistivity Structure
The subsurface resistivity structure of the preferred model is shown in Figures 6-8.When we look at the map view at a depth of 1 km (Figure 6a), the conductive areas (C1, C2, C3, and C4) beneath the sedimentary planes and basins are significant.These conductive anomalies at a depth of 1 km are spatially consistent with the negative or lower Bouguer anomalies (Figure 9).Thus, near-surface conductive anomalies are considered sedimentary layers containing relatively large porosities filled with low-density, conductive materials (aqueous fluid and/or smectite).Sato and Kato (2010) presented the upper crustal seismic cross section and its geological interpretation along line XX' in Figure 2. Figure 8a depicts the resistivity cross section along line XX'.The prominent shallow conductive layers (C1 and C2) are spatially consistent with the layers deemed to be Late Neogene sediments by Sato and Kato (2010).This consistency supports the validity of the proposed model.
In contrast, a high-resistivity area (>1,000 Ωm) was noticeable beneath the Asahi Mountains (R1), as shown in Figure 6.Resistive anomaly R1 continued to depths of several tens of kilometers (Figures 6f-6h).In the Asahi Mountains, Cretaceous basement rocks (granite, granodiorite, and tonalite) dating back to more than 60 Ma are widely exposed (Geological Survey of Japan, 2022).Therefore, the high resistivity appears to have resulted from old basement rocks with relatively low porosity.The resistive anomaly beneath the Asahi Mountains likely contributed to the along-arc direction of the induction arrows at higher periods mentioned above.
At a depth of 10 km, corresponding to the upper crust, the resistivity beneath the Shonai Plain is also higher than 1,000 Ωm (R2).At this depth, conductive areas appear locally beneath Quaternary volcanoes, for example,      study area.In addition, at the depths, conductive anomalies in the east side (C7 and C9) and in the southwest (C8) are prominent (<10 Ωm).The above features are also visible in the cross-sections (Figures 7 and 8).These cross sections are drawn along the profiles shown in Figures 2 and 6.The sounding curves of the observed and calculated response functions are compared in Supporting Information S1 (Figures S6-S29).The preferred model reproduces the respective response functions.The resultant normalized root mean square (NRMS) is 4.36, which was significantly smaller than that of the initial model (39.0).Ichihara et al. (2011) and Asamori et al. (2011) estimated 2-D resistivity structures using the data along lines CC' and FF', respectively.Comparing our result along lines CC' with the 2-D analysis result of Ichihara et al. (2011), the mid-lower-crustal fault-zone conductor estimated by Ichihara et al. (2011) is not prominent in our 3-D resistivity structure (Figure 7c).The conductive areas C5 and C6, which are adjacent to line CC', might have emphasized the faultzone conductor in the 2-D analysis because the apparent resistivity of the Z yx component upwardly deviated from the observed values when we made C5 and C6 resistive.In addition, the upper crust around line FF' in the 3-D structure (Figure 7f) is more resistive than that in Asamori et al. (2011).We compared the 2-D resistivity structures of the previous works and the structure of the present work along the same lines in Supporting Information S1 (Figures S30).
Because the subsurface resistivity structure is not 2-D, as visible in Figure 6, the 3-D resistivity structures obtained by this work should be more accurate.

Sensitivity Test
Sensitivity tests were performed to determine whether the observed data required characteristic electrical resistivity anomalies.Specifically, we performed forward calculations by changing the resistivities of the respective anomalies, and then compared the resultant NRMSs to those of the preferred model.For conductive anomalies, we made the resistivity higher than or equal to the thresholds, that is, if the resistivities of some parts of a target anomaly were lower than a threshold, the resistivities were changed to the threshold while maintaining the resistivities of the other areas.Threshold resistivities were 0.03, 0.1, 0.3, 1, 3, 10, 30, and 100 Ωm.Similarly, for the resistive anomalies, we set the resistivities of the resistive anomalies lower than or equal to the thresholds.Thresholds for the resistive anomalies were 1,000, 3,000, 10,000, and 30,000 Ωm.From the results of the forward calculations, we computed the variance ratios and compared them with a 5% significance level in the same manner as for determining the preferred model.As it is natural for the data of stations distant from a local anomaly to have negligible sensitivity to the resistivity change of the anomaly, we compared the variance ratios to the confidence level station by station, as in Ichiki et al. (2021).If the variance ratio was greater than the confidence level, at least for one station, we judged the increase in the NRMS to be significant.Table 1 lists the boundary resistivities between which the NRMS increases became significant from insignificant.For C5 and C6, we performed sensitivity tests separately for the upper (0-15 km in depth) and lower parts (15-40 km in depth).The ranges listed in Table 1 do not directly indicate the acceptable and unacceptable ranges of resistivity because resistivities other than the target anomaly did not change in the sensitivity tests.However, the results of sensitivity tests suggested that the observed data were sensitive to resistivity anomalies.The sounding curves computed using the preferred model and those changed for the sensitivity tests are compared in Supporting Information S1 (Figures S31-S43).In addition, we confirmed adequate resolutions to the conductive anomalies at 10-40 km in depth, as shown in Supporting Information S1 (Text S1 and Figures S44 and S45).

Relationship With Strain Concentration
As noted in Section 1, the study area contained the following strain concentration zones: SCZGS, SCZEJS, and SCZOBR.Ohzono et al. (2012) estimated the deformation mechanisms of strain concentration zones in Note.From the calculated responses after changing the resistivities of respective anomalies, we computed the variance ratios and compared them with the 5% significance level.We made the resistivity higher than or equal to thresholds (0.03, 0.1, 0.3, 1, 3, 10, 30, and 100 Ωm) for conductive anomalies (C1-C9), whereas we made the resistivity lower than or equal to thresholds (1,000, 3,000, 10,000, and 30,000 Ωm) for resistive anomalies (R1 and R2).
northeastern Japan by comparing the observed coseismic displacements of the 2011 off the Pacific coast of Tohoku Earthquake (Mw 9.0) with the predicted displacements in a half-space elastic medium.Figure 10 depicts the ratios of the observed co-seismic strain changes to the predictions (SC ratio) calculated by Ohzono et al. (2012).SC ratios higher than one, denoted by red, suggest a low elastic modulus in the upper crust, whereas SC ratios lower than one (blue) suggest a low viscosity in the lower crust (Ohzono et al., 2012).Based on the SC ratios, Ohzono et al. (2012) interpreted that the low-viscosity lower crust caused strain concentration in the SCZOBR, and the low-elastic-modulus upper crust caused strain concentration in the southernmost part of the SCZEJS (the northernmost part of the NKTZ).The SC ratio map (Figure 10) shows that shallow conductors (C1, C2, C3, and C4) correspond to areas with higher SC ratios (>1), and deep conductors (C7, C8, and C9) correspond to areas with lower SC ratios (<1).We estimated the strain concentration mechanisms in these zones by combining our resistivity structure with the SC ratio map determined by Ohzono et al. (2012).
The northern part of the SCZEJS contains a shallow conductive layer (C1 in Figure 10b) that corresponds to the thick sedimentary Shonai Plain, as stated in the previous section.C4 is located above the northern tip of the Echigo Plain, which is another thick sedimentary plain along the coast of the Sea of Japan.A conductive sedimentary layer can act as a low-elastic-modulus body responsible for strain concentration (Ichihara et al., 2013).
Because there was no prominent conductive area in the lower crust beneath the SCZEJS in our study area (Figure 10c), it appears that sedimentary plains such as C1 and C4 were responsible for the major portion of the deformation of the SCZEJS.This interpretation is consistent with the results of Ohzono et al. (2012) because the SC ratio is higher than one in the sedimentary plains along the eastern margin of the Sea of Japan (Figure 10a), suggesting a low elastic modulus of the upper crust.
In contrast, under the SCZOBR, a remarkable conductive area exists in the lower crust (C7 in Figure 10c).Previous EM studies have also identified conductive areas beneath the OBR (Ichihara et al., 2014;Mishina, 2009;Ogawa et al., 2014).Based on laboratory velocity measurements of xenoliths derived from the northeast Japan arc crust and a seismic tomography image (Nakajima et al., 2001), Nishimoto et al. (2008) suggested a supra-solidus temperature (>900°C) and the existence of 0.4-8.8vol% melts in the lower crust under the OBR.C7, a lowercrustal conductive anomaly under the OBR, has a resistivity lower than 10 Ωm (Figure 10c).Given that the temperature is 900°C and 6 wt% water content in the melt, the electrical resistivity of the melt is 0.58 Ωm.We calculated the melt resistivity by SIGMELTS (Pommier & Le-Trong, 2011) using the composition of a xenolith sample (hornblende gabbro) measured by Nishimoto et al. (2008) and the pressure (816 MPa) estimated by a onedimensional density structure.Even if we conservatively estimate the melt volume fraction by assuming the bulk resistivity is 10 Ωm and interconnected melt films completely wet the grain boundaries (Waff, 1974), the resultant estimate is 8.7 vol%, which is approximate to the upper bound of the melt fraction estimated by Nishimoto et al. (2008).If the bulk resistivity is lower than 10 Ωm, more melt volume fraction is required.The diffusion creep with 7 vol% of partial melting or higher decreases one order of magnitude in viscosity (∼10 19 Pa s) in the lower crust compared to the melt-free condition (Dimanov et al., 1998;Muto, 2011).
Seismic tomography showed that, under the SCZOBR, there are low-velocity and low-Q areas near the Moho and in the upper mantle (Hasegawa & Nakajima, 2004;Matsubara et al., 2022;Nakajima et al., 2013).Based on the findings of seismological studies and the locally high heat flow on the OBR (Matsumoto et al., 2022;Tanaka & Ishikawa, 2002), Hasegawa et al. (2005) suggested that fluids supplied from the mantle wedge lower the crustal strength beneath the SCZOBR, and the partly anelastic deformation contributes to the concentrated deformation along the SCZOBR.The low viscosity in C7 is consistent with their interpretation.Ohzono et al. (2012) also suggested that the lower crust has low viscosity under the SCZOBR because the SC ratio is smaller than one (Figure 10a).Thus, we consider the lower crustal conductive area under the SCZOBR (C7) as a fluid-rich and low-viscosity area that contributes to the high strain concentration observed on the OBR.
In the SCZGS, shallow conductive areas (eastern half of C1, C2, and C3) underlain by conductive lower crust (<100 Ωm) (Figures 10b and 10c).The shallow conductive areas correspond to thick sedimentary plains and basins.The SC ratio was higher than in this area (Figure 10a), suggesting a low elastic modulus in the upper crust.However, because the SCZGS is the strain concentration zone on a geological time scale (in the past approximately 3 million years), it is not appropriate to investigate its strain concentration mechanism based on geodetic studies.It should be noted that several active faults exist in the SCZGS.The SCZGS was defined based on these faults and folds (Okamura, 2002;Okamura et al., 2007).In other words, deformation caused by large intraplate earthquakes on these faults contributed to the SCZGS.As discussed in the next section, mechanically weak zones in the crust under the SCZGS, imaged as conductive areas, may play a key role in large intraplate earthquakes in the SCZGS, leading to strain concentration on a geological time scale.
Figure 12 shows a schematic illustration of our interpretation of the strain concentration mechanism along line XX' (Figures 2 and 10a), which crosses the boundary between the higher (>1.02) and lower (<0.98)SC ratios.As noted above, shallow sedimentary plains and basins, which were imaged as shallow conductive layers (C1, C2, and C3 in Figures 6a and 10b), are considered to behave macroscopically as low-elastic-modulus bodies to the west of the boundary.In contrast, to the east of the boundary, the hot and fluid-rich lower crust, imaged as an outstanding lower-crustal conductor (C7 in Figures 6f-6h and 10c), is considered to play a major role in strain concentration as a low-viscosity body.For the SCZGS, deformation by intraplate earthquakes on active faults contributes to the strain concentration zone on a geological timescale.Resistivity structure modeling allows us to identify the 3-D position of the heterogeneity causing the strain concentration, whereas the strain field investigated by geodetic and geological studies reflects the structures of whole depths.Thus, to elucidate the 3-D stressstrain field in the crust, it would be essential to image the high-resolution subsurface structure in addition to performing dense geodetic measurements.

Relationship With Large Intraplate Earthquakes
Next, we discuss the relationship of the resistivity structure to large intraplate earthquakes.It would be better to mention that we do not discuss interplate nor intraslab earthquakes because the survey area is more than 300 km away from the Japan trench and the top of the Pacific slab is at a depth deeper than 100 km beneath the survey area (the locations of our survey area and the Pacific slab surface are compared in Figure S46 in Supporting Information S1).
Figure 10c shows the epicenters of historical inland earthquakes (M > 5) (Usami et al., 2013) and the inland earthquakes (M > 5) determined by Japan Meteorological Agency (1997-2021).Most of these large intraplate earthquakes occurred around the rim of the SCZGS because the SCZGS in our study area is partly bounded by active faults (Figures 1, 10b and 10c).These active faults were originally as normal faults owing to extensional back-arc rifting in the Early to Middle Miocene, associated with the opening of the Japan Sea (Sato, 1994).Under the present-day compressional stress field since the Late Pliocene, these faults were reactivated as reverse faults (Sato, 1994).In the Shonai Plain fault zone (Figure 10c), which is one of the most active fault zones in the area, two M 7.0 earthquakes occurred: the 1804 Kisakata and 1894 Shonai earthquakes.Kato (2014) suggested the importance of the ancient rift system for nucleating large intraplate earthquakes along the eastern margin of the Japan Sea based on seismic studies of recent large intraplate earthquakes (2004Niigataken Chuetsu, 2007Niigata-ken Chuetsu-oki, and 2007 Noto-Hanto earthquakes).His suggestion is presumably applicable to our study area because the active faults in the study area belong to a Miocene rift system, as is the case for the earthquake source faults of the 2004 Niigata-ken Chuetsu and 2007 Chuetsu-oki earthquakes (Kato et al., 2009(Kato et al., , 2010)).Figures 10b and 10c illustrates that the shallow conductive layers (C1, C2, and C3) and lowercrustal conductors (C5, C6, and C7) are located in the SCZGS.Beneath the gaps between the shallow conductors in the upper crust, there are deep conductors in the lower crust.As we discussed above, both of them can act as mechanically weak zones.Hence, those shallow and lower-crustal conductors beneath the SCZGS are considered to cause stress loading on the brittle parts of the pre-existing faults within the ancient rift system.The resultant stress accumulation presumably leads to large intraplate earthquakes in the which also contribute to the strain concentration over geological time.
Figure 10c shows that the epicenters of large earthquakes (M > 5) are located near the edges of conductors in the middle to lower crust.Similar spatial relationships have also been reported in other regions (Cai et al., 2017;Ichihara et al., 2014;Sun et al., 2020).It should be noted that pore fluids can also work more directly to generate large earthquakes, in addition to local stress accumulation around conductive areas.Pore fluid can reduce the frictional strength of faults (e.g., Scholz, 2019).The migration of pore fluids also contributes to earthquakes (Ogawa et al., 2001).Recently, Aizawa et al. (2021) suggested that ruptures of large earthquakes selectively initiate near the outer edges of conductors, based on the high-resolution resistivity structure and hypocenters of the 2016 Kumamoto earthquake sequence (M 6.5, M 7.3, and aftershocks).They proposed that the large pressuretemperature gradient of pore fluids near the edge of the conductors contributes to the high probability of large earthquakes.Because large earthquakes occur near the edges of conductors (Figure 10c), the results of this study are consistent with their suggestion.

Relationship With Volcanic Activities
Finally, we discuss the relationship between the resistivity structure and volcanic activity.Under the Quaternary volcanos below which deep low-frequency earthquakes (LFEs) occur, the conductive areas are rooted in the deep crust (Figures 6-8), as noted in Section 3. In this subsection, we focus on the conductors beneath the Chokaisan and Gassan volcanoes (C5 and C6).Figures 7a and 7e show that the roots of the conductive areas exist east of the respective volcanoes.Okada et al. (2015) imaged a low-seismic velocity anomaly east of the Gassan, and the island-arc-scale velocity map of Matsubara et al. (2022) shows that the lower crust east of the Chokaisan volcano has a low velocity.C5 and C6 revealed in the present work are spatially consistent with the low-velocity anomalies found by seismic tomography.
Figures 7a and 7e show that the hypocenters of the LFEs, denoted by white-filled circles, are seemingly located within the conductors (C5 and C6).However, the LFE hypocenters were placed near the edges of the conductors, as evident from the map views (Figures 6g and 6h). Figure 8b depicts the cross section along line YY', which crosses the eastern foot of Chokaisan volcano along a nearly N-S direction.LFEs occurred around the northern edges of the conductor (C5).These adjacencies of LFEs to the edges of conductive anomalies have also been reported in other MT studies (e.g., Aizawa et al., 2004Aizawa et al., , 2022)).In northeast Japan, by combining the seismic velocity structure and relocated LFEs, Niu et al. (2018) suggested that most LFEs are located at the edges of lowvelocity and high-Poisson's ratio areas.They interpreted that the LFE hypocenters were located at the boundary of magma transport channels or reservoirs.The presence of melts can lower bulk rock resistivity because the melt resistivities are significantly lower than those of dry rocks (Figure 11).Thus, it would appear that C5 and C6 represent the trans-crustal magmatic systems (e.g., Cashman et al., 2017) under the Chokaisan and Gassan volcanoes, respectively.However, more dense observations, particularly in the along-arc direction, are required to elucidate detailed structures of the magmatic systems, including magma transport paths just beneath the volcanoes.
Because high 3 He/ 4 He values (>5 × 10 6 ) have been reported over the entire study area (Asamori et al., 2011;Horiguchi et al., 2010;Sano & Wakita, 1985) Sato and Kato (2010).A linear pattern shows the sedimentary plains and basins inferred by Sato and Kato (2010).It would appear that the high strain rates of the SCZEJS and SCZOBR are ascribable to thick sedimentary layers and lowviscosity bodies in the lower crust, respectively.As for the SCZGS, the intraplate earthquakes along the active faults is considered to contribute to the strain concentration in the geological time scale.(Kariya & Shankland, 1983).The red, pink, and orange lines indicate the resistivities of NaCl-H 2 O fluids with salinities of 0.58, 3.38, and 9.52 wt %, respectively (Sakuma & Ichiki, 2016).Blue dashed lines denote the resistivities of dacitic melts (Laumonier et al., 2015), whereas green dashed lines denote the resistivities of andesitic melts (Laumonier et al., 2017).
frequency response functions, respectively.ϕ c (m) denotes the regularization term for the galvanic distortion (Equation 33 of Usui, 2015).Positive real values α and β are the trade-off parameters that control the influence of each term.After linearizing Equation A2, the model parameters were calculated using the data-space Gauss-Newton method, as described by Usui et al. (2017).

Figure A1.
Relationship of the electric field of the adjacent elements sharing a face with hanging nodes.As for the redcolored edges, the same electrical field value as the coarser element is assigned to the finer element.For each edge connecting two hanging nodes, denoted by a blue line, the average electric field of the two parallel edges of the coarser element is given.
Hanging nodes are located only on horizontal faces.

Figure 2 .
Figure 2. Map of the study area with elevation gradation.Observation stations used in this study are shown as colored circles and inverted triangles.Along the respective dash lines in this figure, vertical cross sections of the resultant resistivity structure are depicted (Figures 6-8).The observation stations along line FF' were deployed by Japan Atomic Energy Agency(Asamori et al., 2011), whereas the other stations were deployed by the "Multidisciplinary Research Project for High Strain Rate Zones"(2008)(2009)(2010)(2011).Active faults(Nakata & Imaizumi, 2002) are denoted by solid red lines.Pink diamonds indicate the hypocenters of the low-frequency earthquakes from 1997 to 2021 determined by the Japan Meteorological Agency.Most low-frequency earthquakes occurred at depths ranging 20-40 km.
Figure 2. Map of the study area with elevation gradation.Observation stations used in this study are shown as colored circles and inverted triangles.Along the respective dash lines in this figure, vertical cross sections of the resultant resistivity structure are depicted (Figures 6-8).The observation stations along line FF' were deployed by Japan Atomic Energy Agency(Asamori et al., 2011), whereas the other stations were deployed by the "Multidisciplinary Research Project for High Strain Rate Zones"(2008)(2009)(2010)(2011).Active faults(Nakata & Imaizumi, 2002) are denoted by solid red lines.Pink diamonds indicate the hypocenters of the low-frequency earthquakes from 1997 to 2021 determined by the Japan Meteorological Agency.Most low-frequency earthquakes occurred at depths ranging 20-40 km.

Figure 3 .
Figure 3. Phase tensor ellipses and induction vectors (Parkinson convention).The color of the ellipse in (a), (b), and (c) depends on Φ min , Φ max , and β (Caldwellet al., 2004), respectively.Black arrows in (c) denote the real part of the Parkinson vectors.Because the periods at which the response functions were estimated were not the same between the MTU and the ADU stations, except for the long-period responses (≥42.67 s), periods slightly differ among the two measurement systems in the two left columns.Red triangles indicate Quaternary volcanos.Shaded areas in the leftmost columns are the schematic illustration of sedimentary plains and basins (Geological Survey of Japan, 2022).The response functions that appeared to be heavily affected by cultural noises were not plotted.

Figure 4 .
Figure 4. Computational mesh used for the 3-D MT inversion.The positive x-direction corresponds to N20°E.(a) Map view of the surface mesh around the survey area.(b) Surface mesh around the Chokaisan and Gassan volcanoes viewed from an obliquely upward part.

Figure 5 .
Figure 5. L-curve obtained from the results of the 3-D MT inversion with various trade-off parameter α values.The horizontal and vertical axes represent the model roughness term and data misfit term of the objective function, respectively.Numbers adjacent to the points are the corresponding α values.The preferred model selected in this study is the model for α = 5.6.

Figure 6 .
Figure 6.Map views of the estimated electrical resistivity structure.Small gray squares indicate MT stations.Active faults(Nakata & Imaizumi, 2002) are denoted by solid gray lines.Open circles and white diamond marks indicate the hypocenters of regular and low-frequency earthquakes, respectively, determined by the Japan Meteorological Agency.Hypocenters located within ±5 km of the respective depth are shown.In the study area, the Moho depth ranges 30-35 km(Matsubara et al., 2017), and the bottom of the upper crust is approximately 15 km(Iwasaki et al., 2001).

Figure 7 .
Figure 7. Vertical cross sections of the estimated electrical resistivity structure.The locations of those profiles are shown in Figures 2 and 6.Black inverted triangles indicate MT stations.Open circles and white diamond marks indicate the hypocenters of regular and low-frequency earthquakes (1997-2021), respectively, determined by the Japan Meteorological Agency.The hypocenters located within ±5 km of respective profiles are shown.In the study area, the Moho depth ranges 30-35 km(Matsubara et al., 2017), and the bottom of the upper crust is approximately 15 km(Iwasaki et al., 2001).

Figure 8 .
Figure 8. Vertical cross sections of the estimated electrical resistivity structure.(a) Cross section along the line where Sato and Kato (2010) showed a geological and seismic section.(b) Nearly N-S cross-section through the eastern foot of Chokaisan volcano.The locations of those profiles are shown in Figures2 and 6.Open circles and white diamonds indicate the hypocenters of regular and low-frequency earthquakes, respectively, determined by Japan Meteorological Agency.The hypocenters located within 5 km of respective profiles are shown.In the study area, the Moho depth ranges 30-35 km(Matsubara et al., 2017), and the bottom of the upper crust is approximately 15 km(Iwasaki et al., 2001).

Figure 9 .
Figure 9.Comparison of the gravity anomaly with the map view of the near-surface resistivity structure.(a) Bouguer anomaly assuming a terrain density of 2.67 g/cm 3 .The Bouguer anomaly map was drawn using a database edited by the Geological Survey of Japan, AIST (2013).(b) Estimated resistivity structure at a depth of 1 km.

Figure 10 .
Figure 10.(a) Map of the ratio of the observed coseismic strain changes to the predictions (SC ratio) calculated by Ohzono et al. (2012).(b), (c) Electrical resistivity structure at depths of 1 and 20 km, respectively.The thick broken lines are the boundaries of the SCZEJS and SCZOBR, whereas the horizontal-line pattern denotes the SCZGS.The thin dash line in (a) indicates line XX', along which the vertical cross sections in Figures 8a and 12 are depicted.Active faults (Nakata & Imaizumi, 2002) are denoted by solid gray lines in (b) and (c).Blue stars indicate the epicenters of the historical inland earthquakes (M > 5)(Usami et al., 2013) and the inland earthquakes (M > 5) determined byJapan Meteorological Agency  (1997-2021).Years and magnitudes of respective earthquakes are written adjacent to the epicenters.
Figure 12.(a) Schematic image of the relationship between the strain concentration zones and the subsurface structure along line XX' of Figures 2 and 10.(b) Vertical cross sections of the estimated electrical resistivity structure along line XX' (same as Figure 8a).Open circles and white diamonds represent the hypocenters of regular and low-frequency earthquakes (1997-2021) (Japan Meteorological Agency) located within ±5 km of the profile.Black lines indicate the fault configurations shown bySato and Kato (2010).A linear pattern shows the sedimentary plains and basins inferred bySato and Kato (2010).It would appear that the high strain rates of the SCZEJS and SCZOBR are ascribable to thick sedimentary layers and lowviscosity bodies in the lower crust, respectively.As for the SCZGS, the intraplate earthquakes along the active faults is considered to contribute to the strain concentration in the geological time scale.

Table 1
Boundary Resistivities Between Which the NRMS Increases Became Significant From Insignificant