Magma reservoir beneath Azumayama Volcano, NE Japan, as inferred from a three-dimensional electrical resistivity model explored by means of magnetotelluric method

An electrical resistivity model beneath Azumayama Volcano, NE Japan, is explored using magnetotelluric method to probe the magma/hydrothermal fluid distribution. Azumayama is one of the most concerning active volcanoes capable of producing a potential eruption triggered by the 2011 Tohoku-Oki Earthquake. The three-dimensional resistivity model reveals a conductive magma reservoir (< 3 Ωm) at depths of 3–15 km below sea level (bsl). The 67% and 90% confidence intervals of resistivity are 0.2–5 Ωm and 0.02–70 Ωm, respectively, for the magma reservoir. We assumed dacitic melt + rock at a shallow depth of 4 km bsl and andesitic melt + rock at a greater depth of 9 km bsl. The confidence interval of resistivity cannot be explained by using dacitic melt + rock condition at a depth of 4 km bsl. This suggests that very conductive hydrothermal fluids coexist with dacitic melt and rock in the shallow part of the magma reservoir. For the depth of 9 km bsl, the 67% confidence interval of resistivity is interpreted as water-saturated (8.0 weight %) andesitic melt–mafic rock complex with melt volume fractions greater than 4 volume %, while the shear wave velocity requires the fluid and/or melt volume fraction of 6–7 volume % at that depth. Considering the fluid and/or melt volume fraction of 6–7 volume %, the conductive hydrous phase is likewise required to explain the wide range of the 67% confidence interval of resistivity. The Mogi inflation source determined from geodetic data lies on the resistive side near the top boundary of the conductive magma reservoir at a depth of 2.7 or 3.7 km bsl. Assuming that the resistivity of the inflation source region is above the upper bound of the confidence interval of resistivity for the conductive magma reservoir and that the source region is composed of hydrothermal fluid + rock, the resistivity of the source region is explained by a hydrothermal fluid volume fraction below 5 volume %, which is the percolation threshold porosity in an effusive eruption. This indicates that the percolation threshold characterizes the inflation source region.


Introduction
Great earthquakes have often triggered eruptions in neighboring volcanoes (Bebbington and Marzocchi 2011;Koyama 2015;Nishimura 2017). After the 2011 Tohoku-Oki Earthquake, seismicities of deep low-frequency tremor and of volcano-tectonic earthquake have turned active beneath a few volcanoes in northeastern (NE) Japan (e.g., Yamamoto 2014; Ikegaya and Yamamoto 2020). While the potential eruptions of these volcanoes have raised concerns, Azumayama Volcano is one of the most concerning of the active volcanoes in NE Japan triggered by the 2011 Tohoku-Oki Earthquake. Because this volcano is close (~ 250 km) to the hypocenter of the 2011 Tohoku-Oki Earthquake, and sporadic seismic swarms have been observed around this volcano since the year 2000. Furthermore, the Global Navigation Satellite System (GNSS) baseline extensions and tilt changes have been reported accompanied those seismic swarms (Yoshida et al. 2012;Kawanabe and Ueki 2013;Japan Meteorological Agency 2019). Although outstanding seismic swarms were observed during the years 2014-2015 and 2018-2019 (Japan Meteorological Agency 2019; Yoshigai et al. 2019), no eruption has taken place so far. Azumayama is located in the Nasu volcanic chain of the Quaternary volcanic front in NE Japan (Yamamoto 2005;Takahashi et al. 2013), and is composed of several andesitic edifices that have emerged since the Pleistocene ( Fig. 1). Azumayama has active fumaroles in Oana Crater in the eastern part with recent maximum effusion rates of 600 tons/day. The recent phreatic eruptions occurred in Oana Crater in the year 1977, 1952, and 1950(Kawanabe and Ueki 2013. For the eruption forecast, shallow depths beneath Azumayama, from the surface to sea level, are well investigated (e.g., Yoshida et al. 2012) (Altitudes around the study area are 600-2000 m, and those of the edifice summits are found in Fig. 1 and its caption. The depths below the surface are shown using the altitude of Oana Crater, 1720 m, as a reference altitude, hereafter). However, the transport and distribution of magma/hydrothermal fluid in the depth range from 0 to 20 km below sea level (bsl), or 1.7 to 21.7 km below surface (bsf ), are unknown. This depth range coincides with a seismicity gap: the volcanotectonic earthquakes are primarily restricted to depths less than 1 km bsl (2.7 km bsf ), and deep low-frequency tremors occur at depths greater than 20 km bsl (21.7 km bsf ) (Kawanabe and Ueki 2013). A magma reservoir beneath active volcanic fields often lies at depths of a few to 20 km bsf (e.g., Chaussard  . The contour intervals on the surface topography are labeled in meters above sea level Ichiki et al. Earth, Planets and Space (2021) Cashman et al. 2017;Pritchard et al. 2018), and a magma reservoir beneath Azumayama is expected to lie within this depth range. Takada and Fukushima (2013) model the low-viscosity ellipse in a depth range from 1.26 to 14.4 km bsl (2.98 to 16.1 km bsf ) beneath Azumayama from analysis of the interferometric synthetic aperture radar (InSAR) data recording the subsidence of the edifice caused by the 2011 Tohoku-Oki Earthquake. This is the first geophysical study that represents existing magma reservoir within this depth range beneath Azumayama. The slow shear wave velocity (V s ) region in a depth range from 0 to 20 km bsl (1.7 to 21.7 km bsf ) around Azumayama is imaged by using seismic body wave tomography (Okada et al. 2015) and ambient noise tomography (Chen et al. 2018). While these studies provide broad and blurred images of the magma reservoir (the resolution scale is 10 to 20 km), the inflation source location presumably gives a more precise boundary of the magma reservoir's outline. Japan Meteorological Agency (2019) and Seki et al. (submitted to Quart. J. Seismology) show superposition of two-or three-Mogi inflation sources explains the ground deformations around Azumayama observed by GNSS and the tiltmeters (see Fig. 15). One Mogi inflation source solution of the two-or three-Mogi source model is located at a depth of 3.7 or 2.7 km bsl (5.4 or 4.4 km bsf ) beneath Oana Crater, and the location may indicate a margin of the magma reservoir. Electrical resistivity (the reciprocal of conductivity) models beneath volcanic fields using magnetotellurics (MT) have succeeded in exposing clear images of magma reservoirs (Hill et al. 2009;Heise et al. 2010;Bertrand et al. 2012;Aizawa et al. 2014;Comeau et al. 2016;Hata et al. 2018;Cordell et al. 2018). The primary purpose of this study is to model the three-dimensional resistivity distribution down to a depth of 20 km bsl, in order to provide a clear image of the magma reservoir beneath Azumayama. Along with defining the magma reservoir, we discuss the petrological interpretation for the magma reservoir and the threshold of a phreatic eruption of Azumayama.
Melt volume fraction and melt water content of magma reservoirs have been estimated from resistivity models based on laboratory data of melt and rock resistivity (McGary et al. 2014;Ichiki et al. 2015;Heise et al. 2016;Laumonier et al. 2015Laumonier et al. , 2017Guo et al. 2017). To gain a more reliable quantitative interpretation of the magma reservoir, we perform checker board resolution tests (CRTs) following seismic tomographic studies (e.g., Inoue et al. 1990), and estimate a confidence interval of resistivity (CIR) of the magma reservoir. This contributes to a more rigorous and quantitative interpretation of the magma reservoir in a petrological context based on laboratory data.
Ascent of hydrothermal fluids exsolved from melts in a magma reservoir could trigger a phreatic eruption. If we hypothesize such a concept, a key parameter controlling the eruption triggering is permeability along the hydrothermal fluid pathway. At shallow depths less than several hundred meters bsf, clay minerals likely constitute the primary impermeable layers suppressing the ascent of hydrothermal fluids. The impermeable clay layers are imaged by MT as conductive bodies (Komori et al. 2013;Tsukamoto et al. 2018;Yoshimura et al. 2018;Matsunaga et al. 2020;Tseng et al. 2020). On the other hand, densified magma on the pathway probably plays the impermeable role at higher temperatures and greater depths (Okumura and Sasaki 2014;Heap et al. 2015;Kushnir et al. 2016;Colombier et al. 2017). This study assumes the latter concept in which the densified magma interrupts the ascent of hydrothermal fluids. The reduction of porosity and the low hydrothermal fluid volume fraction resulting from the densification process is expected to result in a more resistive image than the conductive hydrothermal fluids underlying the densified magma. Although the theoretical relationship between the Kozeny-Carman equation (permeability) and Archie's Law (resistivity) in porous media suggests that resistivity can be converted into permeability, simple applications of the relationship are usually misleading (Wright et al. 2009). Thus, we estimate porosity (hydrothermal fluid volume fraction), rather than permeability, of the impermeable region imaged as being resistive in accordance with another relationship between resistivity and porosity, Hashin-Shtrikman bound (Hashin and Shtrikman 1962). On the basis of the estimated porosity of the impermeable region, we test which index characterizes the impermeable region, percolation threshold (e.g., Blower 2001;Burgisser et al. 2017) or "change point" in porosity-permeability relationship (Heap et al. 2014Farquharson et al. 2015;Kushnir et al. 2016). The index is expected to be used as a potential diagnostic threshold for the forecast of a phreatic eruption of Azumayama.

Methods, data, and analyses
Basic concept of method, data and data processing The details of MT and geomagnetic depth sounding (GDS) used in this study are found in Chave and Jones (2012), and we give only a brief introduction to those methods. The MT and GDS are resistivity exploration methods that use the frequency response of electromagnetic (EM) field variation observed at the Earth's surface. The frequency response varies with EM induction in the Earth induced by external EM field variation. The amplitude and phase of the frequency response reflect the subsurface electrical resistivity. Depth of sensitivity is related to period or frequency of EM data, with data at longer periods providing information about greater depths. MT and GDS frequency responses, which are impedance tensor and geomagnetic transfer functions, are defined as follows: is the Fourier spectrum of electric and magnetic field variation; subscript x, y, z is north-south (N-S), eastwest (E-W), down-up component; (f) is the function of frequency, superscript obs means observed data. Z obs ij , A, and B are complex number. Impedance tensor has following nature depending on the dimensionality of subsurface resistivity distribution: One-dimensional (1-D): Z obs xx f = Z obs yy f = 0 and Z obs xy f = − Z obs yx f . Two-dimensional (2-D): Z obs xx f = Z obs yy f = 0 , when x or y axis coincides with the strike.
Three-dimensional (3-D): All four elements are usually non zero.
The 1-D situation means that the subsurface resistivity distribution varies only in the z-direction. The 2-D situation means the resistivity distribution does not vary along one horizontal direction, the strike direction. The 3-D situation means the subsurface resistivity is heterogeneous in any direction. However, the impedance tensor generally contains galvanic distortion, a kind of observational site effect, by which direct interpretation of raw observed impedance is misleading (e.g., Jiracek 1990). The galvanic distortion is expressed as: is the galvanic distortion tensor. Note that C ij is real and independent of frequency. Z ind ij (i, j = x, y) means distortion-free or "regional" impedance element. The distortion tensor cannot be fully corrected without performing additional observations. Some researchers impose the correction by assuming spatial smoothness or other constraints on the distortion tensor (Avdeeva et al. 2015). To cope with the galvanic distortion, Caldwell et al. (2004) Z obs xx f Z obs xy f Z obs yx f Z obs yy f where superscript −1 is the inverse matrix; i is the imaginary number; X ij (i, j = x, y), Y ij (i, j = x, y) are real and imaginary part of each impedance tensor element, respectively. The phase tensor indicates the degree of subsurface resistivity heterogeneity (dimensionality). The phase tensor is diagonalized to determine its principal axis using the horizontal axis rotation matrix, R: Phase tensor ellipticity, Λ, is defined as: The phase tensor is graphically represented by an ellipse (Heise et al. 2006). The tensor principal axis represents eigenvalue and eigenvector, that is, the tensor character changes only in amplitude along the principal axis direction. Thus, the principal axis of the phase tensor ellipse means that phase tensor element information in this direction is not contaminated with information from other directions. Accordingly, the phase tensor ellipse is represented by a circle ( Φ max = Φ min ) in a 1-D situation, an elliptical shape ( Φ max ≠ Φ min ) and β = 0 in a 2-D situation, and an ellipse and β ≠ 0 in a 3-D situation (Caldwell et al. 2004). The arc tangent of the square root of the phase tensor determinant is considered to be the distortion-free equivalent 1-D MT phase response around an observation site (Heise et al. 2008). The square root of phase tensor determinant is denoted as Φ 2 : Z obs xx f Z obs xy f Z obs yx f Z obs yy f GDS data are graphically represented by an induction vector. Parkinson vector, which is the most commonly used induction vector, points towards the conductor: where e N , e E represent northward and eastward unit vector, respectively (e.g., Hobbs 1992).
The MT and GDS data were acquired from 18 sites at Azumayama (Fig. 1). The study area is about 18 km N-S by 14 km E-W, with an average spacing between the sites of 3.00 km. The EM field variation data were acquired with an MTU-5 system (Phoenix Geophysics) at two sites during four days in September 2015, and at the other 16 sites in September 2015 and 2017 using an ADU-07e system (Metronix Geophysics). All sites recorded fivecomponent EM field variations. The two-component horizontal telluric field was observed using Pb-PbCl 2 electrodes with 40 m dipole lengths. The three-component magnetic field for remote-reference processing was recorded at Okura village, Yamagata prefecture (about 97 km north of Mt. Issaikyo) using the ADU-07 system, and at Nishiwaga town, Iwate prefecture (about 210 km NNE of Mt. Issaikyo) using the MTU-5A system. The MT and GDS frequency responses of the two sites were calculated from 0.00055 to 320 Hz using SSMT2K and MTEDITOR (Phoenix Geophysics), and for the other 16 sites from 0.00012 to 384 Hz using BIRRP code (Chave and Thomson 2004).

Three-dimensional resistivity modeling
The WSINV3DMT inversion code (Siripunvaraporn et al. 2005;Siripunvaraporn and Egbert 2009) is used to model the distribution of 3-D resistivity beneath Azumayama. The WSINV3DMT code seeks an optimal model by minimizing the following objective function, W, using the Occam inversion algorithm (Constable et al. 1987): where m is the model parameter vector, m 0 is the prior model parameter vector (fixed in iteration), C m is the model covariance matrix, d is the data vector, C d is the data covariance matrix, F[] is the 3-D forward operator, λ is the Lagrange multiplier. The first and second terms are regularization and data misfit, respectively. The m is sought iteratively, and the model is updated in each iteration by determining the stationary point of Eq. 13. The formulation is a function of λ, and the λ is determined so the normalized root-mean-squared (11) Quadrature-phase Parkinson vector: Im(A)e N + Im(B)e E , (RMS) misfit (Eq. 14 as described later) takes a minimum value (Siripunvaraporn et al. 2005). The inversion data are real and imaginary parts of the MT impedance and GDS response at 18 periods, which have octaves at 1024, 512, 256, 128, …, 1.5625 × 10 -2 , and 1.0417 × 10 -2 s. The number of data is 3888 (18 sites × 18 periods × 12 responses). We give 5% of each MT impedance component to the error floor of respective MT impedance component, and 10% of each GDS response to the error floor of respective GDS response. The modeling space is 420 km × 417 km × 150 km, which is divided into 57 × 59 × 64 blocks in N-S, E-W and vertical directions. Resistivity of blocks in the marine area is fixed at 0.3333 Ωm and the SRTM30_PLUS model (Becker et al. 2009) is used to distinguish seafloor. Land topography is incorporated into the model, and resistivity of air blocks created by the land topography is fixed at 10 12 Ωm. Each element of the model covariance matrix, C m , is imposed to follow a 1-D diffusion equation solution in each N-S, E-W and vertical direction (Siripunvaraporn et al. 2005). The relaxation parameter 4τδ (Siripunvaraporn and Egbert 2000) is fixed at 2.0 in all three directions. One inversion operation is over when the RMS misfit achieves a given threshold, or when the model iteration repeats a given maximum number of iterations. The RMS misfit of all data is defined as: where N is number of data. The RMS misfit threshold and the maximum number of iterations are set to 1.0 and 10, respectively.
We start inversion using the initial prior model, where each marine, air, and subsurface region is uniform, respectively. We test the four initial prior models, of which the subsurface resistivities are 100, 300, 1000, and 3000 Ωm. For each initial prior model, we perform the inversion in the following procedure: (a) 1. When the RMS misfit achieves the RMS misfit threshold, 1.0, at some number of iterations (≤ the maximum number of iterations, 10), the inversion operation is ended, and the solution is accepted as a final model.
2. When the inversion achieves the maximum number of iterations, 10, without achieving the RMS misfit threshold, 1.0, we replace the initial prior model into the best solution (giving a minimum RMS misfit) from the 10 iteration solutions . We denote this operation as prior model update. After the prior model update, we proceed to (b). (14) We continue the inversion operation using the new prior model, and then 1. If the RMS misfit achieves 1.0, we deal with like (a) 1. 2. When the inversion achieves the maximum number of iterations, 10, without achieving the RMS misfit threshold, 1.0, we assess whether the best solution selected from the 10 iteration solutions significantly decreases RMS misfit compared with the RMS misfit of the best model in (a) 2, viz. the new prior model. When the RMS misfit is not significantly improved, the inversion operation is ended, and the solution is accepted as a final model. When the RMS misfit is significantly improved, we again perform a prior model update, and return to the top of (b).
To judge whether the RMS misfit significantly decreases or not, we employ the F-test. The essence of the F-test is the test of variance inequality between two populations. The test statistic F, defined below, follows an F-distribution with the degree of freedom (DOF) N and N (e.g., Chase and Brown 1986, p. 309, 330): where subscript i and i + 1 means ith and (i + 1)th prior model update. The χ 2 -misfit in Eq. 15 is defined as: Here, we use sample variance for simplicity. The F-statistic threshold is set to 1.054, which is the 95% confidence level of the left-tailed (or 5% significance level of the right-tailed) F-test with DOF 3888 and 3888 (Chase and Brown 1986, p. 333).
The best final model is selected from the four final models, each of which is originated from each of the four initial prior models.
To show the data misfit, we calculate phase tensor residual representing differences between two phase tensors. Two kinds of phase tensor residual have been proposed: one by Peacock et al. (2012) and the other by Heise et al. (2007Heise et al. ( , 2008, and this study uses the phase tensor residual of Peacock et al (2012). The phase tensor residual of Peacock et al. (2012), ΔΦ Peacock , and the misfit of the arithmetic mean of two principal axis magnitudes, ΔΦ principal are: where I is the 2 × 2 identity matrix, Φ obs is the observed phase tensor, Φ cal is the phase tensor calculated from the final resistivity model, Φ obs max , Φ obs min : two principal axis magnitudes of Φ obs , Φ cal max , Φ cal min : two principal axis magnitudes of Φ cal , δΦ max , δΦ min : standard errors of two principal axis magnitudes of Φ obs . Phase tensor residual is graphically drawn in the same manner as a phase tensor. However, when an observed phase tensor has a large ellipticity, a large maximum principal axis often remains in both the Peacock et al. (2012)'s and Heise et al. (2007Heise et al. ( , 2008)'s phase tensor residual ellipses, even if the data fit is fairly good. The reason we use Peacock et al. (2012)'s phase tensor residual ellipse is the straightforward nature, that is, the maximum principal axis azimuth approaches α obs + β obs + π/2 in a large ellipticity case as data fit gets better, while Heise et al. (2007Heise et al. ( , 2008)'s phase tensor residual ellipse has a rather complex nature. We explain this in Additional file 1. Figure 2 shows the observed tan −1 (Φ 2 ) at various frequencies or periods. The tan −1 (Φ 2 ) in the frequency range higher than 4 Hz has higher values in the area east of Oana Crater. The high tan −1 (Φ 2 ) area assumes E-W elongation towards Oana Crater. The tan −1 (Φ 2 ) distribution typically shows N-S trends in the period range longer than 4.0 s, and the tan −1 (Φ 2 ) greater than 45° dominates in the western part of the area. The tan −1 (Φ 2 ) falls below 45° across the study area at a period of 512 s. These findings are consistent with the final resistivity model. Figure 3 shows the RMS misfit in terms of prior model updates. The RMS misfits in all inversion operations do not achieve 1.0. We obtain the final resistivity model at the third prior model update starting with the initial prior model of 1000 Ωm subsurface. The best RMS misfit is 2.42. Figure 4 shows the RMS misfit of each site calculated from the final resistivity model. Figures 5, 6, 7 represent the results of the final resistivity model in plan or map view, as a series of vertical profiles in the E-W direction, and as a similar series of vertical profiles in the N-S direction. Figure 5 indicates no conductor at the depth of 0.05 km bsl beneath Oana Crater. Note that the altitude of Oana Crater is 1720 m. At the depth of 3.5 km bsl, the E-W elongated conductor

Results
of each site at various frequencies or periods (colored circles). The white contour lines represent with 2.5° width. The black and red triangles are the same as in Fig. 1 (label A) appears beneath Oana Crater. The E-W elongated conductive feature extends down to the depth of 5.5 km bsl. The structure is delineated by the convex shape of the conductor in the depth range from 3.5 to 5.5 km bsl in the C-C' profile of Fig. 6 (label A). The Φ 2 response distribution in the high-frequency range (> 4 Hz) of the final resistivity model also reveals this feature (Fig. 2). The conductor beneath Oana Crater expands in the N-S direction at depths greater than 6 km bsl. The plan view at the depth of 7.5 km bsl in Fig. 5 and the H-H' profile in Fig. 7 provide clear images of the expansion (label A). The Φ 2 distribution also shows this N-S trend at the periods of 4.0 and 8.0 s (Fig. 2). The area in the depth range from 7 to 9 km bsl beneath Oana Crater has the lowest resistivity value (~ 0.1 Ωm) in the final model, which also shows a second conductive block in the depth range from 5 to 10 km bsl beneath Naka-azuma. The twin conductor image can be found in the plan views at the depth of 7.5 and 9.5 km bsl of Fig. 5 (labels A and B) and in the D-D' profile of Fig. 6 (labels A and B). No apparent conductor is found below the observation sites at depths greater than 20 km bsl. Figure 8 shows the observed and calculated phase tensor, phase tensor residual ellipses. Figure 9 shows the observed, calculated and difference in-phase Parkinson vector. The phase tensor ellipses across the study area have a large ellipticity in the observed period bands, and the in-phase Parkinson vectors have scattered directions at periods less than 8 s. These features indicate the 3-D heterogeneous subsurface in the study area. Although both the phase tensor ellipses and in-phase Parkinson vectors are complicated, we can see their relationship with the resistivity model. Viewing from long to short periods, the ellipse directions of some northern and southern sites are aligned in the N-S direction at the periods of 8 and 4 s. This is due to the current pointing to the conductor A at the depths of 9.5 and 12.5 km bsl (Fig. 5). The in-phase Parkinson vectors at the periods of 8 and 4 s also generally point towards the conductor A at depths of 9.5 and 12 km bsl. As depth decreases, the conductor A moves to the center (beneath Oana Crater) at the depths of 5.5 and 7.5 km bsl. Corresponding to this, the in-phase Parkinson vectors at the periods or frequencies of 2 s to 4 Hz point towards the central part more than those at the periods of 8 and 4 s. The direction of the phase tensor ellipses at 2 s to 4 Hz assume more radial pattern.
Because the phase tensor ellipticity is so great, some large maximum principal axes remain for the phase tensor residuals. However, over half of these align to the azimuth of α obs + β obs + π/2. The data fits are reasonable in the middle period range between 64 s and 32 Hz, while the phase tensor residual ellipses show a poor data fit at the 512 s period or longer. The in-phase Parkinson vectors show generally good fittings. All data fittings of MT apparent resistivity, phase responses, and Parkinson vectors are presented in Additional file 2.

Resolution and confidence interval of resistivity
The better data misfit in the middle period or frequency range from 8 Hz to 8 s indicates that the conductive region in a depth range from 3 to 15 km bsl is well resolved. Focusing on the conductor, we perform CRTs to investigate the resolution of the resistivity model, and estimate the CIR of the conductor.  3.5, 5.5, 7.5, 9.5, 12.5 km bsl. The background color represents the logarithmic resistivity value. The black circles, white lines, red character, and black triangles, are observation sites, topographic contours, Oana Crater, and edifice summits, respectively (see Fig. 1). The A and B for the conductors are referred to in the text Ichiki et al. Earth, Planets and Space (2021)  Figures 10, 11, 12 show the two checkerboard models used to test whether the final model has 3-or 5-km resolution. The checkerboard is a single layer, and background and checkerbox resistivity is 1000 and 1 Ωm, respectively. The marine and air blocks are fixed, as described in "Three-dimensional resistivity modeling". Synthetic The observed and calculated phase tensor ellipse magnitudes are normalized with each observed maximum principal axis magnitude at each site and period/frequency, and are filled with colors corresponding to tan −1 (Φ 2 ) designated by the color scale. The phase tensor residual uses the definition by Peacock et al. (2012), and these are filled with colors corresponding to Δ Φ principal (Eq. 18) designated by the color scale. Azimuth of the white bars represents α obs + β obs + π/2 direction (length is arbitrary) datasets are generated of the respective components at each site that are composed of MT impedance and geomagnetic transfer functions with 5% amplitude Gaussian noises. Because a substantial amount of time and effort are required to demonstrate that the CRTs replicate the precise Occam process of the final model, the Lagrange multiplier, λ in Eq. 13, is fixed constant in the CRTs. The λ takes 0.1 at iteration 9 on the initial prior model, 0.1 at iteration 8 on the first update prior model, 0.032 at iteration 5 on the second update, and 10 at iteration 2 on the third update, until we achieve the final model (viz. blue plots in Fig. 3). Here, we fix λ at 0.1 in the CRTs. Figure 10 shows the result of CRTs using checkerboxes with dimensions of about 3 km to a side. Outputs using each prior model show blurred checkerboxes, and the number of conductive checkerboxes cannot be properly determined. The northeastern most output conductive checkerbox moves horizontally from the original position, and conductive checkerbox resistivities are overestimated to be 5-10 Ωm or more. Figures 11 and 12 show the results from using checkerboxes with dimensions of about 5 km to a side. While the conductive checkerboxes at depth inside the observation network are blurred, the five conductive checkerboxes could be distinguished. Moreover, the conductive checkerbox resistivities are properly estimated to be 1 Ωm at the third update of the prior model. The eight conductive checkerboxes at depth outside the observation network are hardly recognizable. These two CRTs reveal that the conductor in the depth range from 5 to 10 km bsl has no resolution with the scale of less than 3 km, but reasonable resolution with the scale of over 5 km.
Next, we estimate the CIR of the conductor. We assess the CIR of the area surrounded by a 1-Ωm isosurface in the depth range from 3 to 15 km bsl in the final model, and that surrounded by a 3-Ωm isosurface (Fig. 13). The labels A and B in Fig. 13 correspond to the conductive features in Figs. 5, 6, 7. The horizontal dimensions of the defined area are 10 × 5 km and 20 × 15 km, respectively. The target area's resistivity is changed to a variety of constant values between 0.01 and 100 Ωm, and the data misfit change of each site is obtained by forward modeling using those modified models. Plan view of CRT as in Fig. 10, but using checkerboxes about 5 km to a side. The checkerboxes are embedded in the depth range from 5 to 10 km bsl, and the output models are shown at the depth of 7.5 km bsl Ichiki et al. Earth, Planets and Space (2021) 73:150 The data misfit change of each site is shown by using the F-statistics defined as Eq. 15. The RMS misfit in the denominator of Eq. 15 is fixed at the RMS misfit of each site of the final resistivity model (Fig. 4). Figure 14 shows the F-statistics for each of the 18 sites. Note that total RMS misfit using all sites show no significant change in terms of the target area's resistivity. The 95, 90, and 67% confidence levels of left-tailed F-test with DOF 216 and 216 (18 periods × 12 responses per site) are 1.252, 1.191, and 1.062, respectively. The CIR of the area under 1 Ωm cannot be constrained using any of the 95, 90 and 67% confidence levels (Fig. 14a). This is presumably due to Vertical profile in the E-W direction of CRT input and output models. The left and right columns represent the models using checkerboxes with one side 3 and 5 km, respectively the insufficient volume of the target area. The integral equation formulation for the 3-D EM induction problem states that the product of conductivity perturbation and volume, and distance between conductivity perturbation (source) and a site (receiver) primarily controls the EM field and impedance response changes at the site (Hohmann 1975;Wannamaker et al. 1984;Zhdanov et al. 2006). On the other hand, the 90 and 67% CIR of the under 3 Ωm area are estimated to be 0.02 ≦ (90% CIR) ≦ 70 Ωm and 0.3 ≦ (67% CIR) ≦ 5 Ωm (Fig. 14b). The RMS misfits of some sites get significantly worse as the target area's resistivity decreases. The same behavior also appears as the target area's resistivity increases. This nature enables both the upper and lower bounds of CIR to be constrained.

Comparisons with other geophysical results and other active volcanoes
The conductive region in the depth range from 3 to 15 km bsl is consistent with the slow compressive (V p ) and shear (V s ) wave velocity regions by seismic body wave and ambient noise tomographies (Okada et al. 2015;Chen et al. 2018). Figures 15, 16, 17 show the V s distribution by Okada et al. (2015) superimposed on the resistivity model. The V p and V p /V s comparisons with the resistivity model are shown in Additional file 3. The CRT of the seismic tomography indicates the seismic structure resolution is 10 to 20 km (Okada et al. 2015). Although the resolutions of seismic and resistivity models are different, the minimal V s , corresponding to the conductive area beneath Oana Crater, is found to be 3.2 to 3.25 km/s at the depth of 6 km bsl (7.7 km bsf ), and 3.35 to 3.45 km/s at the depth of 12 km bsl (13.7 km bsf ). Figure 18 shows a comparison of the resistivity model and the low-viscosity ellipsoid body inferred from the InSAR data recording   Fig. 13. a F-statistic change obtained by changing resistivity inside isosurface in Fig. 13a. b Same with (a), but obtained by changing resistivity inside Fig. 13b

Fig. 15
Plan or map view of the V s model by Okada et al. (2015) superimposed on the resistivity model (Fig. 5). The white lines are V s contours with the interval of 0.05 km/s. The background color represents the logarithmic resistivity value. The black and white circles, red character, and black triangles, are observation sites, sampling point of V s model, Oana Crater, and edifice summits, respectively. The labels A and B for the conductors are referred to in "Results" section Ichiki et al. Earth, Planets and Space (2021)

Fig. 16
Vertical profiles of the V s model by Okada et al. (2015) superimposed on the resistivity model (Fig. 6) in the E-W direction. The white lines and circles are V s contours and sampling point of the V s model. Other features are the same as those in Fig. 6 subsidence around Azumayama caused by the 2011 Tohoku-Oki Earthquake (Takada and Fukushima 2013). The low-viscosity ellipsoid encompasses almost all the conductive region less than 3 Ωm. The resistivity model gives a more focused image of the magma reservoir than the slow velocity and low-viscosity models.  Figure 18 shows the Mogi model inflation source models superimposed on the final resistivity model. The inflation sources of the three-source model at a depth of 2.7 km bsl (4.4 km bsf ) and of the two-source model at a depth of 3.7 km bsl (5.4 km bsf ) lie on the resistive side near the top boundary of the conductor. The resistivity around the inflation source is 30-300 Ωm. Although it is hard to evaluate the CIR of the inflation source region, the lower bound of the CIR could be taken as the upper bound of the CIR of the conductor estimated in the previous section. On the other hand, the source at the depth of 9 km bsl (10.7 km bsf ) from the Mogi inflation source model is located inside the conductor. However, this Mogi source location is determined with an accuracy of a few km. The depth range of the conductor beneath Oana Crater is about 4.7 to 16.7 km bsf (3-15 km bsl). Similar conductors beneath other active volcanic fields are found in the same depth range from 3 to 20 km bsf [e.g.,  (2018)]. Melt inclusion samples from active volcanoes imply that a magma reservoir or pre-eruptive magma storage could lie in the depth range from a few to 13 km bsf (Cashman et al., 2017). According to the compilation by Chaussard and Amelung (2014), magma reservoir depths in arc volcanoes are typically a few to 10 km bsf, and correlate with crustal thickness and age. Thus, the conductor beneath Oana Crater can be interpreted as the magma reservoir from the viewpoint of magma reservoir depth. However, MT explorations have also revealed magma reservoirs or pre-eruptive magma storage lying outside of this depth range. Volcán Uturuncu and the Bolivian Altiplano have conductors corresponding to pre-eruptive magma storage at depths less than 5 km bsf (0 km bsl), and outstanding large conductors at depths greater than 20 km bsf (15 km bsl) (Comeau et al. 2016). The Ethiopian rift resistivity model (Hübert et al. 2018) indicates an enormous conductor in a depth range from 10 to 30 km bsf (a width of about 20 km). Another extreme example is the resistivity model beneath the base of Mt. Etna (Siniscalchi et al. 2012), which has a broadly resistive zone (> 100 Ωm) at depths greater than 4 km bsf. The variety of these conductors in volcanic and active fields can be explained by the trans-crustal magmatic system (Cashman et al. 2017).

Melt volume fraction of the magma reservoir
We interpret that the shallower part of the conductor or the magma reservoir is composed of dacitic meltrock-hydrothermal fluid complex. The phenocrysts in the andesite and basaltic andesite samples of Oana unit are divided into mafic, mixed magma-derived, and felsic silicates (Ban et al. 2016). The felsic phenocrysts have the shallowest origin, which is dacitic melt with water contents of 2.75-3.25 weight % (wt %) at the depth of about 4 km bsl (5.7 km bsf ) by Rhyolite-MELTS analysis, and at 880-890 ℃ by pyroxene geothermometry (Ban et al. 2016). The depth of 4 km bsl (5.7 km bsf ) is near the top of the conductor. Hence, we first assume that the shallower part of the conductor is composed of two phases of dacitic melt and rock. The melt water content can range up to the water-saturated condition of 5.5 wt % (Sisson and Bacon 1999;Wallace 2005). Based on recent laboratory data of dacitic melt resistivity and modified Archie's Law with m = 1.05 (Laumonier et al. 2015), we calculate bulk resistivity for the dacitic melt and rock complex at 900 °C and 4 km bsl (5.7 km bsf; 0.15 GPa). The pressure is calculated to be the lithostatic pressure assuming 2750 kg/m 3 as the average density, and 9.8 m/s 2 as the gravity acceleration. Figure 19 compares the bulk resistivity and the 67% CIR (0.3 to 5 Ωm) at various water contents. While the upper bound of 67% CIR corresponds to the complex with a melt volume fraction greater than 0.22 at a water content less than 5.5 wt %, the complex with any reasonable melt and water content cannot explain the lower bound of 67% CIR. Hydrothermal fluids and/ or mafic melts are required to explain the lower bound of 67% CIR (see also Laumonier et al. 2015). Regarding the shallower part of the magma reservoir where the magma differentiation is mature, hydrothermal fluids rather than mafic melts are a probable interpretation for this discrepancy. The minimal V s is found to be 3.2-3.25 km/s corresponding to the shallower part of the conductor (Fig. 15). V s of rock + fluid and/or melt complex is expressed as (Dvorkin 2008): where φ is the fluid and/or melt volume fraction, V rock is the V s of rock velocity, ρ rock is the density of rock, ρ f is the density of fluid and/or melt. Using ρ rock = 2750 kg/m 3 , V rock = 3.5 km/s (see Comeau et al. 2016 (Hack and Thompson 2011), the fluid volume fraction is estimated to be the same, that is 0.03 from a V s of 3.25 km/s. Therefore, the resistivity and V s estimates can be reconciled if the hydrothermal fluid volume fraction is more dominant than the melt volume fraction. Magma mixing and differentiation probably occur in the magma reservoir (e.g., Toya et al. 2005;Kimura and Yoshida 2006;Tatsumi et al. 2008;Takahashi et al. 2013;Yoshida et al. 2014), and result in the magma reservoir becoming increasingly more mafic with depth. Thus, we consider more mafic conditions, an andesitic melt and rock complex at a greater depth inside the conductor. Figure 20 shows the bulk resistivity of the complex at 900 °C, and 9 km bsl (10.7 km bsf; 0.29 GPa), with the melt water contents up to 8.0 wt % water saturation (Sisson and Bacon 1999;Wallace 2005). The modified Archie's Law with m = 1.05 is used with andesitic melt resistivity models by Guo et al. (2017) and Laumonier et al. (2017). The upper bound of 67% CIR is explained by the complex with a melt volume fraction greater than 0.04 at a water content less than 8.0 wt %. The lower boundary corresponds to the complex with a melt volume fraction of 0.61 (Guo et al. (2017)'s model) or 0.94 (Laumonier et al. (2017)'s model) at almost the water saturation, 8.0 wt %. On the other hand, assuming V rock = 4.0, and using ρ rock = 2750 kg/m 3 and ρ f =2400 kg/ m 3 for andesitic melt (Comeau et al. 2016), the melt volume fraction is estimated to be 0.06 to 0.07 from V s of 3.35 to 3.45 km/s (Fig. 15) Fig. 14b, and the solid lines show water content of dacitic melt in weight percent. Note that the upper and lower bound of 67% CIR is the bottom and top of the gray shaded area, respectively. Melt resistivity is calculated from Laumonier et al. (2015). Solid phase resistivity is fixed at 1000 Ωm [e.g., granite data of Kariya and Shankland (1983)] probably required even for the deep part of the magma reservoir (cf., Laumonier et al. 2017).

Interruption of hydrothermal fluid ascent
The shallower part of the magma reservoir is likely composed of the dacitic melt-rock-hydrothermal fluid complex. We interpret that an impermeable roof prevents the hydrothermal fluids from ascending and causing a phreatic eruption. The Mogi inflation source at a depth of 2.7 or 3.7 km bsl (4.4 or 5.4 km bsf ) near the top of the conductor indicates the potential location of the interruption. The impermeable roof results from porosity reduction through magma densification (Okumura and Sasaki 2014;Heap et al. 2015). The densification has presumably reached the percolation threshold in porosity (e.g., Blower 2001;Burgisser et al. 2017) or the "change point" porosity in the permeability-porosity relationship (Heap et al. 2014Farquharson et al. 2015;Kushnir et al. 2016). We test which index gives a consistent explanation for the resistivity of the impermeable region. The change point porosity is reported to be 10 to 15 volume % (vol %) (Heap et al. 2014Farquharson et al. 2015;Kushnir et al. 2016). On the other hand, the percolation threshold varies depending on the conditions such as the transition between explosive and effusive eruptions, the inherent hysteresis between emerging and converging processes in eruptions, and others (e.g., Rust and Cashman 2004;Mueller et al. 2005;Michaut et al. 2009;Rust and Cashman 2011;Burgisser et al. 2017;Colombier et al. 2017). We consider the general percolation threshold for effusive eruption, 5 vol %, and for explosive eruption, 30 vol % (Mueller et al. 2005;Colombier et al. 2017). As mentioned in "Comparisons with other geophysical results and other active volcanoes" section, we take the lower bound of the CIR for the Mogi inflation source region as the upper bound of the CIR for the magma reservoir's conductor. Figure 21 shows the bulk resistivity of two phases (hydrothermal fluid-rock) as a function of the hydrothermal fluid volume fraction. The Hashin-Shtrikman upper bound in conductivity (Hashin and Shtrikman 1962 Guo et al. (2017) are shown by the broken line, and from Laumonier et al. (2017) are shown by the solid line. Solid phase resistivity is fixed at 500 Ωm [e.g., rhyolite data of Kariya and Shankland (1983)]. Others are the same as Fig. 19 H y d .F lu id @ 0 .  Hydrothermal fluid resistivity is assumed to be 0.1 or 0.3 Ωm. Rock resistivity is fixed at 1000 Ωm (see Fig. 19). The green, orange, and pink columns show percolation threshold porosity in effusive and explosive eruptions, and change point of permeability-porosity relationship, respectively. The red and green lines show the upper bounds of 67% and 90% CIR of the conductor (see Fig. 14b). The oblique pattern represents likely resistivity range at the Mogi inflation source location, assuming that the upper bound of CIR of the conductor restricts the lower bound of resistivity at the Mogi inflation source in modified Archie's Law on this two phase complex [m should be determined depending on phase combination, see Glover et al. (2000), Glover (2010)]. The hydrothermal fluid is assumed to be supercritical with a resistivity of 0.3 to 0.1 Ωm at conditions of 750-900 °C, 0.13-0.15 GPa, and nearly saturated with CO 2 exsolution (Marshall and Frantz 1987;Nesbitt 1993). Considering the lower bound of 90% CIR, that is, the Mogi inflation source region has resistivities greater than 70 Ωm, the hydrothermal fluid volume fraction around the Mogi inflation source is estimated to be less than 5 vol %, the percolation threshold porosity for effusive eruption. Even if we consider the lower bound of 67% CIR, 5 Ωm, the hydrothermal fluid volume fraction is estimated to be less than 5 vol % using 0.1 Ωm for the hydrothermal fluid resistivity. If using 0.3 Ωm for the hydrothermal fluid resistivity, the lower bound of 67% CIR results in a hydrothermal fluid volume fraction between the percolation threshold for effusive eruption and the change point. However, even if using 0.3 Ωm for the hydrothermal fluid resistivity, a hydrothermal fluid volume fraction less than 5 vol % can explain a wide range of the 67% CIR. As a result, the percolation threshold porosity for effusive eruption likely characterizes the impermeable roof associated with the Mogi inflation source at the depth of 2.7 or 3.7 km bsl (4.4 or 5.4 km bsf ). This suggests that we can use the percolation threshold and resistivity around the Mogi inflation source region as a potential diagnostic threshold for the forecast of a phreatic eruption at Azumayama, if we can monitor the precise resistivity around the source region.

Conclusions
A schematic diagram summarizing the main results of this study is shown in Fig. 22 horizontal dimension of the less than 3-Ωm conductor is 15 to 20 km. The 67% and 90% CIR of the conductor is 0.3 to 5 and 0.02 to 70 Ωm, respectively. The confidence interval is estimated by using the F-test to evaluate the RMS misfit change of each site with changes in the conductor's resistivity. The CRTs suggest the dataset can resolve a conductive body over 5 km in scale with unbiased resistivity at depths of 5 to 10 km bsl (6.7 to 11.7 km bsf ). The magma reservoir imaged from the final resistivity model overlaps with the slow V s region (Okada et al. 2015;Chen et al. 2018) and the low-viscosity model from InSAR data (Takada and Fukushima 2013). The Mogi inflation sources derived from GNSS and tilt data are located near the top boundary of the magma reservoir's conductor at a depth of 2.7 or 3.7 km bsl (4.4 or 5.4 km bsf ). The shallower part of the magma reservoir is interpreted as a water-saturated (5.5 wt %), dacitic melt-silicic rock-hydrothermal fluid complex. The hydrothermal fluids appear necessary, because the lower boundary of the 67% CIR (0.3 Ωm) cannot be explained by a water-saturated dacitic melt alone. Simultaneous consideration of the V s requires that hydrothermal fluid volume fraction is dominant to the dacitic melt volume fraction. The deeper part of the conductor is interpreted as a water-saturated (8.0 wt %) andesitic melt-mafic rock complex. The V s at a depth of 9 km bsl (10.7 km bsf ) suggests a fluid and/or melt volume fraction of 6-7 vol %. Taking prior consideration of the fluid and/or melt volume fraction of 6-7 vol %, the conductive hydrous phase is likewise required to explain the whole range of the 67% CIR. We interpret that the Mogi inflation source at a depth of 2.7 or 3.7 km bsl (4.4 or 5.4 km bsf ) indicates the interruption of hydrothermal fluid ascent. Assuming the Mogi source region is composed of rock-hydrothermal fluid, the hydrothermal fluid volume fraction (viz. porosity) at the Mogi inflation source is estimated to be below 5 vol %, the percolation threshold for effusive eruption. This suggests that we can use the percolation threshold and resistivity around the Mogi inflation source region as a potential diagnostic threshold for the forecast of a phreatic eruption at Azumayama. This may be applied to a phreatic eruption forecast of other active volcanoes.
The model shows no lower crust (depth range from 20 to 35 km bsl) conductor, and no correlation between the conductor beneath Oana Crater and deep, low-frequency tremors. This may be due to limited data beyond a period of 1024 s, which can be remedied by acquiring longerperiod data, or due to conductive near-surface features or the magma reservoir's conductor shielding the signals and concealing deeper features, which may be remedied with a larger survey footprint and regional array.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s40623-021-01451-y. Additional file 3. Comparisons between the electrical resistivity and seismic velocity models. Figure S3.1: Plan or map view of the V p model by Okada et al. (2015) superimposed on the resistivity model (Fig. 5). The white lines are V p contours with the interval of 0.1 km/s. The background color represents the logarithmic resistivity value. The black and white circles, red character, and black triangles, are observation sites, sampling point of V p model, Oana Crater, and edifice summits, respectively. The labels A and B for the conductors are referred to in the main text. Figure S3.2: Vertical profiles of the V p model by Okada et al. (2015) superimposed on the resistivity model (Fig. 6) in the E-W direction. The white lines and circles are V p contours and sampling point of the V p model. Other features are the same as those in Fig. 6. Figure S3.3: Vertical profiles of the V p model by Okada et al. (2015) superimposed on the resistivity model (Fig. 7) in the N-S direction. Other features are the same as those in Figs. 7 and S3.2. Figure S3.4: Plan or map view of the V p /V s distribution superimposed on the resistivity model (Fig. 5). The white lines are V p /V s contours with the interval of 0.025. The background color represents the logarithmic resistivity value. The black and white circles, red character, and black triangles, are observation sites, sampling point of V p and Vs models, Oana Crater, and edifice summits, respectively. The label A and B for the conductors are referred to in the main text. Figure S3.5: Vertical profiles of the V p /V s distribution superimposed on the resistivity model (Fig. 6) in the E-W direction. The white lines and circles are V p /V s contours and sampling point of the V p and V s model. Other features are the same as those in Figures 6. Figure S3.6: Vertical profiles of the V p /V s distribution superimposed on the resistivity model (Fig. 7) in the N-S direction. Other features are the same as those in Figures 7 and S3.4.