論文 / 著書情報
Article / Book Information
Title

Three-dimensional resistivity structure of Asama Volcano revealed bydata-space magnetotelluric inversion using unstructured tetrahedral
elements

Authors

Yoshiya Usui, Yasuo Ogawa, Koki Aizawa, Wataru Kanda, Takeshi
Hashimoto, Takao Koyama, Yusuke Yamaya, Tsuneomi Kagiyama

Citation

Geophysical Journal International, Volume 208, 3, 1359‒1372

Pub. date2017, 3

DOI

https://doi.org/10.1093/gji/ggw459

Copyright

This article has been accepted for publication in Geophysical JournalInternational (c): 2017Yoshiya Usui,Yasuo Ogawa,Koki Aizawa,Wataru
Kanda,Takeshi Hashimoto,Takao Koyama,Yusuke Yamaya,Tsuneomi
Kagiyama Published by Oxford University Press on behalf of Royal
Astronomical Society. All rights reserved.

Powered by T2R2 (Science Tokyo Research Repository)

Geophysical Journal International
Geophys. J. Int. (2017) 208, 1359–1372
Advance Access publication 2016 December 9
GJI Geomagnetism, rock magnetism and palaeomagnetism

doi: 10.1093/gji/ggw459

Three-dimensional resistivity structure of Asama Volcano
revealed by data-space magnetotelluric inversion using unstructured
tetrahedral elements
Yoshiya Usui,1,2 Yasuo Ogawa,1,2 Koki Aizawa,3 Wataru Kanda,1,2 Takeshi Hashimoto,4Takao Koyama,5 Yusuke Yamaya6 and Tsuneomi Kagiyama7
1 Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Tokyo 152-8550, Japan. E-mail: usui.y.ad@m.titech.ac.jp
2 Volcanic Fluid Research Center, Tokyo Institute of Technology, Tokyo 152-8550, Japan
3 Institute of Seismology and Volcanology, Faculty of Sciences, Kyushu University, Nagasaki 855-0843, Japan
4 Institute of Seismology and Volcanology, Faculty of Science, Hokkaido University, Hokkaido 060-0810, Japan
5 Volcano Research Center, Earthquake Research Institute, University of Tokyo, Tokyo 113-0032, Japan
6 Fukushima Renewable Energy Institute, National Institute of Advanced Industrial Science and Technology, Fukushima 963-0298, Japan
7 Aso Volcanological Laboratory, Institute for Geothermal Sciences, Graduate School of Science, Kyoto University, Kumamoto 869-1233, Japan

Accepted 2016 December 1. Received 2016 November 30; in original form 2016 July 24

SUMMARY
Asama Volcano is an andesitic composite volcano and one of the most active volcanoes
in Japan. In order to reveal electrical resistivity structure beneath the volcano accurately,
we performed a 3-D inversion of dense magnetotelluric survey data. In order to prevent
misinterpretation of the subsurface resistivity due to the steep topography around Asama
Volcano, we used an unstructured tetrahedral mesh to represent the topography. Furthermore,
we reduced the calculation time by transforming the inverse problem from the model space into
the data space. Comparison of the new data-space method to the original model-space method
showed that the calculation time required to update the model parameters was reduced as a resultof the transformation, whereas the resistivity structure obtained remained unchanged. In thesubsurface resistivity structure around Asama Volcano that was estimated from the inversion,
resistive bodies were discovered to be located under the old eruption centres. In particular,
under the 24 ka collapse caldera to the west of the presently active crater, a spherical resistive
body was found to exist in isolation. In addition, there was a widespread conductive layer
below the resistive surface layer. By comparison with previous hydrological and geochemical
studies, the conductive layer was interpreted as being a high-water-content layer and an
overlying layer rich in altered clay minerals. Because the western part of the volcanic conduit
was considered to be the resistive area, which is inferred to consist of unfractured rocks
with lower permeability than their surroundings, it would appear that the area obstructs thewestward flow of the hydrothermal fluid beneath the summit, thereby contributing to higher
concentrations of SO4 2− and Cl− in the spring water at the northern and eastern feet as wellas the uneven location of a diffuse CO2 anomaly.
Key words: Inverse theory; Electrical properties; Magnetotellurics; Volcanic arc processes.

1 I N T RO D U C T I O N
Asama Volcano is an andesitic composite volcano and one of themost active volcanoes in Japan. It is composed of three consecutive volcanic edifices (Kurofu Volcano, Hotoke-iwa Volcano andMaekake Volcano in chronological order) and three lava domes(Sekison-zan, Ko-Asama and Hanareyama lava domes) (Fig. 1).
The highest point of Asama Volcano is 2568 m above sea level(ASL) and corresponds to the summit of Maekake Volcano, where
the present active crater is located. Asama Volcano constitutes a series of WNW–ESE aligned volcanoes, known as the Asama-Eboshi

Volcanoes, along with the Eboshi Volcanoes located to the west.
The volcanic activity of the Asama-Eboshi Volcanoes started fromthe westernmost part around 1 Ma, after which the active centre
migrated from west to east during its history (Nishiki et al. 2013).
Therefore, the topography around Asama Volcano has an east–west
trend (Fig. 1). Prior to the formation of the latest edifice (Maekake
Volcano), a large conical-shaped stratovolcano at 2800–2900 m
ASL, Kurofu Volcano, was formed approximately 2 km west of thepresent active crater (Aramaki 1963). However, around 24 ka, thestratovolcano collapsed eastward and a horseshoe-shaped caldera
remained (Nishiki et al. 2013). After the activity of Kurofu Volcano


C The Authors 2016. Published by Oxford University Press on behalf of The Royal Astronomical Society.

1359

1360

Y. Usui et al.
40˚

135˚

140˚

138˚24'

138˚30'138˚36'5 km
36˚30'

36˚30'

35˚

100

0
1000

15
00

1500

20

Koasamalava dome

36˚24'

1500

36˚24'

00

Maekake
Volcano

2000

15

00

Eboshi Volcanos

Hotokeiwa
VolcanoKurofuVolcano

1500

Tomi fault
1000

Hanareyama
lava dome

AMT Stations
MT Stations

Sekizonsan
lava dome

1000

MT Stations (E-field only)

138˚24'

138˚30'138˚36'

Figure 1. Spatial distribution of the magnetotelluric observation sites used in this study. Blue triangles and green circles indicate AMT and MT sites,
respectively, at which both electric and magnetic fields were observed. Green squares indicate MT sites where only an electric field was observed. The altitudes
and locations of rivers were extracted from the Fundamental Geospatial Data developed by GSI-Japan (Geospatial Information Authority of Japan 2011). In
addition, the locations of the normal faults around the volcano (Takahashi et al. 2013) are also shown. Red filled circle indicates the location of the Jigokudani
fumarolic area.

ceased, from around 17–11 ka, Hotokeiwa Volcano had been active
slightly east of the present eruption crater (Nishiki et al. 2013),
which was almost buried by the volcanic ejecta of the following activities of Maekake Volcano. After approximately 8.5 ka (Nishiki etal. 2013), volcanic eruptions occurred at Maekake Volcano, which
is currently active. Some normal faults are recognized in the summit area (Takahashi et al. 2013; Fig. 1), which were formed by thedepression of the ground after the pumice eruptions during the laststage of the Hotokeiwa activity (Aramaki 1963).
Around Asama Volcano, continuous seismic and GPS observations have been carried out to obtain information about its subsurface structure and the mechanism of magma transportation (Aokiet al. 2013). Takeo et al. (2006) used the measurements recordedby the GPS stations to suggest the occurrence of a dyke intrusion
several kilometres to the west of the summit preceding the 2004
eruption. They also revealed the sinuous hypocentre distribution
under Asama Volcano and suggested that it represents the magma
supply path beneath the volcano. Aoki et al. (2009) attempted toreveal the relationship between volcanic activities and the subsurface structure by conducting an active source seismic survey andfound an area with high P-wave velocity to the west of the summit. They also showed that the high-velocity area corresponds tothe high-density area inferred from Bouguer gravity anomaly. Theyinterpreted that the area to the west of the summit contains solidified intrusive magma responsible for constraining magma migration under the volcano. On the other hand, Nagaoka et al. (2012)
estimated the S-wave velocity structure of the upper crust around
Asama Volcano by using ambient noise seismic tomography. They
revealed the existence of a low-velocity area at depths of 5–10 km
under the area about 8 km to the west of the summit. They interpreted that it represents a deep magma chamber beneath thedyke.
Preceding the seismic and geodetic studies, Aizawa et al. (2008)
conducted a dense magnetotelluric (MT) survey around Asama
Volcano and has already estimated the 2-D resistivity structuresof three north–south profiles and one east–west profile, and theyfound the resistive anomalies surrounded by the conductive area.
As the resistive bodies are located below the old eruption centres,
the previous work suggested that the resistive anomalies are oldsolidified magma and that the surrounding conductive area corresponds to a hydrothermal system driven by the remaining heat of themagma. They also suggested that the old solidified magma prevents
the ascent of magma. However, as the resistivity structures of theprevious study were estimated from 2-D inversion, these results

Resistivity structure of Asama Volcano
could have been subject to 3-D effects including the ones caused bysteep topography around Asama Volcano.
The topography of a mountainous area can significantly influence
the observational data of an MT survey (e.g. Müller & Haak 2004;
Usui 2015) such that these topographies can cause 3-D distortions to
the observed data and prevent us from accurately estimating a resistivity structure from 2-D analysis. In recent years, 3-D MT inversion
has become practical, and it has been applied to mountainous areas.
However, in almost all the cases, the finite-difference method has
been used and topographies were represented by rectangular grids.
Therefore, in these calculations, mountain slopes were modelled asa sequence of steps, which resulted in incorrect discontinuities in
the response functions between narrow periods unless two blocks ofgrid were added to each direction of the observation sites (Müller& Haak 2004). The addition of grids to the areas around all of
the sites increases the number of required grids and computationalcost. Recently, Kordy et al. (2016b) developed an inversion scheme
using unstructured hexahedral elements and applied it to the fielddata observed on Mount St. Helens Volcano. By using unstructured
elements, they succeeded in representing the smooth dipping surface precisely in the computational mesh. However, when a small
element was created around observation sites in a hexahedral meshto represent precise topography, the total number of elements significantly increases since the edge lengths of the outermost area are
required to be equal to those around sites.
Usui (2015) overcame these problems by developing an inversion scheme using unstructured tetrahedral elements and showed
its applicability to synthetic data affected by topography. Therefore,
in order to obtain accurate 3-D resistivity structure around Asama
Volcano by incorporating 3-D topographic effects into the inversion
model, we applied the inversion scheme of Usui (2015) to the data
set of Aizawa et al. (2008). It is the first application of 3-D MT
inversion using unstructured tetrahedral elements to field-observed
data.
Application of the inversion scheme of Usui (2015) to an actual MT data set, however, is problematic because this approach
uses the model-space Gauss–Newton method to obtain new model
parameters. Recovery of a resistivity distribution with higher resolution around the observation sites requires hundreds of thousands
or more parameter cells. In such cases, the model-space Gauss–
Newton method is not practical (Kordy et al. 2016b). To overcome
the problem and reduce the calculation time, we transformed theinverse problem from the model space into the data space according
to Siripunvaraporn et al. (2005) and Kordy et al. (2016b).
In the following sections, a brief review of the original modelspace method of Usui (2015) is first presented, followed by a detailed
description of the newly developed data-space method. Next, theobservation data and calculation conditions of the inversion targeted
for Asama Volcano are presented. Subsequently, the results of theinversion are described, after which the subsurface structure around
Asama Volcano is inferred from the result.

2 M O D I F I C AT I O N O F I N V E R S I O NSCHEME
2.1 Description of the modified inversion scheme
The 3-D inversion scheme proposed by Usui (2015) uses the modelspace Gauss–Newton method to obtain new model parameters. Assuming that these model parameters consist only of log-resistivitiesof parameter cells for ease of explanation, and the number of unknown log-resistivities and those of the observed data are Mρ and1361

N, respectively, it is necessary to solve the following linear equation
at the (k + 1)th iteration of the method(W J k )T (W J k ) + α 2 R T R δmρ,k+1 = dˆ
(1)



dˆ = −R T Rmρ,k + (W J k )T W d − F mρ,k ,

(2)

N ×N

is the diagonal matrix consisting of the recipwhere W ∈ R
rocals of the standard deviations of observed data, J k ∈ R N ×Mρ
is the sensitivity matrix, R ∈ R Mρ ×Mρ is the roughening matrix,δmρ,k+1 ∈ R Mρ consists of the increments of the model parametersat the (k + 1 )th iteration and α ∈ R is a trade-off parameter. Furthermore, mρ,k ∈ R Mρ is a vector composed of the model parameters at
the k th iteration; d ∈ R N consists of the real and imaginary partsof observed response functions, whereas F(mρ,k ) ∈ R N consists of
those of the calculated response functions. Although the solution
vector is δmρ,k+1 in eq. (1), it is equivalent to eq. (34) of Usui(2015), in which the solution vector is the model parameters at the(k + 1)th iteration. The roughening matrix R ∈ R Mρ ×Mρ gives thedifferences of log-resistivities of adjacent parameter cells r k ∈ R Mρ
as
T

(3)
Rmρ,k = r k = rk1 · · · rki · · · rkMρirki = N Fi m iρ,k −

NF


j,i

m ρ,k

(4)

j=1

where N Fi is the number of faces of the ith parameter cell, m iρ,k ∈ R
is the log-resistivity of the ith parameter cell at the kth iterationj,iand m ρ,k ∈ R is the log-resistivity of the parameter cell adjacentto the ith parameter cell through its jth face. Notably, N Fi does
not include the number of faces located on the outer boundary of
the computational mesh, which are not shared by any adjacent cell.
The coefficient matrix of eq. (1) is a dense matrix, and its dimension is Mρ . Therefore, the more model parameters are estimated ininversion, the more CPU time is required.
To overcome the problem, we transformed the computational
space from model space to data space in the same manner as
Siripunvaraporn et al. (2005) and Kordy et al. (2016b). If theroughening matrix R is non-singular, it is possible to apply the
Sherman–Morrison–Woodbury formula (Golub & Van Loan 2013)
to the inverse of the coefficient matrix of eq. (1) as
−1
(W J k )T (W J k ) + α 2 R T R−1  2 T −1
−1= α2 RT R− α R R (W J k )T  k −1 W J k α 2 R T R
(5)
−1

 k = I + (W J k ) α 2 R T R (W J k )T .

(6)

Although the matrix  k is dense, its dimension is N , which
represents the number of data. Therefore, even if the number of
unknown log-resistivities Mρ is large, it does not increase the timerequired to calculate the inverse of  k . Although the rank of theroughening matrix is less than Mρ (Uchida 1993), the roughening
matrix R is regularized by adding a small positive value  ∈ R to its
diagonals as in Kordy et al. (2016b). Kordy et al. (2016b) showedthat this procedure makes the roughening matrix regularized and,
when  is very small, its influence on the result of the inversion is
negligibly small. When R is regularized in this way, from eqs (1),(2), (5) and (6), the increments of the model parameters at the(k + 1)th iteration δmρ,k+1 can be calculated as


−1
−1
dˆ − α 2 R T R (W J k )T  k −1
δmρ,k+1 = α 2 R T R
−1

×W J k α 2 R T R dˆ
(7)

1362

Y. Usui et al.

In our inversion code, the inverse of  k ∈ R N ×N is not directly
calculated; instead, it is calculated by solving the linear equation
whose coefficient matrix is  k with the aid of the direct solver of
LAPACK (Anderson et al. 2000). As described above, the computer resources required to solve the linear equation are independent of the number of unknown log-resistivities Mρ . Furthermore,
we parallelized the construction of the coefficient matrix  k byhybrid programming with MPI and OpenMP. This hybrid parallelization becomes possible by virtue of the transformation to thedata space. Since, in the inversion code of Usui (2015), each of
the MPI processes performs calculations of different frequencies,
the construction of  k can be parallelized with the multiple process as in Siripunvaraporn & Egbert (2009) while it is difficult toparallelize the construction of the coefficient matrix of eq. (1) with
MPI in an efficient way. Matrix (α 2 R T R)−1 is equal to the square of
the inverse of α R since the roughening matrix R is symmetric, andthe inverse of the matrix is calculated by solving the linear equation whose coefficient matrix is α R. Although the dimension of thelinear equation is Mρ , it is a real sparse matrix and can be solved
quickly by the MKL PARDISO library as confirmed by Kordy et al.
(2016b).
When distortion matrices of observation sites are also estimatedin inversion, the regularization term for them is added to the objective function as in Avdeeva et al. (2015). If the number of unknown
components of distortion matrices is Mc , the increments of model
parameters at the (k + 1)th iteration δmk+1 ∈ R Mρ +Mc is calculatedby the equation T −1  T −1
T −1 d̄ − R̄ R̄
W J k ¯ k
W Jk
δmk+1 = R̄ R̄ T −1
× R̄ R̄
d̄
(8)
T 



d̄ = − R̄ T R̄mρ,k + W J k W d − F mρ,k(9)

  T −1 
T

W Jk
¯ k = I + W J k R̄ R̄

(10)

R̄ =

αR
0

0
,
βI

(11)

where J k ∈ R N ×(Mρ +Mc ) is the sensitivity matrix including thederivatives of data with respect to the distortion matrices and β ∈ R
is the trade-off parameter for the regularization term of the distortion matrices. This equation is calculated in the same way assolving eq. (7). The modified method is equivalent to the originalmodel-space method.

2.2 Synthetic data test
We confirmed that the modified inversion method gives a result
equivalent to that obtained by the original method in much lesscomputational time by performing an inversion for the synthetic
data using both methods. The target model was the model used in
the synthetic case II of Siripunvaraporn et al. (2005). It consists ofa conductive (1  m) block and resistive (100  m) block buried
in a 10  m horizontal layer and the underlying 100  m layer(Fig. 2). The 64 observation sites were set on the surface of themodel, which are marked as solid dots in Fig. 2. The full components of the impedance tensors for 16 periods from 0.1 to 1000 s
were used as the data in the calculation, and the total number of data
(real and imaginary part of the response function) was 8192. Gaussian random noises were added to all impedance tensor componentsof the synthetic data. The standard deviation of the noise was set
to 5 per cent of the maximum amplitude of the off-diagonals, thatis, 0.05 × max(|Z x y |, |Z yx |), where Z x y and Z yx denote the xy andyx components of the impedance tensor, respectively. The size ofthe computational mesh used in the inversion was 1000 × 1000 ×1100 km, of which the uppermost 100 km was the air layer and themodel was discretized with 619 070 tetrahedral elements. All theelements, except for those of the air layer, were grouped into 59 027
parameter cells by the algorithm written in Usui (2015). It should
be noted that the computational mesh used to calculate the synthetic data was different from that used for the inversion, in which
there were no pre-defined boundaries of the resistivity anomalies.
Galvanic distortion was not considered in the calculation, and thetrade-off parameter α 2 was set as 10.0. The initial resistivities ofthe parameter cells were 100  m. The small value  added to the
diagonals of the roughening matrix was 0.0001. The calculation
was performed using the computer systems of the Earthquake Information Center of the Earthquake Research Institute, Universityof Tokyo (Xeon E5-4627 v2@3.3 GHz, 512 GB), and four MPI
processes and eight OpenMP threads per process were used in eachcalculation.
After the seventh iteration of the Gauss–Newton method, convergent solutions were obtained by both the original and the modified
method. Both the rms misfit of the result of the original method
and the one of the modified method were 1.14. Fig. 2 shows that
there was not any notable difference between the resistivity structures obtained by the two methods. Therefore, we confirmed thatthe modified method gives the same result with the original method.
However, the new method was computationally more efficient thanthe original method in the phase in which new model parameters were calculated by solving the normal equation: the modified
method required 73 s per iteration, whereas the original method
required 2176 s per iteration. This confirmed that the calculation
time can be reduced significantly by transformation from the modelspace method to the data-space method.
In that model-update phase, two procedures, the dense matrix–
matrix multiplication needed to construct the coefficient matrix ofthe linear equation ([(W J k )T (W J k ) + α 2 R T R] of eq. (1);  k ofeq. (6)) and the numerical factorization of the matrix, took up mostof the calculation time. Due to the transformation to the data-spacemethod, the calculation time per iteration was reduced from 1425
to 65 s for the former procedure, while 748 to 3 s for the latter.
The speed-up of the former procedure was mainly attributed to
the reduction of the number of multiplications and summations offloating-point numbers. Since both of the coefficient matrices of themodel-space method and the data-space method are symmetric, itis necessary to calculate only the upper triangle part of them. The
model-space method, then, requires M × (M + 1) × N floatingpoint operations for the multiplication of (W J k )T and (W J k ),where M and N indicate the number of model parameters and data,
respectively. On the other hand, the data-space method requires
N × (N + 1) × M floating-point operations for the multiplicationof (W J k ) and (α 2 R T R)−1 (W J k )T . Hence, when N is less thanM, the data-space method needs less floating-point operations. In
addition, in the procedure, the use of multiple MPI processes inthe data-space method also contributed the speed-up (it took 217 s
when only one process was used). On the other hand, the speed-upof the latter procedure was owing to the decrease of the dimensionof the coefficient matrix to be factorized. The dimension of thecoefficient matrix of the data-space method is N while that of themodel-space method is M. Since the number of the floating-point
operations required by the numerical factorization with LAPACK

True model

Original method(Model-space method)

1363

Modified method(Data-space method)0z (km)

Cross section at x = 0 (km)Resistivity structure of Asama Volcano

1020

10 Ωm 1 Ωm 100 Ωm 10 Ωm
100 Ωm-40 -30 -20 -10 0 10 20 30 40 -40 -30 -20 -10 0 10 20 30 40 -40 -30 -20 -10 0 10 20 30 40
y (km)y (km)y (km)

30
20
10
x (km)

Plan view at z = 3 (km)40

01 Ωm 100 Ωm

-10
-20
-30

10 Ωm

-40
-40 -30 -20 -10 0 10 20 30 40 -40 -30 -20 -10 0 10 20 30 40 -40 -30 -20 -10 0 10 20 30 40
y (km)y (km)y (km)

0
1
2
Resistivity [log(Ωm)]
Figure 2. Comparison of the resistivity structures obtained by the original (model-space) method with those obtained by the modified (data-space) method.
The upper figures are vertical cross-sections through the centre, and the lower figures are map views at 3 km in depth. The leftmost figures indicate the true
resistivity structure. Figures adjacent to them (the middle panel) were drawn from the result obtained by the original method, whereas the rightmost panel
indicates the results obtained by the modified method. The observation points are shown by dots in all sections.

(Anderson et al. 2000) is proportional to the cube of the dimension of the matrix, the calculation time was reduced significantly
by the transformation to the data-space method when N is lessthan M.
In addition, we compared the computation time of the forward
calculation part of our code to that of a 3-D code using a hexahedral mesh. Kordy et al. (2016a) reported the computational times
required for their inversion code using unstructured hexahedral elements, in which MKL PARDISO was used in the forward calculation as in our inversion code. The computational times for thereordering and for the factorization in a forward calculation for ahexahedral mesh with 734 820 degrees of freedom were 17 and 60 s,
respectively, and the solving time for 100 right-hand sides was 11 s,
when the calculation was done with 24 cores. When we performed
the forward calculation for the unstructured tetrahedral mesh of thesynthetic test described above, which had 712 922 degrees of freedom, with 24 cores following to Kordy et al. (2016a), the times forreordering and factorization were 9 and 37 s, respectively, and thesolving time for 100 right-hand sides was 15 s. Thus, there was nodifference of more than two times between the calculation timesof the forward part for the unstructured tetrahedral mesh and forthe unstructured hexahedral mesh reported in Kordy et al. (2016a).
However, further researches are needed to conclude that the required
calculation time for an unstructured tetrahedral mesh is comparable
in a general way to that for an unstructured hexahedral mesh withsimilar degrees of freedom because it depends on the performanceof the computer system, computational mesh and other calculation
conditions.
3 METHOD
3.1 Observational data
In this study, the response functions estimated by Aizawa et al.
(2008) were used although the data recorded at the three observationsites were not used since they were noisy (one audiomagnetotelluric(AMT) site and one MT site at the north of the volcano, whichwere more than 5 km away from the summit, and the southernmost
AMT site). The measurement sites used in this study consisted of
36 MT sites and 37 AMT sites (Fig. 1). Of the former sites, 15
sites were two-channel MT (2E) sites and measured only electricfields. Except for the 2E sites, both impedance tensors and verticalmagnetic transfer functions (VMTFs) were estimated. On the other
hand, at the 2E sites, only impedance tensors were estimated byusing the magnetic field of another MT site. At every site, response
functions at about four periods per decade were used.
For the MT sites, the frequency range used in the inversion wasfrom 320 to about 0.094 Hz, whereas the frequencies ranged from1300 to 0.59 Hz for the AMT sites. In addition, VMTFs lower than
10 Hz were not used at AMT sites as in Aizawa et al. (2008). The
number of data (real and imaginary part of response function) usedin the inversion was 8472. The error floor of each impedance tensor

1364

Y. Usui et al.

Figure 3. Computational mesh used in the 3-D inversion. The horizontal centre of the mesh was the summit of Asama Volcano, and the positive x-axis
corresponds to the northern direction. (a) Overall view. Elements of the air layer are omitted to improve visibility of the Earth’s surface within the mesh. (b)Surface mesh around the survey area. Edges of elements shown in the lower figure are omitted from the upper figure. Around the observing sites, mesh sizes
are locally refined and edge lengths are shortened to approximately 10 m around the sites. Height differences around the mountainous area are emphasized byenlarging the vertical coordinates 1.5 times.

component was 5 per cent of the largest absolute value of the offdiagonals, that is, if the error of an impedance tensor component
was less than 0.05 × max(|Z x y |, |Z yx |), it was force to be 0.05 ×
max(|Z x y |, |Z yx |). The error floor of each VMTF component was
5 per cent of the maximum absolute values of its components, thatis 0.05 × max(|Tzx |, |Tzy | ) where Tzx and Tzy denote the x and ycomponents of the VMTF, respectively. In addition, when the errorof a VMTF component was smaller than 0.05, it was set to 0.05.
3.2 Computational mesh with precise topography
In this study, the influences of the topography are incorporated in
the inversion model with the aid of unstructured tetrahedral elements. The computational mesh is shown in Fig. 3. The size of themesh was 100 km × 100 km × 150 km, of which the uppermost
50 km is the air layer. The horizontal centre of the mesh was thesummit of Asama Volcano. Whereas the lengths were approximately
5 km at the outermost region, the edge lengths were shortened whenapproaching the observation sites and were scaled down to about
500 m within 15 km radius from the summit of Asama Volcano
and about 200 m within 5 km radius from the summit. Furthermore, within 50 m from each observation site, the edge lengths areapproximately 10 m.
The total number of elements was 1 140 072. All the elements,
except those of the air layer, were distributed into 195 607 parameter
cells by the algorithm used in Usui (2015). Since parameter cells
located in the uppermost 3 km around the survey area consist of
at most only a few elements, the resolution of the resistivity distribution around the area is considered to be sufficiently high. In the
deeper part, more elements constitute a parameter cell, however, we

think it is reasonable because the electromagnetic field diffuses andattenuates as the depth increases, which causes the lower sensitivity
to resistivity structure.
The mesh was created by the following two steps. First, the surface mesh covering the computational domain was generated. In
the surface mesh, the top, bottom and sides of the computational
domain were defined as well as the Earth’s surface. The last ofthese surfaces was the inner surface, where precise topography was
represented. All of the surfaces consist of triangles and they werecreated by the Delaunay triangulation method (Cheng et al. 2012).
The surfaces were produced by first generating a mesh consisting of 2-D triangles, and then the heights (z-coordinates) of thepoints of the 2-D mesh were interpolated from the 10 m grid digital elevation model of Fundamental Geospatial Data developed by
GSI-Japan (Geospatial Information Authority of Japan 2011) andETOPO1 (Amante & Eakins 2009); the former was used in the areaaround Asama Volcano (Longitude: E138.25◦ –E138.75◦ , Latitude:
N36.25◦ –N36.58◦ ), whereas the latter was used in the other areas.
Next, the 3-D mesh was generated with the aid of a tetrahedral meshgenerator TetGen (Si 2007) by using the surface mesh as piecewise
linear complexes. When generating the mesh, the inputted surface
mesh was preserved, that is, the Steiner points (automatically generated points) were only inserted in the interior space of the surface
mesh.
3.3 Calculation conditions of inversion
To calculate response functions in inversion, tangential componentsof electric fields were used while horizontal components were used
for magnetic fields. The different locations of electric and magnetic

Resistivity structure of Asama Volcano

4 T H R E E - D I M E N S I O NA L R E S I S T I V I T Y
S T RU C T U R E
Converged solutions were obtained for all of the six different tradeoff parameters α 2 . The calculation time taken to update the model
parameters in each iterative cycle was about 400 s using eight MPI
processes and eight OpenMP threads per process for the same computer system that was used in the synthetic test. As it took about
60 000 s when the original model-space inversion (Usui 2015) wasused, the calculation time was dramatically reduced by the transformation to the data-space method. Fig. 4 shows the L-curve drawn
from the data misfit and model roughness of these results. From thefigure, the point corresponding to the model of α 2 = 10.0 appears
to be the ‘knee’ of the curve, that is, the most preferable model.
Therefore, we selected the model of α 2 = 10.0, whose rms was
1.108, as the best model. Figs 5 and 6 compare the observed response functions with the ones calculated from the best model. The
apparent resistivity and the phase of the determinant impedance tensor for the representative frequency ranges are compared in Fig. 5,
whereas the induction arrows are compared in Fig. 6. These figures
illustrate that, overall, the best model explains the observed data.
In the following, the characteristic features of the best model are
shown.
Fig. 7 shows the vertical cross-sections of the 3-D resistivity
structure along the profiles of the 2-D inversions that were previously determined by Aizawa et al. (2008). As with the result of theprevious study, the obtained resistivity structure of the 3-D inversion
has a resistive (>100  m) surface layer and an underlying conductive layer (<5  m). At the westernmost profile (Line-A) and the
profile through the horseshoe-shaped caldera (Line-B), the obtained
resistivity distributions are similar to those of Aizawa et al. (2008)
although the conductors of the 2-D model appear to have slightly

2.5E+04316
2.0E+04

100
Data Misfit

fields of the 2E sites were taken into account in the inversion. The
initial subsurface resistivity was 10  m. Even if the initial subsurface resistivity was set to 100  m, this had little influence on the
result. The small value added to the diagonals of the roughening matrix was set to 0.0001 as in the synthetic test of the previous section.
In the inversion, the two-step approach was used to estimate both thesubsurface resistivity values and the components of the distortion
matrices as in Kordy et al. (2016b). First, only subsurface resistivity values were estimated (the first step) until the Gauss–Newton
iteration converged. Next, adding the components of distortion matrices to unknown model parameters, the Gauss–Newton iteration
was restarted (the second step).
The trade-off parameter α 2 was selected by the L-curve criteria from the results obtained by six different α 2 s (1.00, 3.16, 10,
31.6, 100 and 316) as in previous work (e.g. Patro et al. 2005;
Matsuno et al. 2014). The L-curve and the optimal trade-off parameter that were obtained are presented in the next section. On
the other hand, the value of the other trade-off parameter β 2 was
set to be 0.01 since Avdeeva et al. (2015) and Kordy et al. (2016b)
showed that the resistivity structure can be estimated properly withsuch a small trade-off parameter. Furthermore, it was confirmedthat when β 2 was changed to be 0.1, the estimated resistivity structures were rarely different. Even if the components of the distortion
matrix were not estimated in the inversion, the rms misfit was only
10 per cent larger than for the case in which these components were
estimated as model parameters. Therefore, the data set used in thisstudy was not considered to be significantly affected by galvanic
distortion.

1365

1.5E+0431.6
1.0E+04

10.0
1.00

3.16

5.0E+03

0.0E+00
0.0E+00

5.0E+02

1.0E+03

1.5E+03

Model Roughness
Figure 4. L-curve obtained from the results of the 3-D inversions with various trade-off parameters. The horizontal and vertical axes represent themodel roughness term and data misfit term of the objective function, respectively. The numbers adjacent to the points indicate the value of parameterα 2 giving the corresponding result. The point surrounded by a dotted circle
corresponds to the best model selected in the present study (α 2 = 10.0).

lower resistivities. At Line-A, the resistive area (R1) is found near
the crossover point to the E–W profile and is sandwiched between
more conductive layers (C1 and C2 in Fig. 7a). At about 1 km below the ground surface of the crossover point of Line-B, there isa circular resistor of about 30  m (R2) under the 24 ka collapse
caldera (Fig. 7b).
Nam et al. (2007) calculated the MT response functions along
the centre line of a 3-D trapezoidal hill using a computational meshwith the 3-D topography, and compared them with those obtained by
the 2-D forward calculation in which only the topography parallelto the profile was included in the mesh (Wannamaker et al. 1986).
They revealed that the TM-mode apparent resistivity obtained by
the 2-D calculation was not much different from that obtained by
the 3-D calculation while, in TE mode, the 2-D calculation did
not reproduce the decrease of the TE-mode apparent resistivity on
the hill obtained by the 3-D calculation. Since Line-A and Line-B
are approximately perpendicular to the level line of Asama-Eboshi
Volcanoes, the topographic effect to the TM-mode responses might
have been adequately modelled by the mesh with the 2-D topography. On the other hand, it can be presumed that the topographic
effect to the TE-mode responses, which could not be represented by
the 2-D mesh, made the resistivity values around the summit area
to be underestimated than the real values. Furthermore, the 2-D
resistivity structures of Aizawa et al. (2008) were estimated using
the 2-D inversion code developed by Ogawa & Uchida (1996). The
2-D code utilizes the finite-element method as our inversion code
does, and the regularization functional of the 2-D code is similar tothat of our 3-D code: the roughening matrix of the 2-D code gives
the difference of the log-resistivity of each parameter cell from theaverage log-resistivity of its surroundings as that of our code andboth codes do not use prior models for the regularization. Similarity to the resistivity structures of Aizawa et al. (2008) seen in our

Y. Usui et al.

Apparent resistivity
1000

1000

1000

10
00

10
00

10
00

1000

1000

1000

10
00

1000

200 - 400 Hz Cal.

200 - 400 Hz Obs.

200 - 400 Hz Cal.

1000

200 - 400 Hz Obs.

10

1000

12

Phase
1000

1366

1000

1000

6
4

15

15

00

15
00

00

15
0

0

2
2000

−4

00
10

1000

00

00

1500

1500

1500

20

00
20

1500

00

0

00

20

0
20

20

−2

2000

20

20

00

0

2000

00

2000

20

Northing (km)

8

1500

1500

00
10

1000

1500

1500

00
10

1000

00
10

1000

−6
00
10

00

20 - 40 Hz Cal.

20 - 40 Hz Obs.

1000

1000

1000

00

20 - 40 Hz Cal.

20 - 40 Hz Obs.

1000

1000
10
00

10
00

1000

1000

10
00

10

00
10

10

10

−8
12

1000

1000

6
4

15
00

2

15

2000

00

20
00

00
20

0

00
20

0
20

−2

00

2000

20

0

15

00

2000

20
00

Northing (km)

8

1500

1500

−4

1500

00
10

1000

1500

1500

1500

00
10

00
10

1000

1000

−6

1000

1000

1000

1000

10
00

1000

10
00

10
00

10
00

2 - 4 Hz Cal.

1000

1000

1000

2 - 4 Hz Obs.

1000

1000

00

00

0

2 - 4 Hz Cal.

2 - 4 Hz Obs.

10

10

10

0
10

−8
12

1000

1000

6
15

00

00

00
20

0.2 - 0.4 Hz Cal.

1000

1000

1000

10
00

10
00

1000

1000

10
00

00

20

00

00

1000

10

10

00
10

0.2 - 0.4 Hz Obs.
1000

1000

10
00

00

00
10

1000

1000

0.2 - 0.4 Hz Cal.

0.2 - 0.4 Hz Obs.
1000

1500

1500

00
10

00
10

1000

00
10

−8
12

20

00
10

1000

1500

1500

1500

−6

10

0

1500

1500

1500

2000

0
20

00

−4

2000

20

00
20

−2

15
00

00

2000

2000

0

15

00

1000

15

2

20
00

4

20

Northing (km)

8

1000

1000

6
15
00

15
00

00
10

1000

00
20

00

1500

1500

00
10

1000

00
10

1000

2

4 6

Easting (km)1
2
3Resistivity [log(Ωm)]

8 10 −8 −6 −4 −2 000

8 10 −8 −6 −4 −2 0

10

00

4 6

00
10

1000

00
10

10

00

2

Easting (km)

1500

1500

1500

10

−6
−8
−8 −6 −4 −2 0

20

0

−4

1500

1500

1500

2000

00
20

0
20

00
20

−2

15
00

0

2000

20
00

0

15
0

2000

2000

00

2

20

4

20
00

Northing (km)

8

2

4 6

Easting (km)8 10 −8 −6 −4 −2 0

2

4 6

8 10

Easting (km)0 15 30 45 60 75 90
Phase [deg]Figure 5. Maps for comparison of observed data (Obs.) with calculated responses (Cal.). The apparent resistivity and phase of the determinant impedance
tensor are shown. Frequencies at which response functions were estimated are not the same at all observation sites; thus, representative frequencies are indicatedas ranges instead of fixed values.

results could have been obtained also by this similarity of modelling
technique.
On the other hand, at the easternmost N–S profile (Line-C) and theE–W survey line (Line-EW), there are notable differences between

the resistivity structures of this study and those of Aizawa et al.
(2008). At Line-C (Fig. 7c), the resistivity structure of this study
is smoother than that of Aizawa et al. (2008). Therefore, it is presumed that some of the 3-D effects were responsible for artificially

Resistivity structure of Asama Volcano
20 - 40 Hz Obs.
1000

10
00

1000

10
00

1000

1000

10
00

1000

1000

1000

20 - 40 Hz Cal.

1000

200 - 400 Hz Cal.

1000

200 - 400 Hz Obs.

10

10
00

12

1367

1000

1000

1000

6
4

15

15

00

15
00

15
00

00

2
2000

00
10

1000

1500

1500

1500

00
10

1000

1500

1500

1500

−4

20
00

00

20

0

1500

1500

00
20

00
20

0
20

00
20

−2

2000

20

00

00

0

2000

2000

20

Northing (km)

8

00
10

1000

00
10

1000

00

1000

1000

1000

1000

10
00

10
00

1000

10
00

1000

1000

10
00

0.2 - 0.4 Hz Cal.

0.2 - 0.4 Hz Obs.

2 - 4 Hz Cal.

1000

t

1000

2 - 4 Hz Obs.

10

10

00

12

10

00
10

−8

00
10

−6

1000

1000

1000

6
4

15
00

2

15

2000

00

2000

2000

2000

00

0

1500

00
10

1000

00
10

2

4 6

Easting (km)8 10 −8 −6 −4 −2 0

2

4 6

Easting (km)8 10 −8 −6 −4 −2 0

00
10

1000

00
10

00
10

00

00

−8
−8 −6 −4 −2 0

00
10

1000

10

10

−6

1500

1500

1000

1500

1500

1500

1500

1500

20

20
00

00

00
20

20

20

00
20

0
20

−2
−4

15

00

00

0

15
00

00
20

Northing (km)

8

2

4 6

Easting (km)8 10 −8 −6 −4 −2 0

2

4 6

8 10

Easting (km)1 UnitFigure 6. Induction vectors (Parkinson convention) for the comparison of observed data (Obs.) with calculated response functions (Cal.). Red and blue arrows
indicate real and imaginary vectors, respectively.

creating a rougher structure inferred from the 2-D inversion along
this profile. In the resistivity structure obtained by this study, theconductive layer (C1) is widely spread along Line-C, particularly
in its northern part. The conductive area found by the 2-D study
along Line-EW is located several hundred metres shallower thanthe corresponding conductor inferred from this study (Fig. 7d), andshows the resistivity value of almost a half of the 3-D conductor.
In addition, in the 2-D analysis, the resistivity structure of LineEW was not consistent with those of Line-A and Line-B aroundthe intersections. We guess 3-D topography and/or 3-D resistivity
anomalies caused the inconsistency in the 2-D study. It should be
noted that, along Line-EW, we can see a difference in the resistivity
values between the eastern side of the presently active crater and its
western side: the resistivity on the eastern side at about 1 km below
sea level (BSL) is lower than that on the western side.
The horizontal section views of the 3-D resistivity structure areshown in Fig. 8. At altitudes from 0.5 to 1.5 km ASL, a circular
resistive area (R2) can be found under the caldera, and this area
corresponds to the resistive body found in the vertical cross-section
along Line-B and Line-EW. Aizawa et al. (2008) inferred that theresistive body was isolated because a low-phase zone under the collapse caldera of Kurofu Volcano, that appeared not to be connectedto the other surrounding low-phase zone, was found in the phasemaps for several frequencies calculated from the determinant ofimpedance tensor. The 3-D inversion we performed enabled us toconfirm that the resistive body is isolated. At depths ranging from
0.5 km BSL to 0.5 km ASL (Figs 8d–f), a relatively high-resistivity
area (R1) exists below the eastern end of the Eboshi Volcanoes, andthis area corresponds to the sandwiched resistive area of Line-A(Fig. 7a). Fig. 8 shows that R1 seems to extend to beneath R2, although Aizawa et al. (2008) suggested that the resistive bodies ofLine-A and Line-B exist in isolation. Therefore, we attempted toverify whether R1 and R2 are connected below sea level by performing a sensitivity test. Specifically, the resistivity of the area to thesoutheast of R1 (marked by the blue solid line in Fig. 8) was forced
to be higher than or equal to 3  m at a depth between 1.0 km BSL
and 1.0 km ASL (to be more precise, if the resistivity of a parameter
cell whose gravity centre lies in that range was lower than 3  m, it
was forced to be 3  m), and a forward calculation was performed
for the model. The rms obtained by the forward calculation was
1.111; thus, it was quite similar to the original rms (1.108). Since
only AMT sites were located over the area and there was shallower
conductor (<3  m) above the area, the data in this study had lesssensitivity to the area. Therefore, it is unclear whether R1 and R2
are connected below sea level from the data in this study.
By performing the 3-D inversion using the mesh with precise
topography, we succeeded in obtaining more plausible resistivity
structure of Asama Volcano than that proposed by Aizawa et al.
(2008), which was a first-order structure as described in that paper.
As was pointed out above, the resistivity structures along LineC and Line-EW showed some differences between the previous
2-D study and this study. Although some of the features seen in
the resistivity structure had already been suggested by the previous
study, this study gave credibility to these features. Furthermore, thisstudy provided the information about the horizontal distributions ofresistivity anomalies, and revealed that conductive areas (C1 and
C2) are widely distributed under the mountainous area.

5 I N T E R P R E TAT I O N O F R E S I S T I V I T Y
S T RU C T U R E
Fig. 9 shows precisely relocated hypocentres of the earthquakes
from 2004 January to 2015 January (M. Takeo 2015, unpublisheddata). These hypocentres are plotted over the resistivity structure.
The distribution of the hypocentres is considered to indicate themagma pathway beneath Asama Volcano (Takeo et al. 2006; Aoki
et al. 2013). The hypocentres are located around the spherical

1368

Y. Usui et al.
(a) Line-A
Line-EW

S30W

N30E

Depth(km)

-2
-1

R1

C2

0
-6

-4

0

-2

C1

2
4
Distance(km)

6

10

8

(b) Line-B

Depth(km)-3Collapse Caldera
Line-EW

S20W

N20E
AMT

-2

MT

-1

C2

0
-4

MT (E-field only)

R2
-2

C1
0

4

0
2
Distance(km)

1

2

3

4

Resistivity [log(Ωm)](c) Line-C
Line-EW

South

North
Boring hole

Depth(km)-2-1

1.6 km

C1

0
-8

-6

-4

-2

0
2
Distance(km)

4

6

8

10

(d) Line-EW
East

Depth(km)

West

-3

Line-A

Line-B

Line-C

-2
-1

C1

R2

0

-6

-4

-2

0
2
Distance(km)

4

6

8

Figure 7. Vertical cross-sections of the 3-D resistivity structure along the profiles of Aizawa et al. (2008). (a) Cross-section along Line-A. (b) Cross-section
along Line-B. (c) Cross-section along Line-C. (d) Cross-section along the E–W survey line (Line-EW). The locations of these profiles are shown in Fig. 8(a).
Red arrow in (b) indicates the location of the Jigokudani fumarolic area, which corresponds approximately to the intersection of Line-B and the Tomi fault(Fig. 1). In addition, the location of a borehole for a commercial hot spring is indicated in (c).

resistive body below the collapse caldera (R2): horizontally distributed earthquakes are located below the resistive body and vertically distributed earthquakes are lined up to the east of the resistor.
As with Aizawa et al. (2008), we interpret that the resistive bodyconsists of old and solidified magma and it forces the ascending
magma to move eastward to the present conduit. Such structuralcontrols were also suggested for Iwate Volcano, Japan (Aizawa
et al. 2009b), Aso Volcano, Japan (Kanda et al. 2008), Kirishima
Volcanoes, Japan (Aizawa et al. 2014) and El Hierro Island (CanaryIslands), Spain (Garcia-Yeguas et al. 2014). Because the resistive
body is located beneath the centre of the old stratovolcano, repeatedintrusions of magma under the old eruption centre might have reduced the connected porosity of the area, causing it to become highly
resistive. On the other hand, resistive body R1 is located in the areaof high P-wave velocity and high density revealed by Aoki et al.
(2009), which coincides with the area of dyke intrusion during the2004 eruptions (Takeo et al. 2006). Aoki et al. (2009) suggested thatthe high velocity is due to the solidification of repeatedly intruded

magma since the slowly cooled magma becomes unfractured rockwith high density and high seismic velocity. Thus, R1 is consideredto also probably have lowered the connected porosity compared tothe surrounding area, and is then imaged as a resistive area.
Next, we focus on the resistive layer directly beneath the land
surface and the underlying conductive layer commonly found in
the obtained resistivity structure. We interpret the resistive layer
as an undersaturated zone. Since volcanoclastic rocks often havelarge connected porosity, the resistivity depends heavily on the water content. Therefore, for the layer to be resistive, its water content
necessarily needs to be lower. On the other hand, the underlying
conductive layers (C1 and C2) are considered to have high watercontent and to be rich in altered clay minerals. In volcanic terrain,layer rich in altered clay minerals, especially smectite, is commonly
considered to be a major cause of conductive anomalies (e.g. Ogawa
et al. 1998; Ussher et al. 2000; Nurhasan et al. 2006; Aizawa et al.
2009a; Kanda et al. 2010; Yamaya et al. 2013). Morita et al. (2016)
inferred from the diffuse CO2 mapping that the conductive area

Resistivity structure of Asama Volcano(b) 1.5 km ASL(c) 1.0 km ASL

12

12

12

4

8

6

6

4

Line-EW

2000

20

00

0

20

Line-B

1500

R2

1500

00

00

00

20

20

−2

−2

−4

00
10

R2
00
10

1000

−6

−6

10

00

00

10

00
10

−8
−8 −6 −4 −2

1500

1500

−4

00
10

1000

−6
0
4
2
Easting (km)

6

8

10

−8
−8 −6 −4 −20
42
Easting (km)

6

8

10

−8
−8 −6 −4 −20
42
Easting (km)(f) 0.5 km BSL

12

12
1000

10
00

10
00

10

1000

1000

1000

8

8

8

6

6

6

2

4

R1

15

00

C1

2

4

R1

15

00

C1

C1

20

20

−2

−2

1500

C2

−4

00
10

1500

1000

00

6

8

10

−8
−8 −6 −4 −2

00
10

00
10

10

00
10

0
4
2
Easting (km)

C2

−6

−6

−6

1500

−4

00
10

1000

1500

0

R2

−8
−8 −6 −4 −2

00

0
20

C2
1000

0
00

00

1500

1500

−4

15

20

20

−2

R1

1000

2000

00

0

00

0

2

2000

2000

10

20
00

4

8

1000

10
00

10

1000

6

1000

(e) Sea level

12
1000

(d) 0.5 km ASL
10

Northing (km)

10
00

2000

2000

15

1000

00

2

0

−4

15

00

2

0

1000

4
15

00

2

−2

1000

8

00

15

1000

1000

20

10
00

Line-C

6

10

1000

1000

Line-A

8

Northing (km)

10

1000

10
00

10

1000

1000

(a) 2.0 km ASL

1369

0
4
2
Easting (km)

6

8

10

−8
−8 −6 −4 −20
42
Easting (km)

6

8

10

AMT Stations

0
1
2
3
4
Resistivity [log(Ωm)]

MT Stations
MT Stations (E-field only)

Figure 8. Horizontal cross-sections of the resistivity structure obtained by 3-D inversion. In (a), profiles of Aizawa et al. (2008) are also shown. The area of
which the resistivity was changed to 3  m for a sensitivity test is marked as a blue solid line in (d)–(f).

beneath the summit of Asama Volcano is constituted of to theuppermost clay-rich layer and the underlying hydrothermal fluidlayer. Furthermore, there are several studies indicating the presenceof high-water-content layer under Asama Volcano. Kazama et al.
(2015) calculated the steady distribution of the unconfined groundwater level around Asama Volcano based on hydrological physics.
They revealed that the groundwater level around the summit area isabout 1.5 km ASL and becomes higher to the west, which is consistent with the top of the conductive layer. In addition, at a distanceof 7 km along Line-C (Fig. 7c), a commercial hot spring with waterat about 73◦ C is drawn from 1.6 km below ground (Aizawa et al.
2008). Fig. 7(c) shows that the bottom of the borehole corresponds
to the conductive layer. In particular, around the summit area, it
was suggested that there exist a hydrothermal system driven by theheat of magma (Aizawa et al. 2008; Morita et al. 2016). Near theintersection of Line-B and the Tomi fault, which is a normal fault
cutting the southern part of the collapse caldera of Kurofu Volcano, a fumarolic area called Jigokudani is located (Figs 1 and 7b).

The volcanic gases rich in SO2 and CO2 and the acid spring water (pH 3–6) rich in SO4 2− are discharged from the fumarolic area
(Aramaki 1963; Suzuki 2002). The resistivity just beneath the fumarolic area (about 10  m) is lower than the surroundings, andthis low-resistivity area is connected to the conductive layer C2(Fig. 7b). It appears that the volcanic fluids supplied from the underlying hydrothermal system to the fumarolic area through thefractures around the Tomi fault.
In addition to geophysical studies, geochemical studies werealso performed around Asama Volcano. Suzuki (2002) analysed thechemical and isotopic compositions of 80 samples of spring wateraround the volcano and found that the spring water on the northern and eastern slopes have higher SO4 2− and Cl− concentrations
than on the southern and western slopes (Fig. 10). Suzuki (2002)
interpreted these results to mean that the spring water discharged
at the northern and eastern feet originates from the groundwater
affected by the sulphate-chloride thermal water formed beneaththe presently active crater, and that the groundwater flows to the1370

Y. Usui et al.
WestLine-C

Line-B−3
Depth (km)

East

−2
−1
R2

0
1
-4

-2

0
Distance(km)

2

AMT

0

1

2

3

4

MT

Resistivity [log(Ωm)]

MT (E-field only)

Figure 9. Comparison of the resistivity structure to the hypocentre distribution around Asama Volcano. On the vertical cross-section along Line-EW,
precisely relocated hypocentres of earthquakes detected between 2004 January and 2015 January (M. Takeo 2015, unpublished data), within ±500 mof the section, are plotted. The spherical resistive body beneath the collapsed
caldera (R2) is marked by a broken circle.

northern and eastern feet along the surface of the collapsed calderaof Kurofu Volcano. Kazama et al. (2015) also suggested that groundwater flows from the western side of Asama Volcano to the eastern
flank through the volcanic conduit below the presently active crater.
Furthermore, the helium isotope ratio of the hot spring water to thenorth of Line-C (Fig. 7c) suggests the contribution of magmatic fluids to the hot spring (Aizawa et al. 2008). Therefore, we suggested
that the conductive layer extending from beneath the presently active crater to the northern and eastern area (C1) represents the highwater-content layer with the water affected by the hydrothermal
water beneath the presently active crater. To the west of the volcanic
conduit, there were resistive areas (R1 and R2), both of which wereinferred to have low connected porosities, that is, low permeability(Farquharson et al. 2015). These resistive areas are considered toprevent groundwater from flowing to the west of the conduit andmake the spring water rich in SO4 2− and Cl− eccentrically locatesat the eastern and northern slopes.
In addition, Morita et al. (2016) observed the diffused CO2 andsoil temperature around the summit area and revealed that highly
diffuse CO2 anomalies were found only to the east of the present
vent. They ascribed the heterogeneous distribution of diffuse CO2
to the different depths of the fluid layer. However, as in the case of
the anion compositions discussed above, the differences among theconnected porosities beneath the summit area may also contribute to
the heterogeneous diffuse CO2 flux. That is, it would appear that themagmatic gases of the conduit flow preferentially into the eastern(b) Cl- concentrations

12

12
1000

1000

(a) SO42- concentrations
10

100

10

0

100

0
1000

1000

8

8

6

6
4

4
15

15

00

00

2

2
2000

20

20

20

00
20

00

−2

00

0

00

0

2000

1500

−2

1500

−4

−4

00
10

1000

1500

1500

00
10

1000

−6

−6
−6 −4 −2

mg/L

0
4
2
Easting (km)

6

8

10

0 30 010 ~ 2 ~ ~ 4 40
< 10 20 30
>

−8
−8

00

00

10

10

−8
−8

−6 −4 −2042Easting (km)

6

8

10

0 30
10 2
0
5
< 5 ~ 0 ~ 20 ~ > 3
1mg/L
AMT Stations

0
1
2
3
4
Resistivity [log(Ωm)]

MT Stations
MT Stations (E -field only)

Figure 10. Comparison of the resistivity structure at the sea level to the SO4 2− and Cl− concentrations of the spring water around Asama Volcano (Suzuki
2002). (a) SO4 2− concentrations. (b) Cl− concentrations.

Resistivity structure of Asama Volcano

12High-water-content layer under thenorthern and eastern feet

10

4

15

AC K N OW L E D G E M E N T S
We are grateful to the members of the Volcanic Fluid ResearchCenter of the Tokyo Institute of Technology for their many helpful suggestions. We would like to sincerely thank Prof M. Takeo(the University of Tokyo) for providing us with the data of the
precisely relocated hypocentres around Asama Volcano. We would
like to express our gratitude to Dr H. Suzuki (Komazawa Univer-

00

2
2000

Diffuse CO2 anomaly

00

0

00

The resistivity structure around Asama Volcano was able to be revealed with the aid of an unstructured tetrahedral mesh with precise
topography. In order to reduce the calculation time required to update unknown resistivities in inversion, we transformed the inverse
problem from the model space into the data space. A synthetic test
showed that the new method can significantly reduce the calculation time required to update model parameters while giving thesame result with the original method. The resistivity structure determined around Asama Volcano was commonly characterized by aresistive layer directly beneath the land surface and the underlying
conductors. The former was considered to consist of volcanoclastic
rocks with low water content, whereas the latter was interpreted ascomprising high-water-content and cray-rich layers. In addition, thefeature of the resistivity structure changed in the east–west direction beyond the volcanic conduit below the crater that is presently
active. To the west of the conduit, a spherical resistive body was
found to exist beneath the old eruption centres, which would consistof unfractured rocks and prevent the volcanic gases and hydrothermal water beneath the active crater from flowing to the southwestof the conduit. Furthermore, it would contribute to the higher concentrations of SO4 2− and Cl− in the spring water discharged at
the northern and eastern feet and would cause a diffuse CO2 flux
anomaly to appear only to the east of the present crater. Therefore,by estimating the 3-D resistivity structure, we were able to infer
the relationship between the unfractured area formed by magma
intrusion, as found by geophysical studies, and the characteristic
distribution of the anion concentrations of the spring water and the
diffuse CO2 flux, as revealed by geochemical studies. The present
work was the first application of 3-D MT inversion using unstructured tetrahedral elements to field-observed data, and demonstrated
that it can provide sensible results with actual MT data.

Subsurface water affected bythermal water formed beneaththe present active crater

20

Northing (km)

6

−2

1500

1500

Resistive area consisting of
unfractured rocks

−4

00
10

1000

−6
6 C O N C LU S I O N S

1000

8

20

part, and may be responsible for the highly diffuse CO2 anomalies
that are only found in the eastern part of the summit area.
Fig. 11 shows a schematic view of the relations between the resistive structure and the subsurface fluids. In the east–west direction,
the feature of the subsurface resistivity structure changes beyond
the volcanic conduit. To the west of the conduit, resistive areas (R1and R2) exist below old eruption centres, which may have lower
connected porosity than their surroundings due to the repeated intrusion of the magma from the underlying magma chamber. These
resistive areas are likely to influence on the anion concentrationsof the ground water and the heterogeneity of the diffuse CO2 flux
around the summit area, which also varies beyond the conduit, by
forcing magmatic gases and hydrothermal water to flow preferentially to the east and north of the active crater. By performing the3-D inversion using the mesh with precise topography, we wereable to infer the relationship between the resistive area formed bymagma intrusion and the characteristic distributions of the SO4 2−
and Cl− concentrations of the spring water and diffuse CO2 flux.

1371

High-water-content layer
under the southern foot−8
−8 −6 −4 −20
42
Easting (km)

6

8

10

Figure 11. Schematic of the relation between the resistive area to the westof the present crater and subsurface fluids. The resistive area is inferred
to contribute to the heterogeneity of the SO4 2− and Cl− concentrationsof the spring water and the diffuse CO2 flux around the summit area by
preferentially directing the flow of subsurface fluids to the east of the summit.

sity) for kindly sending us a copy of his doctoral thesis. In addition,
we would like to thank P. Wannamaker, an anonymous reviewer
and the editor in charge (U. Weckmann) for their valuable comments that significantly improved the manuscript. This study was
supported by the 7th National Project for Volcano Eruption Prediction. For this study, we have used the computer systems of the
Earthquake and Volcano Information Center of the Earthquake Research Institute, the University of Tokyo. Some figures were produced with the aid of Generic Mapping Tools (GMT) software(Wessel & Smith 1998).

REFERENCES
Aizawa, K., Ogawa, Y., Hashimoto, T., Koyama, T., Kanda, W., Yamaya,
Y., Mishina, M. & Kagiyama, T., 2008. Shallow resistivity structure of
Asama volcano and its implications for magma ascent process in the 2004
eruption, J. Volc. Geotherm. Res., 173, 165–177.
Aizawa, K., Ogawa, Y. & Ishido, T., 2009a. Groundwater flow and hydrothermal systems within volcanic edifices: delineation by electric self-potential
and magnetotellurics, J. geophys. Res., 114(1), 1–12.
Aizawa, K., Ogawa, Y., Mishina, M., Takahashi, K., Nagaoka, S., Takagi,N., Sakanaka, S. & Miura, T., 2009b. Structural controls on the 1998 volcanic unrest at Iwate volcano: relationship between a shallow, electricallyresistive body and the possible ascent route of magmatic fluid, J. Volc.
Geotherm. Res., 187(1–2), 131–139.
Aizawa, K. et al., 2014. Three-dimensional resistivity structure and magma
plumbing system of the Kirishima Volcanoes as inferred from broadband
magnetotelluric data, J. geophys. Res., 119, 198–215.
Amante, C. & Eakins, B.W., 2009. ‘ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis’. NOAA Technical Memorandum NESDIS NGDC-24, 19 pp. Available at:
http://www.ngdc.noaa.gov/mgg/global/global.html, last accessed 31 May
2015.
Anderson, E. et al., 2000. LAPACK User’s Guide, 3rd edn, SIAM.

1372

Y. Usui et al.

Aoki, Y. et al., 2009. P-wave velocity structure beneath Asama Volcano,Japan, inferred from active source seismic experiment, J. Volc. Geotherm.
Res., 187, 272–277.
Aoki, Y., Takeo, M., Ohminato, M., Nagaoka, Y. & Nishida, K., 2013.
Magma Pathway and Its Structural Controls of Asama Volcano, Japan,
Geol. Soc. London, Spec. Publ., 380 pp.
Aramaki, S., 1963. Geology of Asama volcano, J. Fac. Sci. Univ. Tokyo,
14(2), 229–443.
Avdeeva, A., Moorkamp, M., Avdeev, D.B., Jegen, M. & Miensopust, M.P.,
2015. Three-dimensional inversion of magnetotelluric impedance tensor
data and full distortion matrix, Geophys. J. Int., 202(1), 464–481.
Cheng, S.-W, Dey, T.K. & Shewchuk, J.R., 2012. Delaunay Mesh Generation, CRC Press.
Farquharson, J., Heap, M.J., Varley, N.R., Baud, P. & Reuschlé, T., 2015.
Permeability and porosity relationships of edifice-forming andesites:
a combined field and laboratory study, J. Volc. Geotherm. Res., 297,
52–68.
Garcia-Yeguas, A., Ibanez, J.M., Koulakov, I., Jakovlev, A., Romero-Ruiz,
M.C. & Prudencio, J., 2014. Seismic tomography model reveals mantle
magma sources of recent volcanic activity at El Hierro Island (CanaryIslands, Spain), Geophys. J. Int., 199(3), 1739–1750.
Geospatial Information Authority of Japan, 2011. ‘Fundamental geospatial
data’. Available at: http://fgd.gsi.go.jp/download/, last accessed 7 June2015 (in Japanese).
Golub, G.H. & Van Loan, C.F., 2013. Matrix Computations, 4th edn, Johns
Hopkins Univ. Press.
Kanda, W., Tanaka, Y., Utsugi, M., Takakura, S., Hashimoto, T. & Inoue,H., 2008. A preparation zone for volcanic explosions beneath Naka-dakecrater, Aso volcano, as inferred from magnetotelluric surveys, J. Volc.
Geotherm. Res., 178(1), 32–45.
Kanda, W., Utsugi, M., Tanaka, Y., Hashimoto, T., Fujii, I., Hasenaka, T. &Shigeno, N., 2010. A heating process of Kuchi-erabu-jima volcano, Japan,
as inferred from geomagnetic field variations and electrical structure, J.
Volc. Geotherm. Res., 189(1–2), 158–171.
Kazama, T., Okubo, S., Sugano, T., Matsumoto, S., Sun, W., Tanaka, Y. &Koyama, E., 2015. Absolute gravity change associated with magma mass
movement in the conduit of Asama Volcano (Central Japan), revealed byphysical modeling of hydrological gravity disturbances, J. geophys. Res.,
120(2), 1263–1287.
Kordy, M., Wannamaker, P., Maris, V., Cherkaev, E. & Hill, G., 2016a.
3-D magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on SMP
computers—Part I: Forward problem and parameter Jacobians, Geophys.
J. Int., 204(1), 74–93.
Kordy, M., Wannamaker, P., Maris, V., Cherkaev, E. & Hill, G., 2016b.
3-dimensional magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized onsymmetric multiprocessor computers—Part II: Direct data-space inverse
solution, Geophys. J. Int., 204(1), 94–110.
Matsuno, T., Chave, A.D., Jones, A.G., Muller, M.R. & Evans, R.L.,
2014. Robust magnetotelluric inversion, Geophys. J. Int., 196(3),1365–1374.
Morita, M., Mori, T., Kazahaya, R. & Tsuji, H., 2016. Diffuse carbon dioxide
emissions from hidden subsurface structures at Asama volcano, Japan,
Bull. Volcanol., 78(3), 17, doi:10.1007/s00445-016-1008-5.
Müller, A. & Haak, V., 2004. 3-D modeling of the deep electrical conductivity of Merapi volcano (Central Java): integrating magnetotellurics,
induction vectors and the effects of steep topography, J. Volc. Geotherm.
Res., 138(3–4), 205–222.
Nagaoka, Y., Nishida, K., Aoki, Y., Takeo, M. & Ohminato, T., 2012. Seismic
imaging of magma chamber beneath an active volcano, Earth planet. Sci.
Lett., 333–334, 1–8.

Nam, M.J., Kim, H.J., Song, Y., Lee, T.J., Son, J.S. & Suh, J.H., 2007.
3D magnetotelluric modelling including surface topography, Geophys.
Prospect., 55(2), 277–287.
Nishiki, K., Takahashi, K. & Matsumoto, A., 2013. Spatial and temporal change of volcanic field in Asama-Eboshi volcanoes, central
Japan, J. geol. Soc. Japan, 119(7), 474–487 (in Japanese with English
abstract).
Nurhasan, Ogawa, Y., Ujihara, N., Tank, S.B., Honkura, Y., Onizawa, S.,
Mori, T. & Makino, M., 2006. Two electrical conductors beneath KusatsuShirane volcano, Japan, imaged by audiomagnetotellurics, and their implications for the hydrothermal system, Earth Planets Space, 58(8), 1053–
1059.
Ogawa, Y. & Uchida, T., 1996. A two-dimensional magnetotelluric inversion
assuming Gaussian static shift, Geophys. J. Int., 126, 69–76.
Ogawa, Y., Matsushima, N., Oshima, H., Takakura, S., Utsugi, M., Hirano,K., Igarashi, M. & Doi, T., 1998. A resistivity cross section of Usu volcano, Hokkaido, Japan, by audiomagnetotelluric soundings, Earth Planets
Space, 50(4), 339–346.
Patro, B.P.K., Brasse, H., Sarma, S.V.S. & Harinarayana, T., 2005. Electrical structure of the crust below the Deccan Flood Basalts (India), inferred from magnetotelluric soundings, Geophys. J. Int., 163(3), 931–
943.
Si, H., 2007. ‘TetGen: a quality tetrahedral mesh generator and 3D Delaunay triangulator’. Available at: http://tetgen.berlios.de, last accessed 22
February 2014.
Siripunvaraporn, W. & Egbert, G., 2009. WSINV3DMT: vertical magnetic
field transfer function inversion and parallel implementation, Phys. Earthplanet. Inter., 173(3–4), 317–329.
Siripunvaraporn, W., Egbert, G., Lenbury, Y. & Uyeshima, M., 2005. Threedimensional magnetotelluric inversion: data-space method, Phys. Earthplanet. Inter., 150(1–3), 3–14.
Suzuki, H., 2002. Chemical evolution of spring water and groundwaterflow system in Asama volcano, Central Japan, PhD thesis, University ofTsukuba, (A), p. 2991.
Takahashi, M., Otsuka, T., Hirakawa, T., Nagai, M., Yasui, M. & Aramaki, S.,
2013. Whole-rock major element chemistry for eruptive products of theEboshi Volcano Group, Central Japan: summary of 222 analytical data,
in Proceedings of The Institute of Natural Sciences, Nihon University,
Vol. 48, pp. 111–140 (in Japanese with English abstract)
Takeo, M., Aoki, Y., Ohminato, T. & Yamamoto, M., 2006. Magma supply
path beneath Mt. Asama volcano, Japan, Geophys. Res. Lett., 33(15),
1–5.
Uchida, T., 1993. Smooth 2-D inversion for magnetotelluric data basedon statistical criterion ABIC, J. Geomag. Geoelectr., 45(1990), 841–
858.
Ussher, G., Harvey, C., Johnstone, R., Anderson, E. & Zealand, N., 2000.
Understanding the resistivities observed in geothermal systems, in Proceedings, World Geothermal Congress 2000, Kyushu-Tohoku, Japan,
pp. 1915–1920.
Usui, Y., 2015. 3-D inversion of magnetotelluric data using unstructured
tetrahedral elements: applicability to data affected by topography, Geophys. J. Int., 202(2), 828–849.
Wannamaker, P., Stodt, J.A. & Rijot, L., 1986. Two-dimensional topographic
responses in magnetotellurics modeled using finite elements, Geophysics,
51(11), 2131–2144.
Wessel, P. & Smith, W.H.F., 1998. New, improved version of generic mapping
tools released, EOS, Trans. Am. geophys. Un., 79(47), 579–579.
Yamaya, Y., Alanis, P.K.B., Takeuchi, A., Cordon, J.M., Jr., Mogi, T.,
Hashimoto, T., Sasai, Y. & Nagao, T., 2013. A large hydrothermal
reservoir beneath Taal Volcano (Philippines) revealed by magnetotelluric resistivity survey: 2D resistivity modeling, Bull. Volcanol., 75(729),
doi:10.1007/s00445-013-0729-y.