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

Electrical image of subduction zone beneath northeastern Japan

Authors

Masahiro Ichiki, Yasuo Ogawa, Toshiki Kaida, Takao Koyama, Makoto
Uyeshima, Tomotsugu Demachi, Satoshi Hirahara, Yoshimori Honkura,
Wataru Kanda, Toshio Kono, Masaki Matsushima, Takashi Nakayama,
Syuichi Suzuki, Hiroaki Toh

Citation

Journal of Geophysical Research, , ,Pub. date2015, 12

Copyright

Copyright (c) 2015 American Geophysical Union (AGU)

Creative Commons

See next page.

Powered by T2R2 (Science Tokyo Research Repository)

License

Creative Commons : CC BY-NC-ND

PUBLICATIONSJournal of Geophysical Research: Solid EarthRESEARCH ARTICLE
10.1002/2015JB012028
Key Points:• Flow overturn toward the back arc in
the mantle wedge is imaged byresistivity• Fluid pathways from the slab top
to lower crust are revealed beneath
NE Japan• The lower crust conductor requires atleast 0.7 vol. % of ﬂuid fraction

Supporting Information:• Texts S2, S3, and S5, and Figures S1and S3–S5
Correspondence to:
M. Ichiki,ichiki@m.tohoku.ac.jp

Citation:
Ichiki, M., et al. (2015), Electrical image ofsubduction zone beneath northeasternJapan, J. Geophys. Res. Solid Earth, 120,
doi:10.1002/2015JB012028.
Received 11 MAR 2015
Accepted 1 OCT 2015
Accepted article online 8 OCT 2015

Electrical image of subduction zone beneath
northeastern Japan
Masahiro Ichiki1, Yasuo Ogawa2, Toshiki Kaida1, Takao Koyama3, Makoto Uyeshima3,
Tomotsugu Demachi1, Satoshi Hirahara1, Yoshimori Honkura2,4, Wataru Kanda2, Toshio Kono1,
Masaki Matsushima4, Takashi Nakayama1, Syuichi Suzuki1, and Hiroaki Toh5
1

Graduate School of Science, Tohoku University, Sendai, Japan, 2Volcanic Fluid Research Center, Tokyo Institute ofTechnology, Tokyo, Japan, 3Earthquake Research Institute, University of Tokyo, Tokyo, Japan, 4Graduate School of Scienceand Engineering, Tokyo Institute of Technology, Tokyo, Japan, 5Graduate School of Science, Kyoto University, Kyoto, Japan

Abstract We conducted long-period magnetotelluric observations in northeastern Japan from 2010
to 2013 to investigate the three-dimensional electrical resistivity distribution of the subduction zone.
Incorporating prior information of the subducting slab into the inversion scheme, we obtained a three-dimensional
resistivity model in which a vertically continuous conductive zone is imaged from the subducting slab surface tothe lower crust beneath the Ou Backbone Range. The conductive body indicates a saline ﬂuid and/or melt
pathway from the subducting slab surface to the lower crust. The lower crust conductor is less than 10 Ω m, andwe estimate a saline ﬂuid and/or melt fraction of at least 0.7 vol. %. Other resistivity proﬁles in the across-arc
direction reveal that the conductive body segregates from the subducting slab surface at 80–100 km depth andtakes an overturned form toward the back arc. The head of the conducting body reaches the lower crust justbeneath Mt. Gassan, one of the prominent back-arc volcanoes in the system.

1. Introduction
Water in the deep earth accelerates the strain rate of mantle silicates and diminishes the solidus temperature[e.g., Mei and Kohlstedt, 2000a, 2000b; Green, 1973]. In this sense, water plays a key role in mantle dynamics
and magma generation. Presumably, water inside the earth is supplied mostly from subduction zones, whichis estimated to be 870–1000 Tg/yr of water [Peacock, 1990; van Keken et al., 2011]. Many earthquakes andvolcanic eruptions develop around subduction zones, which are affected by ﬂuids released from the slabs.
Investigations to determine the distribution of ﬂuids in the mantle wedge are pertinent for revealing ﬂuid
transport after dehydration of the subducting slab and for accounting for earthquakes and volcanic activity.
Northeastern Japan (NEJ) contains a well-studied subduction system and is an ideal place to use magnetotelluric (MT) methods to delineate melt and ﬂuid transport. To delineate the ﬂuid distribution in the mantle
wedge beneath NEJ, we performed a three-dimensional (3-D) electrical resistivity (or conductivity) study.
Long-period (longer than 103 s) geomagnetic transfer functions and magnetotellurics (MT) in NEJ have been
widely observed during the 1970s and 1980s [see Kanda and Ogawa, 2014, and references therein].
Incorporating the ocean bathymetry and sea water into the models, Ogawa et al. [1986] and Utada [1987]
proposed two-dimensional (2-D) resistivity models in the crust and upper mantle beneath NEJ for the ﬁrst
time. Since the 1990s, Network-MT observation targeting the upper mantle has been expanded to NEJ[Uyeshima et al., 2002]. However, these previous studies analyzed the data over a limited number of periods
and could only identify uniform resistivity in the mantle wedge beneath NEJ. This study collected MT data in
the wide period range (20–61,440 s) and provides a more detailed view of the resistivity distribution in the
mantle wedge beneath NEJ.

©2015. The Authors.
This is an open access article under theterms of the Creative Commons
Attribution-NonCommercial-NoDerivs
License, which permits use and distribution in any medium, provided theoriginal work is properly cited, the use isnon-commercial and no modiﬁcations
or adaptations are made.

ICHIKI ET AL.

It is well known that ﬂuids in the upper mantle have various types of conﬁgurations, such as supercritical ﬂuid
at grain boundaries, hydrogen in nominally anhydrous minerals, and ﬂuid inclusions. Electrical conductivity
measurements of mantle silicates under wet and dry conditions in laboratory experiments have revealed thatelectrical conductivity is consistently sensitive to the presence of water irrespective of conﬁguration [e.g.,Karato, 1990; Jödicke, 1992; Yoshino et al., 2006]. Recent 3-D MT, particularly of the USArray project, haverevealed rippled conductors around the Moho depth in map view [e.g., Patro and Egbert, 2008; Kelbert
et al., 2012; Meqbel et al., 2014]. The conductors are interpreted as saline ﬂuids and/or melt (hereafter denoted
as ﬂuid/melt) associated with underplating [e.g., Wannamaker et al., 2008]. In this paper, we show some
ELECTRICAL IMAGE BENEATH NE JAPAN1

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

conductive plumes in the crust and upper mantle below NEJ, which correspond well with active volcanoes atthe surface. The conductive plumes indicate ﬂuid/melt pathways from the surface of the subducting slab tolower crust conductors. One of the pathways appears very speciﬁc and seems to be associated with the
back-arc volcanism. We discuss this result and the implications of this speciﬁc pathway later in the paper.
In the following section, we ﬁrst give an overview of the tectonic settings of NEJ and summarize thelong-period MT experiments, data processing, and distortion analysis of the MT. Section 3 describes how
we adapted the 3-D inversion framework to the subduction zone. Section 4 describes the obtained resistivity
structure and its key features. Section 5 describes the sensitivity tests for the key features using the Fréchet
derivative and boxcar test, and then we discuss all the results and conclude this paper in sections 6 and 7.

2. Experiment2.1. Outline of Tectonic Settings
NEJ is located on the Okhotsuk plate in the MOVEL model [DeMets et al., 2010]. Beneath the Okhotsuk plate,the Paciﬁc plate subducts at a dip angle of about 30° and at a rate of 8–9 cm/yr [e.g., Cruciani et al., 2005].
Below NEJ, the Paciﬁc plate is very cold and old (~130 Ma), and its thickness has been estimated to be about
90 km by seismic receiver function [Kawakatsu et al., 2009]. A prominent orogenic belt called the Ou
Backbone Range (OBR) runs through central NEJ in the along-arc direction (Figure 1). The OBR likely
developed by crustal softening owing to ﬂuids squeezed from the magma chamber in the lower crust andby upthrust caused by the principal regional stress axis that is parallel to the across-arc direction [e.g.,
Hasegawa et al., 2005]. In this paper, the back arc and fore arc are deﬁned as the Japan Sea side andPaciﬁc Ocean side of the OBR. Active Quaternary volcanoes in NEJ are concentrated along OBR. However,there exists some enigmatic Quaternary volcanism in the back-arc region [e.g., Tatsumi, 1989]. Tamura
et al. [2002] pointed out that the back-arc volcanism is located on some zonal clusters of volcanoes in theacross-arc direction (Figure 1). They called the zonal clusters hot ﬁngers. Volcanism on each hot ﬁnger track
likely developed with a 70–80 km eastward migration in the last 5 Ma [Kondo et al., 2004; Honda andYoshida, 2005].
2.2. Observation and MT Response
Long-period MT data were acquired at 59 land sites in NEJ (Figure 1) by the GEOFLUID project in Japan(http://www.geoﬂuids.titech.ac.jp/eng/index_e.html). The dimensions of the survey area are 150 × 200 km
in the north-south (N-S) and east-west (E-W) directions, respectively, and the typical spacing between sites
is 20 km. Either a LEMI-417M or NIMS system, using a ﬂuxgate magnetometer, was installed at each land site
during a 2-to-4 month time period. The sampling rate was 1 Hz. Telluric ﬁeld variation was recorded with adipole length of 40 m using Pb-PbCl2 electrodes. The installation of the electrodes was performed by usinga magnetic compass, and so the obtained electric ﬁeld data are assumed to have been oriented to local
geomagnetic coordinates. We carefully maintained the level of the ﬂuxgate magnetometers and kept themagnetic y component within plus or minus a few hundred nanoteslas when installing the sites. The
magnetic ﬁeld data at eight sites were acquired using LEMI-417M gimballed ﬂuxgate magnetometers, but
the ﬂuxgate magnetometers at 51 sites were not gimballed.
To remove tilt changes, baseline steps, and drifts of the ﬂuxgate magnetometers, we ﬁrst subtracted
magnetic ﬁeld variations to which a median ﬁlter was applied from raw data (Figure S1 of the supporting
information). The time-window width of the median ﬁlter was 1 day. To delineate the response of the medianﬁlter, we calculated the power spectrum density (PSD) of magnetic ﬁeld variation before and after thedemedian operation. The test magnetic ﬁeld variation was acquired using a gimballed ﬂuxgate magnetometer. The PSD revealed that the demedian operation eliminated natural magnetic ﬁeld variations atperiods longer than 1 day as well as tilt changes, baseline steps, and drifts, whereas it did not modify thenature of the magnetic ﬁeld variation in periods shorter than 1 day (Figure S1). Henceforth, period ranges
longer than 1 day are not dealt with in this paper. The horizontal coordinate of the ﬂuxgate magnetometerat each site should be somewhat rotated from the local geomagnetic coordinate depending on the local timeof installation and on the installation itself. We rotated the horizontal coordinate so that the declination of theobservation coordinate coincided with that of the reference local geomagnetic coordinate. The declination of
the observation coordinate at each site was calculated from the median ﬁltered data at the beginning of
ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN2

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 1. Observation site superimposed on morphological and geophysical background of NEJ. Circle shows the MT
observation sites. The triangle shows volcanoes and the numbering volcanoes are mentioned in the text: 1. Ichinomegata
Maar, 2. Mt. Yakeishi, 3. Mt. Kurikoma, 4. Mt. Naruko, 5. Mt. Chokai, and 6. Mt. Gassan. Most arc volcanoes are aligned north tosouth at the center of the northeastern Japan arc, which is called the Ou Backbone Range. Mts. Chokai and Gassan are locatedin the back-arc side from the Ou Backbone Range. These are deﬁned as back-arc volcanism in this paper. Back-arc volcanism islocated along some zonal volcano clusters delineated by white line. The zonal volcano distribution is called hot ﬁngers[Tamura et al., 2002]. Yellow lines show isodepths of subducting Paciﬁc plate surface [Kita et al., 2010; Nakajima and Hasegawa,2006]. Solid lines are permanent major rivers. Child map: Plate boundaries around study area. Dashed line shows the plateboundary. Rectangle surrounded by thick solid line represents the region of the main ﬁgure. Paciﬁc plate subducts beneathNorth American plate.

observation because the beginning of each observation was expected to have been performed on leveled
coordinates. The reference local geomagnetic declination at each site was obtained from interpolation of grid
data of Geomagnetic Chart 2010 observed by the Geospatial Information Authority of Japan (http://vldb.gsi.
go.jp/sokuchi/geomag/menu_04/index-e.html). The grid width is 3 min.
The obtained electromagnetic ﬁeld variation data were processed into MT responses with the bounded inﬂuence robust remote reference processing code (BIRRP) [Chave and Thomson, 2004]. We calculated the MT
responses for all sites at 20 and 30 s and their overtone periods up to 40,960 and 61,440 s (24 periods in total).
Each site had three to seven adjacent sites that simultaneously observed the electromagnetic ﬁeld. We
selected the site that had the best quality data among those sites as a remote reference. A 20 km separation
between target and remote reference site is sufﬁcient to obtain good quality MT responses in the period
range of 30 to 20,000 s because of the high ﬁdelity electromagnetic signals measured in NEJ. Figure 2 shows
the spatial distribution of the MT nondiagonal phase response at periods of 320, 1920, 5120, and 15,360 s. The
Pyx represents the characteristics of the E-W telluric current that is nearly perpendicular to the Japan Trench.
The Pyx distribution at a period of 320 s shows several isolated patches corresponding to local crustal heterogeneity. At periods of 1920 and 5120 s, a low-Pyx domain of less than 135° occurs to the east and a high-Pyx
domain of more than 135° occurs to the west of OBR. A few patches that correspond to complex
heterogeneity around Mts. Yakeishi, Kurikoma, Naruko, and Gassan are superimposed on the background.

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN3

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 2. Spatial distribution of MT nondiagonal phase response at the periods of 320, 1920, 5120, and 15,360 s. Left andright columns show Pxy and Pyx distribution, respectively. The Pxy is MT phase response between electric ﬁeld in the northsouth component and magnetic ﬁeld in the east-west component, and Pyx is that between electric ﬁeld in the east-westcomponent and magnetic ﬁeld in the north-south component. Cross represents observation site. The spatial distribution isdrawn using cosine arch ﬁltering with 20 km width (use of “grdﬁlter” in GMT) [Wessel and Smith, 1998].

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN4

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 3. (a) Spatial distribution of the skew angle (β angle) and reference axis direction of the phase tensor at the periodsof 320, 1920, 5120, and 15,360 s. The background color shows the β angle distribution drawn by the same means of
Figure 2. White bar direction shows the reference axis direction at each site and period. The length is arbitral. (b) Rose
diagram of the reference axis direction in each period range.

Although the phase tensor analysis revealed that the conductivity distribution is entirely 3-D in this study
region (described later), this spatial background distribution of Pyx indicates a regional two-dimensional
image in which a cold and electrically resistive subducting slab lies to the east and a warm and conductive
mantle wedge is located in the west. At 15,360 s, the Pyx distribution no longer assumes a two-dimensional
pattern. In contrast to the Pyx distribution, the Pxy spatial distribution is dominantly affected by the 3-D
character inherent in the study region.
2.3. Dimensionality Analysis
The MT response is distorted by the static electric ﬁeld caused by charge buildup at boundaries of local
surﬁcial conductivity heterogeneity around the observation site, which are smaller than the shortest skindepth [e.g., Jiracek, 1990]. Hence, the phase tensor, which is calculated from the MT response but which is
not inﬂuenced by galvanic effect, is often used to highlight details of the conductivity distribution as afunction of period (or depth) [Caldwell et al., 2004; Bibby et al., 2005]. The skew angle β and reference axis
α of the phase tensor allow us to probe the dimensionality of the conductivity distribution [Bibby et al.,2005]. Figure 3a shows the phase tensor reference axis of each site superimposed on the spatial distributionof skew angle magnitude at the same periods used in Figure 2. The regions around Mts. Yakeishi, Kurikoma,
Naruko, and Gassan and Kitakami Low Land show larger skew angle magnitudes of more than 10° throughout the period range of 10 to 105 s. This suggests a 3-D resistivity distribution in the crust and upper mantle
around these active volcanoes and the Kitakami River. In other areas, small skew angle magnitudes of less

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN5

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

than 5° are dominant in the period
range of 10 to 1000 s. Hence, the crusts
outside of the volcanic regions and the
Kitakami River are consistent with a
two-dimensional resistivity distribution.
The skew angle magnitude increases
over the entire region at periods over1000 s. The rose diagram of reference
axis direction (Figure 3b) provides aclearer image of these trends. Note that
there is a 90° ambiguity in strike direction
to the reference axis. In the period rangeof 10 to 1000 s, the statistics of the reference axes have a distinctive modal direction of west-northwest. This indicates a
two-dimensional resistivity distribution
in correspondence with the skew angle
result. A small amount of scattering of
the reference axis direction is seen in
the period range of 10 to 1000 s. This
Figure 4. Horizontal mesh in the central region of model space. Black dot
scattering stems mainly from the dataand open triangle represent observation site and volcano, respectively.
in the 3-D conductivity region of theactive volcanoes and Kitakami River area,
as described above. An obvious modal direction is observed in the period range of 1000 to 104 s. However, it isrotated about 30° clockwise from the modal direction in the period range of 10 to 1000 s. Furthermore, a large
skew angle magnitude region covers the area in the period range of 1000 to 104 s (see the 5120 s period
illustration in Figure 3a). At periods of more than 104 s, the reference axis directions are completely random
statistics. Consequently, no consistent reference axis direction exists throughout the period range of 10 to
105 s. Moreover, the spatially wide and periodic distribution of large skew angle magnitude reveals the 3-Dfeature of the study area. Hence, 3-D analysis and interpretation are appropriate for this data set.

3. Inversion
3.1. Outline of the Inversion
The input data set for the subsequent 3-D inversion was arranged depending on the data quality. We excluded
the MT impedance data for which the standard error (SE) was inﬁnity or was “not a number” in the output of
BIRRP processing from the input data set. We also eliminated extreme outlier MT impedances. Here extreme
outlier means that the two-sided twice SE ranges lie away from the main trend of the MT impedance response
curve as a function of period. The error ﬂoors were all assigned as 10% of each individual impedance component throughout the inversion process. Because of the limitation of our computer resources, we used the MT
impedances at 30 s and its harmonics up to 61,440 s (12 periods). Note that we conﬁrmed that the use of the
MT impedances at the alternative 12 periods (20 s and its harmonics) gave the same conductivity models
and the root-mean-squared misﬁts normalized by the standard errors (hereafter denoted as RMS misﬁts).
The 3-D inversion was carried out by using WSINV3D-MPI-1.2 [Siripunvaraporn et al., 2005; Siripunvaraporn andEgbert, 2009]. The WSINV3D seeks a likelihood model giving a stationary point of the following functional U:hiUðm; λÞ ¼ ðm  m0 ÞT Cm 1 ðm  m0 Þ þλ1 ðd  F ½mÞT Cd 1 ðd  F ½mÞ ;(1)where m is the conductivity model, m0 is the reference or prior model, Cm is the model covariance matrix[Siripunvaraporn and Egbert, 2000], d is the observed data, F[ ] is the forward operator, Cd is the data
covariance matrix, and λ1 is the Lagrange multiplier. The ﬁrst term is called the regularization or penalty
term. We describe more details of WSINV3D and its algorithms in the supporting information.
Figure 4 shows the mesh distribution in the central model space. The total model space domain is ±700 km in the
N-S direction, ±744 km in the E-W direction, and 1200 km in depth. The ﬁnest width of the lateral mesh is 3.8 km.
ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN

6

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

The surface vertical grids down to 1 km
are 100, 200, 400, and 600 m. The total
number of meshes is 47 × 66 × 57
(with seven air layers) = 176,814.
Ocean bathymetry (ETOPO1) andseawater resistivity (0.33 Ω m) were
ﬁxed and incorporated beforehand
into all models (Figure 5). The total
number of unknown model parameters was 148,938.
3.2. Search for Reference Modeland Model Covariance Parameters
We took the following three steps inthe inversion to obtain the highest likelihood conductivity model: (a) search
for the best reference model by trial
and error using the default model covariance parameters, (b) search for the
Figure 5. Bird′s eye view of ocean bathymetry employed in the inversion.
optimal model covariance parameters
using the best reference model obtained in (a), and (c) complement the search for the reference model by using
model covariance parameters other than the default ones in (a).
In the ﬁrst step, we gave τ = 1, δi = 0.5 (i = x, y, z) to the default model covariance parameters. Here 4τδi (i = x, y, z)
gives a model covariance width or model smoothness in x, y, and z directions (see supporting information).
The reference model, m0, in the regularization term of equation (1) is requisite in the MT inversion targetingthe subduction zone. The regularization term in the form of λ1 Lm without the reference model is the mostpopular [e.g., Constable et al., 1987; Uchida, 1993]. Here L is a ﬁnite difference matrix of the model parameters
between laterally and vertically adjacent blocks (smoothness constraint). This regularization form without thereference model has successfully overcome unstable and/or unrealistic solutions in many cases; however, when
we use the regularization term without the reference model for a subducting slab and mantle wedge, the inversion usually provides only a less-informative model in which the electrically resistive region corresponding to
the subducting slab is widely blurred [e.g., Toh et al., 2006; Matsuno et al., 2010]. This is ascribed to an MT
response insensitive to the surface boundary of the resistive block [e.g., Torres-Verdin and Bostick, 1992]. To
avoid recovering such models, we experimented with assigning the surface boundary and thickness of the
subducting slab to the reference model. Accordingly, we attempted three types of reference models in operating the 3-D inversion. The model types were uniform models (ﬁve models), layered models (ﬁve models), andmodels incorporating the subducting slab into the layered model (hereafter denoted as slab models; fourmodels). The reference models attempted are shown in Figure 6 and Table 1. The resistivities of the attempted
uniform reference models were 5000, 1000, 500, 100, and 10 Ω m. The boundaries of the layered models were
ﬁxed at depths of 40 (crust), 200 (shallow upper mantle), and 400 (deep upper mantle) km. The resistivities ofthe deep upper mantle and mantle deeper than 400 km were ﬁxed at 10 and 1 Ω m, respectively. The slab
models were composed of the best layered model and the subducting slab. The slab geometry was
assigned by referring to the upper boundary of the subducting Paciﬁc plate below NEJ deﬁned by the seismicity [Kita et al., 2010; Nakajima and Hasegawa, 2006] and to the Paciﬁc plate thickness of 90 km[Kawakatsu et al., 2009]. We experimented with 1000, 5000, 10 k, and 20 k Ω m of slab resistivity. The bestreference model was selected using the following rules. (1) The reference model giving a minimum RMS
misﬁt was selected as optimal when no inversion attained stage II of the Occam process. (2) When there
were multiple types of reference models that progressed inversion to stage II, we adopted the model inwhich more prior information was incorporated as the best reference model.
The next step was to determine the optimal model covariance parameters using the best reference modelobtained in the ﬁrst step. We adopted the model covariance parameters that gave the minimal RMS misﬁt.
Here the attempted values of 4τδi were 0.5, 1.0, 2.0 (default), 4.0, and 8.0. In this step, we removed the schemeof stage II of the Occam process. Generally, there likely exists a minimal RMS misﬁt as a function of model

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN7

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 6. Schematic view of vertical proﬁle of (a) layered and (b) slab reference model in east-west direction. White crossindicates grid point. White line box in Figure 6a corresponds to the domain of Figure 6b.

covariance parameters in the WSINV3D inversion. We conﬁrmed this fact by using two synthetic data sets; theexample data given by Siripunvaraporn et al. [2005] and the other were calculated from a model in which asingle conductive cube was embedded in the best reference model (see supporting information). Aizawa
et al. [2014] also indicated that there was a minimal RMS misﬁt as a function of model covariance parametersin their data set.
In the third step, we considered the dependence of the best reference model on the model covarianceparameters. In the ﬁrst step, the best reference model was evaluated by using 4τδi = 2.0. If we started
the three steps using different model covariance parameters, another reference model might be adopted,
and it might possibly give a smaller RMS misﬁt. Hence, we looked for a likelihood model in the model
domain that was far from the search range in steps 1 and 2. Speciﬁcally, we repeated step 1 using4τδi = 0.5 and 8.0, the two cases of extreme model covariance parameter, respectively. Then we conﬁrmed
whether the best reference model was changed and whether the RMS misﬁt fell below the smallest RMS
misﬁt obtained in the second step.

a

Table 1. Attempted Reference Resistivity Models
Model
ID

Crust
(0–40 km Depth)

Shallow Upper Mantle(40–200 km Depth)

Deep Upper Mantle(200–400 km Depth)

400–1200 km
Depth

Slab

Uniform

1
2
3
4
5

5000
1000
500
100
10

5000
1000
500
100
10

5000
1000
500
100
10

5000
1000
500
100
10

–
–
–
–
–

Layered

6
7
8
9
10

1000
500
1000
500
100

100
100
10
10
10

10
10
10
10
10

1
1
1
1
1

–
–
–
–
–

Slab

11
12
13
14

1000
1000
1000
1000

100
100
100
100

10
10
10
10

1
1
1
1

1000
5000
10 k
20 k

Model Typea

ICHIKI ET AL.

Unit in Ωm. See text and Figure 6 for detail.

ELECTRICAL IMAGE BENEATH NE JAPAN8

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

aTable 2. The Best RMS Misﬁt of the Inversion Using Various Kinds of Reference Models
4τδi = 2.0

4τδi = 0.5

4τδi = 8.0

Model ID

Iteration

λRMS Misﬁt

Iteration

λRMS Misﬁt

Iteration

λRMS Misﬁt

Uniform

1
2
3
4
5

5
7
4*
3
1

1.0
0.5
0.5*
0.5
1.0

3.75
2.81
3.01*
3.18
3.75

4
3
3*
8
1

0.5
0.0
0.0*
0.0
0.5

6.05
3.32
3.28*
3.08
3.65

7
2
5
4
3

1.5
3.0
1.0
1.5
1.5

4.17
4.77
3.66
3.49
3.74

Layered

6
7
8
9
10

8
9
9*
7
6

0.5
0.5
0.0*
0.5
0.5

2.72
2.85
2.95*
3.11
3.48

7*
3
3*
3
2

0.0*
0.0
0.5*
0.5
0.0

3.41*
3.33
3.61*
3.89
3.47

7
6
6
6
3

1.5
1.5
1.5
1.5
1.5

3.26
3.35
3.34
3.32
3.53

Slab

11
12
13
14

7
5
5
5

0.5
0.5
0.5
0.5

2.66
2.48
2.47
2.47

3
3
3
3

0.0
0.5
0.0
0.0

3.26
2.88
3.08
3.07

4*
4
4*
4*

1.0*
1.0
1.0*
1.0*

3.27*
3.34
3.24*
3.24*

Model TypeaThe λ (unit in common logarithm) and Iteration represent Lagrange multiplier and iteration times at which the best RMS misﬁt was achieved, respectively. The *means manually determined model (see text).

We found that the Occam process occasionally failed. This occurred because an extremely small Lagrange
multiplier gave rise to a large-contrast conductivity model, and thereby, an inversion process did notconverge. In this case, we repeated the inversion with the Lagrange multiplier ﬁxed at several values anddetermined the Lagrange multiplier and conductivity model that gave a minimal RMS misﬁt. This isrepresented as the manually determined model in Table 1, shown in the next section.

4. Results
Table 2 shows the RMS misﬁts that resulted from the various kinds of reference models (Table 1) and thedefault model covariance parameters, 4τδi = 2.0. No inversion attained stage II of the Occam process. The
RMS misﬁt decreased as we incorporated more prior information into the reference model. In other words,
the slab reference models provided smaller RMS misﬁts than any other layered or uniform reference model.
The 5000 Ω m slab resistivity (model 12) diminished the RMS misﬁt by 0.18 compared to the 1000 Ω m
resistivity (model 11), while an increase of slab resistivity to more than 5000 Ω m led to little improvement in
RMS misﬁt. This result suggests that the gross resistivity of the subducting Paciﬁc plate is at least 5000 Ω m
in the shallow upper mantle below NEJ. We adopted the slab model with 5000 Ω m slab resistivity (model
12) as the best reference model. Figure 7 shows RMS misﬁt as a function of model covariance parameters,
4τδi, and the resultant Lagrange multiplier in terms of the best reference model. We determined the optimal
4τδi to be 2.0, which gave the minimal RMS misﬁt. Table 1 shows the RMS misﬁts obtained by using themodel covariance parameters of 4τδi = 0.5 and 8.0, respectively. Having carried out inversions in the wider
model space, we could not ﬁnd a superior model to the result of the best reference model and 4τδi = 2.0. Accordingly,
we took the model obtained by usingreference model 12 and 4τδi = 2.0 as
the ﬁnal result.

Figure 7. RMS misﬁt (circle) and Lagrange multiplier, λ, (triangle) as afunction of model covariance parameter, 4τδi, in terms of the bestreference model (model 12 in Table 1).

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN

Figure 8a shows map views of the resistivity distribution at the depths of 22.5,
32.5, 52.5, and 72.5 km. Figure 8b
depicts them as resistivity perturbations
relative to the reference model.
Prominent conductive anomalies, C1 to
C3, in the mid-crust (22.5 km) are
observed beneath and in the vicinity of9Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 8. (a) Map views of resistivity distribution at depths of 22.5, 32.5, 52.5, and 72.5 km. C1 to C3 and R1 to R2 areprominent conductors and resistors discussed in text, respectively. Black dot and open triangle represent observationsites and volcanoes, respectively. (b) Same as Figure 8a but showing in the resistivity perturbation in terms of the bestreference model (model 12 in Table 1).

Mts. Yakeishi, Gassan, and Chokai. The C1 and C2 conductors below Mt. Yakeishi and Mt. Gassan, respectively,
extend to deeper levels, while the C3 conductor below Mt. Chokai disappears in the lower crust (32.5 km). The
horizontal location of the C1 conductor follows OBR, whereas the C2 conductor retreats about 50 km to theback-arc side from OBR. These two aspects are clearly delineated in vertical proﬁles perpendicular to
the trench axis (Figure 9). The B–B′ proﬁle in Figure 9 runs through the C1 conductor below Mt. Yakeishi.
The proﬁle reveals that a conductive zone adhering to the subducting slab in the upper mantle detaches
from the slab at a depth of 80–100 km and rises to the C1 conductor in the lower crust and uppermost mantle.
On the other hand, the C–C′ proﬁle passing through Mt. Gassan shows that the conductive zone rises toward
the back-arc side after detaching from the slab. As a result, the conductive zone assumes the shape of an
overturned plume, and the overturn likely gives rise to Mt. Gassan. The conductive zone adhering to the subducting slab at depths greater than 100 km is common to all the vertical proﬁles. Figure 9 also demonstrates
that low-frequency earthquake swarms (represented by pink circles) occur in and around the C1 to C3
conductors in the lower crust and uppermost mantle, which is consistent with previous MT results targeting

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN10

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 8. (continued)

the crust [e.g., Mitsuhata et al., 2001; Ogawa et al., 2001, 2014]. These low-frequency earthquake swarms areinterpreted as volcanic or magmatic tremors because of their depths greater than the brittle-ductile transition zone and the nondouble-couple component of focal mechanisms [Hasegawa and Yamamoto, 1994].
We observed a couple of signiﬁcant resistive anomalies, R1 and R2, in the crust below the fore arc. The R1
resistor is convincingly related to Tono batholiths, early Cretaceous granitic plutons in the Kitakami
Mountains [Tsuchiya and Kanisawa, 1994, and references therein]: it is exposed at the earth’s surface in theA–A′ proﬁle, and the location of R1 in the map view at a depth of 22.5 km coincides with the surface distribution of Tono batholiths. The R2 resistor in proﬁle C–C′ is not exposed at the earth’s surface and has no correlation with surface geology. The presence of the R2 resistor is reliable, but in the next section, we demonstrate
the limited resolution to the resistivity of this feature.
Figure 10a shows the spatial distribution of RMS misﬁt of the ﬁnal model at eight periods. Each distribution is
calculated from all four complex components of MT impedance tensor. The total RMS misﬁt is 2.48. The small
RMS misﬁt at the period of 30,720 is due to the large standard errors. Although a few sites in the northeasternand southwestern marginal regions indicate poor RMS misﬁts throughout the period range, the ﬁnal modelaccomplishes good spatially uniform ﬁts of the data. Figure 10b shows the RMS distribution calculated fromICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN11

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 9. (a) Vertical proﬁle of resistivity distribution perpendicular to the trench axis. Black dots and pink circles in resistivity proﬁle represent normal microearthquakes and low-frequency tremors. Red triangles and thick horizontal solid lines show volcanoes and land area. C1 to C3 and R1 to R2 are prominent conductors andresistors discussed in text, respectively. (b) Same as Figure 9a but showing in the resistivity perturbation in terms of the best reference model.

the 3-D model obtained by using reference model 6. The total RMS misﬁt is 2.72. The reference model 6 does
not contain a slab. An improvement in RMS misﬁt (Figure 10a) is found in the region surrounding the yellow
line at periods of 960 and 1920 s. Incorporating the subducting slab into the reference model results in the
inclusion of an intricate mantle wedge toe into the model. This improves the RMS misﬁts over the mantle
wedge toe region at periods of 960 and 1920 s.
So far, we have not mentioned the geomagnetic transfer function. Because of the limitation of our computer resources and calculation time, geomagnetic transfer function data were not used simultaneously
in the inversion. In this paper, we only compared the observed induction arrows with those calculated
from the ﬁnal resistivity model. The observed and calculated induction arrows at periods of 320, 1920,
5120, and 15,360 s are shown in the supporting information. The observed and calculated induction
arrows ﬁt fairly well, although the ﬁts of induction arrows as well as magnetotellurics at a period of
15,360 s are poor.

5. Sensitivity Test for Key Features
To test the sensitivity and resolvability of key features in the obtained electrical conductivity model, we
derived a modiﬁed Fréchet derivative calculation and a boxcar test simulating the pseudo-impulse
response of the inversion process. The key features treated here are (1) the C1 conductor near the Moho
below Mt. Yakeishi, (2) the R2 resistor below the fore arc, (3) the C2 conductor and the overturn of the

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN12

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 9. (continued)

conductive body on proﬁle C–C′, and (4) the conductive zone adhering to the subducting slab at depths
greater than 100 km.
Fréchet derivative permits the calculation of a locally linearized change in response due to a change in areference model [e.g., Boerner and West, 1989]. In magnetotellurics, the Fréchet derivative is deﬁned asδZ ij ðrk ; x; T Þ ¼ K ij ðrk ; r’; m; T Þδmðr’Þdr’ þ O kδmðr’Þk2 ;(2)∫

where Kij is the Fréchet derivative (i, j = x, y components). Other parameters are Zij: MT impedance, rk:
the location vector of each observation station k, r′: the location of the perturbed model parameter, x: acertain representative location of the perturbed region (x ∈ r′), m: a model parameter as a function of space,and Τ: the period. Next, we treat the summation of the above equation in terms of the observation station(hereafter denoted as Δij):
Δij ðx; T Þ ¼

sites
X
k

∫K ij ðrk ; r’; σ0 ðr’Þ; T Þδlogρðr’Þdr’;(3)

where ρ and σ 0 indicate resistivity and reference conductivity, respectively. We took the perturbation oflogarithmic resistivity as δm(r′). In the one-dimensional case, Oldenburg [1979] derived the following
analytical form of equation (2):E ðz; T Þ 2
δZ ðT Þ ¼ μ0 σ 0 ðzÞδlnρðz Þdz;(4)
Bð0; T Þ∫

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN13

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 10. (a) Spatial distribution of RMS misﬁt at eight periods in terms of the ﬁnal model (obtained by using reference
model 12 and 4τδi = 2.0). Cross and close triangle show observation site and volcano, respectively. Yellow circles represent
the areas in which RMS misﬁts are signiﬁcantly improved by using the reference model 12 incorporating the subducting slab
into. (b) Spatial distribution of RMS misﬁts of the 3-D model obtained by using reference model 6 and 4τδi = 2.0. The reference
model 6 is a slab-free model from reference model 12. RMS misﬁts in the area as shown by yellow circles are inferior to those in
Figure 10a. See text for details.

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN

14

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 10. (continued)

where E(z,Τ ) and B(0,Τ ) are the electric ﬁeld at depth z and magnetic ﬂux density at the earth’s surface,respectively. If we assume that the 3-D Fréchet derivative takes a form similar to that of the one-dimensionalcase, equation (3) is reduced tosites
XΔij ðx; T Þ ¼ μ0 σ 0 ðxÞδ log ρðxÞδV ðxÞ
κij ðrk ; x; T Þ;
(5)k

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN15

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 11. (top) Total modiﬁed Fréchet derivative and (bottom) upper bound of total phase change as a function of depth
below (a) back arc, (b) Ou Backbone Range, and (c) fore arc, at the periods of 60 (solid line), 480 (dotted line), 3840 (dashedline), and 30,720 (chain line) s, respectively. Vertical axis corresponds to the depth of the center of perturbed cube. Unit of2 22total modiﬁed Fréchet derivative is Ω m log (S/m). The perturbed cube of 20 km on side is located along B–B′ proﬁle,and the assigned resistivity perturbation is three orders of magnitude more conductive than the best reference model.
Locations of B–B′ proﬁle (solid line) and perturbed cubes (shaded areas) are shown in each child map in the plain view.

where x is a representative location of perturbed region r′. In equation (5), we assumed the following
conditions to simplify the expression: (1) the perturbed region r′ is comparatively small and (2) Kij[rk,r′,σ 0(r′),T ],sitesXκij ðrk ; x; T Þ indicates the relative
σ 0, and δlogρ are constant within the perturbed region, respectively.
ksensitivity in terms of perturbed location x on the condition of observation network, rk. We hereafter callthis the total modiﬁed Fréchet derivative (TMFD). We do not divide TMFD by the number of sites to discuss
the resolvability of the observation network—when the MT response change is locally restricted in an
observation network, the average is affected by the number of insensitive observation sites. If we compare aresolvability of an observation network to that of a narrower one with a same observation density, theaverage operation results in an unfair evaluation of the wider observation network. In the following forward
calculation for estimating TMFD, we allotted 2.0 and an approximate cube with length on a side of 20 km toδlogρ(x) (two orders of magnitude more conductive) and the perturbed region, respectively. Note that we
obtained almost the same results using the resistivity contrast 1.0. To display the ﬁnal TMFD, we show thesites X
κij ðrk ; x; T Þ . The TMFD provides a kind of amplitude responsesummation of absolute amplitude, viz.
ksensitivity. To delineate phase response sensitivity, we show the upper bound on the total phase change interms of observation sites. The phase change at each site, δφij(rk,x,Τ ), is expressed asδZ ij ðrk ; x; T Þ(6)δφij ðrk ; x; T Þ≾ 
Z ij ðrk ; x; T ÞICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN16

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 11. (continued)

and the upper bound on the total phase change, Φij(x,Τ ), is
sites X
δZ ij ðrk ; x; T ÞΦij ðx; T Þ ¼
Z ij ðrk ; T Þ ;(7)kwhere Zij(rk,T) in equations (6) and (7) is the impedance response at site k and period T to the reference model.
The TMFD and Φij(x,T) down to 150 km below the back arc, OBR, and fore arc in the B–B′ proﬁle are plotted in
Figure 11. The site locations are identical to those for the real data, and the reference model is ﬁxed on thebest reference model (model 12 in Table 1). The TMFD deeper than 100 km below back arc is an order ofmagnitude smaller than that at 50 km depth and about two orders of magnitude less than that of the crust,
respectively. The Φij(x,T) of each component deeper than 100 km below the back arc is less than 10°. These
results suggest that feature (4) is poorly resolved with our MT data. On the other hand, the TMFD and Φij(x,T )of the entire period range have sufﬁcient sensitivity down to 40 km below the back arc, OBR, and fore arc.
Hence, the C1 conductor is well constrained by our MT data. With regard to feature (3), the lateral sensitivity contrast in proﬁle C–C′ is important, as is the sensitivity at crustal depths. If the lateral sensitivity
contrast of the crust is signiﬁcant in the C–C′ proﬁle, the C2 conductor might actually be located in an
insensitive region in the crust and uppermost mantle, which could cause the overturn image in proﬁleC–C′ to disappear. In other words, if the location below OBR has a poor sensitivity in proﬁle C–C′ andthe C2 conductor relocates below OBR, the overturn image would immediately vanish. Figure 12 shows
the lateral sensitivity contrast of the crust in proﬁle C–C′. To estimate Figure 12, we allocated a 10 Ω m
cube with sides of 20 km to between 20 and 40 km depth and moved the cube laterally along proﬁleC–C′. Our MT site distribution provides the best sensitivity at the C2 conductor location in proﬁle C–C′.

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN17

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 11. (continued)

For example, the Φyy[x(C2),60s] is at most one order of magnitude greater than the least sensitive region.
This strongly suggests that the C2 conductor below the back arc in proﬁle C–C′ is a required feature.
However, the TMFD and Φij(x,T) diminish below OBR because of the lack of observation sites, and the conductive region could expand to the low resolution region below OBR. Assuming that the conductive
region expands below OBR, the conductor in the C–C′ proﬁle may be a funnel shape rather than an overturn shape. To clarify this conductive image, we estimated the TMFD and Φij(x,T) in another proﬁle (Z–Z′)
across the C2 conductor between the depths of 20 and 40 km (Figure 13). The TMFD and Φij(x,T) in Z–Z′
show almost the same lateral sensitivity. None of the conductive bodies in the crust below OBR and fore
arc is imaged, and the conductive region also shows the overturn in Z–Z′. As a result, the middle to lower
crust in the eastern vicinity of the C2 conductor is not as conductive as C2, and the conductive overturnshape is convincingly delineated with our MT data.
The boxcar test is used to reinvert synthetic data calculated from a model in which a conductive or resistive
cube is embedded in a certain reference model. The reinverted model indicates a kind of impulse responseof the inversion process in terms of the reference model, and it allows us to evaluate how much theresistivity image obtained by the inversion blurs and the bias of the resistivity contrast. In other words,
the inversion process always outputs more or less blurred (smoothing) image than an actual contrast.
The boxcar test gives the degree of the blurring. We checked the resolution of the C1 conductor and theR2 resistor in Figures 8 and 9 by use of this boxcar test. The reference model and site location are the same
as for the Fréchet derivative calculation, and the employed periods of synthetic data are identical to the 12
used in the inversion. To make the synthetic data for the C1 conductor and R2 resistor, we embedded a
1 Ω m (three orders of magnitude conductive) and a 100 k Ω m (20 factors of magnitude resistive) cube withsides of 20 km into the reference model, respectively (Figure 14). We appended Gaussian random noisewith a 5% standard deviation to the synthetic impedance data. The Lagrange multiplier (100.5) and model

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN18

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 12. Total modiﬁed Fréchet derivative and upper bound of total phase change as a function of lateral location along
C–C′ proﬁle. Solid, dotted, dashed, and chain line represents the results at periods of 60, 480, 3840, and 30,720 s, respectively.
The resistivity perturbation is embedded in 20 to 40 km depth in C–C′ proﬁle as is shown by white rectangles in the bottom
ﬁgures. Thick line and triangle represent location of land and volcanoes.

covariance parameters (4τδi = 2.0) were ﬁxed at those of the ﬁnal model in section 4. The vertical proﬁle of
the reinverted model is shown in Figure 14 for comparison with the input model. The number of iterationsof each inversion was one, and the RMS misﬁt fell below 1.0. Because of the ﬁxed Lagrange multiplier,stage II of the Occam process was not performed. In the boxcar test for the C1 conductor, the outlineof the conductive region is almost comparable to that of the C1 conductor in Figure 9. The recovered
conductive region greater than 0.2 orders of magnitude assumes a width as large as 40 km. The recovered conductivity contrast is only up to 0.6 orders of magnitude more conductive than the reference
model despite the three orders of magnitude conductivity contrast in the input model. This boxcar test
reveals that the absolute resistivity obtained in this method is deﬁnitely biased toward resistive values.
Indeed, the resistivity of the C1 conductor obtained with the inversion (ca. 800 Ω m) is considerably moreresistive compared with lower crustal conductors seen in previous studies [e.g., Ogawa et al., 2001, 2014;
Mishina, 2009; Ichihara et al., 2014]. Figure 15 shows another boxcar test for the C1 conductor using 1,
10, and 100 Ω m cubes with the same volume. The recovered conductive regions using 1 and 10 Ω m cube
has a similar volume and resistivity contrast with that of the C1 conductor in proﬁle B–B′, while the result of
100 Ω m cube obviously shows a weak contrast resistivity image. Consequently, the actual volume of the C1
conductor is comparable with that of a cube with sides of 20 km, but the resistivity is about 10 Ω m or more
conductive. The boxcar test for the R2 resistor resulted in complete disappearance of the resistive block

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN19

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 13. Same as Figure 12 but along Z–Z′ line as is shown in the inset child map. Z–Z′ line is in east-west direction.

(Figure 14). The absence of an observation site above the R2 resistor as well as the boxcar test suggests that
the R2 resistor is poorly resolved.

6. Discussion
Our 3-D electrical resistivity model contains the following important features: (1) A conductive body has an
overturned-plume shape in proﬁle C–C′, with the conductive plume head reaching the lower crust beneaththe back-arc volcano Mt. Gassan. (2) The most prominent lower crust and uppermost mantle conductors
occur below Mt. Yakeishi on OBR. The resistivity of the conductive body in the lower crust is less than
10 Ω m, and the volume is comparable to that of a 20 km cube on a side. The conductor is continuous fromthe lower crust to the subducting slab surface. We next discuss the tectonic implications and geophysicalsigniﬁcance of these key features.
6.1. Overturn of Electrically Conductive Plume and Back-Arc Volcanism
Back-arc volcanism in NEJ has been interpreted as resulting from a mantle diapir originating from the subducting
slab at a depth of 150–200 km [e.g., Hasegawa and Nakajima, 2004]. Petrological studies have suggested thatICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN20

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 14. Boxcar test for (left) C1 conductor in B–B′ vertical proﬁle and for (right) R2 resistor in C–C′ proﬁle. Upper(map view) and middle (vertical view) ﬁgures are input models with conductive or resistive cube. The lower ﬁgure is arecovered model in vertical view.

dehydration of phlogopite [Tatsumi, 1989], phengite, and lawsonite [Pawley, 1994; Poli and Schmidt, 1995] could
yield such a mantle diapir below the back-arc region. A numerical study also demonstrated the onset of the
mantle diapir separating from the subducting slab at 150–200 km depth [Gerya and Yuen, 2003]. Hence, we
had expected a conductor isolated in the upper mantle beneath the back arc. However, our electrical resistivity
model revealed no such conductor isolated in the shallow upper mantle beneath the back arc. In contrast, themodel proﬁle in C–C′ (Figure 9) delineates a conductive plume overturned toward the back-arc direction, thehead of which reaches the lower crust below Mt. Gassan. This indicates that ﬂuids released from the subducting
slab at depths of 80–100 km move upward toward the back-arc direction and that overturning of the ﬂuids yields
the back-arc volcanism. The conductive plume at 40 to 80 km depth shows a weak resistivity contrast in themodel. As we demonstrated in the previous section, this resistivity would be biased toward the resistive values.

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN21

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 15. Boxcar test for C1 conductor resistivity. (a) The input model composed of the cube with 20 km on side in verticalview. The 1, 10, and 100 Ω m is allotted to the resistivity of the cube, respectively. The recovered model calculated from 1,
10, and 100 Ω m cube is shown in Figures 15b–15d, respectively, in vertical view.

The consequential issue is what enables the overturning of the upwelling ﬂuids in the mantle wedge. Currieand Hyndman [2006] advocated a conceptual model in which small-scale (ca. 100 km wavelength) poloidal
convections in the trench-perpendicular direction were depicted [see Figure 16 proposed by Currie andHyndman, 2006]. The small-scale poloidal convections are derived from thermal convections. We have interpreted that the overturn in our model corresponds to lateral ﬂow in the across-arc direction from upwelling
toward the downﬂow of the thermal convections. A numerical study of subduction zone dynamics realized
such an overturn in the mantle wedge [Davies and Stevenson, 1992, denoted as DS92 hereafter]. DS92
employs a two-dimensional Stokes ﬂow formulation with the Boussinesq approximation. When a buoyant
body with a buoyancy number greater than 106 emerges in the mantle wedge close to a subducting slab,
the buoyant body assumes an overturned shape (described as ﬂow reversal by DS92). The buoyancy number
is deﬁned as gd3Δρ/(κη), where d: the speciﬁc length of the buoyant body, Δρ: density perturbation,
g: gravity acceleration, κ: thermal diffusivity, and η: viscosity. DS92 estimated 2 wt. % water content
and 1130–1180°C temperature of melt to account for 106 of buoyancy number. Using these parameters,a laboratory experiment of peridotite melt conductivity yielded 0.9–2.0 S/m [Ni et al., 2011]. If we assume
that the resistivity of the mantle wedge conductive zone is comparable to the C1 and C2 conductors (viz.
less than 10 Ω m in resistivity as is shown in section 5), this estimation is consistent with our model. Notethat conductivity with 100% melt fraction should be used in the argument above because of Stokes ﬂow
condition. Volcanism on the hot ﬁnger track from Mt. Chokai to Mt. Kurikoma has migrated 70–80 km
eastward in the last 5 Ma [Kondo et al., 2004]. This eastward migration of volcanism might result from suppression of the overturning of upwelling ﬂuid/melt owing to a reduction of the buoyancy number of the
upwelling ﬂuid/melt. Honda and his colleagues represented the hot ﬁngers as a result of roll-like convection
ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN22

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

or Rayleigh-Taylor instability [Honda
and Saito, 2003; Zhu et al., 2009], i.e.,a two-dimensional and along-arc poloidal convection image. On the other
hand, the overturn in our resistivity
model indicates that there are smallscale poloidal convections in theacross-arc direction. These two aspects
suggest that the mantle wedge ﬂow
below NEJ assumes Rayleigh-Bénard
convection rather than Rayleigh-Taylor
instability and that the overturn inFigure 16. Conceptual model of mantle wedge ﬂow proposed by Currie
our model corresponds to the lateraland Hyndman [2006]. A slab subducts from the upper left to the lower
ﬂow in the across-arc direction fromright in this ﬁgure. Solid and dashed arrows indicate the viscous couplingof slab and thermal convection, respectively.
upwelling to downﬂow of Bénard cells.
Currie and Hyndman [2006] cited shorthorizontal wavelength (about 100 km)
seismic velocity variation in the western United States as evidence [Gao et al., 2004; Yang and Forsyth,2006] to advocate the presence of small-scale convections in the mantle wedge. Recent high-resolution
electrical resistivity [Meqbel et al., 2014] and seismic velocity studies [Schmandt and Humphreys, 2010;
Gilbert et al., 2012; Tian and Zhao, 2012; Jones et al., 2014] have revealed undulations of resistivity and seismic velocity distribution at the Moho and lithosphere-asthenosphere boundary in the western UnitedStates. Those results indicate the presence of Bénard cells beneath the back arc. The distance between
C2 and C3 is about 70 km, which is comparable to the wavelength of the seismic velocity and resistivity
variations in the western United States.
Allowing for upwelling of partial melt or inﬁltration of melt into grain boundaries or cracks, we should
consider the porous ﬂow governed by Darcy’s law [McKenzie, 1984; Katz et al., 2007; Cagnioncle et al.,2007; Wilson et al., 2014] rather than the Stokes ﬂow. However, it has yet to be studied whether numericalsimulation allowing for the porous ﬂow can cause the overturned ﬂow in the mantle wedge and underwhich condition the overturn can occur. Here we consider the porous ﬂow in the mantle wedge assuming the permeability distribution model proposed by Wada et al. [2011]. They modeled grain size distribution in the mantle wedge and converted the grain size distribution into a permeability distribution
by using the Blake-Kozeny-Carman equation [Bear, 1972]. The permeability distribution has a similar pattern to the overturn conﬁguration. When ﬂuid/melt emerges in a mantle wedge with this permeability
distribution calculated by Wada et al. [2011], Darcy’s law tells us that ﬂuid/melt moves up toward the backarc (see supporting information). Thus, the porous ﬂow could cause the overturned ﬂow. However, this heterogeneous permeability cannot explain the vertical conductive body below OBR (proﬁles A–A′ and B–B′).
Although the overturn of the conductive body is not prominent below Mt. Chokai (or the C3 conductor) in
proﬁle B–B′, the overturn likely takes place in the middle crust in this proﬁle. Seismic tomography studies
in NEJ [e.g., Huang et al., 2011; Nakajima et al., 2001] show that positive Vp and Vs perturbations at 40 km
depth below the back-arc region of our proﬁle B–B′ decrease abruptly to negative perturbations at 25 km
depth. Moreover, their vertical seismic wave velocity proﬁles evidently delineate the overturn image of slow
seismic wave velocity perturbation in the middle crust below the back arc in our proﬁle B–B′ location.
Hasegawa et al. [2008] depicted an overturn image of slow Vs on a vertical proﬁle passing the northern vicinity of Mt. Chokai in the across-arc direction. The inﬂection point of the slow Vs overturn is as deep as 20 km,
which is consistent with our resistivity model. Seismic wave velocity and resistivity images seem to corroborate the presence of overturn of ﬂuid/melt toward the back-arc direction in the middle crust.
6.2. On the Lower Crust Conductor
6.2.1. Quantitative Interpretation of the Lower Crust Resistivity
The most distinctive conductive body, C1, is observed at depths of 20 to 50 km below OBR in proﬁle B–B′. The
boxcar test in section 5 suggests that the resistivity of C1 is less than 10 Ω m. We calculated the volume fraction of ﬂuid/melt by which a resistivity less than 10 Ω m in the lower crust can be explained. In the calculation,
ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN23

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

we used 1 GPa and 800°C of pressure
and temperature, respectively. The mantle xenoliths sampled at Ichinomegata
indicate that the equilibrated temperature in the uppermost mantle is about
800°C [Takahashi, 1980, 1986]. Because
the water-saturated solidus of felsic andmaﬁc rocks at 1 GPa is 600–700°C [e.g.,
Wolf and Wyllie, 1994, and references
therein], melt could be present in thelower crust. Melt conductivity in thelower crust, σ melt, is calculated fromthe equation for amphibolite (obsidian)melt obtained by Gaillard [2004]:
σ melt ¼ σ 0 exp½ðE a þ 20PÞ=ðRT Þ (8)Figure 17. Bulk resistivity of the composites of amphibolite and hydrous
melt at 1 GPa in pressure and 800°C in temperature. Hydrous melt resistivity depends on melt water content. Contour line unit is in log(Ω m). The
amphibolite and melt conductivity data are from Fuji-ta et al. [2011] andGaillard [2004], respectively. Hashin-Shtrikman upper bound is used to
calculate bulk resistivity.

σ 0 ¼ 78:9 lnðW Þ þ 754(9)

E a ¼ 2925 lnðW Þ þ 64132(10)

where P is pressure in MPa, T istemperature in Kelvin, R is the gas
constant in J mole1 K1, and W is thewater content in wt. % (note that
Gaillard’s original paper contains atypographical omission). The conductivityof amphibolite isσ solid ¼ σ 0 exp ½E a =ðkT Þ;(11)where σ 0 = 0.8385 S/m, Ea = 0.4874 eV, andk is the Boltzman constant in eV/K. The
data are from BC-3 shown by Fuji-ta et al.
[2011; Fuji-ta, personal communication].
The Hashin-Shtrikman upper bound[Hashin and Shtrikman, 1962] is used to
calculate bulk resistivity.

Figure 18. Bulk resistivity of the composites of amphibolite, watersaturated melt, and aqueous sodium chloride with 17 wt. % at 1 GPa inpressure and 800°C in temperature. Conductivity of aqueous sodium
chloride with 17 wt. % is from Shimojuku et al. [2014]. Other data andrepresentation are the same with Figure 17.

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPANWe consider two cases. The ﬁrst is thatall saline ﬂuid dissolves in melt, and thebulk conﬁguration is composed of lower
crustal rock and melt. Figure 17 shows
the relation between volume fractionof melt and bulk resistivity of the meltrock composite. The melt fractionneeded to account for the 10 Ω m inbulk resistivity is 1.5 vol. %. This estimate is consistent with 1–2 vol. % at
40 km depth of ﬂuid/melt calculated
from the seismological data [Nakajima
et al., 2005]. A 2 vol. % (upper bound)melt/ﬂuid fraction derived from seismological data corresponds to a resistivity
in the lower crust of 7–8 Ω m. The secondcase is the coexistence of melt and saline
ﬂuid. That is, melt is saturated with water24

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Figure 19. Thermal model beneath NEJ by Iwamori [2007] superimposed on the resistivity proﬁle in B–B′. White line shows
isotherms. Unit is in °C. Slab surface temperature at 100 km depth is about 600°C and is much cooler than Cascadia. See text
for detail.

and excess water exists as saline ﬂuid. The solubility of water in rhyolitic melts is 10 wt. % using the equation of
Tamic et al. [2001]. The conductivity of melt with 10 wt. % water content is 8.658 S/m. Fluid inclusions [e.g.,
Roedder, 1984], ﬂuids dehydrated from the subducting slab [e.g., Manning, 2004], and ﬂuids sampled in the deep
crust [e.g., Bucher and Stober, 2010] are all predominantly composed of aqueous sodium chloride and calcium
chloride solutions. Because calcium cation mobility is negligible [Gaillard, 2004; Shimojuku et al., 2014], the major
charge carriers in aqueous ﬂuid and melt under lower crustal conditions are sodium and chloride ions. Thus, we
assume that the saline ﬂuid is aqueous sodium chloride solution. For simplicity, we do not allow for the partition
coefﬁcients of sodium and chlorine between melt and saline ﬂuid phases, but rather employ the conductivity ofsaline ﬂuids estimated by Shimojuku et al. [2014]. Figure 18 shows the bulk resistivity of water-saturated melt,
saline ﬂuid with 17 wt. % salinity (20 S/m), and amphibolite. The Hashin-Shtrikman upper bound for multiple
phases is reduced toσ bulk ¼ σ 1 þ A=½1  A=ð3σ 1 Þ
(12)
A¼

N
X

fi
;
1=
ð
σ

σ
Þ þ 1=ð3σ 1 Þ
i
1
i¼2

(13)

where σ 1 > … > σ i > σ i + 1 > … > σ N. σ bulk is the bulk conductivity and fi is the volume fraction of i’s medium[Park and Ducea, 2003]. The end member of the second case is that melt volume is approximately zero. In thiscase, the total amount of ﬂuid/melt fraction is also a minimum. To account for 10 Ω m in the lower crust, only
0.7 vol. % of ﬂuid/melt is required. For consistency with the upper bound volume fraction estimated by theseismic data (2 vol. %), the lower crustal conductor would need to be more than 0.3 Ω m.
ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN25

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

6.2.2. Fluid Transport Pathway From Slab Surface to Lower Crust
The conductive zone extends from the lower crust conductor to the slab surface at 80–100 km depth below
OBR. This indicates that the ﬂuid/melt pathway, that is, ﬂuid accompanying the subducting slab is“segregated” from the slab at 80–100 km depth and rises to the lower crust conductor below OBR.
Because our resistivity model has insufﬁcient resolution in the mantle wedge deeper than 100 km, the model canneither resolve the depth down to which the conductive zone extends nor determine the source location of
ﬂuid/melt. Based on a 2-D model of Cascadia, McGary et al. [2014] revealed a ﬂuid/melt pathway similar with thatin proﬁle B–B′ of our model. They interpreted that the primary source existed at the slab surface at about 100 km
depth. However, this interpretation cannot be adapted to the subduction zone beneath NEJ because the age ofthe subducting plate beneath Cascadia (up to 10 Ma) is much younger than that beneath NEJ (about 130 Ma at the
trench). The temperature at the slab surface at 100 km depth is as low as 600°C beneath NEJ [Iwamori, 1998, 2007],
which is much cooler than the 1150°C of the Cascadia geothermal setting (see Figure 19). According to seismic
tomography [e.g., Zhao et al., 1992; Nakajima et al., 2001], laboratory experiments [e.g., Iwamori, 1998, 2007;
Schmidt and Poli, 1998], and ocean bottom magnetotellurics [Toh et al., 2006], the dehydration of dragged hydrous
peridotite and of subducting MORB beneath NEJ likely takes place at depths as deep as 150 to 200 km. The connection between the dehydration ﬁeld at 150 to 200 km depth and the segregation of ﬂuid/melt at 80–100 km
depth remain unresolved because of low resolution below 100 km depth in our model. The resistivity contrastof the plume at 40 to 80 km depth is weak (60–100 Ω m in the ﬁnal model). However, this resistivity would be
biased toward the resistive values, as we demonstrated in the boxcar test.
As it reached the lower crust, the ﬂuid/melt remains in the lower crust. This feature is consistently observed in
subduction zones of the world [South America: Brasse and Eydam, 2008; North America: Patro and Egbert, 2008;
Rippe et al., 2013; McGary et al., 2014; Meqbel et al., 2014; New Zealand: Heise et al., 2007; Wannamaker et al.,2009; see also Shankland and Ander, 1983]. The reason why ﬂuid/melt remains in the lower crust and uppermost
mantle is possibly due to the formation of impermeable cap layers owing to increasing dihedral angle and/or
silica precipitation in the saline ﬂuid. The dihedral angle hypothesis has been proposed by many researchers[e.g., Hyndman and Shearer, 1989]. Indeed, laboratory experiments have revealed that the dihedral angle inanorthite aggregates with aqueous ﬂuid was more than 60° below 1 GPa pressure, whereas it decreased below
60° with increasing pressure [Yoshino et al., 2002]. Silica precipitation also contributes to interrupting the rise of
ﬂuid/melt in the crust. Such silica precipitation is probably stimulated as saline ﬂuid rises to the middle or upper
crust because silica solubility in saline ﬂuid decreases with the exchange of concomitant solute cation fromsodium to calcium [Shumlovich et al., 2006]. The exchange of solute cation could take place when saline ﬂuid
rises to the middle or upper crust with lithological facies change. Furthermore, silica solubility in saline ﬂuid
essentially decreases as pressure and temperature diminish [e.g., Manning, 1994].
The location of the C1 conductor is below Mts. Yakeishi and Kurikoma and the hypocentral region of the 2008
Iwate-Miyagi Nairiku Earthquake (M 7.2) [e.g., Okada et al., 2010, 2012, 2014]. Three isolated conductive
bodies (less than 10 Ω m in resistivity and 10 to 20 km in width and depth) were observed in this region by
Ichihara et al. [2014] and Ogawa et al. [2014]. The C1 conductor of our study corresponds to the merger of
two conductors below Mts. Yakeishi and Kurikoma.

7. Conclusion
Our 3-D model of the subduction zone below NEJ contains a vertically continuous conductive body from thelower crust to the slab surface below Mts. Yakeishi and Kurikoma on OBR. The conductive body indicates aﬂuid/melt pathway from the slab surface to lower crust. The lower crust conductor is 10 Ω m or more conductive in resistivity. The 0.7 vol. % at minimum of ﬂuid/melt is required to account for 10 Ω m resistivity in thelower crust, and the resistivity of the lower crust conductor greater than 0.3 Ω m is in agreement with the
1 to 2 vol. % of ﬂuid/melt volume fraction derived from seismological data [Nakajima et al., 2005].
Other resistivity proﬁles in the across-arc direction show that the conductive body segregates from the subducting slab at 80–100 km depth and takes an overturned form toward the back-arc direction. The head of
the conducting body reaches the lower crust just below Mt. Gassan, which is one of the prominent backarc volcanoes. Hence, the back-arc volcanism is likely not caused by a mantle diapir but rather by
ﬂuid/melt overturn rising toward the back-arc direction. A numerical simulation with Stokes ﬂow formulation
ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN26

Journal of Geophysical Research: Solid Earth

Acknowledgments
Hiroshi Ichihara (Japan Agency for
Marine-Earth Science and Technology),
Hiromi Fukino, Izumi Sato, and Oliver
Hartkorn (Tokyo Institute of Technology),
Sadato Ueki, Yasuo Yabe, Takuro Iida, andKeisuke Yoshida (Tohoku University), andNobutaka Hozumi (Kyoto University)
helped greatly with the ﬁeldwork. Special
thanks are extended to Eiichi Takahashi(Tokyo Institute of Technology) for his
continuous encouragement. Kiyoshi
Fuji-ta (Osaka University) kindly allowed
us to show the result of the laboratory
experiments. Hikaru Iwamori (Japan
Agency for Marine-Earth Science andTechnology), Junichi Nakajima, Tomomi
Okada, Toru Matsuzawa, and MichihikoNakamura (Tohoku University), andHiroshi Sakuma (National Institute forMaterials Science) provided us withfruitful discussions. Magnetotelluric
equipments with the gimballed systemwere provided by the University ofTokyo. The Geomagnetic Chart in 2010
used in this paper is available at
Geospatial Authority of Japan. We usedseismic hypocenter location data
determined by the Japan Meteorological
Agency, the National Research Institutefor Earth Science and Disaster Prevention,
the Geoinformation Authority of Japan,
Hokkaido University, Hirosaki University,
Tohoku University, the University ofTokyo, Nagoya University, Kyoto
University, Kochi University, Kyushu
University, Kagoshima University,
Japan Agency for Marine-Earth
Science and Technology, the National
Institute of Advanced Industrial
Science and Technology, the Tokyo
Metropolitan Government, the Aomori
Prefectural Government, the Shizuoka
Prefectural Government, the Hot
Springs Research Institute of
Kanagawa Prefecture, and Yokohama
City. We partly carried out the numerical calculation using the computer
systems of the Earthquake andVolcano Information Center of the
Earthquake Research Institute, theUniversity of Tokyo, and the TSUBAME
computer service, Tokyo Institute ofTechnology. Alan Chave (Woods Hole
Oceanographic Institution) kindly provided the bounded inﬂuence remotereference processing code. Most partsof the ﬁgures were created with the
Generic Mapping Tool [Wessel andSmith, 1998]. Constructive comments
by the associate editor and two anonymous reviewers greatly improved themanuscript. This work was supported
ﬁnancially by the Japan Society for the
Promotion of Science KAKENHI,
21109003, 20403006, 23403004, and26109006.

ICHIKI ET AL.

10.1002/2015JB012028

showed that the overturn or ﬂow reverse toward the back arc takes place in the case of a buoyancy number
greater than 106 [Davies and Stevenson, 1992]. The conductive overturn distribution indicates the presence ofsmall-scale poloidal convections in the across-arc direction. Allowing for the presence of poloidal convectionsin the along-arc direction as well as those in the across-arc direction, the mantle wedge ﬂow below NEJ likely
assumes Rayleigh-Bénard convection rather than Rayleigh-Taylor instability.

References
Aizawa, K., et al. (2014), Three-dimensional resistivity structure and magma plumbing system of the Kirishima Volcanoes as inferred frombroadband magnetotelluric data, J. Geophys. Res. Solid Earth, 119, 198–215, doi:10.1002/2013JB010682.
Bear, J. (1972), Dynamics of Fluids in Porous Media, 764 pp., American Elsevier Publ. Com., New York.
Bibby, H. M., T. G. Caldwell, and C. Brown (2005), Determinable and non-determinable parameters of galvanic distortion in magnetotellurics,
Geophys. J. Int., 163, 915–930, doi:10.1111/j.1365-246X.2005.02779.x.
Boerner, D., and G. F. West (1989), Fréchet derivatives and single scattering theory, Geophys. J. Int., 98, 385–390, doi:10.1111/
j.1365-246X.1989.tb03359.x.
Brasse, H., and D. Eydam (2008), Electrical conductivity beneath the Bolivian Orocline and its relation to subduction processes at the SouthAmerican continental margin, J. Geophys. Res., 113, B07109, doi:10.1029/2007JB005142.
Bucher, K., and I. Stober (2010), Fluids in the upper continental crust, Geoﬂuids, 10, 241–253, doi:10.1111/j.1468-8123.2010.00279.x.
Cagnioncle, A.-M., E. M. Parmentier, and L. T. Elkins-Tanton (2007), Effect of solid ﬂow above a subducting slab on water distribution andmelting at convergent plate boundaries, J. Geophys. Res., 112, B09402, doi:10.1029/2007JB004934.
Caldwell, T. G., H. M. Bibby, and C. Brown (2004), The magnetotelluric phase tensor, Geophys. J. Int., 158, 457–469, doi:10.1111/
j.1365-246X.2004.02281.x.
Chave, A. D., and D. J. Thomson (2004), Bounded inﬂuence magnetotelluric response function estimation, Geophys. J. Int., 157, 988–1006,
doi:10.1111/j.1365-246X.2004.02203.x.
Constable, S. C., R. L. Parker, and C. G. Constable (1987), Occam’s inversion: A practical algorithm for generating smooth models fromelectromagnetic sounding data, Geophysics, 52, 289–300, doi:10.1190/1.1442303.
Cruciani, C., E. Carminati, and C. Doglioni (2005), Slab dip vs. lithosphere age: No direct function, Earth Planet. Sci. Lett., 238, 298–310,
doi:10.1016/j.epsl.2005.07.025.
Currie, C. A., and R. D. Hyndman (2006), The thermal structure of subduction zone back arcs, J. Geophys. Res., 111, B08404, doi:10.1029/
2005JB004024.
Davies, J. H., and D. J. Stevenson (1992), Physical model of source region of subduction zone volcanics, J. Geophys. Res., 97, 2037–2070,
doi:10.1029/91JB02571.
DeMets, C., R. G. Gordon, and D. F. Argus (2010), Geologically current plate motions, Geophys. J. Int., 181, 1–80, doi:10.1111/j.1365-246X.2009.04491.x.
Fuji-ta, K., T. Katsura, M. Ichiki, T. Matsuzaki, and T. Kobayashi (2011), Variations in electrical conductivity of rocks above metamorphic
conditions, Tectonophysics, 504, 116–121, doi:10.1016/j.tecto.2011.03.008.
Gaillard, F. (2004), Laboratory measurements of electrical conductivity of hydrous and dry silicic melts under pressure, Earth Planet. Sci. Lett.,
218, 215–228, doi:10.1016/S0012-821X(03)00639-3.
Gao, W., S. P. Grand, W. S. Baldridge, D. Wilson, M. West, J. F. Ni, and R. Aster (2004), Upper mantle convection beneath the central Rio Grande
rift imaged by P and S wave tomography, J. Geophys. Res., 109, B03305, doi:10.1029/2003JB002743.
Gerya, T. V., and D. A. Yuen (2003), Rayleigh-Taylor instabilities from hydration and melting propel ‘cold plumes’ at subduction zones, Earth
Planet. Sci. Lett., 212, 47–62, doi:10.1016/S0012-821X(03)00265-6.
Gilbert, H., Y. Yang, D. W. Forsyth, C. H. Jones, T. J. Owens, G. Zandt, and J. C. Stachnik (2012), Imaging lithospheric foundering in the structureof the Sierra Nevada, Geosphere, 8, 1310–1330, doi:10.1130/GES00790.1.
Green, D. H. (1973), Experimental melting studies on a model upper mantle composition at high pressure under water-saturated andwater-undersaturated conditions, Earth Planet. Sci. Lett., 19, 37–53, doi:10.1016/0012-821X(73)90176-3.
Hasegawa, A., and A. Yamamoto (1994), Deep, low-frequency microearthquakes in or around seismic low-velocity zones beneath active
volcanoes in northeastern Japan, Tectonophysics, 233, 233–252, doi:10.1016/0040-1951(94)90243-7.
Hasegawa, A., and J. Nakajima (2004), Geophysical constraints on slab subduction and arc magmatism, in The State of the Planet: Frontiersand Challenges, Geophys. Monogr. Ser., vol. 150, edited by R. S. J. Sparks and C. J. Hawkesworth, pp. 81–94, AGU, Washington, D. C.
Hasegawa, A., J. Nakajima, N. Umino, and S. Miura (2005), Deep structure of the northeastern Japan arc and its implications for crustal
deformation and shallow seismic activity, Tectonophysics, 403, 59–75, doi:10.1016/j.tecto.2005.03.018.
Hasegawa, A., J. Nakajima, S. Kita, Y. Tsuji, K. Nii, T. Okada, T. Matsuzawa, and D. Zhao (2008), Transportation of H2O in the NE Japan
Subduction Zone as inferred from seismic observations [in Japanese with English abstract and ﬁgure caption], J. Geogr. (Chigaku Zasshi),
117, 59–75, doi:10.5026/jgeography.117.59.
Hashin, Z., and S. Shtrikman (1962), A variational approach to theory of effective magnetic permeability of multiphase materials, J. App. Phys.,
33, 3125–3131, doi:10.1063/1.1728579.
Heise, W., H. M. Bibby, T. G. Caldwell, S. C. Bannister, Y. Ogawa, S. Takakura, and T. Uchida (2007), Melt distribution beneath a young
continental rift: The Taupo Volcanic Zone, New Zealand, Geophys. Res. Lett., 34, L14313, doi:10.1029/2007GL029629.
Honda, S., and M. Saito (2003), Small-scale convection under the back-arc occurring in the low viscosity wedge, Earth Planet. Sci. Lett., 216,
703–715, doi:10.1016/S0012-821X(03)00537-5.
Honda, S., and T. Yoshida (2005), Application of the model of small-scale convection under the island arc to the NE Honshu subduction zone,
Geochem. Geophys. Geosyst., 6, Q01002, doi:10.1029/2004GC000785.
Huang, Z., D. Zhao, and L. Wang (2011), Seismic heterogeneity and anisotropy of the Honshu arc from the Japan Trench to the Japan Sea,
Geophys. J. Int., 184, 1428–1444, doi:10.1111/j.1365-246X.2011.04934.x.
Hyndman, R. D., and P. M. Shearer (1989), Water in the lower continental crust: Modeling magnetotelluric and seismic reﬂection results,
Geophys. J. Int., 98, 343–365, doi:10.1111/j.1365-246X.1989.tb03357.x.
Ichihara, H., S. Sakanaka, M. Mishina, M. Uyeshima, T. Nishitani, Y. Ogawa, Y. Yamaya, T. Mogi, K. Amita, and T. Miura (2014), A 3-D electricalresistivity model beneath the focal zone of the 2008 Iwate-Miyagi Nairiku earthquake (M 7.2), Earth Planets Space, 66, 50, doi:10.1186/
1880-5981-66-50.
Iwamori, H. (1998), Transportation of H2O and melting in subduction zones, Earth Planet. Sci. Lett., 160, 65–80, doi:10.1016/S0012-821X(98)00080-6.

ELECTRICAL IMAGE BENEATH NE JAPAN27

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Iwamori, H. (2007), Transportation of H2O beneath the Japan arcs and its implications for global water circulation, Chem. Geol., 239, 182–198,
doi:10.1016/j.chemgeo.2006.08.011.
Jiracek, G. R. (1990), Near-surface and topographic distortions in electromagnetic induction, Surv. Geophys., 11, 163–203, doi:10.1007/
BF01901659.
Jödicke, H. (1992), Water and graphite in the earth’s crust—An approach to interpretation of conductivity models, Surv. Geophys., 13,
381–407, doi:10.1007/BF01903484.
Jones, C. H., H. Reeg, G. Zandt, H. Gilbert, T. J. Owens, and J. Stachnik (2014), P-wave tomography of potential convective downwellings andtheir source regions, Sierra Nevada, California, Geosphere, 10, 505–533, doi:10.1130/GES00961.1.
Kanda, W., and Y. Ogawa (2014), Three-dimensional electromagnetic imaging of ﬂuids and melts beneath the NE Japan arc revisited by using
geomagnetic transfer function data, Earth Planets Space, 66, 39, doi:10.1186/1880-5981-66-39.
Karato, S. (1990), The role of hydrogen in the electrical conductivity of the upper mantle, Nature, 347, 272–273, doi:10.1038/347272a0.
Katz, R. F., M. G. Knepley, B. Smith, M. Spiegelman, and E. T. Coon (2007), Numerical simulation of geodynamic processes with the portable
extensible toolkit for scientiﬁc computation, Phys. Earth Planet. Inter., 163, 52–68, doi:10.1016/j.pepi.2007.04.016.
Kawakatsu, H., P. Kumar, Y. Takei, M. Shinohara, T. Kanazawa, E. Araki, and K. Suyehiro (2009), Seismic evidence for sharp lithosphere-asthenosphereboundaries of oceanic plates, Science, 24(324), 499–502, doi:10.1126/science.1169499.
Kelbert, A., G. D. Egbert, and C. deGroot-Hedlin (2012), Crust and upper mantle electrical conductivity beneath the Yellowstone Hotspot
Track, Geology, 40, 447–450, doi:10.1130/G32655.1.
Kita, S., T. Okada, A. Hasegawa, J. Nakajima, and T. Matsuzawa (2010), Anomalous deepening of a seismic belt in the upper-plane of the
double seismic zone in the Paciﬁc slab beneath the Hokkaido corner: Possible evidence for thermal shielding caused by subducted
forearc crust materials, Earth Planet. Sci. Lett., 290, 415–426, doi:10.1016/j.epsl.2009.12.038.
Kondo, H., K. Tanaka, Y. Mizouchi, and A. Ninomiya (2004), Long-term changes in distribution and chemistry of middle Miocene to
Quaternary volcanism in the Chokai-Kurikoma area across the Northeast Japan Arc, Island Arc, 13, 18–46, doi:10.1111/j.14401738.2003.00417.x.
Manning, C. E. (1994), The solubility of quartz in H2O in the lower crust and upper-mantle, Geochim. Cosmochim. Acta, 58, 4831–4839,
doi:10.1016/0016-7037(94)90214-3.
Manning, C. E. (2004), The chemistry of subduction-zone ﬂuids, Earth Planet. Sci. Lett., 223, 1–16, doi:10.1016/j.epsl.2004.04.030.
Matsuno, T., et al. (2010), Upper mantle electrical resistivity structure beneath the central Mariana subduction system, Geochem. Geophys.
Geosyst., 11, Q09003, doi:10.1029/2010GC003101.
McGary, R. S., R. L. Evans, P. E. Wannamaker, J. Elsenbeck, and S. Rondenay (2014), Pathway from subducting slab to surface for melt andﬂuids beneath Mount Rainier, Nature, 511, 338–341, doi:10.1038/nature13493.
McKenzie, D. (1984), The generation and compaction of partially molten rock, J. Petrol., 25, 713–765.
Mei, S., and K. L. Kohlstedt (2000a), Inﬂuence of water on plastic deformation of olivine aggregates 1. Diffusion creep regime, J. Geophys. Res.,
105, 21,457–21,469, doi:10.1029/2000JB900179.
Mei, S., and K. L. Kohlstedt (2000b), Inﬂuence of water on plastic deformation of olivine aggregates 1. Dislocation creep regime, J. Geophys.
Res., 105, 21,471–21,481, doi:10.1029/2000JB900180.
Meqbel, N. M., G. D. Egbert, P. E. Wannamaker, A. Kelbert, and A. Schultz (2014), Deep electrical resistivity structure of the northwestern U.S.
derived from 3-D inversion of USArray magnetotelluric data, Earth Planet. Sci. Lett., 402, 290–304, doi:10.1016/j.epsl.2013.12.026.
Mishina, M. (2009), Distribution of crustal ﬂuids in Northeast Japan as inferred from resistivity surveys, Gondwana Res., 16, 563–571,
doi:10.1016/j.gr.2009.02.005.
Mitsuhata, Y., Y. Ogawa, M. Mishina, T. Kono, T. Yokokura, and T. Uchida (2001), Electromagnetic heterogeneity of the seismogenic region of1962 M6.5 northern Miyagi earthquake, northeastern Japan, Geophys. Res. Lett., 28, 4371–4374, doi:10.1029/2001GL013079.
Nakajima, J., and A. Hasegawa (2006), Anomalous low-velocity zone and linear alignment of seismicity along it in the subducted Paciﬁc slab
beneath Kanto, Japan: Reactivation of subducted fracture zone?, Geophys. Res. Lett., 33, L16309, doi:10.1029/2006GL026773.
Nakajima, J., T. Matsuzawa, A. Hasegawa, and D. Zhao (2001), Three-dimensional structure of Vp, Vs, and Vp/Vs beneath northeastern Japan:
Implications for arc magmatism and ﬂuids, J. Geophys. Res., 106, 21,843–21,857, doi:10.1029/2000JB000008.
Nakajima, J., Y. Takei, and A. Hasegawa (2005), Quantitative analysis of the inclined low-velocity zone in the mantle wedge of northeasternJapan: A systematic change of melt-ﬁlled pore shapes with depth and its implications for melt migration, Earth Planet. Sci. Lett., 234, 59–70,
doi:10.1016/j.epsl.2005.02.033.
Ni, H., H. Keppler, and H. Behrens (2011), Electrical conductivity of hydrous basaltic melts: Implications for partial melting in the upper
mantle, Contrib. Mineral. Petrol., 162, 637–650, doi:10.1007/s00410-011-0617-4.
Ogawa, Y., T. Yukutake, and H. Utada (1986), Two-dimensional modeling of resistivity structure beneath the Tohoku district, northern Honshu
of Japan, by a ﬁnite element method, J. Geomag. Geoelectr., 38, 45–79, doi:10.5636/jgg.38.45.
Ogawa, Y., et al. (2001), Magnetotelluric imaging of ﬂuids in intraplate earthquake zones, NE Japan back arc, Geophys. Res. Lett., 28,
3741–3744, doi:10.1029/2001GL013269.
Ogawa, Y., M. Ichiki, W. Kanda, M. Mishina, and K. Asamori (2014), Three-dimensional ﬂuid distribution around Naruko volcano, inferred fromwide-band magnetotelluric measurement, Earth Planets Space, 66, 158, doi:10.1186/s40623-014-0158-y.
Okada, T., N. Umino, and A. Hasegawa (2010), Deep structure of the Ou mountain range strain concentration zone and the focal area of the2008 Iwate-Miyagi Nairiku earthquake, NE Japan—Seismogenesis related with magma and crustal ﬂuid, Earth Planets Space, 62, 347–352,
doi:10.5047/eps.2009.11.005.
Okada, T., N. Umino, A. Hasegawa, and Group for the aftershock observations of the Iwate-Miyagi Nairiku Earthquake in 2008 (2012),
Hypocenter distribution and heterogeneous seismic velocity structure in and around the focal area of the 2008 Iwate-Miyagi Nairiku
Earthquake, NE Japan—Possible seismological evidence for a ﬂuid driven compressional inversion earthquake, Earth Planets Space, 64,
717–728, doi:10.5047/eps.2012.03.005.
Okada, T., T. Matsuzawa, J. Nakajima, N. Uchida, M. Yamamoto, S. Hori, T. Kono, T. Nakayama, S. Hirahara, and A. Hasegawa (2014), Seismic
velocity structure in and around the Naruko volcano, NE Japan, and its implications for volcanic and seismic activities, Earth Planets Space,66, 114, doi:10.1186/1880-5981-66-114.
Oldenburg, D. W. (1979), One-dimensional inversion of natural source magnetotelluric observations, Geophysics, 44, 1218–1244, doi:10.1190/
1.1441004.
Park, S. K., and M. N. Ducea (2003), Can in situ measurements of mantle electrical conductivity be used to infer properties of partial melts?,
J. Geophys. Res., 108(B5), 2270, doi:10.1029/2002JB001899.
Patro, P. K., and G. D. Egbert (2008), Regional conductivity structure of Cascadia: Preliminary results from 3D inversion of USArray transportable array magnetotelluric data, Geophys. Res. Lett., 35, L20311, doi:10.1029/2008GL035326.

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN28

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012028

Pawley, A. R. (1994), The pressure and temperature stability limits of lawsonite—Implications for H2O recycling in subduction zones, Contrib.
Mineral. Petrol., 118, 99–108, doi:10.1007/BF00310614.
Peacock, S. A. (1990), Fluid processes in subduction zones, Science, 248, 329–337, doi:10.1126/science.248.4953.329.
Poli, S., and M. W. Schmidt (1995), H2O transport and release in subduction zones—Experimental constraints on basaltic and andesitic
systems, J. Geophys. Res., 100, 22,299–22,314, doi:10.1029/95JB01570.
Rippe, D., M. J. Unsworth, and C. A. Currie (2013), Magnetotelluric constraints on the ﬂuid content in the upper mantle beneath the southern
Canadian Cordillera: Implications for rheology, J. Geophys. Res. Solid Earth, 118, 5601–5624, doi:10.1002/jgrb.50255.
Roedder, E. (1984), Fluid inclusions, in Reviews in Mineralogy, vol. 12, 646 pp., Miner. Soc. Am., Chantilly, Va.
Schmandt, B., and E. Humphreys (2010), Seismic heterogeneity and small-scale convection in the southern California upper mantle,
Geochem. Geophys. Geosyst., 11, Q05004, doi:10.1029/2010GC003042.
Schmidt, M. W., and S. Poli (1998), Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation,
Earth Planet. Sci. Lett., 163, 361–379, doi:10.1016/S0012-821X(98)00142-3.
Shankland, T. J., and M. E. Ander (1983), Electrical conductivity, temperatures and ﬂuids in the lower crust, J. Geophys. Res., 88, 9475–9484,
doi:10.1029/JB088iB11p09475.
Shimojuku, A., T. Yoshino, and D. Yamazaki (2014), Electrical conductivity of brine-bearing quartzite at 1 GPa: Implications for ﬂuid content
and salinity of the crust, Earth Planets Space, 66, 2, doi:10.1186/1880-5981-66-2.
Shumlovich, K. I., B. W. D. Yardley, and C. M. Graham (2006), Solubility of quartz in crustal ﬂuids: Experiments and general equations for salt
solutions and H2O–CO2 mixtures at 400–800°C and 0.1–0.9 GPa, Geoﬂuids, 6, 154–167, doi:10.1111/j.1468-8123.2006.00140.x.
Siripunvaraporn, W., and G. D. Egbert (2000), An efﬁcient data-subspace inversion method for 2-D magnetotelluric data, Geophysics, 65,
791–803, doi:10.1190/1.14444778.
Siripunvaraporn, W., and G. D. Egbert (2009), WSINV3DMT: Vertical magnetic ﬁeld transfer function inversion and parallel implementation,
Phys. Earth Planet. Int., 173, 317–329, doi:10.1016/j.pepi.2009.01.013.
Siripunvaraporn, W., G. D. Egbert, Y. Lenbury, and M. Uyeshima (2005), Three-dimensional magnetotelluric inversion: Data-space method,
Phys. Earth Planet. Inter., 150, 3–14, doi:10.1016/j.pepi.2004.08.023.
Takahashi, E. (1980), Thermal history of lehrzolite zenoliths from the Ichinomegata crater, Oga peninsula, northeast Japan, Geochim.
Cosmochim. Acta, 44, 1643–1658, doi:10.1016/0016-7037(80)90217-3.
Takahashi, E. (1986), Genesis of calk-alcali andesite magma in a hydrous mantle crust boundary—Petrology of lehrzolite xenoliths fromthe Ichinomegata crater, Oga peninsula: Northeast Japan. 2, J. Vol. Geotherm. Res., 29, 355–395, doi:10.1016/0377-0273(86)90051-X.
Tamic, N., H. Behrens, and F. Holtz (2001), The solubility of H2O and CO2 in rhyolitic melts in equilibrium with a mixed CO2-H2O ﬂuid phase,
Chem. Geol., 174, 333–347, doi:10.1016/S0009-2541(00)00324-7.
Tamura, Y., Y. Tatsumi, D. Zhao, Y. Kido, and H. Shukuno (2002), Hot ﬁngers in the mantle wedge: New insights into magma genesis in
subduction zones, Earth Planet. Sci. Lett., 197, 105–116, doi:10.1016/S0012-821X(02)00465-X.
Tatsumi, Y. (1989), Migration of ﬂuid phases and genesis of basalt magmas in subduction zones, J. Geophys. Res., 94, 4697–4707, doi:10.1029/
JB094iB04p04697.
Tian, Y., and D. Zhao (2012), P-wave tomography of the western United States: Insight into the Yellowstone hotspot and the Juan de Fuca
slab, Phys. Earth Planet. Inter., 200–201, 72–84, doi:10.1016/j.pepi.2012.04.004.
Toh, H., K. Baba, M. Ichiki, T. Motobayashi, Y. Ogawa, M. Mishina, and I. Takahashi (2006), Two-dimensional electrical section beneath theeastern margin of Japan Sea, Geophys. Res. Lett., 33, L22309, doi:10.1029/2006GL027435.
Torres-Verdin, C., and F. X. Bostick (1992), Principles of spatial surface electric ﬁeld ﬁltering in magnetotellurics: Electromagnetic array
proﬁling (EMAP), Geophysics, 57, 603–622, doi:10.1190/1.1443273.
Tsuchiya, N., and S. Kanisawa (1994), Early Cretaceous Sr-rich silicic magmatism by slab melting in the Kitakami Mountains, northeast Japan,
J. Geophys. Res., 99, 22,205–22,220, doi:10.1029/94JB00458.
Uchida, T. (1993), Smooth 2-D inversion for magnetotelluric data based on statistical criterion ABIC, J. Geomag. Geoelectr., 45, 841–858,
doi:10.5636/jgg.45.841.
Utada, H. (1987), A direct inversion method for two-dimensional modeling in the geomagnetic induction problem PhD thesis, Graduate
School of Science, Univ. of Tokyo, Japan.
Uyeshima, M., et al. (2002), Network-MT survey in Japan to determine nation-wide deep electrical conductivity structure, in Seismotectonics
in Convergent Plate Boundary, edited by Y. Fujinawa and A. Yoshida, pp. 107–121, Terra Scientiﬁc Comp., Tokyo, Japan.
Van Keken, P. E., B. R. Hacker, E. M. Syracuse, and G. A. Abers (2011), Subduction factory: 4. Depth-dependent ﬂux of H2O from subducting
slabs worldwide, J. Geophys. Res., 116, B01401, doi:10.1029/2010JB007922.
Wada, I., M. D. Behn, and J. He (2011), Grain-size distribution in the mantle wedge of subduction zones, J. Geophys. Res., 116, B10203,
doi:10.1029/2011JB008294.
Wannamaker, P. E., et al. (2008), Lithospheric dismemberment and magmatic processes of the Great Basin–Colorado Plateau transition, Utah,
implied from magnetotellurics, Geochem. Geophys. Geosyst., 9, Q05019, doi:10.1029/2007GC001886.
Wannamaker, P. E., T. G. Caldwell, G. R. Jiracek, V. Maris, G. J. Hill, Y. Ogawa, H. M. Bibby, S. L. Bennie, and W. Heise (2009), Fluid anddeformation regime of an advancing subduction system at Marlborough, New Zealand, Nature, 460, 733–736, doi:10.1038/nature08204.
Wessel, P., and W. H. F. Smith (1998), New, improved version of generic mapping tools released, Eos Trans. AGU, 79(47), 579, doi:10.1029/
98EO00426.
Wilson, C. R., M. Spiegelman, P. E. van Keken, and B. R. Hacker (2014), Fluid ﬂow in subduction zones: The role of solid rheology andcompaction pressure, Earth Planet. Sci. Lett., 401, 261–274, doi:10.1016/j.epsl.2015.05.052.
Wolf, M. B., and P. J. Wyllie (1994), Dehydration-melting of amphibolite at 10 kbar: The effects of temperature and time, Contrib. Mineral.
Petrol., 115, 369–383, doi:10.1007/BF00320972.
Yang, Y., and D. W. Forsyth (2006), Rayleigh wave phase velocities, small-scale convection, and azimuthal anisotropy beneath southern
California, J. Geophys. Res., 111, B07306, doi:10.1029/2005JB004180.
Yoshino, T., K. Mibe, A. Yasuda, and T. Fujii (2002), Wetting properties of anorthite aggregates: Implications for ﬂuid connectivity incontinental lower crust, J. Geophys. Res., 107(B1), 2027, doi:10.1029/2001JB000440.
Yoshino, T., T. Matsuzaki, S. Yamashita, and T. Katsura (2006), Hydrous olivine unable to account for conductivity anomaly at the top of theasthenosphere, Nature, 443, 973–976, doi:10.1038/nature05223.
Zhao, D., A. Hasegawa, and S. Horiuchi (1992), Tomographic imaging of P and S wave velocity structure beneath northeastern Japan,
J. Geophys. Res., 97, 19,909–19,928, doi:10.1029/92JB00603.
Zhu, G., T. V. Gerya, D. A. Yuen, S. Honda, T. Yoshida, and J. A. D. Connolly (2009), Three-dimensional dynamics of hydrous thermal-chemical
plumes in oceanic subduction zones, Geochem. Geophys. Geosyst., 10, Q11006, doi:10.1029/2009GC002625.

ICHIKI ET AL.

ELECTRICAL IMAGE BENEATH NE JAPAN29