Introduction

An interaction of the Pacific, Okhotsk, Eurasian and Philippine Sea lithosphere plates with the deeper mantle around the Japanese islands (Fig. 1) is complicated by active subduction of the plates1,2 and back-arc spreading3, which cannot be understood by the plate kinematics only. The Pacific plate subducts under the Okhotsk and the Philippine Sea plates with the relative speed of about 9 cm yr−1 and 5 cm yr−1, respectively, whereas the Philippine Sea plate subducts under the Eurasian plate with the relative speed of about 5 cm yr−1 (ref. 4). Back arcs of these subduction zones are also known as the site of active spreading in the past and recent as inferred from both the geophysical and geological surveys3. Elucidating a cause of the Japan Sea back-arc opening is one of scientific challenges in geosciences. Although its kinematic description and/or qualitative images of dynamics are fairly well understood3, quantitative interpretations based on sound physical and other principles are yet missing.

Figure 1
figure 1

Topography map of the Japanese Islands and surroundings.

The plate motions and deformations are presented by arrows. The rate of the motions is determined from geodetic data4 and for the Philippine Sea plate from the PB2002 model23. The white star indicates the place of sampling for geochemical analysis47.

Full size image

P-wave seismic tomography of the mantle beneath the subducting Pacific plate near the Japanese islands revealed a low velocity region extending oceanward at depths around the 410-km seismic discontinuity and this low velocity anomaly region was interpreted as a zone with an excess temperature of 200 K and the associated fractional melt of less than 1% (Ref. 5). To clarify the origin of the hot temperature anomaly beneath the Pacific plate and its implication for back-arc basin evolution, we study the mantle evolution beneath the Japanese islands and their surroundings based on the assimilation of temperature inferred from seismic tomography1, the present movements derived from geodetic observations4 and the past plate motion inferred from paleogeographic and paleomagnetic plate reconstructions6,7,8,9.

A quantitative restoration of mantle evolution in the geological past using present geophysical, geodetic and geological data becomes a key to unravelling the Earth's mantle history. Data assimilation methods employed for geodynamic restorations allow constraining the conditions for the mantle temperature and flow velocity in the past from the present temperature and plate velocity using an explicit dynamic model10. The temperature in the past is considered as a control variable to optimise the temperature minimising the discrepancy between the temperature and velocity field at present and those derived from the solution of the dynamic model. Data assimilation applied recently to restore the evolution of the relic slab beneath the southeastern Carpathians11, of the Farallon plate12,13 and of the dynamic topography of Africa14 provides a bridge between the past and present and serves as a tool for the structural geology of the mantle.

The problem of assimilation of regional observations into a mantle convection model to restore quantitatively the thermal structures in the geological past beneath the Japanese islands and their surroundings is reduced here to solving backwards in time the numerical problem described by a set of discrete, regularized, thermal convection equations with the prescribed initial (the present temperature model) and boundary (e.g., the plate motion) conditions in the three-dimensional model domain [0,4000 km] × [0,4000 km] × [0,800 km] (see Methods and the Supplementary Material).

Results

The present temperature model beneath the Japanese islands is developed by using the high-resolution seismic tomography (P-wave velocity anomalies) for the region1,5,15. The seismic tomographic model consists of 16 horizontal layers of different thickness (12 km up to 88 km) starting from the depth of 12 km and extending down to the depth of 800 km. Each layer of 4000 × 4000 km2 is subdivided horizontally into 80 × 80 km2 blocks. To restrict numerical errors in our data assimilation, the velocity anomaly data are smoothed between the layers using a spline interpolation. The temperature anomalies are inferred from the seismic wave anomalies using a non-linear inversion method16 and considering the effects of mantle composition, anelasticity and partial melting on seismic velocities17. The temperature anomalies are then added to the vertical temperature profile (the background temperature; see Supplementary Material) to obtain the present temperature model. At shallow depth (down to 110 km), the background temperature is the solution of the cooling half-space model18 with 48 million years after the start of cooling. This temperature field gives a heat flow of 70 mW m−2, which is close to the mantle heat flux. At deeper level, the background temperature is presented by the adiabatic temperature distribution19 with the potential temperature of 1330°C. The crustal temperature inferred directly from the seismic velocity anomalies is relatively high and the resulting heat flow is unrealistically elevated. Thus, the crustal temperature Tcr is calculated by estimating the geothermal gradient δTm/δx3 from measured regional surface heat flow20,21 (Ts is the surface temperature): .

The present thermal state of the back-arc region so obtained is characterized by shallow hot anomalies reflecting the remnants of the back-arc spreading3 and deep cold anomalies related to the stagnation of the lithospheric slabs1 (Fig. 2). Meanwhile the state of the mantle beneath the Pacific plate is characterized by the shallow cold anomalies reflecting the existence of the old oceanic Pacific plate and the deep broad hot anomaly5 of unknown origin. This temperature model is used as the initial condition for restoration models.

Figure 2
figure 2

Thermal evolution of the mantle beneath the Japanese Islands.

The model covers 4000 km (horizontal) × 4000 km (horizontal) × 800 km (depth). The top panel presents the plate motion prescribed in the model. Horizontal cross-sections present temperature anomalies and the projection of mantle flow, which are obtained by the assimilation of the present temperature to the Middle Eocene times. Yellow closed curves indicate the hot mantle. White boxes (bold and dashed) show the places of mantle flow changes are likely to be associated with the rotation of northeastern and southwestern parts of the Japan Sea in the Miocene times. See the text for characters A, B, C, a, b and c.

Full size image

Present and past plate motions

Although the mantle dynamics is coupled to the lithosphere dynamics, the coupling can be weak or strong depending on the viscosity contrast between the lithosphere and the underlying mantle22. To explore how plate kinematics is related to the dynamics of the mantle, we assume kinematic conditions at the upper model boundary, where the direction and rate of plate motion are prescribed. The plate motion velocity is determined from the Actual Plate Kinematic and Deformation Model (APKIM2005) derived from geodetic data4 and from the PB2002 model for the Philippine Sea and Okinawa Plates23. The reference frame for the velocities is the International Terrestrial Reference Frame (ITRF, http://itrf.ensg.ign.fr), which is realized by a set of stations on the Earth's surface. The reference for the velocities is defined with respect to Earth rotation, that is, there is no net rotation of the frame with respect to the rotating Earth24. The velocities are determined in APKIM2005 such a way to minimise a root mean square sum of all velocities over the Earth surface. The recent MORVEL model25 is similar to APKIM2005, although the velocities in the MORVEL model are taken from different sources and have not the same reference frame.

Plate motions in the geological past are required to determine past mantle dynamics processes. Our model incorporates the history of plate motions as the boundary conditions changing with time. Based on paleogeographic reconstructions of the Philippine Sea6 and Japanese Islands26, Philippine Sea plate motion from paleomagnetic studies8, relative motion of the Pacific plate9 and Cenozoic plate tectonic evolution of the southeastern Asia7, we construct the velocity field to define the plate motion in the model for the past 40 million years, meanwhile recognizing the uncertainties in directions and magnitudes of the regional plate motion in the mentioned studies. To simplify the model, we omit the present boundaries between the Okhotsk and Eurasian plates and between the North American and Pacific plates. Although change in the plate motion is a gradual process, for the sake of simplicity we assume that the plates moved with a constant velocity within four time intervals in the past: from 0 to 10 million years ago (0–10 Ma), 10–20 Ma, 20–30 Ma and 30–40 Ma (Fig. 2).

In numerical models of mantle dynamics, when a surface velocity is prescribed, the flow is driven by a combination of the velocity and internal buoyancy of the fluid. The choice of the prescribed velocities affects the flow velocities in the uppermost mantle. To analyse the influence of the prescribed velocity, we varied the velocities changing the magnitude of the velocities (within about 3%) and direction (within about ±2 degrees) The small changes in the imposed velocities are found not to change significantly the mantle flow and hence the pattern of upwellings. We observed that the additional work introduced by the prescribed velocity at the upper boundary is dissipated within the uppermost part of the model domain and has little influence on the thermal convection below27.

Rheology

Mantle rheology is rather complex and depends on temperature and pressure28, stresses29, water content30, grain size31 and composition32. Slow deformation of minerals in the upper mantle under the influence of stresses and temperature is governed by diffusion and dislocation creeps29. Diffusion creep dominates deformation at cooler temperatures and larger grain sizes, whereas dislocation creep dominates deformation at higher strain rates and is not grain-size dependent33. Diffusion-controlled creep is characterized by the Newtonian rheology, meanwhile dislocation creep by a non-Newtonian rheology (e.g., a power-law relationship between stress and strain rate). If in the upper mantle diffusion and dislocation creep mechanisms act together to accommodate deformation, the lower mantle is likely to deform by diffusion creep as seismic anisotropy studies detect less seismic anisotropy in the deep mantle34. Effects of a non-Newtonian viscosity on steady and time-dependent convection have been studied intensively using 2-D numerical models35,36,37,38 and 3-D numerical models39,40. For a stationary convection of a temperature- and pressure-dependent viscosity fluid it was shown that the use of the Newtonian viscosity with activation enthalpy one-third to half of the experimentally determined value for olivine could mimic the dynamics of the convection with a strongly non-Newtonian power-law viscosity35. Although a non-stationary thermal convection introduces changes in the pattern of viscous flow, the qualitative behavior of the flow is similar and the difference comes in the way of time-dependence.

In the modelling we assume that the Earth's mantle behaves as a temperature- and pressure-dependent (Arrhenius-type) Newtonian fluid with the activation energy of 3 × 105 J mol−1, activation volume of 4 × 10−6 m3 mol−1 (see the Supplementary Material) and the viscosity beneath the lithosphere being around 3 × 1019 Pa s. According to experimental studies of olivine deformation, the values of activation energy and activation volume in a dry upper mantle are 3 × 105 J mol−1 and 6 × 10−6 m3 mol−1 for diffusion creep and 5.4 × 105 J mol−1 and 2 × 10−5 m3 mol−1 for dislocation creep, respectively17. Recent measurements of the activation volume for creep of dry olivine at pressures of 2.7–4.9 GPa and temperatures near 1473° K showed a wide range of the values 9.5 ± 7 × 10−6 m3 mol−1 (ref. 41). For the prescribed model values of the activation energy and the activation volume the model activation enthalpy is less than half of the experimentally determined values and therefore, the Newtonian fluid with the chosen model values can approximate to some extent the mechanical behaviour of the upper mantle. The inclusion of “realistic” rheology of the upper mantle in the model would be preferable. However, because the complex (non-linear) rheology is more susceptible to ambiguities (see Methods) associated with the temperature estimates inferred from seismic tomography and also because more parameters are requited to describe the non-Newtonian rheological model, we prefer to use the Newtonian rheology with reduced activation enthalpy to be conservative.

Reconstructing thermal structures in the past

Using the quasi-reversibility method for data assimilation42 the set of the discrete equations governing the backward mantle convection flow are solved together with the initial and boundary conditions (see Supplementary Material) to restore the mantle evolution in the region of the Japan islands. In backward sense, the high-temperature patchy anomaly beneath the back-arc Japan Sea basin splits into two prominent anomalies showing two small-scale upwellings beneath the southwestern and northern part of the Japan Sea (Fig. 2). The present hot anomalies in the back-arc region marked by A and B move down to the spots marked by a and b at 38.9 Ma (Fig. 2). Meanwhile the broad hot anomaly under the Pacific plate moves slowly down westward as depicted by C at present and c at 38.9 Ma in Fig. 2. The hot anomalies (spot b) and (spot c) tend to merge at 38.9 Ma and below 560 km. The model shows the link between spot b in the back-arc region and spot c in the sub-slab mantle at depths of about 440–560 km in the Middle to Late Eocene times (see the movie in the Supplementary Material). The upwellings a and b are likely to be generated in the sub-slab hot mantle (Figs. 2 and 3) and penetrated through breaches/tears of the subducting Pacific plate into the mantle wedge toward spots A and B. Hence, we propose that the present hot anomalies in the back-arc and sub-slab mantle had a single origin located in the sub-lithospheric mantle.

Figure 3
figure 3

Schematic representation of the thermal evolution beneath the Japanese islands and their surroundings.

See the text for explanation.

Full size image

It is also noteworthy that the shallow mantle velocity field in the back arc, which potentially affects the surface tectonics there, shows rather complex pattern with time compared to that in other area. The southwestward uppermost mantle flow (the boxed area at depth 80 km of Fig. 2) beneath the northeastern (NE) part of the Japan Sea region in the Middle Eocene – Oligocene times, 38.9–25.7 Ma, changed to the south-southeastward flow in the Miocene times, 12.6 Ma (Fig. 2). This flow change could contribute to counterclockwise rotation of the NE Japan Islands, which resulted in the back-arc basin opening3,43,44. In the southwestern (SW) Japan, clockwise rotation was likely to be instantaneous at about 15 Ma44, which led to opening of this part of the Japan Sea. Although the model cannot explicitly predict the changes in the mantle flow in the SW Japan Sea region, the flow pattern in the uppermost mantle (the dash boxed area at 80 km depth of Fig. 2) indicates that the northwestward flow established in Oligocene becomes much weaker in Miocene.

Discussion

An analysis of the onshore and offshore borehole data from South Sakhalin, the Oga Peninsula and Dolgorae shows that the regional tectonic subsidence commenced in the Early Oligocene and continued until the Pliocene with a few breaks for a tectonic uplift (Fig. 4; ref. 45). The rapid subsidence of South Sakhalin (in the Early Miocene, some 23–22 My ago) was followed by the slow subsidence and uplift until about 16 My ago, which was interrupted by the second phase of rapid subsidence in the Middle Miocene (15 My ago) followed again by slow subsidence and uplift. The alternation of subsidence and uplift phases in the basin is indicative of rifting style changes from rollback-induced passive rifting/extension to upwelling-induced (upwelling B, Figs. 3 and 4) active rifting/extension46. Thus the small-scale upwelling beneath the northern part of the Japan Sea (predicted by the data assimilation modelling) could contribute to the rifting and back arc basin opening. The evolution of the southwestern Japan Sea recorded in the sediments (drilled by the offshore Dolgarae-1 borehole) shows a persistent subsidence of the basin until about 11 My ago followed by rapid uplift. This uplift is likely to be associated with the small-scale upwelling A (Figs. 3 and 4). Meanwhile preceded by a long phase (about 15 My) of a slow subsidence, the relatively rapid subsidence of the Oga Peninsula (which is located in-between two upwellings A and B) started in the Middle Miocene (about 16 My ago) followed by slow subsidence for about 13 My and then rapid uplift in Quarternary (about 2 My ago).

Figure 4
figure 4

Subsidence history of the Japan Sea: tectonic subsidence curves45 for the southern Sakhalin section representing the northern margin of the Japan Sea (dotted curve), the Oga Peninsula section representing the inner arc area of north-western Honshu (solid curve) and Dolgorae-1 well representing the southern Tsushima Basin and the southern margin of the Japan Sea (dashed curve).

The location of the wells is marked in (b). (b) Three-dimensional view of snapshots of the iso-surfaces of positive (5%) temperature anomalies; colours mark the depth.

Full size image

If our hypothesis, that the hot mantle upwellings under the Japan Sea originated beneath the Pacific plate and penetrated through the subducting lithosphere contributing to the back-arc spreading, is valid, a signature of such a phenomenon may be found in magmatic rock samples. The change in the mantle source of magmatic rocks related to the opening of the Japan Sea from enriched to depleted with time43 suggests a possibility for a physical replacement of the enriched subcontinental upper mantle with the depleted asthenospheric mantle by an injection during back-arc development. The sample showing the mixture of asthenospheric and lithospheric slab material is found47 near the triple junction of the Pacific, Okhotsk and Philippine Sea plates (shown by star in Fig. 1). The age of this sample is about 23 Ma and corresponds to early stages of the Japan Sea opening. These observations can be explained by our hypothesis of hot upwellings coming from the Pacific side.

Presently, there is no consensus on the origin of hot anomaly under the Pacific plate. It can be a relic of the hot plume originated at the thermal boundary layers, although its closeness to the cold subduction zone seems to reject the nearby origin of hot upwellings. It can be a relic of large plume head conveyed horizontally by the plate movement or the hot thermal anomaly close to the sinking plume typical to the internally heated convection. The previous study shows that the latter scenario is more likely48. There is also a suggestion that the high temperature anomaly adjacent to the cold plume in internally heated convection may even reverse the overall direction of flow49 implying that the dynamical effect of such a high temperature anomaly is fairly significant.

It is known that the opening of the Japan Sea is not a simple homogeneous and instantaneous spreading, but has temporal and spatial variations and shows several stages of deformation3. Inhomogeneous spreading is consistent with the patchy character of hot materials as evident in our results. Non-instantaneous deformation may imply that the hot materials have penetrated through, or affected, the overlying subducting Pacific lithosphere several times.

The slab break-off, detachment and tear have been intensively investigated for the last decades. A lithospheric slab tear has been proposed in the Indonesian arc50 and in the Mediterranean region51,52. A horizontal tear in the subducting slab has been recognized in the Izu-Bonin-Mariana arc2,5,53. There are numerical studies54 related to this phenomenon, but they usually do not assume the existence of hot anomalies under the subducting lithosphere. A result of recent study48 shows a possibility of slab break-off by hot materials in the sub-lithospheric mantle. However, we note that our results do not show the slab break-off but the penetration via breaches/tears in the subducting lithosphere, which is probably more likely than the total break-off.

The geometry of the thermal structures in the mantle changes with time due to heat advection, which deforms the structures and heat conduction, which smoothes the complex shapes of the structures. This creates difficulties in understanding the evolution of the mantle structures in the past. A quantitative assimilation of the present mantle temperature and flow into the geological past provides a tool for restoration of thermal shapes of prominent structures in the past from their diffusive shapes at present. Using the present seismic tomography data as well as some other geological, geophysical and geodetic data and constraints, we were able to reconstruct prominent features of hot upwelling beneath the Japan Sea and Japanese islands.

There are several principal qualitative models proposed to explain back-arc spreading: mantle diapirism55, convection induced by the subducting slab56, movements of adjacent plates57 and the injection of hot asthenospheric material43. Our model (Fig. 3) as the latter one supports the hypothesis of an asthenospheric injection. Meanwhile the major difference of the model compared to our model is related to where the hot material came from. While the previous model43 assumes that it may have come from the back-arc side, the results of our dynamic restoration of the thermal state of the mantle in the past show that it came from the opposite side, that is, from the mantle under the Pacific plate.

High-resolution experiments on seismic wave attenuation, improved knowledge of crustal and mantle mineral composition and enhanced paleo-reconstruction models of plate movements in the region will assist to refine the present model of the mantle thermal evolutions beneath the western Pacific region. The basic knowledge we have gained from the study is the past mantle dynamics beneath the Japanese islands, which could lead to its present dynamics. The results of our study may have an important implication for understanding back-arc spreading, tearing of subducting slabs and geochemical mixing of mantle rocks. Although our principal finding is related to mantle dynamics beneath Japan and its surroundings, such a phenomenon (a linkage between the mantle beneath the slab and above it) may be rather ubiquitous in other subduction zones.

Methods

Numerical methods

The model domain is a bottom-heated Cartesian box 800 km deep extending 4000 km in both horizontal directions. The equation of energy conservation (or the heat balance equation) with a small additional term responsible for heat relaxation is integrated backward in time using the quasi-reversibility method42. The additional term represents mathematically a higher order differential term and regularizes the ill-posed backward heat balance equation. The numerical solution is inferred by optimisation of the misfit between observations and model solution. The equations of mass and momentum conservation with an infinite Prandtl number at the extended Boussinesq approximation are coupled to the energy balance equation and the full system of equations together with the prescribed boundary and initial conditions are solved numerically using the finite-volume method10 (see the Supplementary Material for more detail).

Analysis of data uncertainties

Data assimilation can be influenced primarily by uncertainties of a present temperature model (initial condition) and uncertainties in boundary conditions. Uncertainties in a present temperature model come from seismic tomography models, mantle composition, seismic attenuation models and poor knowledge of the presence of water or melt content at mantle depths10. Seismic tomography models may introduce some uncertainties in data assimilation due to a resolution power. A qualitative comparison of the images related to three seismic tomography models (WEPP25, MIT58 and GAP_P215; Fig. S11 in the Supplementary Material) shows that their resolution power is almost similar. Although the resolution of the MIT model is a bit higher than WEPP2 and GAP_P2, the image of the slab and that of the sub-slab hot anomaly in the MIT model are less pronounced. The temperature anomaly under the Pacific plate estimated from WEPP2 (the present model) is nearly the same as that estimated by using the different approach5. Thus, we used WEPP2 in this study.

The temperature at the lower boundary of the model domain used in the numerical modelling is an approximation to the real temperature, which is unknown and may change over time at this boundary. The conditions at lateral boundaries are prescribed to satisfy certain properties of the model (e.g., conservation of mass, conservation of momentum), meanwhile the true conditions at the boundaries are unknown and hence contribute to uncertainties of the model. The velocities prescribed at the upper boundary of the model domain comes from recent geodetic observations and paleogeographic/paleomagnetic reconstructions of plate motion and hence the uncertainties in the prescribed velocity field are associated with an accuracy of geodetic measurement, viscosity contrast between the lithosphere and the underlying mantle and uncertainties in paleo-reconstructions. In numerical modelling sensitivity analysis assists in understanding the stability of model solutions to small perturbations in input variables or parameters.

We performed a sensitivity analysis to understand how stable is the numerical solution to small perturbations of input data (the present temperature). The model of the present temperature has been perturbed randomly by 0.5 to 1% and then assimilated to the past to find the initial temperature. A misfit between the initial temperatures related to the perturbed and unperturbed present temperature is about 3-5%, which proves the stability of the solution.

For a given temperature, mantle dynamics depends on various factors including rheology33, phase changes59,60 and boundary conditions10. We conduct a search over the ranges of uncertain parameters in the temperature- and pressure-dependent viscosity (activation energy and activation volume) to achieve ‘plate-like' behaviour of the colder material (see Supplementary Material for numerical tests on viscosity variations). When restoring the mantle dynamics, we observed that a highly viscous descending slab (5 × 1021 Pa s about 100 times higher than that of the surrounding asthenosphere) causes a slow rise of the slab (in a backward sense); alternatively, a less viscous slab (5 × 1020 Pa s about 10 times higher than that of surrounding asthenosphere) generates a faster and unrealistic vertical rise. Nevertheless, the small-scale upwellings in the models show persistency despite the variations in the slab viscosity.

Also we tested the influence of phase changes, model depth variations and boundary conditions on the model results (see the Supplementary Material). In numerical models of mantle dynamics the choice of boundary conditions and the size of the model domain influence the pattern of flow and slab dynamics10,33. If the depth of the model domain is significantly smaller than the horizontal dimensions of the domain, the thermo-convective flow in the model with a low viscosity upper mantle and higher viscosity lower mantle will generate the return flow focused in the upper mantle. Increasing the model domain's depth removes the artificial lateral return flow in the upper mantle. To study a sensitivity of the model results to the conditions imposed on the lower boundary, we varied the depth of the model domain from 800 km to 2000 km (see Supplementary Material for numerical tests on model depth variation). The sensitivity analysis related to the presence of phase transformations and to changes in boundary conditions show that the model is robust and the principal results of the model do not change.