Early magnetic studies of the Japan Trench showed that seafloor spreading magnetic anomalies progressively fade away and disappear during subduction, reflecting the increasing distance to magnetized sources and the removal of their remanent magnetization with alteration and increasing temperature. An improved magnetic anomaly map derived from both scalar and vector magnetic anomaly data, coupled with a better knowledge of the slab geometry in one hand, of the magnetic structure of the oceanic crust on the other hand, allow us to constrain the thermal structure of the subducting slab. We, for the first time, identify two steps in the anomaly disappearance: first the magnetization of extrusive basalt is rapidly erased between 9 and 12 km, where titanomagnetite reaches its blocking temperature between 150 °C and 350 °C, then the magnetization of deeper crustal layers slowly decreases down to ∼20 km, reflecting the progressive slab heating toward the Curie temperature of magnetite, 580 °C. The resulting slab temperatures are higher than predicted by most thermal models. Recent observations and models suggest rejuvenated hydrothermal activity triggered by lithospheric flexure before subduction that may significantly heat up the subducting oceanic crust through thermal blanketing and possibly serpentinization, with consequences on the depth of the seismogenic zone.

Subduction is a major geodynamic process that leads to the consumption of oceanic lithosphere, the creation of volcanic arcs and the largest volcanic edifices on Earth, and the generation of natural disasters, such as earthquakes and tsunamis, resulting in major damage and casualties. Many geophysical studies have been conducted over the Japan Trench to understand the structure of the subducting plate. Detailed slab geometry models Slab 1.0 and Slab 2 (Hayes et al., 2012, 2018) have been developed based on active seismic profiles and seismicity studies. Thermal structure models have been proposed based on heat flow surveys (Hyndman and Peacock, 2003; van Keken et al., 2012; Kawada et al., 2014; Wada et al., 2015). A pioneer study of magnetic anomalies on a subducting plate has shown that, from the Japan Trench landward, the magnetic anomalies decrease in amplitude and their short wavelength content attenuates (Okubo et al., 1991), as a result of increasing distance between the magnetized sources and the observation points, thermal demagnetization, and continuous oxidation of magnetic minerals within the extrusive oceanic crust (Okubo et al., 1991; Kido and Fujiwara, 2004). Today, the structure of the slab is well constrained (Hayes et al., 2012, 2018); moreover, the magnetic structure of the oceanic crust is better understood (e.g., Dyment and Arkani-Hamed, 1995; Gee and Kent, 2007).

Based on this existing information, the magnetic anomalies over the Japan Trench offer an independent means to access the slab thermal structure. To address this problem, we take advantage of a unique scalar and vector marine magnetic data set to build an improved high-resolution magnetic anomaly map. We identify two steps of thermal demagnetization, corresponding to the two major magnetic minerals of the oceanic crust layers, which constrain the thermal structure at shallow depths within the subduction system.

The oceanic crust in the northwestern Pacific plate off Japan was formed at ca. 125–140 Ma at the Pacific-Izanagi plate boundary (Nakanishi et al., 1989). The area is covered by pelagic sediments 1.6 km thick, and subducts beneath the Japan islands at a speed of 70–90 km/m.y. (Seno, 2017). A thicker and wider accretionary prism is observed in the northern Japan Trench (Kodaira et al., 2017). The free-air gravity anomaly locally delineates seamounts and fracture zones (Fig. 1A). In the entire area, it displays a negative anomaly associated with the trench and a positive anomaly on the flexural bulge ∼150 km seaward. This bulge reflects the bending of the subducting plate, and induces horsts and grabens as the plate approaches the trench (Tsuru et al., 2002; Kodaira et al., 2017).

We gathered two sets of marine magnetic data. Scalar magnetic data (i.e., total field vector intensity) acquired by proton precession magnetometer (PPM) were obtained from the Data and Sample Research System for Whole Cruise Information (DARWIN) program of the Japan Agency for Marine-Earth Sciences (JAMSTEC) (http://www.godac.jamstec.go.jp/darwin/e), the Geophysical Data System (GEODAS, www.ngdc.noaa.gov/mgg/geodas) of the U.S. National Oceanic and Atmospheric Administration’s National Center for Environmental Information, and the Nautilus data base of the Institut Français de Recherche pour l’Exploitation de la MER (IFREMER, http://donnees-campagnes.flotteoceanographique.fr). Vector magnetic data (i.e., total field vector components) acquired by a shipboard three-component magnetometer (STCM), mostly on Japanese research vessels, were obtained from JAMSTEC (the DARWIN database). The International Geomagnetic Reference Field (IGRF) model (Thébault et al., 2015) was subtracted from the PPM data. The STCM data were corrected for the ship magnetic effect and motion using the method of Isezaki (1986). To minimize the misfit at cross-over points, a modified cross-over error analysis technique was applied, both to the PPM data and to the STCM data (see the GSA Data Repository1). We reduced the corrected magnetic anomaly grid to the pole (RTP) assuming a 53.6° inclination and − 7.6° declination (IGRF averaged over 20 years), 33° paleoinclination and 11° paleoazimuth (average value for the study area from the global grids of Dyment and Arkani-Hamed [1998]).

High-resolution bathymetric data were obtained from Global Multi-Resolution Topography (GMRT; Ryan et al., 2009). The top of the magnetic source, i.e., of the extrusive basalt layer, for the subducting plate was computed by merging and correcting different grids. The World sediment thickness grid (Divins, 2003) was subtracted from the bathymetry grid, and the resulting grid was merged with Slab 1.0 (Hayes. et al., 2012; for a discussion on the choice of the Slab 1.0 model, see the Data Repository) validated by published seismic profiles (e.g., Tsuru et al., 2002; Kodaira et al., 2017). The geometrical misfit between the two grids along the trench boundary was erased and re-interpolated using the Partial Differential Equation (PDE) surface method (D’Errico, 2005).

We carefully selected magnetic anomalies showing no tectonic or volcanic local complexities, such as fracture zones, propagators, or seamounts. Anomalies older than M10r (130.8 Ma) and younger than M8 (129.0 Ma) are therefore discarded (Fig. 1). We extracted profiles from the magnetic grid across the selected anomalies and considered separately the anomaly profiles before and after subduction; i.e., located east and west of the Japan Trench. The anomaly profiles before subduction were inverted to equivalent magnetization assuming a 500-m-thick magnetized source layer with no vertical variation of magnetization (Parker and Huestis, 1974). The equivalent magnetization shows little variation among profiles and has been averaged (Fig. 2A). We use this average equivalent magnetization and the inferred top of the magnetic source layer along the anomaly profiles after subduction (Fig. 2B) to compute synthetic magnetic anomalies along these profiles (Fig. 2C). These modeled anomalies represent the contribution of the subducting plate at the sea surface if the magnetic structure of the plate remains unchanged. Comparison of these synthetic anomalies with the observed ones gives us the opportunity to estimate how the magnetic structure of the plate has been changed (Fig. 2C).

The ratio of peak to trough anomaly amplitudes of the observed and synthetic anomalies, hereafter named RAM (Remaining Amount of Magnetization), provides an estimate of the remaining fraction of magnetization in the subducting plate: when it is close to 1, demagnetization remains negligible, whereas when it tends to 0, demagnetization is almost complete. We expect demagnetization to progress as a function of depth, and display RAM as a function of the depth to the top of the magnetized layer under the maximum (respectively minimum) of the observed positive (respectively negative) anomalies (Fig. 3).

We compiled both scalar and vector marine magnetic data available in the study area and merged them into a unique scalar magnetic anomaly map (Fig. 1B). Many gaps in the scalar anomaly coverage could be filled with the vector data. The resulting map shows two types of anomalies. Near the Japanese Islands, a strong positive anomaly (FIMAB in Fig. 1) is caused by the induced magnetization of serpentinite in the forearc mantle (Okubo and Matsunaga, 1994; Blakely et al., 2005) and, possibly, of the volcanic arc. On the Pacific plate and subducting slab, alternating positive and negative lineated anomalies are caused by the remanent magnetization of the oceanic crust. Magnetic anomalies M5 to M17 (ca. 124.6–139.7 Ma; Malinverno et al., 2012) are identified between two NNW-SSE–trending fracture zones depicted on the free-air gravity anomaly (Fig. 1A), confirming the interpretation of Nakanishi et al. (1989).

In our study area, only anomalies M8r–M10r (129.0–130.8 Ma) are suitable for a detailed analysis, as they are linear and only disrupted by a few isolated seamounts. They display a clear seafloor-spreading magnetic signal before and after subduction (Fig. 1A). Conversely, anomaly M5 is still recognizable but totally subducted, while anomalies M6–M7 are only partially subducted beneath the southern Kuril Trench, roughly parallel to the magnetic anomaly trend. The presence of propagating rifts (Nakanishi, 2011) and the Joban seamount chain affect anomalies M10n1–M15 south of 38°40′N. During this period, the (half) spreading rate was 70–80 km/m.y. and the studied oceanic crust was formed at a fast spreading center (Nakanishi et al., 1989).

The amplitude of the slab magnetic anomalies continuously decreases landward, and the anomalies disappear beyond 20 km depth below sea level (bsl) (Fig. 1C). The increasing depth of the magnetized source preferentially attenuates shorter anomaly wavelengths, as does demagnetization of the oceanic crust induced by fluids (alteration) and increasing temperature (thermal demagnetization). To separate the effects of increasing depth and demagnetization, we compute synthetic magnetic anomalies assuming an unchanged magnetic structure of the subducting plate. The observed anomalies decrease faster than the synthetic ones (Fig. 2C). Indeed, the RAM decays rapidly, by 20% per km, between 9 and 12 km bsl, and more slowly, by 2% per km, beyond (Fig. 3A).

The observation of two distinct depth intervals with contrasted RAM decay is in good agreement with our knowledge of the magnetic structure model of oceanic crust formed at fast spreading centers, where two magnetic layers are distinguished (Dyment and Arkani-Hamed, 1995). The shallower one, <1 km thick (e.g., Karson, 2002), is made of extrusive basalt of which the magnetic mineral is titanomagnetite with various Ti content and oxidation state (Curie temperature Tc of 100–350 °C; Zhou et al., 2001). The deeper one, ∼5 km thick, is made of dolerite, gabbro, and serpentinized peridotite in which magnetite (Tc 580 °C) is the dominant magnetic phase. The natural remanent magnetization (NRM) carried by these layers varies: the extrusive layer bears a strong NRM (>10 A/m) at the ridge axis, which decays to ∼3 A/m for old oceanic crust. It is the main contributor to the lineated marine magnetic anomalies observed at the sea surface. Conversely, the dolerite and gabbro bear a weaker NRM (∼1–1.5 A/m) and the serpentinized peridotite has a variable NRM (0–6 A/m) depending the degree of serpentinization (Harrison, 1987). We, therefore, suggest that the sudden decay of RAM between 9 and 12 km bsl corresponds to the thermal demagnetization of the extrusive basalt layer, and the slow decrease beyond 12 km bsl to that of the deeper layers (Fig. 3B). Observed magnetic anomalies beyond 18 km bsl are very low and their shape does not match that of the synthetic ones, precluding the calculation of RAM. Their short wavelength content is not consistent with the expected depth of subducting slab, suggesting the presence of shallower sources in the upper plate continental crust.

A global compilation over subduction zones suggests that, apart from local effects, heat flow over oceanic lithosphere entering subduction does not significantly deviate from that of normal oceanic lithosphere (Stein, 2003), whereas a more recent study supports higher heat flow over oceanic lithosphere approaching subduction (Harris et al., 2017). Conversely, heat flow over the forearc basin is generally low (Stein, 2003; Yamano et al., 2014). Heat flow measurements on the Pacific plate off the Japan Trench range between 50 and 100 mW/m2, significantly higher than the ∼50 mW/m2 expected for oceanic lithosphere of this age (Yamano et al., 2014).

Estimating the slab thermal structure is difficult because the hydrothermal circulation in the accretionary prism and overriding continental crust is hard to quantify. Proposed models consider constraints such as the age of oceanic lithosphere, rate of convergence, shear heating, and dip of the slab (van Keken et al., 2012; Wada et al., 2015). However, strong uncertainties remain in the shallow subduction zones on the effectiveness of hydrothermal circulation, the permeability of the igneous crust, and its evolution with depth. Consequently, models propose a wide variety of temperature ranges for the shallow subducting slab. For instance, the obtained temperatures are as low as ∼100–200 °C from the deformation front to 20 km depth (Hyndman and Peacock, 2003; van Keken et al., 2012; Wada et al., 2015), a consequence of the rapid convergence of old oceanic crust (van Keken et al., 2012). Our study suggests that the Tc of titanomagnetite, within the range 150–350 °C, is reached by the extrusive basalt layer at 9–12 km bsl, and the Tc of magnetite, 580 °C, is reached by the deeper crustal layers when the slab surface is at 20 km bsl (i.e., at ∼22–26 km bsl, considering the initial depth of these layers). These temperatures are significantly higher than those predicted by most published models. They agree better with the background thermal gradient of 26.29° ± 0.13 °C/km measured on the overriding plate near the trench (Fulton et al., 2013), which predicts 170 °C at 10 km bsl (6.5 km below seafloor) and 460 °C at 20 km bsl (17.5 km bsf).

A recent study explains the high heat flow before subduction by rejuvenated hydrothermal activity induced by the flexural bending of the lithosphere and associated normal faults (Kawada et al., 2014). After entering subduction, the convecting fluid mines heat from depth but is prevented to transfer this heat upward by the thick sediments of the accretionary prism, the subducted pelagic sediments, and possibly an impermeable décollement surface, warming up the subducting oceanic crust through thermal blanketing (e.g., Granot and Dyment, 2019). Additional heat may be provided by serpentinization, i.e., the hydration of mantle rocks by the percolating fluid, an exothermic transformation most active between 100 and 400 °C (Macdonald and Fyfe, 1985) reported in various active tectonic contexts involving the oceanic crust—slow seafloor spreading (Kelley et al., 2001), compressive deformation (Delescluse and Chamot-Rooke, 2008), and subduction (Grevemeyer et al., 2018). Conversely, the frictional heat generated by recurring megathrust earthquakes, although reaching high temperatures (Yang et al., 2016) for a short time at a very local scale, remains negligible if the estimate of Fulton et al. (2013) is considered.

A warmer subducting oceanic crust has implications on the rheology of the subducting slab. The brittle-ductile transition in basalt occurs at 550 ± 100 °C (e.g., Violay et al., 2012). If our inference on the Curie temperatures are correct, this depth is reached at ∼20 km, explaining why most earthquakes of magnitudes stronger than 5 are observed at depths <20 km, as Figure 3C suggests, despite the poor vertical resolution of epicenter determinations.

We thank Jeff Gee, Ingo Grevemeyer, and an anonymous reviewer for their helpful comments. We also thank Koichiro Obana for providing precise earthquake catalogs even though we could not use them due to their limitation in deeper areas. Choe has been supported by an Ecole Doctorale STEP’UP Fellowship. We thank all scientists and crews who collected the marine geophysical data used in this study. This is Institut de Physique du Globe de Paris contribution 4089.

1GSA Data Repository item 2020076, supplementary information on the method, and discussions of (1) crossover correction for scalar and vector magnetic anomalies, (2) uncertainties in the slab geometry, and (3) the effect of a dipping slab on the inclination and declination of magnetization, is available online at http://www.geosociety.org/datarepository/2020/, or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY license.