On a large magmatic fluid reservoir oblique to the volcanic front in the southern part of NE Japan revealed by the magnetotelluric survey

Many active volcanoes and various types of seismic activities exist in the southern part of the Northeast Japan subduction zone. One of the geologically most interesting features in this area is the sequential explosive eruptions of a group of volcanoes. The group consists of Mt. Azuma and Mt. Adatara on the volcanic front line, Mt. Bandai west of the volcanic front, and Mt. Numazawa on the back-arc side. A previous petrological study on the eruption products regarded Mt. Numazawa as an anomalous back-arc volcano because its lavas are similar to those of volcanoes on the volcanic front. The reason behind this unique connection was unclear, and hence, this study was intended to understand the deep fluids distribution beneath the area. For this purpose, a 3-D regional electrical resistivity structure was estimated from a series of wide-band magnetotelluric surveys, with 45 observation points deployed from the fore-arc to the back-arc sides. The most important feature of the resistivity structure is a large conductive zone in the central part of the area, spanning from the upper mantle to the lower crust. Interestingly, the lateral elongation of the conductor is oblique to the volcanic front line and consistent with the spatial distribution of the group of volcanoes and the low-frequency earthquake clusters around them. Therefore, the conductor most likely represents a large, elongated magmatic fluid reservoir beneath the volcanoes. Hydrous partial melting might be the cause of the enhanced conductivity.


Introduction
The southern part of Northeast (NE) Japan is a region full of in-land crustal activities due to the ongoing subduction of the Pacific plate toward the Eurasian plate.One of the geologically most interesting features in this area is the distribution and history of quaternary volcanoes in the central part, including Mt. Azuma and Mt.Adatara on the volcanic front, Mt.Bandai just westward of the volcanic front, and Mt.Numazawa on the back-arc side (see Fig. 1).The spatial distribution of deep, low-frequency earthquakes (DLFE) is consistent with the position of volcanoes (Niu et al. 2018;Kurihara and Obara 2021;Fig. 1).A tephrochronology study by Kimura (1996) revealed that explosive eruptions often occurred sequentially from the volcanic front to the back-arc side, starting with Mt.Adatara and continuing to Mt. Azuma, Mt.Bandai, and lastly, Mt.Numazawa.One set of sequential eruptions took only several tens of thousands of years (Kimura 1996).Because such a period is geologically very short, Kimura (1996) speculated that movements of magmatic fluids from the upper mantle occurred almost simultaneously.Kimura and Yoshida (2006) found that Mt.Numazawa lavas preserved the characteristics of mantle-driven basaltic melts like those of volcanoes on the volcanic front, especially Mt.Azuma and Mt.Adatara, and not like those of most backarc volcanoes in NE Japan.Because of this, Kimura and Yoshida (2006) regarded Mt.Numazawa as an anomalous back-arc volcano.However, before this study, the reason behind this interesting connection of volcanoes remained an open scientific question.
The magnetotelluric (MT) method is a passive geophysical exploration method that probes the electrical resistivity structure of the Earth based on the measurement of natural electromagnetic fields on the surface (Simpson and Bahr 2005; Chave and Jones 2012).One of the most useful applications of the MT method is the imaging of subsurface fluid distribution.In NE Japan, MT has been very useful in revealing the role of fluids in different scenarios, such as the upwelling flow of fluids in the upper mantle (e.g., Asamori et al. 2011;Ichiki et al. 2015;Motoyama et al. 2020), localized magmatic fluids beneath active volcanoes (e.g., Ogawa et al. 2014;Ichiki et al. 2021), and fluids responsible for seismic activities (e.g., Ogawa et al. 2001;Uyeshima et al. 2005;Umeda et al. 2015;Ichihara et al. 2014Ichihara et al. , 2023)).However, MT studies in the southern part of NE Japan were still limited.Motoyama et al. ( 2020) performed a wide-band MT study with stations deployed across the arc, but their southernmost profile was just north of Mt.Azuma (Fig. 1b).Another study focused on the south of Mt.Nasu, the fore-arc side, to understand the earthquake swarm activity on the fore-arc side induced after the 2011 Great Tohoku Earthquake (Umeda et al. 2015: Fig. 1b).On the back-arc side, Uyeshima et al. (2005) investigated the 2004 M6.8 Mid-Niigata earthquake, close to the swarm of earthquakes south of Echigo Plain (Fig. 1b).None of these studies have investigated the fluid distribution beneath the group of volcanoes in the central part.
This study estimated a three-dimensional (3-D) regional electrical resistivity structure across the southern part of NE Japan.This area has yet to be investigated by previous MT studies.The main objective was to reveal fluidrich conductive zones associated with the volcanism and DLFEs in the central part of the area.The following section describes the data acquisition, estimation of response functions, and 3-D inversion.Then, we discuss the resulting resistivity structure and its interpretation.

Data acquisition
The array of electromagnetic recording stations is shown in Fig. 1b.We deployed the stations along three parallel NW-SE survey lines, namely L01 for the northern line, L02 for the middle line, and L03 for the southern line.There are 14, 15, and 16 stations at L01, L02, and L03, respectively, with about 10 km station intervals.At stations denoted by blue inverted triangles in Fig. 1b, two horizontal components of the electric field ( E x and E y ), two horizontal components of the magnetic field ( H x and H y ), and a vertical component of the magnetic field ( H z ) were recorded using the Metro- nix ADU-07 or ADU-07e data logger with induction coils MFS-06, MFS-06e, and/or MFS-07e.At stations denoted by green inverted triangles in Fig. 1b, only E x and E y were recorded using the ELOG1K data loggers from NT System Design Inc, Japan.There were three survey periods.The first was in 2015 for L02 stations, the second was in 2016 for L01 stations, and the third was in 2017 for L03 stations.At each period, the stations were deployed simultaneously for about 20 days.The data were sampled at 32 Hz (continuously) and 1024 Hz (intermittently).

Response functions
We estimated the MT impedance, the vertical magnetic field transfer function (VMTF), and, in addition, the inter-station horizontal magnetic field transfer function (HMTF).
The MT impedance tensor Z is defined as the transfer function between the horizontal electric field components ( E x and E y ) and horizontal magnetic field compo- nents ( H x and H y ).It is a 2 × 2 complex tensor and can be expressed as For ELOG1K stations, we used magnetic field data at neighboring stations to calculate MT impedance.Apparent resistivity ρ a and phase φ are then defined as (1) where µ is the magnetic permeability of the medium and ω is the angular frequency of the EM fields.From the MT impedance, we calculated the phase tensor to evaluate the dimensionality of the structure (Caldwell et al. 2004).We found that skew angles of longer period data (> 10 s) often exceed 3°, a three-dimensionality threshold recommended by Caldwell et al. (2004), suggesting that the data must be modeled using 3-D inversion.See Additional file 1 for phase tensor maps at various periods.
The vertical component of the magnetic field is sensitive to lateral variations in resistivity.Its relationship to the horizontal magnetic field can be expressed as where T = [ T zx T zy ] T is a complex vector of VMTF.
The inter-station HMTF tensor H relates the horizontal components of the magnetic field at a local station and that at a reference station.It can be expressed as Superscripts i and ref denote the local and reference sta- tions, respectively.Due to the excellent quality of magnetic field data in the first survey period, we decided to use the easternmost station on L02 as the reference station for HMTF analysis (see Fig. 1b).We then continued to deploy this station during the second and third survey periods.
HMTF could be easily calculated under the standard MT field operation provided simultaneous measurement at the local and the reference stations.Previous studies have demonstrated the advantage of using HMTF in the inversion.Numerical experiments by Campanyà et al. (2016) showed that joint inversion with HMTF data improves the accuracy of the inverted 3-D model, especially high resistivity anomalies that are difficult to resolve using MT impedance.Habibian and Oskooi (2014) compared synthetic inversions using HMTF and VMTF data.They showed that HMTF inversion recovered better structure when there is a conductive overburden, whereas VMTF inversion is more sensitive to lateral resistivity changes.
The time series was processed into frequency-domain response functions using the BIRRP code (Chave and Thomson 2004).To obtain better response functions, magnetic fields at a station in the Yamagata Prefecture, more than 100 km north of the survey line, were used as remote (2) (3) reference data (Gamble et al. 1979;Chave and Thomson 2004).

3-D inversion
FEMTIC inversion code (Usui 2015(Usui , 2021;;Usui et al. 2017) was used to obtain the 3-D resistivity structure as it can jointly invert all the response functions used in this study.The Earth was discretized into 94 × 134 × 82 hexahedral elements.We used the N40°E(x)-E40°S(y) coordinate system so that the y-axis is approximately parallel to the survey lines.In the target area with the array of stations, the lateral spacing of the grid is 2 km, which is the finest mesh size.The thickness of the first layer is 0.1 km.Land topography and seafloor bathymetry based on Tozer et al. (2019) were incorporated in the model for the inversion.The starting model of the inversion was a uniform 100 Ωm Earth, including the seawater (0.25 Ωm) and air layer ( 10 10 Ωm) that were kept unchanged during the inversion.The input data of the inversion were the real and imaginary parts of all the components of the MT impedance tensor, VMTF, and HMTF, calculated at the coordinate system described above.We inverted the data at 15 representative periods from 0.1 to 10,000 s, evenly spaced in the logarithmic scale.The error floor was set at 5% for all the response functions.
In the inversion, the following objective function was minimized, where �(m) is the total objective function, � d (m) is the data misfit term, � m (m) is the model roughness term, and α 2 is a trade-off parameter.The optimum α was determined objectively using the L-curve criteria (Hansen 1992).Among different inversion results obtained with different α , we found that the model at α = 3 corresponds to the knee point of the L-curve (see Additional file 1 for details).Thus, we selected this model as the preferred 3-D resistivity structure.

Result
Cross-sections of the preferred 3-D resistivity structure are shown in Fig. 2. A large, low-resistivity zone of about 20 Ωm exists in the central part of the area (C1).The upper boundary of C1 is about 20 km in depth.The horizontal cross-sections at 40 and 60 km depth show that C1 is laterally elongated in the NE-SW direction, nearly orthogonal to the survey lines (Fig. 2f, g).Also, at these depths, another moderately low-resistivity zone (C2) is found on the fore-arc side north of L01.The resistivity of C2 is about 40 Ωm.South of C2, the structure is very resistive, with a resistivity value exceeding 10,000 Ωm.This alongarc contrast in the fore-arc side differs from the back-arc ( 5) side, where the resistivity structure is relatively uniform.At a depth of 10 km, a 10 Ωm conductor exists on the fore-arc side south of L03, and a 20 Ωm conductor exists on the back-arc side of the survey lines (Fig. 2d).
The model explains the data with an RMS misfit of 3.66, calculated without the error floor.Figure 3 shows the observed and calculated response functions at three different periods.MT impedance tensor is shown in terms of the apparent resistivity and phase of the effective impedance, which is one of the rotational invariants of MT impedance (e.g., Pedersen and Engels 2005).The magnitude of the HMTF tensor and the induction vector are also shown.A complete comparison of sounding curves between the observed and calculated data is available in Additional file 1.We checked the robustness of the inversion result by running the inversion using different starting resistivity values: 30, 100, and 300 Ωm.The resultant structures and RMS misfits were very similar.
Among the features revealed in the structure, we focused more on C1 because it is distributed in the target area with many observation stations and is geologically most interesting.To investigate C1 further, a simple sensitivity test was performed.We replaced the resistivity value of C1 (the area enclosed by the red dashed line in Fig. 2f ) from about 20 Ωm to 30, 50, 60, 70, and, finally, 100 Ωm, which was the starting resistivity value and calculated the response functions, respectively.We found that the modification mainly affected ρ a and φ data in longer periods (> 10 s) at the stations above the conductor, i.e., S05-S08 at L01, L02, and L03.Difference of the sounding curves, respectively, between the preferred and test models are shown in the Additional file 1.We then performed a statistical F-test to quantitatively assess the significance of the RMS difference between models, as was done in some previous MT studies (e.g., Ichiki et al. 2021;Ichihara et al. 2023).For the F-test, we considered the RMS misfit without the error floor.Because, unlike the statistical data errors calculated by the BIRRP code (Chave and Thomson 2004), the error floor had no statistical meaning and mainly needed to stabilize the inversion, the use of the error floor was not appropriate for the statistical F-test.The RMS of the preferred inversion model where C1 was about 20 Ωm was 3.66.When C1 was replaced by 30, 50, 60, 70, and 100 Ωm, the RMS misfits increased to 3.68, 3.74, 3.79, 3.82, and 3.96, respectively, as shown in Fig. 4a.According to a left-tailed F-test at a 95% confidence level (degree of freedom = 2464), the threshold RMS was 3.78.So, the F-test demonstrates that C1 must be less resis- tive than 60 Ωm.In Fig. 2f, a red dotted line divides C1 laterally into northeastern and southwestern parts.The northeastern part is outside the station array and has a similarly low resistivity structure as in the southwestern part beneath the station array; hence, C1 seems to elongate in the NE-SW direction, consistent with the spatial distribution of active volcanoes.To confirm the lateral extent of C1, we performed a sensitivity test for only the northeastern part by replacing the resistivity values from around 20 Ωm to 30, 50, 70, and 100 Ωm.As in Fig. 4b, the RMS misfits increased from 3.66 to 3.67, 3.72, 3.79, and 3.89, respectively, while the threshold RMS based on a left-tailed F-test at a 95% confidence level was 3.78.RMS misfits of the 70 and 100 Ωm test  The blue line denotes the RMS misfit of the final inversion model shown in Fig. 2, whereas the red line denotes the threshold RMS based on a 95% confidence limit F-test.Here, the RMS misfits were calculated without the error floor models surpassed the threshold RMS, indicating that the northeastern part of C1 needs to be conductive (< 70 Ωm).This test indicates that the laterally elongated C1 is required to explain the data.C2 is laterally outside of the array of stations, yet we also examined its sensitivity in the same manner.The resistivity value of C2, which was about 40 Ωm in the preferred inversion model, was replaced by 50, 70, and 100 Ωm.It turned out that ρ a and φ at S13-S15 on L01 and L02 were slightly affected by the modification.However, adopting the same F-test criteria used in the previ- ous C1 tests, the RMS misfits of the test models did not exceed the threshold RMS.The RMS of the preferred inversion model where C2 was 40 Ωm was 3.66.Even when C2 was replaced by 50, 70, and 100 Ωm, the RMS misfits changed only to 3.66, 3.67, and 3.71, respectively, while the threshold RMS was 3.78 (Fig. 4c).This test indicates that the data poorly constrain the low resistivity of C2.Thus, we refrained from interpreting C2 further.

Discussion
We superimposed hypocenters, volcanoes, and the Moho depth on the preferred resistivity structure in Fig. 5. Considering the Moho depths, the central conductor (C1) spans from the uppermost mantle to the lower crust.The conductor lies beneath the group of volcanoes, and DLFEs occurred within the upper body of the conductor.According to Niu et al. (2018), who studied various types of DLFE in NE Japan, DLFEs occurring in volcanic areas are often associated with deep magmatic fluid activity.Previous MT studies in different areas in NE Japan have also revealed conductive fluid-rich areas spatially consistent with DLFE clusters (e.g., Asamori et al. 2011;Ichiki et al. 2015;Motoyama 2020;Motoyama et al. 2020).Thus, it is possible that the conductive area C1 is rich in magmatic fluids.Katz et al. (2003) stated that hydrous partial melting at 1GPa requires a temperature greater than 900 ℃.The ground temperature structure of Matsumoto et al. (2022) showed that the temperature around the top of C1 is roughly 800 ℃.Although Matsumoto et al. (2022) did not show the temperatures at depths greater than 20 km, assuming increasing temperature with depth, the temperature would reach the solidus at some depths greater than 20 km.Therefore, we suggest that the conductor represents a deep magmatic fluid reservoir containing interconnected electrically conductive hydrous melts.
The lateral elongation of the conductor in the NE-SW direction is consistent with the spatial distribution of active volcanoes consisting of Mt.Azuma, Mt.Adatara, Mt.Bandai, and Mt.Numazawa, and the occurrence of DLFEs around them (dashed line in Fig. 5d).North of our survey lines, Motoyama et al. (2020) revealed that the uppermost mantle structure beneath Mt.Azuma is conductive.This result supports the lateral elongation of the conductor.As mentioned in the introduction, petrological studies found that the group of volcanoes erupted sequentially and periodically (Kimura 1996) and that the lavas of Mt.Numazawa have similar characteristics to those of volcanic front volcanoes (Kimura and Yoshida 2006).Considering this evidence, our study seems to confirm that the deep source of fluids of those volcanoes is connected, which is imaged as the elongated conductor.Although fluids might have taken different paths to the surface, a common source in the upper mantle may account for the connection of eruptions.
The reason why the fluid-rich area deviates from the volcanic front to the back-arc side requires further studies, but there are previous findings that might be related to it.The seismic S-wave anisotropy study by Nakajima et al. (2006) indicated a relatively complex mantle flow under the study area.The complex mantle flow might be associated with the fluid-rich area deviating from the volcanic front line.Another potential candidate for the deviation is early slab dehydration.To the south of our survey lines, Kato et al. (2013) and Umeda et al. (2015) revealed a zone rich in aqueous fluids in the lower crust on the fore-arc side.If early slab dehydration under the fore-arc side caused the lower-crustal fluids, it might have also caused a relatively dry upper mantle under the volcanic front.

Conclusion
In order to understand volcanic and seismic activities in the southern part of NE Japan, we estimated a 3-D regional resistivity structure beneath the area from a series of wide-band MT surveys.A total of 45 MT stations were deployed from the fore-arc to the back-arc sides.MT impedance, VMTF, and HMTF data at periods from 0.1 to 10,000 s were used as the input of the inversion.The most important feature of the resistivity structure is a large conductive zone in the central part of the area, spanning from the upper mantle to the lower crust.Sensitivity tests suggest that a resistivity value less than 70 Ωm is necessary to explain the data satisfactorily.The conductor is laterally oblique to the volcanic front line and consistent with the spatial distribution of active volcanoes such as Mt.Azuma, Mt.Adatara, Mt.Bandai, and Mt.Numazawa and some DLFE clusters around them.Considering previous petrological findings about the connection of those volcanoes, we infer that the conductor represents a large and elongated magmatic fluid reservoir supplying the volcanoes.Deviation from the volcanic front line might be attributed to the apparently complex mantle wedge flow or unusually early slab dehydration, but it requires further investigation.

Fig. 1 a
Fig. 1 a Japan subduction system.Dashed lines reflect trenches, and arrows indicate the direction of motion of the respective plates relative to the Eurasia Plate.Red triangles denote active volcanoes based on the Japan Meteorological Agency (JMA) database.The red band highlights the volcanic front line in NE Japan.b Map of the study area.Inverted triangles denote MT stations.At blue stations, we measured two horizontal electric field components, two horizontal magnetic field components, and the vertical magnetic field component.At green stations, we only measured two horizontal electric field components.Brown lines are active faults taken from Nakata and Imaizumi (2002).Normal-frequency earthquakes (black dots) and low-frequency earthquakes (orange circles) are based on the JMA Unified Catalog from 1997 to 2019, shallower than 60 km.AZ Azuma, AT Adatara, BD Bandai, NZ Numazawa, NS Nasu

Fig. 2
Fig. 2 The preferred 3-D resistivity structure.Vertical cross-sections beneath L01, L02, and L03 are shown in a, b, and c, respectively.Horizontal cross-sections at 10, 20, 40, and 60 km depths are shown in d, e, f, and g, respectively.Inverted triangles denote the MT stations, and the purple one denotes the HMTF reference station.Red triangles denote active volcanoes.In the first sensitivity test of C1, the whole area within the red dashed line in f was replaced by testing resistivity values, whereas in the second sensitivity test of C1, only the northeast area of the red dotted line was tested.In the sensitivity test of C2, the area enclosed by the blue dashed line was replaced by testing resistivity values

Fig. 3 Fig. 4
Fig.3Map of the observed and calculated response functions at three different periods.The coordinate system is adjusted so that the y-axis is horizontal.In a and b, the color of the circles denotes the apparent resistivity and phase of the effective impedance, respectively.In c, the color of the circles denotes the magnitude of the HMTF tensor, and the arrows denote the Parkinson vector(Parkinson 1962).The circle with a red outline is the HMTF reference station.Data with large errors were excluded

Fig. 5
Fig. 5 Interpreted 3-D resistivity structure.In a, b, and c, grey dots and orange circles are 1997-2019 JMA hypocenters of normal-frequency earthquakes and DLFEs, respectively, within 10 km from the profile.The white dashed line is the Moho (crust-mantle) boundary based on Matsubara et al. (2017).In d, earthquake hypocenters within a 6 km distance from the 40 km depth are plotted.The red band denotes the volcanic front line.The grey dashed line enclosed the estimated area of a deep magmatic fluid reservoir under the cluster of volcanoes and DLFEs.Inverted triangles denote the observation stations, and the purple one is the HMTF reference station