Thermal Structure and Intermediate-depth Seismicity in the Tohoku-hokkaido Subduction Zones

The cause of intermediate-depth (> 40 km) seis-micity in subduction zones is not well understood. The viability of proposed mechanisms, which include dehydration embrittlement, shear instabilities and the presence of fluids in general, depends significantly on local conditions, including pressure, temperature and composition. The well-instrumented and well-studied subduction zone below Northern Japan (Tohoku and Hokkaido) provides an excellent testing ground to study the conditions under which intermediate-depth seismicity occurs. This study combines new finite element models that predict the dynamics and thermal structure of the Japan subduction system with a high-precision hypocenter data base. The upper plane of seismicity is principally contained in the crustal portion of the subducting slab and appears to thin and deepen within the crust at depths > 80 km. The disappearance of seismicity overlaps in most of the region with the predicted phase change of blueschist to hydrous eclogite, which forms a major dehydration front in the crust. The correlation between the thermally predicted blueschist-out boundary and the disappearance of seismicity breaks down in the transition from the northern Japan to Kurile arc below western Hokkaido. Adjusted models that take into account the seismically imaged modified upper mantle structure in this region fail to adequately recover the correlation that is seen below Tohoku and eastern Hokkaido. We conclude that the thermal structure below Western Hokkaido is significantly affected by time-dependent, 3-D dynamics of the slab. This study generally supports the role of fluids in the generation of intermediate-depth seismicity.


Introduction
The subduction of the Pacific plate below central and northern Japan (Fig. 1) provides a classic cold end-member in the thermal structure of subduction zones, due to the high convergence rate and the old age of the oceanic lithosphere (e.g., Peacock and Wang, 1999).The available high quality seismic catalogues provide excellent opportunities to determine the seismic structure of the subduction zone below the Japanese islands and the possible relation between intermediate-depth seismicity and dehydration of the oceanic crust.The low velocity of the hydrated oceanic crust provides a sharp velocity contrast between crust and mantle that is imaged well in receiver function studies (e.g., Kawakatsu and Watada, 2007) and tomographic studies (Zhang et al., 2004;Tsuji et al., 2008;Nakajima et al., 2009).Intriguingly, the upper plane of seismicity seems to be confined to the oceanic crust (Fig. 2; Kita et al., 2006Kita et al., , 2010a)), with a marked reduction in number and deepening of the earthquakes with respect to the top of the crust at greater depth.A similar pattern, and correlation with low seismic velocities in the crust, is seen in the central Alaska subduction zone (Abers et al., 2006).Thermalpetrological modeling following van Keken et al. (2011) strongly suggests that earthquakes in a particular cross section below northern Tohoku are limited by the blueschist-tohydrous-eclogite phase boundary (Fig. 2), which is the first main deep dehydration front in cold subduction zones (van Keken et al., 2011).It is important to note that the location of the blueschist-out boundary in Fig. 2 is based on independent thermal-petrological modeling and is not a fit to the observed seismicity.This independent correlation suggests the hypothesis that, at least in cold subduction zones, the upper Published by Copernicus Publications on behalf of the European Geosciences Union.Zhao et al., 1997;Kita et al., 2010b).Heat flow (Tanaka et al., 2004) is shown with the colored circles.Volcanoes are shown in black triangles.The down-dip limit of low-angle thrust events from Kita et al. (2010a), which is taken to be the down-dip limit of the region of strong seismic coupling between upper and lower plate, is shown in the bold red line.Cross sections are labeled following Kita (2009) with P2-P23 over Hokkaido and T2-T25 over Tohoku.Arrow shows Pacific (PA) to North America (NA) convergence.
plane of seismicity is not only confined to the descending oceanic crust, but also related to the presence of fluids that are liberated by the transformation of blueschist to eclogite.Support for this idea is provided in the qualitative sense by the deepening of the seismic belt below Kanto and western Hokkaido (Hasegawa et al., 2007;Kita et al., 2010a).Thermal shielding of the slab by an anomalous overriding plate structure is supported by heat flow data (Fig. 1; Tanaka et al., 2004).The correlation between the predicted blueschistout boundary and the cessation of seismicity is even more intriguing since the eclogite that forms still contains signifi- cant H 2 O (∼ 2-3 wt %).The crust is generally assumed to be dry after the lawsonite-out reaction, which does not happen until significantly higher temperature (∼ 650 • C at 3 GPa, and at higher temperature at higher pressure; Schmidt and Poli, 1998;Forneris and Holloway, 2003).The blueschist-out boundary provides a relative sharp dehydration front, where over a short distance in pressure and temperature space the H 2 O content drops by a few wt % (Hacker, 2008;van Keken et al., 2011).At this stage we do not understand the fundamental reasons for the good correlation between blueschistout and the cessation of seismicity, but we will assume that this is related to a significant pulse of fluids upon blueschist dehydration.We will use this hypothesis to test whether significant along-arc variations exist in seismicity that indicate changes in thermal structure.
The main goal of this paper is to develop new models to determine the thermal structure of the subducting Pacific slab below Tohoku and Hokkaido.We study along-arc variations in the correlation between the blueschist-eclogite boundary and test whether simple thermal shielding proposed by Kita et al. (2010a) leads to the observed deepening of the seismic belt.

Modeling approach
We use two-dimensional finite element models with strong grid refinement to determine the thermal structure in the Fig. 3. General model geometry (top) and partial mesh (bottom) used in the finite element modeling.The velocity is prescribed in the slab and set to zero in the crust of the overriding plate (with the Moho indicated by M).The wedge is modeled as a fluid with a rheology appropriate for dry olivine (as in van Keken et al., 2008).The wedge flow is driven by coupling with the slab.As in Syracuse et al. (2010) and Wada and Wang (2009) we assume in general that this depth is 80 km.Between the down-dip limit of the earthquake thrust zone (bold red line in Fig. 1) and 80 km depth we assume partial coupling.The governing equations are discretized on a finite element mesh with a finest resolution of 1 km along the slab surface and in the crust.The Stokes equations are solved in the dynamic wedge only, with boundary conditions indicated in red.The heat equation is solved on the entire domain, with boundary conditions indicated in blue.In one particular case we assume the presence of a subducted fore-arc sliver as indicated by the grey lines.T s is surface temperature (assumed to be 293 K here), T GDH1 is the temperature of the half-space cooling model of Stein and Stein (1992).u is velocity, σ n is normal stress, σ t is tangential stress, u slab is the imposed slab velocity.
Tohoku and Hokkaido subduction systems (Fig. 3).We use a selection of cross sections from Kita et al. (2010a) (Fig. 1).For each cross section we determine the slab surface geometry from the 3-D slab surface geometry provided by Zhao et al. (1997) and Kita et al. (2010b) (Fig. 1).Details on the interpolation are provided in Appendix A. We use a simplified Moho representation that is adapted from Katsumata (2010).We show the Moho in the cross sections below by a solid black line indicated with the symbol M.
The finite element models follow the approach used by Syracuse et al. (2010) andvan Keken et al. (2011).We use a kinematic-dynamic approach where the slab is prescribed as a kinematic entity.Slab velocities are parallel to the slab surface.The crust over the overriding plate is assumed to be rigid.The wedge is modeled as a fluid with a dislocation creep rheology that is appropriate for dry olivine follow-ing Karato and Wu (1993), that is, η(T ,II) = A exp(E/nRT ) II (1−n)/n , where η is viscosity, II is the second invariant of the strainrate tensor, A = 28968.6Pa • s 1/n , E = 335 kJ mol −1 , n = 3.5 and R is the gas constant.We use a maximum viscosity cut-off of 10 25 Pa s.For the continental crust we assume density 2700 kg m −3 , conductivity 2.5 W mK −1 , and a twolayer heat production of 1.31 µW m −3 above 15 km depth and 0.27 µW m −3 below this depth.For the rest of the model (including oceanic crust) we assume a density of 3300 kg m −3 , conductivity 3.1 W mK −1 and zero heat production.Since we ignore thermal buoyancy, the thermal expansivity is effectively zero.The full details on the governing equations are provided in Syracuse et al. (2010).The wedge fluid is entrained by the subduction of the slab and the resulting corner flow brings high temperatures to the slab surface (van Keken, 2003).The speed of the slab is the projection of the plate velocity onto the cross section, leading to an effective reduction in plate convergence speed.As is shown in Bengtson and van Keken (2012) this approach of taking cross sections perpendicular to the trench and using the projected velocity component is preferred over the approach where the cross section is taken parallel to convergence velocity.As in Wada and Wang (2009) and Syracuse et al. (2010) we assume in general full coupling between slab and wedge at depths below 80 km.We assume full decoupling between slab and overriding plate along the thrust zone, which has a down-dip limit as shown in Fig. 1.In between we assume partial coupling, where 5 % of the slab velocity is transmitted to the overlying wedge.The effective subduction speed and depth of the decoupling zone for each cross section are provided in Table 1.Note that our use of the terms "coupling" and "decoupling" is opposite to that in the seismological community.At geological time-scales relevant to this study the seismologically "coupled" zone acts as a fully decoupled zone due to large underthrusting events.
The equations that we solve are provided in van Keken et al. (2002) and Syracuse et al. (2010).We solve the Stokes equations (which satisfy the conservation of mass and momentum) in the mantle wedge.We solve the time-dependent heat advection-diffusion equation (satisfying the conservation of energy) on the entire domain.As in Stein and Stein (1992) we use a potential mantle temperature of 1421.5 • C. At the trench inflow we assume the thermal structure of a 130 Ma oceanic lithosphere following the GDH1 plate model (Stein and Stein, 1992).At the arc inflow boundary we assume the thermal structure of 30 Ma old oceanic lithosphere, which we deem appropriate for the back-arc setting of the Japan Sea.We use the two-layer continental heat production (as in Syracuse et al., 2010) but ignore shear heating along the thrust zone.We assume that the sedimentary layer has a thickness of 0.6 km at the trench with a compaction to 0.3 km below 15 km depth.We evolve the model for 20 Myr, which is sufficient for the slab to reach a steady-state thermal structure.Further time integration would only change the solution by the very slow growth of the overriding lithosphere.
Table 1.Starting longitude (E) and latitude (N), strike, effective convergence speed, thermal parameter, and down-dip limit of the thrust zone for each cross section.The convergence speed is based on the projection of the convergence velocity onto the cross section.We use the Pacific-North America convergence velocity from NUVEL1-A (DeMets et al., 1994)  We solve the governing equations on a finite element mesh with strong grid refinement in the boundary layers (Fig. 3).The mesh is composed of Taylor-Hood triangles with linear velocity and linear pressure for the Stokes equations.The same geometry is used for solution of the heat equation with streamline upwinding.The mesh has the highest resolution (1 km spacing) along the top of the slab at shallow and intermediate depths (which is the main focus of our study).The coarsening of the mesh away from the boundary layers allows for efficient yet accurate solution of the equations.We performed convergence tests using meshes that had a finest spacing of 750 m and almost twice the number of elements.We found differences of less than 1 • C in the temperature prediction for the oceanic crust, which is consistent with the performance of our models in a recent benchmark for subduction zones (van Keken et al., 2008).
The projection of seismicity onto each cross section requires some special care, particularly in regions where the slab surface geometry is dipping significantly in the direction normal to the cross section (e.g., P9 in central Hokkaido).The method we use to map slab surface and earthquake hypocenters in a 22 km wide swath around each cross section is described in Appendix A. If we were to use a direct projection of the earthquake data from 3-D onto the cross section, we would incur spatial errors of several kilometers, which are above the precision of the hypocenter data.

A typical 2-D subduction zone model for Tohoku
An example of the new thermal models is shown in Fig. 4 for cross section T18 over northern Tohoku (approximately parallel but slightly to the south of the cross section in Fig. 2).Seismicity is projected onto the cross section from within 11 km on either side of the cross section.The upper-plane seismicity (defined here loosely as to within 10 km from the top of the slab) is shown in with the red circles.The model follows the general characteristics of the "D80" models of Syracuse et al. (2010) closely.The coupling of the slab to the overlying wedge (at 80 km depth) leads to strong thermal gradients at the top of the slab below this depth and draws in hot mantle material, with the hottest part of the wedge below the volcanic arc (which sits here at ∼ 100 km above the Wadati-Benioff seismicity).Heat flow observations are sparse in this region but the predicted heat flow follows the large scale pattern, with low fore-arc values and higher but scattered values in the arc and back arc.

Comparison between the predicted blueschist-out boundary and upper-plane seismicity
We use the thermal models to predict the blueschist-out boundary.Following Hacker (2008) we constrain this boundary by a temperature T = 617 − 52P , where T is in • C and pressure P is in GPa.The hydrostatic pressure is computed assuming a crustal density of 2700 kg m −3 and a mantle density of 3300 kg m −3 .We assume mantle density below the variable thickness Moho and therefore do not take into account the variable density in the subducting crust.This introduces an uncertainty of less than 1 km in depth on the position of the blueschist-out boundary.We plot the predicted boundary and the seismicity for a number of cross sections below Tohoku and eastern Hokkaido in Fig. 5.The blueschist-out boundary is shown for the oceanic crust with an assumed thickness of 7 km for cross sections.As in Fig. 2, the seismicity of the upper plane is significantly reduced above this boundary.
In most cases there is an increase in seismicity just above where the blueschist-eclogite boundary starts, suggesting that the heat from the mantle wedge in this region causes a spike in seismicity.There is a scattering of low magnitude seismicity above this boundary and deeper in the slab.It is tempting to consider the line of seismicity that extends almost linearly down from the blueschist-out boundary in P2 and P5 (indicated by blue question-mark), and to a lesser extent in T18 and T25, to be related to the same phase boundary.If this correlation is correct, the local thickness of the oceanic crust in sections P2/P5, and perhaps in T18/T25, may be slightly larger than the 7 km thickness assumed here.In P2 there is a somewhat more significant cluster of seismicity (at a depth of 120 km) that is in the region where the thermal models predict eclogite.It is intriguing that this region seems connected to the second plane of seismicity through a number of larger events and could potentially signal the release of fluids from the deeper slab.
While the seismicity below Tohoku (T2-T25) is somewhat sparser than that below central and eastern Hokkaido (P2-P9), we can discern a trend where seismicity tends to be a bit further below the blueschist-out boundary to the south (compare for example T25 with T2 or P9).This may indicate regional variations in thermal structure and perhaps slight changes in dynamics.For example, a shift of the full coupling point (assumed at 80 km here) to 75 km depth would bring the blueschist-out boundary closer to the top of the seismicity.This suggests we could use along-strike variations in seismicity to develop a better understanding of regional variations in subduction zone thermal structure and dynamics.

Seismicity and thermal structure at the trench junction
At the junction between the Japan arc and Kurile arcs (Fig. 1; cross sections P15-P23) the correlation between predicted blueschist-out and upper-plane seismicity seems to be less clear (Fig. 6).In particular in P15 the geometry of upperplane seismicity seems nearly independent of the slope and position of the blueschist-out boundary.Intriguingly, seismicity is well below the top of the slab at shallow depths (80-100 km) with an increase in seismicity at shallower levels in the crust below 100 km depth.There are several indications that the thermal structure of the crust and mantle below P15 are anomalous.Low heat flow values (Fig. 1; Tanaka et al., 2004) are significantly farther inland compared to elsewhere between 37-45 • N. Seismicity in the fore-arc wedge, which suggests cold conditions, also extends significantly farther from the trench (Fig. 5).Kita et al. (2010a) used tomographic techniques to image anomalously low P-and S-wave speeds below this section and demonstrated a significant deepening of the seismic belt here (their Fig. 4).Kita et al. (2010aKita et al. ( , 2012) ) suggested the region with low seismic velocity corresponds to a part of the fore-arc that is subducted in this region.Local geology supports the suggestion for crustal thickening.Kimura (1994) demonstrated that the Kurile fore-arc has been colliding with the NE Japan arc since the Middle Miocene due to the migration to the southwest of the Kuril forearc sliver, forming an arc-arc type collision zone.High temperature metamorphic rocks of Tertiary age are exposed and include materials from the upper and lower crust (e.g., Shimura et al., 2004).Ito (2000) and Iwasaki et al. (2004) suggested this region is characterized by delamination of the lower crust, with partial exhumation of the lower crust, and partial subduction along the Pacific slab.The combination of inferences from local geology and tomographic imaging suggests that the presence of subducted fore-arc materials cause thermal shielding of the slab to greater depths than elsewhere.We tested this hypothesis by modifying the crustal structure of the overriding plate in model P15.Specifically, we extended the crust down to include the suggested fore-arc sliver (using the geometry interpreted in figure 14 of Kita et al., 2010a).The resulting thermal structure (along with heat flow) is shown in Fig. 7c.The fore-arc sliver is imaged to a depth of just over 80 km.We assumed that the slab coupled fully to the overriding wedge at 90 km or 100 km depth.The modified predictions for the blueschist-out boundary are shown in Fig. 7c, corresponding (1) to the unmodified section (as in Fig. 5), (2) to the modified section with full coupling at 90 km depth; and (3) to the modified section with coupling at 100 km depth.While the blueschist-out boundary shifts to greater depth (essentially following the shift in depth of the coupling point), the match with seismicity does not improve significantly.At greater depths in the crust the three different predictions for blueschist-out tend to merge and the greatest differences are in the shallow portion of the crust, where seismicity is already sparse.Whereas the heat flow predicted from this model is appropriate in the back-arc, the predicted fore-arc heat flow is significantly higher than that observed, suggesting the fore-arc crust and upper mantle structure are significantly colder than the conditions we model here.

Discussion
The intermediate depth seismicity that we focused on in this study occurs at depths where lithostatic pressures are too high for normal brittle failure to occur.A number of processes have been suggested for causing local weakening.These include shear instabilities (Ogawa, 1987;Kelemen and Hirth, 2007), thermal runaways (John et al., 2010), and local weakening and embrittlement of the crust due to dehydration (e.g., Raleigh and Paterson, 1965;Kirby et al., 1996;Seno and Yamanaka, 1996;Hacker et al., 2003;Yamasaki and Seno, 2003).Possible consequences of dehydration include stress changes to the change in volume of the oceanic crust (Kirby et al., 1996), pore pressure changes (Wong et al., 1997) or just the release and transport of water (e.g., Jung et al., 2004;John and Schenk, 2006).
In general, our new thermal models provide strong quantitative support for the hypothesis that intermediate-depth seismicity is caused by the presence of fluids that are generated by metamorphic dehydration reactions.The predicted blueschist-out boundary tends to delimit the intermediatedepth seismicity in the upper plane below Tohoku and central/eastern Hokkaido.The simplest explanation for this correlation is that intermediate-depth seismicity is not possible at low water content in the downgoing oceanic crust.A similar correspondence between upper-plane seismicity and the blueschist-out dehydration boundary is suggested in Alaska (Abers et al., 2006;Abers et al., 2012) although intriguingly in this subduction zone the seismicity tends to align with the predicted boundary, rather than just delimit it.
The wide-scale presence of seismicity in the blueschistfacies crust in Northern Japan favors the suggestion that the presence of fluids (presumably released by the blueschisteclogite phase change and at least partially transported back up the slab) is more important than dehydration embrittlement itself.We also cannot directly demonstrate that a negative Clapeyron slope is necessary for seismicity, as suggested by e.g., Wong et al. (1997) and Hacker (1997).The blueschist-out boundary runs subparallel to the temperature contours and could, for the purposes of delimiting the upperplane seismicity, be approximated with a constant T ∼ 400-450 • C boundary, with seismicity also occurring at significantly lower temperatures.
We ignored shear heating in the seismogenic zone.While with this assumption our models generally satisfy the heat flow data shown in Fig. 1, it has been suggested that minor shear heating is necessary in this arc (Peacock and Wang, 1999) and the fore-arc thermal structure predicted here might  4b and (c) close-up on earthquake distribution and predicted blueschist-out curves as in Fig. 6.Here, (1) is the prediction for the general model parameters (Fig. 6); ( 2) is for the modified crustal structure; and (3) is like (2) but with deepening of the decoupled zone to 100 km depth.be slightly on the cold side.The thermal structure of slab below 80 km is principally controlled by the speed of the slab and wedge convection, and our assumption of ignoring shear heating does not affect the conclusions of this paper.
The thermal models provide a less convincing explanation for the seismicity (and heat flow) in the junction between the northern Japan and Kurile arcs.We tested the hypothesis that a subducted sliver of the Kurile fore-arc crust could provide sufficient shielding (as in Kita et al., 2010a) but failed to find a significantly improved model.We conclude that our 2-D models, while reasonable for most of the subduction zone geometry considered here, fail to adequately predict the thermal structure at this junction.It is likely that 3-D and time-dependent effects that are not taken into account in the present set of models have a strong influence on the local temperature conditions in this slab.The junction between the Kurile and Japan arcs forms a transition from a moderately steep subduction zone in the north to much less steeply dipping subduction in the south.This has been ascribed to faster trench retreat in the south (and potentially linked to slab stagnation below the Japan Sea and further west; see Fukao et al., 2009, for a review).Slab stagnation appears absent to the north suggesting that the junction is an actively deforming and rapidly modifying corner in this subduction zone.Morishige and Honda (2012) provide a fully dynamical model for this region that takes trench retreat and slab stagnation into account.While they focused on processes somewhat deeper than in the present study, they demonstrated significant along-arc deformation.Their results are in broad strokes consistent with the projected southeastward translation of the Kurile fore-arc.Combined these effects lead to significant out-of-plane advection of material, specifically at the Hokkaido-Tohoku junctions, which is not represented in our 2-D models.

Appendix A
The change in slab shape in parts of this region causes significant dip of the slab in the direction perpendicular to the cross section, and application of a standard projection algorithm of the 3-D data onto the 2-D cross section leads to significant errors, particularly at the junction between the Japan and Kurile arcs.We therefore adopt the following procedures to map the 3-D slab surface and earthquake data onto the plane of the cross section.The cross section-specific coordinates are depth z, horizontal distance from the trench p, and distance from the cross section q (Fig. A1a).To find the intersection of the slab surface we use the GMT project tool (gmt.soest.hawaii.edu) to find the points in (p, q, z) space that are to within a kilometer from the cross section (red and grey symbols).We then use linear interpolation between nearest neighbors to q = 0.The coordinates of the points on the slab surface that are used in the mesh are determined by linear interpolation in (p, z) from the piecewise linear curve spanned by the yellow symbols.This approach makes sure we have the same smoothness to the slab surface in 2-D as in 3-D.
We interpolate earthquakes from a distance of 10 km on each side of the main cross section (Fig. A1b).In this case the position of the earthquakes with respect to the top of the slab is most important.To make this accurate, we set up a series of auxiliary cross sections at a spacing of 2 km from the main cross section.We determine the local slab surface in (p, z) as

PFig. 1 .
Fig. 1.Location of cross sections used in the thermal modeling of the Tohoku-Hokkaido subduction system.The thin red lines show the depth contours (in km) to the Pacific (PA) slab (fromZhao et al., 1997;Kita et al., 2010b).Heat flow(Tanaka et al., 2004) is shown with the colored circles.Volcanoes are shown in black triangles.The down-dip limit of low-angle thrust events fromKita et al. (2010a), which is taken to be the down-dip limit of the region of strong seismic coupling between upper and lower plate, is shown in the bold red line.Cross sections are labeled followingKita (2009) with P2-P23 over Hokkaido and T2-T25 over Tohoku.Arrow shows Pacific (PA) to North America (NA) convergence.

Fig. 2 .
Fig.2.Seismicity in the Pacific slab below Tohoku.Modified afterKita et al. (2006) andAbers et al. (2012).The earthquakes from 120 km wide box shown in the insert are projected on the center cross section, normal to the trench.This region is approximately bordered by profiles T9 and T18 in Fig.1.The upper plane of seismicity is within 10 km (in depth) from the top of the slab and tends to be deeper in the slab as the slab descends.The blueschist-out boundary (predicted for this profile following the methods by vanKeken et al., 2011) seems to delimit the deepest seismicity in the upper plane.

Fig. 4 .
Fig. 4. Thermal model for cross section T18 over northern Tohoku.(Top:) Predicted heat flow (solid line) compared to observations (from Tanaka et al., 2004; plotted here due to the sparsity of the data from 100 km distance from the profile).(Bottom:) Temperature distribution.Temperature contours are in • C. The top of the slab in this model is shown in green.Position of the volcanic arc is shown by the red triangle.The Moho is indicated by M. The seismicity within 10 km on each side from this profile is shown in circles colored by distance to the slab surface to approximately indicate the upper plane (depths from 0-10 km in red), and the lower plane of the double Wadati-Benioff zone (> 25 km depth in green).Scattered seismicity at a depth of 10-25 km is shown in grey.

Fig. 5 .
Fig. 5. Seismicity within 10 km from cross section (orange circles) compared to the location of the blueschist-eclogite transition (blue line) from the thermal models of Tohoku (T2-T25) and eastern Hokkaido (P2-P9).The top of the slab is shown in the solid green line.Volcanoes (within 30 km from cross section) are shown in red triangles.The furthest extent from coast line to coast line is indicated by the solid black bar at the top of each cross section.Dashed blue line in P2 and P5 is speculative extension of the blueschist-out boundary to below Moho depths.

Fig. 6 .
Fig.6.Same as Fig.5, but now for the cross sections in western Hokkaido above the juncture between the Tohoku and Hokkaido subduction systems.

Fig. 7 .
Fig. 7. Modified thermal model for cross section P15 of western Hokkaido, using modified crustal structure suggested by Kita et al. (2010a).(a) heat flow as in Fig. 4a; (b) temperature distribution as in Fig.4b and (c) close-up on earthquake distribution and predicted blueschist-out curves as in Fig.6.Here, (1) is the prediction for the general model parameters (Fig.6); (2) is for the modified crustal structure; and (3) is like (2) but with deepening of the decoupled zone to 100 km depth.