Imaging of a serpentinite complex in the Kamuikotan Zone, northern Japan, from magnetotelluric soundings

We conducted magnetotelluric measurements to investigate a large serpentinite complex in the northern Kamuikotan Zone that intruded a Cretaceous–Paleocene forearc sedimentary sequence. The resistivity model we derived by three-dimensional inversion clearly shows a low-resistivity zone beneath the outcrop of the serpentinite complex. We interpret the low-resistivity zone to represent aqueous pore fluid within a serpentinite mélange derived from the subducting Pacific plate or mantle wedge. Previous geological studies in the area have shown that the serpentinite mélange had uplifted during the early Pleistocene. They indicate that the ultramafic rocks and aqueous fluids have continued to rise in the area. The uplifting serpentinite body might have formed a zone enriched in pore fluid that promoted the occurrence of a previously identified intra-plate slow slip event. These results demonstrate the important role of fluid transport during tectonic processes related to uplift in subduction zones.


Introduction
Serpentinites in subduction zones contain large volumes of aqueous fluid in the form of hydrated minerals and pore water (e.g., Hyndman and Peacock 2003), and play essential roles for aqueous fluid transport, fault rupturing, and tectonic processes. The serpentinites are formed by the hydration of ultramafic rocks in the mantle wedge, the subducting oceanic plate, or in normal fault zones at the outer-rise off the oceanic trench (e.g., Guillot et al. 2015). In outcrop, such serpentinites are commonly associated with uplifted high-P/T metamorphic rocks. However, serpentinites in outcrop are not well understood because the broad structural features of deep subsurface serpentinite bodies have not been directly imaged.
The area of the Shirikomadake serpentinite complex in northern Hokkaido (Fig. 1) is important in studies on the role of serpentinites in subduction tectonics. The Shirikomadake complex consists of serpentinized ultramafic rocks and gabbro-diorite dykes in a serpentinite mélange (Katoh et al. 1979) and is one of the largest ultramafic bodies in the Japan arc. Faulting of the serpentinite complex against Cretaceous to Quaternary sedimentary rocks implies that the complex experienced large-scale uplift and intrusion (e.g., Igi 1959). In addition, Ohzono et al. (2015) detected a slow slip event (SSE) near the edge of the serpentinite body as the first discovery of slow earthquake in intra-plate (Fig. 1b). Because slow earthquakes have been attributed to the presence of pore fluids (e.g., Shelly et al. 2006;Obara 2020), the intra-plate SSE may also be associated with pore fluid derived from the serpentinites. Because of the uniqueness, the Shirikomadake serpentinite complex provides an ideal opportunity to undertake a detailed investigation of this association.

Tectonic and geologic setting
The Shirikomadake serpentinite complex is in northern Hokkaido Island, which lies at Kamuikotan Zone in the Northeast Japan arc. The following summary of the geology of the island is based on the work of Ueda (2016). The Sorachi-Yezo belt (Fig. 1a) is characterized by a Jurassic ophiolite overlain in turn by a Late Jurassic to Early Cretaceous ophiolite and siliceous sedimentary sequence (Horokanai Ophiolite and Sorachi Group) and a Late Cretaceous to Paleocene forearc basin sequence (Yezo Group). The Kamuikotan Zone represents the cores of anticlines within the Sorachi-Yezo belt (Fig. 1a) and characterized by high-P/T metamorphic rocks (Kamuikotan metamorphic rocks) and serpentinite mélange. Serpentinite bodies are dispersed through the Kamuikotan Zone. The metamorphic rocks of the Kamuikotan Zone originated from the subducting slab and accreted rocks, and  Ogura and Kamon (1992) and the cyan rectangle marks the location of the borehole and resistivity survey of Okazaki et al. (2011). The red dashed rectangle encloses the rupture zone of the SSE that occurred between July 2012 and January 2013  were metamorphosed during the early Cretaceous to Paleocene (e.g., Sakakibara and Ota 1994). The source rocks of the Shirikomadake serpentinite complex were dunite, harzburgite and orthopyroxenite and the complex contains antigorite that formed after lizardite, chrysotile, and brucite (Katoh et al. 1979). Late Cretaceous (Yezo Group) and Paleogene-Neogene sedimentary rocks distributed around the complex (Fig. 1b) thicken westward, extending to a depth of more than 4500 m in the western part of the study area (borehole SK-1 in Fig. 1b) (Ogura and Kamon 1992), and are deformed by active folds and thrusts (e.g., Oka 1985). Miocene igneous rocks are distributed in the middle and eastern parts of the study area. Rocks of a Cretaceous accretionary complex in the eastern part of the study area represent the Idonnappu Zone (Suzuki et al. 1997).
Geodetic studies indicate that the study area is presently in a zone of east-west regional compression (e.g., Sagiya et al. 2000), which is consistent with the geology described above. Crustal seismic activity is high in the western and central parts of the study area (e.g., Tamura et al. 2003) (Fig. 1b). Ichiyanagi et al. (2015) identified earthquake swarms that occurred in the south of the study area between July 2012 and January 2013. Based on Global Navigation Satellite System data, Ohzono et al. (2015) inferred a SSE (Mw 5.4) that occurred in the eastern part of the zone of earthquake swarms during the same period (Fig. 1b). The fault mechanisms of both the earthquake swarms and the SSE are consistent with the regional east-west compressional stress field. In the eastern part of the study area, seismicity is low and no active faults or folds have been recognized.

Magnetotelluric measurements and impedance tensors
We conducted broadband MT measurements in 2001MT measurements in , 2002MT measurements in , and 2018, and 3 sites, respectively (Fig. 1b). We measured two horizontal components of the electric field by using Pb-PbCl 2 electrodes and three orthogonal components of the magnetic field by using induction coils. These data were recorded by MTU-5 systems (Phoenix Geophysics, Ltd., Toronto, Canada) for 2-6 days. MT impedance tensors were estimated from the obtained electromagnetic fields during 0.00325 to 1820 s by using SSMT2000 software (Phoenix Geophysics, Ltd.) and remote reference processing (Gamble et al. 1979). For remote reference processing, we used horizontal magnetic field data obtained at the following sites: Esashi station (operated by the Geospatial Information Authority of Japan) for the 2001-2002 MT data, andSawauchi station (operated by Nittetsu Mining Consultants) for the 2018 MT data. Both of the remote stations are about 400 km south-southwest of the study area.
Estimated MT impedances are shown in Fig. 2 and Additional file 1. To visualize the spatial trend of impedance, we estimated sounding curves for the sum of the squared elements' invariant impedances (SSQ; e.g., Rung-Arunwan et al. 2016; see Additional file 2), which show trends of low and high resistivity in the western and eastern parts of the study area, respectively.
In the western area (red symbols in Fig. 1b and Additional file 2), phase angles are high and apparent resistivity decreases with increasing period in the shortmiddle-period band (< 30 s). These results imply the presence of a near-surface conductive layer. In the central area (green and black symbols in Fig. 1b and Additional file 2), apparent resistivity is low and phase angles are neutral in the short-middle-period band (< 30 s). These results imply that the conductive anomaly is distributed in the subsurface. The splitting of off-diagonal components at longer periods implies the presence of deeper 2-D or 3-D structures.
In the eastern area, the sounding curves show considerable variations and large diagonal components (blue symbols in Fig. 1b and Additional file 2). These results are indicative of a strong 2-D or 3-D structure.

3-D electrical resistivity inversion
We constructed a regional 1-D electrical resistivity model to use as the initial model of the 3-D inversion to minimize the local minima problem inherent in nonlinear inversions. The 1-D modeling was based on the 1-D Occam's inversion procedure developed by Constable et al. (1987). We used the averaged SSQ invariant impedances as input data (Additional file 3a) and set the target RMS misfit of the inversion at 1.00. The modeled 1-D resistivity profile explained the data reasonably well (RMS misfit: 1.00) and showed an increase of resistivity with depth (Additional file 3).
The 3-D resistivity distribution was then estimated using the inversion code of Tada et al. (2012), which is a modified version of the WSINV3DMT 3-D inversion code (Siripunvaraporn et al. 2005) that adopts a data-space variant of Occam's inversion procedure. The modified WSINV3DMT code can incorporate boundaries between land and seawater where high-resistivity contrasts can seriously distort MT responses. In this inversion procedure, the following objective function is minimized: where m and m p are model parameters and prior model vectors, respectively. C m is the model covariance matrix that characterizes the expected magnitude and smoothness of resistivity variations relative to m p , and is controlled by hyperparameters δx, δy, δz, and τ. d is the data parameter vector consisting of the observed MT impedances. F[m] is the vector of the forward response to m. C d is a data covariance matrix representing the observation errors. λ is a hyperparameter that balances data misfit and model misfit terms, including model roughness, and is determined in each iteration to minimize the RMS misfit when the misfit does not reach a target value. Once the misfit reaches the target value, the model based on the largest value of λ that maintains the misfit at the target value is selected as the final model following Occam's inversion procedure.
MT impedances at 16 periods between 0.0667 and 1820s were used as data parameters (d) in the inversion. We added error floors to avoid overfitting the calculated responses to observed data for which errors were small. The error floor was defined as 5% of the SSQ impedance and was applied to all components of the MT impedance. Because the geomagnetic transfer function (tipper) was not available in the inversion code we used, it was not used in this study. The resistivity model space defined a volume of 3268 (x-axis) × 3268 (y-axis) × 1049 km (z-axis, without air layers) discretized into 58 × 58 × 37 blocks. The x and y dimensions of the blocks were 2 km within the survey area but were widened outside the study area. The z dimensions of the blocks increased with depth and were the same as the those of the 1-D inversion (Additional file 3b). The inversion procedure began with the 1-D resistivity model (Additional file 3b), except for the area of seawater where we set the resistivity at 0.3 Ω m (RMS misfit = 10.07). The prior model (m p ) was the same as the initial model in this procedure. We used the default values (δx = δy = δz = 0.1 and τ = 5) of the WSINV3DMT 3-D inversion code (Siripunvaraporn et al. 2005) to control the values of C m . The initial value of λ was 1.0 and the target RMS misfit was 1.0. We iterated the inversion five times and obtained a minimum RMS misfit model (RMS misfit = 2.23, λ = 0.0316) at the third iteration. We then updated m p as the minimum RMS misfit model and again iterated the inversion five times. Because these iterations did not reach the target RMS misfit (1.0), we adopted the minimum RMS misfit model obtained at the fifth iteration (RMS misfit = 1.47, λ = 0.0316) as the inverted model. The inverted model shows a shallow conductive layer (C-1, Fig. 3) at 0-5 km depth in the western part of the study area. This layer is also clearly indicated by low apparent resistivities of the off-diagonal components of MT impedance in the western part of the study area (Additional file 2). Also in the western part of the study area, a layer of moderate resistivity (10-100 Ω m) partly extends from the surface to a maximum depth of about 1 km as shown in the section "X = -5 km" in Fig. 3, as indicated by the moderate apparent resistivity and high phase angle in the shortperiod band (e.g., site D13 in Fig. 2). The C-1 layer thins eastward and disappears approaching the eastern edge of the study area, consistent with the west-east trend of the sounding curves (Additional file 2). A near-vertical anomalous conductive zone (C-2, Fig. 3) was modeled under the area where the Shirikomadake serpentinite complex crops out in the south-central part of the study area. The inverted model also shows a zone of high resistivity that extends from the surface in the east of the study area. This modeled high resistivity is consistent with the high apparent resistivity in this area.

Sensitivity test for resistivity anomalies
Because the C-2 anomaly was not apparent in our original sounding curves, we performed the following sensitivity test to examine the anomaly further. We replaced modeled resistivity values within the C-2 anomaly that were < 100 Ω m (yellow polygons in Additional file 4) with an approximation of the resistivity of the surrounding rocks (100 Ω m) and recalculated the model response using the modified WSINV3DMT code. The sounding curves of this model changed considerably from those of the inverted model and did not explain the observed impedances at sites above the C-2 anomaly (e.g., Sites M18 and S18 in Fig. 2). In particular, the low apparent resistivities at periods > 3 s and moderate phase angles (45° and -135°) between periods of 1 and 30 s of the off-diagonal components at Site M18 were not explained. These results indicate that the observed data require the presence of conductive anomaly C-2.
Next, we examined the reliability of the deep extension of the C-2 anomaly. We varied the depth of the upper limit of the zone of replacement resistivities in the yellow polygons in Additional file 4, whereas for other parameters we used the values of the original inverted model. RMS misfits of the replaced models that were significantly larger than that of the original inverted model indicate model degradation, thus requiring the upper limit of the conductive anomaly to be deeper. We adopted the one-sided F-test at the 95% confidence level to indicate statistical significance (i.e., sensitivity models that fail the F-test are statistically inferior to the original inverted model). The results of this analysis show that the sensitivity models with upper limits at depths ≤ 10.6 km failed the F-test and were thus significantly inferior to the inverted model (Fig. 4a). Therefore, our original inversion of the C-2 anomaly is reliable to 10.6 km depth.
To further examine the reliability of the resistivity values of the C-2 anomaly, we replaced modeled resistivity values within the C-2 anomaly that were < 10 Ω m ( Fig. 3b and red polygons in Additional file 4) with a range of replacement resistivities and then recalculated the model responses. F-tests at the 95% confidence level indicate that the models derived using replacement resistivities > 12 Ω m and < 1.2 Ω m were significantly inferior to the inverted model (Fig. 4b). If it is assumed that the conductive area is uniformly distributed within the anomalous zone, these results indicate that the resistivity range of the C-2 anomaly lies between 1.2 and 12 Ω m. Note that the real conductive zone is possibly narrower and of lower resistivity because of the smoothness constraint of the inversion.
We also examined the zone of high resistivities in the eastern part of the study area, where some resistivities exceeded 10,000 Ω m. We calculated responses of the model that the resistive zone (delimited by x = -16 to 16 km, y = 0 to 30 km, and z = 0.85 to 29 km) were replaced with lower values of 600, 800, and 1000 Ω m. One sided F-tests at the 95% confidence level indicated that the models derived by using replacement resistivities less than 800 Ω m were significantly inferior to the inverted model, thus indicating that resistivities are higher than 800 Ω m within the zone of high resistivity in the eastern area.  Ichihara et al. Earth, Planets and Space (2021) 73:154

Interpretation of the resistivity model
The C-2 anomaly lies directly below the outcrop of the Shirikomadake serpentinite complex (Fig. 3). Similar conductive anomalies have been identified under a serpentinite body in the southern part of the Kamuikotan Zone (Ogawa et al. 1994;Ichihara et al. 2019) and in an area around the San Andreas fault (Unsworth and Bedrosian 2004;Becken and Ritter 2012). In a general sense, these relationships can be considered to imply an association between serpentinite bodies and underlying conductive anomalies. Previous geological studies have shown vertically inclined sedimentary rocks near the contact of the Shirikomadake serpentinite complex that imply the intrusion of the complex (Igi 1959;Oka 1985). Other regional geological studies of the vertical distribution of serpentinite complexes in the Kamuikotan Zone (e.g., Ueda 2005) have illustrated that they extend to considerable depths. However, the resistivity (3-100 Ω m) in the shallow part of our 3-D model of the outcropping serpentinite complex (Fig. 3b) is not particularly low. Near-surface resistivity surveys (0-200 m depth) conducted in the south of the Shirikomadake serpentinite complex (Fig. 1b) by Okazaki et al. (2011) detected a range of resistivities consistent with our results. They also showed that the resistivities of their massive, foliated, and clay-rich serpentinite core samples were 250-1,500, 200, and 20-40 Ω m, respectively. Therefore, the resistivities of serpentine group minerals are not the sole control of conductive anomalies in serpentinite, such as our C-2 anomaly.
Candidates for the source of conductive anomaly C-2 include the presence of (1) interconnected pore fluids, (2) interconnected conductive minerals (e.g., magnetite) that are commonly found in serpentinite, and (3) a high-temperature zone or melt that decreases resistivity.
The first of these is the most likely candidate as previous studies of conductive anomalies in subduction zones (e.g., Aizawa et al. 2021;Ichihara et al. 2016;Evans et al. 2014) and an experimental study for conductive antigorite (Reynard et al. 2011) have suggested. The formation of serpentinite in mantle wedges requires an abundant supply of fluids from the subducting plate or possibly from repeated serpentinite dehydration (e.g., Hyndman and Peacock 2003). Large volumes of interconnected fluids are common in mélange zones because of their numerous fractures and high permeability. In addition, the resistivity of dehydrated saline fluid is low at depth (Sakuma and Ichiki 2016;Sinmyo and Keppler 2017). Therefore, serpentinite mélange in subduction zones is a favorable environment for the presence of interconnected conductive aqueous pore fluids resulting in a low-resistivity anomaly. Note that a conductive zone could also be formed by interconnected conductive fluids in a fractured area free of serpentinites, as indicated by the presence of conductors within tectonic boundaries (e.g., Ichihara et al., 2011;Ikeda et al., 2013). If this were the case in our study area, a considerable area of fractured rock, Resistivity (Ωm) Fig. 4 Results of the sensitivity tests for the C-2 anomaly. The green lines are the RMS misfit of the inverted model and the blue lines show the limit of the 95% confidence interval of the F-test. a Upper limits of the zone of replacement of anomalous C-2 resistivities versus RMS misfit. b Resistivities used to replace anomalous C-2 resistivities versus RMS misfit such as a fault zone, would be required within the C-2 anomaly. Although few geological studies have reported on faulting in this area, Igi (1959) identified a gap in the Cretaceous sedimentary sequence between the western and eastern parts of the Shirikomadake serpentinite complex.
There are two possible explanations for the relatively high resistivity zone between the serpentinite outcrop and the C-2 anomaly (Figs. 3, 4). It might represent resistive rocks that contain little or no pore fluid, such as, massive serpentinite, basaltic lavas (Katoh et al. 1979), or metamorphic rocks of the Kamuikotan zone that do not outcrop in the study area. Alternatively, the resistivity of pore fluid remains high in the shallow sequence either because the resistivity of saline water increases with decreasing temperature (i.e., with decreasing depth; Sakuma and Ichiki 2016;Sinmyo and Keppler 2017) or because the pore water in the shallow sequence is highly resistive meteoric water.
The second candidate (interconnected conductive minerals) has been considered for mid-crustal conductive areas (e.g., Myer et al. 2013) because magnetite contained in serpentinites is highly conductive and can produce a conductive anomaly (Stesky and Brace 1973). However, it is difficult to attribute the C-2 anomaly to magnetite alone because both the conductivity of borehole samples collected by Okazaki et al. (2011) and the surface-measured resistivity of the Shirikomadake serpentinites do not indicate a sufficient anomaly, although the presence of magnetite may contribute to the low resistivity of interconnected pore fluids. The third candidate (high-temperature zone or melt) is unlikely to be a direct cause of the C-2 anomaly because the study area is a considerable distance from the nearest area of volcanic activity and no high-temperature phenomena, such as hot springs, are known in the study area.
Borehole SK-1 near the western margin of the study area (Fig. 1b) initially penetrated Paleogene-Neogene sedimentary rocks at 3030 m depth and, was within late Cretaceous sedimentary rocks at its total depth (4505 m) (Ogura and Kamon 1992). Borehole logs indicate that resistivities in the Paleogene-Neogene and late Cretaceous sequences were 1-4 and 4-10 Ω m, respectively (Kanekiyo 1999). The modeled resistivity within anomaly C-1 (Fig. 5) is consistent with that logging data. It is noteworthy that the surface extent of the C-1 anomaly is consistent with the mapped extent of the Cretaceous (Yezo Group) and Paleogene-Neogene sedimentary rocks (Figs. 1b, 3a). Hence, we interpret the C-1 anomaly to represent Cretaceous and Paleogene-Neogene conductive sedimentary rocks. Similar associations between resistivity anomalies and outcrops have been recognized in other parts of middle-western Hokkaido (Ichihara et al. 2008(Ichihara et al. , 2016(Ichihara et al. , 2019Yamaya et al. 2017).
The C-1 anomaly deepens toward the western boundary of the study area (Fig. 5) and the steepest dip of its lower boundary coincides with thrust faults (Fig. 5) that uplifted the eastern part of the study area. The C-1 anomaly is partly overlain by a moderate resistivity (about 100 Ω m) layer from the surface to a maximum depth of about 1 km that is similar to a near-surface layer of moderate resistivity identified by Ueda et al. (2014) near sites B13 and B14. On the basis of seismic reflection survey data, they interpreted this layer to be the Quaternary Sarabetsu Formation, which contains relatively fresh groundwater. Thus, we interpret the layer of moderate resistivity that we modeled overlying anomaly C-1 to represent Quaternary sediments.

Uplift of serpentinite and its implications for intra-plate SSEs
Previous geological studies have identified faults at the boundaries between the Shirikomadake serpentinite complex and surrounding sedimentary rocks (Igi 1959;Katoh et al. 1979;Oka 1985). The early Pleistocene Sarabetsu Formation near the center of the study area (Oka and Igarashi 1993;Niwa et al. 2020) contacts with the serpentinite complex by fault and is vertically inclined near the contact, indicating that uplift of the complex occurred during the Pleistocene or later. In addition, tectonic fault blocks within the serpentinite complex contain Cretaceous sedimentary rocks (Katoh et al. 1979), suggesting that Cretaceous rocks from the surrounding area were dragged into the serpentinite complex during the Quaternary (or earlier) uplift. The extension of the C-2 anomaly to a depth of 10.6 km or more possibly indicates that the uplifting serpentinite complex (Fig. 5) or perhaps the fault zone. Previous researchers (e.g., Guillot et al. 2015) have proposed buoyancy-driven diapiric ascent of serpentinites from depth. Therefore, serpentinite complexes in the Kamuikotan zone may have been uplifted in the surrounding rocks and thus contribute to the transport of aqueous fluids derived from the subducting plate.
Recent seismic studies have proposed that intra-plate SSEs occur in subduction zones when the presence of high-pressure pore fluids weakens the plate interface (e.g., Shelly et al. 2006;Kato et al. 2010). In our study area, an intra-plate SSE was inferred near the western edge of the Shirikomadake serpentinite complex during 2012-2013   (Fig. 1). Because we interpreted the C-2 anomaly to represent upwelling aqueous pore fluid, such fluid may have been available to the area where that SSE occurred (Fig. 5). Ohzono et al. (2015) indicated that the SSE might have occurred in a detachment fault at the base of the Cretaceous-Neogene sedimentary rocks. Moreover, if the sedimentary rocks contained impermeable clay minerals and thus acted as cap rocks, they may have trapped (Fig. 5) and elevated the pressure of the pore fluids within the detachment fault, thus causing an intra-plate SSE. Therefore, the plateinterface SSE reported by Ohzono et al. (2015) might well have been caused by high-pressure fluids associated with the Shirikomadake serpentinite complex. However, our study area does not fully cover the fault zone and the earthquake swarm that accompanied the SSE occurred mainly beyond the southern margin of the MT array . Thus, additional MT soundings south of our study area are required to better understand the fault rupture processes in the study area.

Conclusions
We conducted broadband magnetotelluric soundings at 48 sites around the Shirikomadake serpentinite complex in the Kamuikotan Zone in Hokkaido, northeastern Japan arc. Our 3-D resistivity inversion of the estimated MT impedance revealed two anomalies: a near-surface conductive layer (C-1) and a conductive zone beneath the outcrop of the serpentinite complex (C-2). We interpreted the C-1 anomaly to represent Cretaceous-Neogene sedimentary rocks that showed evidence of thrust faulting (Fig. 5) and the C-2 anomaly to represent interconnected pore fluids, possibly within the serpentinite complex (Fig. 5). Our results and previous geological studies indicate long-term uplift of the serpentinite from the mantle wedge or plate interface that might have continued until the early Pleistocene or later. These results demonstrate the important role of tectonic processes related to the uplift of metamorphic belts for fluid transport in subduction zones and indicate that such fluids may have contributed to the 2012-2013 intra-plate SSE that occurred between conductive anomalies C-1 and C-2.  Ohzono et al. (2015). The other symbols are the same as in Fig. 3. b Geological interpretation of the inverted model at X = -5 km. The pink area with red arrows represents the uplifting serpentinite complex and associated fluids, and includes the C-2 anomaly, the surface conductive area, and the high resistivity area sandwiched between them