Dominant role of temperature on deep earthquake mechanics for the Tonga slab near the bottom of the upper mantle

Abstract

The effects of temperature on the mechanics of deep earthquakes are investigated based on detailed seismic images of a horizontal portion of the Tonga slab (19°S to 22°S). The hypocenter distribution of deep earthquakes and tomographic models have shown that the Tonga slab in this region stagnates laterally around the upper- and lower-mantle boundary. We analyze data from seismic networks in the United States and Japan from the deepest earthquakes that occurred in the stagnating part of the slab (focal depths 658 to 678 km). We observe a clear arrival approximately 8 s after the direct P wave and find that it is an S-to-P converted wave emitted upward from the deepest earthquake foci and reflected downward by a horizontal interface located nearly 30 km above the foci. The S wave speed drops upward across the interface by a few percent, so the interface is best interpreted as the former oceanic Moho in the stagnating Tonga slab. A numerical model of the thermal structure of the Tonga slab is constructed based upon the slab geometry, with the dynamics of the slab motion being incorporated. We find that if the deepest earthquakes occur in the coldest part of the slab, the vertical distance between the foci and the Moho should be approximately 30 km, which agrees well with the observations. The temperature around the foci is determined to be 750°C. We suggest that deep earthquakes are stringently restricted below this temperature, near the bottom of the upper mantle.

1Findings

1.1 Introduction

The overall features of slab subduction are now well known, owing to advances in seismic tomography and dynamic modeling. However, the detailed mechanism of slab descent across the upper- and lower-mantle boundary still needs to be constrained, as important geophysical issues such as deep earthquake mechanics and the fate of oceanic crust hinge upon it.

The thermal control on deep earthquake mechanics has been suggested by previous researchers (e.g., Gorbatov and Kostoglodov [1997]; Wiens [2001]). Wiens ([2001]) shows that the colder a slab is, the more deep earthquake aftershocks tend to occur in the slab. Gorbatov and Kostoglodov ([1997]) show positive correlation between the maximum depths of deep seismicity and the thermal parameters of slabs (the product of the age and downward velocity of the slab, called ‘ Φ’ hereafter). We note that it is by no means straightforward to reveal the accurate thermal structures of deep slabs around deep earthquake foci and that previous studies relying mostly on global seismicity data sets and utilizing Φ of deep slabs have a weakness in this respect. A more recent study by Emmerson and McKenzie ([2007]) re-examines the thermal structure of individual slabs and compares the temperature structure with the deep seismicity for each slab. These authors then estimate the limiting temperature, below which deep earthquakes can occur. However, it is implicitly assumed that the hypocenters of the deepest earthquakes occur in the coldest part of the slabs, which needs to be verified. It, thus, remains a challenge to accurately constrain the locations of hypocenters of the deepest earthquakes inside the slabs for a better understanding of the thermal control on deep earthquake mechanics. An advantage of this study over previous global studies is that accurate observational constraints are available for the hypocenters within the Tonga slab. A recent study by Kubo et al. ([2009]) has determined the temperature near the deep earthquake foci inside the Mariana slab, around 600 km depth, to be about 550°C based upon the observation of a metastable olivine wedge (Kaneshima et al. [2007]). However, experimental studies on the condition below which metastable olivine persists still carry substantial uncertainty. In this study, we investigate the role of thermal structure with unprecedented resolution on the basis of the geometry of the Tonga slab and the hypocenter distribution in it. We shall focus on the problem of whether deep earthquakes can occur under certain temperature and stress conditions and shall not discuss specific processes in which slab rocks rupture during deep earthquakes. The results show a very stringent thermal constraint on deep earthquake mechanics near the bottom of the upper mantle.

Our arguments rely critically on seismological detection of former oceanic crust. It has been debated whether the crust is scraped off near the upper- and lower-mantle boundary and remains in the upper mantle or if it continues subducting into the lower mantle, being attached to the slab (e.g., Karato [1997]), while studies of scattered seismic waves have confirmed the ubiquitous presence of compositional heterogeneities that are most likely to represent rocks with basaltic composition in the mid-lower mantle (e.g., Kaneshima and Helffrich [1999], [2010]). As for the shallower upper mantle (depths of 50 to 350 km), a thin oceanic crust with lower seismic speeds has been detected by ScSp analyses (e.g., Nakanishi [1980]; Helffrich et al. [1989]). Intact subduction of oceanic crust deeper in the upper mantle has not been confirmed by seismic observations but is often assumed in many geophysical and geochemical studies. This is partly because there is no obvious mechanism in the upper mantle to strip the crust from the slab. Thus, it has been an important challenge for seismologists to detect former oceanic crust attached to the slab, if any, deeper in the mantle. In this study, we first present seismological evidence that shows the attachment of former oceanic crust to the underlying peridotitic slab at the bottom of the upper mantle. This observation suggests that the boundary (which shall be called the former oceanic Moho or simply Moho) remains sharp during the descent through the upper mantle.

1.2 Event locations and regional deep seismicity

We focus on a small area in the Tonga subduction zone at latitudes from 19°S to 22°S, where detailed structure of the slab has been imaged by the accumulation of long-term array data. The data we used are short-period seismograms from the University of Washington network (UW), J-array (JA), and Hi-net (HI) in Japan for very deep earthquakes (658 km ≤h≤ 678 km) that occurred at the Tonga subduction zone (Figure 1a) from 1995 to 2009 (Table 1). We relocate the earthquakes using arrival times of P, PKP, and S waves reported by the ISC using the method described in Kaneshima et al. ([2012]). In order to constrain the focal depths, we use pP-P times measured for the UW array by slant stacking the seismograms with the slowness predicted for the pP wave.

Figure 1

Event locations and geometry of the slab. (a) Map of the Tonga region showing the relocated epicenters of the events (circles, Table 1), with focal depths shown with a color scale. Thick dotted line is approximately parallel to the subduction direction, along which the vertical cross section of hypocenters and reflection points is shown. Contours of the depths of the Tonga slab are shown every 50 km from 250 to 700 km with thin dotted line (Gudmundsson and Sambridge [1998]). Plotted here and in the cross section (b) are essentially all the events large enough to be relocated accurately by incorporating array measurements of pP. (Inset) PDE epicenters of the deep and intermediate-depth earthquakes from 1970 to 2013 (deeper than 600 km, all magnitudes) are shown with red dots in a map of the larger area of the Tonga subduction zone. The study area is shown in the box in the middle. Gray irregular shapes are the islands of Fiji. (b) Vertical cross section of the relocated hypocenters and conversion (reflection) points along the thick dotted line in (a). The horizontal axis indicates the distance along the dotted line from the origin taken at 20.5°, 181.5°E, WNW to ESE from left to right. Red circles show the hypocenters of the deep events, while blue circles show the locations of the reflection points corresponding to the X phase (interpreted as the Moho). (Inset) Schematic figure of the Tonga slab showing hypocenters (red circles), former Moho (blue line), crust, and ray path of s-to-P waves (dotted line for S, and solid line with arrowhead for P). Dashed line shows the boundary between the former oceanic crust and the surrounding mantle.

Table 1 Event locations

The relocated hypocenter distribution of the events, combined with those of shallower events (Figure 1b), shows that the Tonga slab descends from ESE to WNW. The events near the northwestern edge of the study area with focal depths 668 ± 10 km are hereafter called either ‘deepest earthquakes’ or ‘deepest events’ in this study. There are 14 deepest events (m _b > 5) since 1995, whose hypocenters are well constrained (Table 1). The cluster of deepest events extends nearly laterally for about 100 km (Figure 1b). The frequent occurrence of deepest earthquakes sharply terminates beyond the northwestern edge of the study area (Figure 1a, inset). If the horizontally lying slab really exists further to the northwest, as suggested by tomography models, the rapid truncation of seismicity should reflect a change in the thermal and/or material conditions of the slab. It is believed that deep earthquakes do not occur after the post-spinel transformation because of significant grain size reduction due to the phase transformation. In the present case, the transition from seismic to aseismic conditions inside the slab does occur not vertically as usual, but laterally, and the post-spinel transformation front is probably located somewhat below the foci (Kaneshima et al. [2012]). Therefore, the transition is most likely to reflect thermal effects.

1.3 Waveform data and detection of former oceanic Moho by array analyses

We use array techniques to detect later arrivals in the P coda, which are converted and/or reflected at an approximately horizontal interface around the focal depth. The data processing methods are the same as in Kaneshima et al. ([2012]). We find anomalous arrivals that are not predicted by the standard seismological Earth models such as ak135 (Kennett et al. [1995]), approximately 8 s after the direct P waves (hereafter called ‘X phase’). The measured slowness of these arrivals relative to P are mostly positive (Δ p ∼ +0.02 s/°) (Figure 2, Table 2). This means that the X phase is not an S-to-P converted wave (emitted downward from the foci as an S wave) generated at the ‘660 km discontinuity’ that is depressed below the foci down to about 680-km depth but is an s-to-P wave that is emitted upward as an S wave and reflected at an interface above the foci at a depth of about 640 km (Table 2). This inference is supported by a positive correlation between the focal depth and the arrival times of the X phase relative to P (Δ t) among the deepest events observed at Hi-net (Table 2, Figure 3). The reflection points are determined using the measured Δ t and the hypocentral locations, based on the assumption that the s-to-P reflection occurs at a horizontal interface. The averaged vertical distance of the conversion point from the foci are about 30 km when ak135 is used (Figure 1b, Table 2), and the correction for a 4% higher speed along the path inside the slab to the reflection points (e.g., Deal et al. [1999]) decreases the interface depth by less than 2 km.

Figure 2

Examples of slant stack vespagrams for selected events. (Top) EV9 (Table 1) recorded at the Hi-net. Seismograms are phase-weight stacked with n=1 (Schimmel and Paulssen [1997]), and the color shows the stacked intensity (the strongest signals are red). The horizontal axis is the delay time (Δ t) in seconds after the first peak of the direct P wave. The vertical axis is the slowness relative to the direct P wave (Δ p) in seconds per degree. Arrival of the X phase (s-to-P reflected wave at horizontal interface above the foci) is seen near 8 s after P (s641P, Table 1). The predictions of s-to-P based on ak135 are marked with blue dots for every 5 km of interface depth. The stack outputs between 3 to 10.5 s after P are magnified 20 times, while those after 10.5 s are magnified 40 times. (Bottom) Same as the top panel, for EV8, showing the arrival of s639P (Table 1). The outputs after 4.5 s after P are magnified 30 times.

Figure 3

Examples of stacked waveforms for the direct P waves and s-to-P waves at Moho observed at Hi-net. Stacked waveforms for the direct P waves (left traces) and the s-to-P waves at Moho (‘X phase’, right traces) observed at Hi-net. Stacking is performed over the entire Hi-net with corresponding slowness values. The seismograms for the X phases are linearly stacked (dashed lines) or phase-weighted stacked (solid lines). Each trace is normalized to its maximum peak-to-peak amplitude. The stacked waveforms of the X phase are magnified relative to the stacked P wave, and the magnification factor for each trace is shown below the event name. When phase-weighted stacked, the amplitudes of weaker anomalous arrivals are reduced compared to those of the linearly stacked waveforms. Note that the polarities of the X phase of EV9 and EV10 are flipped relative to the P waves that are initially positively polarized, while the X phase polarity of EV8 is the same as the P wave that is negatively polarized.

Table 2 Measurements for X and Y phases

For events EV8, EV9, and EV10 (Table 1), the stacked waveforms of the arrivals at Hi-net have high S/N ratios, and their waveforms are quite similar to those of direct P arrivals. However, the polarities relative to P are different between EV8 and the other two events (Figure 3). According to their focal mechanisms and the P wave polarities observed at Hi-net, if the S wave speed (Vs) decreases upward across the interface, an s-to-P wave reflected at a horizontal interface immediately above the foci should have the same (EV8) or flipped (EV9 and EV10) polarity compared with the direct P. This is consistent with what is observed for the three events (Figure 3) and thus strengthens the interpretation that the arrivals are s-to-P reflected waves. Based on the observed amplitudes of the arrivals relative to direct P (Table 2), we estimate the jump of seismic speeds across the interface. The data set is most sensitive to the Vs contrast, so we assume that the Vp anomaly scales with the Vs anomaly. We find that the Vs drops upward from the inside to the outside of the horizontally lying slab by a few percent (approximately 3%) when the interface is a first-order discontinuity, although the estimate carries substantial uncertainty due to the limited accuracy of the focal mechanisms. The observations of the interface location and the sense of anomaly indicate that the interface represents the Moho of former oceanic crust attached to the subducting Tonga slab.

There is also a weaker later arrival nearly 4 s after the s-to-P reflection at the Moho (called ‘Y phase’), which also tends to show positive slowness relative to P (Table 2; Figure 2). We interpret this also to be s-to-P reflections from another interface above the Moho, which would represent the surface of the former oceanic crust. Although the polarity of this weaker arrival has not been constrained, the vertical distance between the two interfaces (14 ± 2 km) is determined using Δ t (Table 2).

1.4 Thermal structure of the Tonga slab

The Tonga slab is one of the coldest slabs, with a Φ of about 9,000 km around the study area. In order to determine the detailed thermal structure of the Tonga slab, we calculate the stream function and temperature fields associated with the subduction of the Pacific Plate beneath the Australian Plate by numerically simulating temperature and mantle flow in a 2-D box model (Yoshioka and Sanshadokoro [2002]) (Figure 4a). As the initial condition concerning mantle flow in the model, no flow condition is given, and mantle flow is entirely driven by the subduction of the Pacific Plate. As the initial condition for the temperature distribution in the model space, we use a layered thermal state of a cooling oceanic plate described by an error function. Heating due to adiabatic compression is taken into account in computing temperatures for depths below 150 km, as is commonly performed (e.g., Torii and Yoshioka [2007]). As the boundary condition for temperature, we assume that materials with this temperature profile for the oceanic Pacific Plate flow into the model domain from the side boundary (Figure 4a). Temperature at the model surface is fixed to 0°C, and an adiabatic condition is given for the left and bottom boundaries (Figure 4a). Friction at the interface between the slab and the overriding continental plate is not included, but its effect on the slab temperature near the bottom of the upper mantle is negligible. Latent heat absorption associated with the phase transformations is not incorporated into the model. Latent heat absorption of post-spinel transformation reduces the temperature only slightly (Yoshioka et al. [1997]). On the other hand, when the transformation from metastable olivine to either wadsleyite or ringwoodite occurs, latent heat release causes a temperature rise. The temperature increase is primarily dependent on where the phase transformation occurs. It would exceed 200°C when the transformation occurs immediately upstream of the foci (Kirby et al. [1996]; Kubo et al. [2009]) or would be less than 50°C for the transformation at a 450-km depth, where the coldest temperature of the Tonga slab is approximately 600°C. According to recent experimental studies, the latter case seems more plausible (e.g., Perrillat et al. [2013]). As there is no observational constraint on the metastable olivine wedge for the study region, further quantitative arguments about this effect seem inadequate.

Figure 4

Models of the temperature and the thermal shear instability of the Tonga slab. (a) Schematic figure of the simulation of the subduction of the Tonga slab. The symbol σ _xx indicates the normal stress acting on the side boundaries. (b) Contour plot for the temperature model of the Tonga slab in the study area. Contour interval is 50°C. Model parameters used in the temperature calculation are as follows: C_p=1.046×10³ (J·kg⁻¹ K⁻¹), k=3.138 (J·m⁻¹ K⁻¹s⁻¹), and α=2×10⁻⁵ (K⁻¹), where C_p is specific heat at constant pressure, k is thermal conductivity, and α is thermal expansivity. The instability parameter is evaluated along the vertical line annotated with ‘A’ and ‘B’. The red circle shows the coldest point, therefore the representative deepest foci, along A-B. Blue and green circles are 30 and 44 km above the coldest point, respectively. (c) Properties associated with instability criterion across the stagnating part of the slab along A-B in (a). Horizontal axis is the distance across the slab from A to B. Top, middle, and bottom images show temperature (°C), stress (MPa), and logarithm of Θ (log10Θ), respectively. Red, blue, and green circles in the top image correspond to those in (a). Solid lines are for diffusion creep for the case of 1- μ m grain sizes when diffusion creep dominates over dislocation. Dashed lines are for dislocation creep for the grain size of 5 μ m when dislocation creep dominates over diffusion creep. The ratio between the dislocation glide length and climb length (l_g/l_c) is 200 (Shimojuku et al. [2009]). Other rheological parameters of ringwoodite are taken from Shimojuku et al. ([2009]).

For the region outside the slab, we solve the momentum equation expressed by the stream function and the energy equation simultaneously as a coupled problem using a 2-D finite difference method. The slab materials inside a prescribed guide are forced to descend kinematically, with a uniform velocity of 10 cm/year at the age about 100 Myr. The trench is kinematically rolled back by 5 cm/year (e.g., Heuret and Lallemand [2005]; van der Hilst [1995]), using the method described in Torii and Yoshioka ([2007]). The Clapeyron slope of the post-spinel transformation is −2.0 MPa/K (Kaneshima et al. [2012]). We assume that mantle viscosity depends on temperature and depth for the calculation of the thermal structure (Christensen [1996]; Kaneshima et al. [2012]). The viscosity of the shallowest 300 km of the lower mantle is reduced to a quarter of that of the transition zone. A guide descending 60° down to 500-km depth is prescribed in order to model the flow which resembles the Tonga slab. Other details of the simulation are the same as those in Kaneshima et al. ([2012]).

We seek a transient temperature structure that resembles the overall shape of the slab imaged by tomography studies (Obayashi et al. [2006]; van der Hilst [1995]) and the deep seismicity in the study area (Figure 1b).

Figure 4b shows the calculated temperature distribution at a time 20 Myr after the initiation of the subduction of the Pacific Plate, which adequately matches the shape of the Tonga slab. The model parameters used for obtaining the model are listed in the caption of Figure 4. It is possible to find different sets of dynamical parameters that fit the observations equally well, but our aim in performing the simulation is to present a model of a temperature field that is physically and tectonically feasible.

1.5 Location of the former oceanic Moho relative to the deep earthquakes

Deep earthquakes are often considered to occur around the coldest part of subducting slabs (Gorbatov and Kostoglodov [1997]; Karato et al. [2001]). The coldest slab temperature at 660-km depth is approximately 750°C for the thermal parameter Φ of about 7,000 km (Karato et al. [2001]), approximately 680°C for Φ approximately 14,000 km (Kubo et al. [2009]), and approximately 700°C for Φ approximately 8,600 km (Bina and Kawakatsu [2010]). Emmerson and McKenzie ([2007]) estimate the coldest temperature of the Tonga slab to be 681°C to 710°C using temperature-dependent thermal conductivity models. The temperature of the coldest portion at depths of 660 to 680 km in our model for the Tonga slab is slightly above 750°C (Figure 4) and is generally consistent with previous estimates.

Next, we evaluate the permissible range of temperature estimates. Thermal conductivity (k) is obviously an influential parameter and is fixed at 3.138 J·m⁻¹ K⁻¹ s⁻¹ in our simulations (Marton et al. [2005]). As k may decrease with temperature (Emmerson and McKenzie [2007]), we investigate a range of k from 2.9 to 3.2 J·m⁻¹ K⁻¹ s⁻¹. We also vary the age of the Pacific Plate within a range from 90 to 120 Myr. The results show that the temperature at the coldest part of the Tonga slab could differ from the current estimate by about 20°C. The value of k is constant in our simulations. If its temperature dependency is considered, the difference could be larger (Emmerson and McKenzie [2007]).

If we assume that the deepest earthquakes in Tonga occur near the coldest portion of the slab, the location of the seismically detected Moho relative to the events (30 km above the foci) indicates that the temperature of the base of the former crust is nearly 1,050°C, according to the calculated thermal model (Figure 4b). Thus, the temperature 14 km below the slab surface is approximately 400°C colder than the surrounding mantle (Figure 4b), which generally agrees with previous studies (e.g., Karato et al. [2001]; Bina and Kawakatsu [2010]). We thus claim that the location of the Moho is quite reasonable. Combined with the consistency among previous tomography images, the deep seismicity (Figure 1b), the dynamically plausible temperature structure (Figure 4b), and the relative location of the deepest earthquakes and the Moho provide strong support for the model that deep earthquakes occur only below a certain temperature that is near 750°C.

The estimated maximum temperature for deep seismicity falls within the range obtained by Gorbatov and Kostoglodov ([1997]) (approximately 600°C to 800°C) but is higher than for the Marianas (550°C; Kubo et al. [2009]) and the higher estimate for Tonga (710°C; Emmerson and McKenzie [2007]). The large difference between the temperatures for the Marianas and Tonga may be caused by the fact that the temperature for the Marianas indicates the temperature at which metastable olivine disappears and does not strictly represent the limiting temperature of deep seismicity.

The upper bound temperature we estimated may carry an uncertainty on the order of 50°C, since it depends on the tectonic parameters of the slab and the properties of the slab materials as described above. However, the narrow range of focal depths of the deepest earthquakes (approximately 20 km) indicates that the temperature range that encompasses the seismogenic zone is as narrow as 50°C. Therefore, the presence of a stringent upper limit of temperature for the deepest seismogenesis is robust.

1.6 Discussion

The mechanism of deep earthquakes is still controversial. The postulated mechanisms include transformational faulting of olivine (e.g., Kirby et al. [1996]), adiabatic shear (plastic) instability (e.g., Karato et al. [2001]), and embrittlement due to dehydration of hydrous minerals (Meade and Jeanloz [1991]). As for dehydration embrittlement, the deepest earthquakes at the Marianas occur mainly inside the metastable olivine wedge that is inferred to be very dry (Kubo et al. [2009]). Green et al. ([2010]) also claim that subducting slabs are almost dry below 400 km. We have shown in the last section that the deepest seismicity is strictly controlled by temperature. There is no dehydration reaction that is sensitive to temperature under the low-temperature conditions relevant to a cold slab like the Tonga slab, according to the phase diagram of hydrated minerals associated with the slab (Ohtani [2005]; Komabayashi et al. [2004]). Dehydration embrittlement can thus be eliminated.

For both the transformational faulting and shear instability mechanisms, faulting is most likely to occur near the coldest places inside slabs. Therefore, the thermal control we have shown seems compatible with these two mechanisms.

Transformational faulting requires the presence of metastable olivine around the deepest foci, which has not been proved for the study region. It also means that the observed seismicity reflects the disappearance of metastable olivine at a distance of about 1,000 km from the trench along the slab. According to recent experimental studies (e.g., Perrillat et al. [2013]; Du Frane et al. [2013]), the estimated temperature at the foci (approximately 750°C) seems too high for a substantial volume of metastable olivine to persist down to this distance from the trench, unless the slab is extremely dry (the H ₂O content is much less than 70 ppm).

In the absence of metastable olivine, the rheology of the slab near the foci is controlled by the rheology of ringwoodite, for which the expected deformation style would be either diffusion creep or dislocation creep, depending primarily on the grain size (e.g., Shimojuku et al. [2009]). We consider whether thermal shear instability readily occurs under the temperature, pressure, and strain rate of the slab model near the deepest foci (Figure 4b). A catastrophic plastic slip occurs when the strain softening effect due to temperature rise following the deformation exceeds the strain hardening effect. The more the stress is exerted and the cooler the temperature, the more likely the instability is to occur. The criterion for instability is expressed as $\Theta \equiv -{\left(\frac{\mathrm{\partial \sigma }}{\mathrm{\partial T}}\right)}_{\epsilon ,\dot{\epsilon }}\frac{\kappa }{\mathrm{aC\rho }}>1$, where σ, C, ρ, κ, and a, are differential stress, specific heat, density, fraction of heat retained in a slipped region during the slip, and strain hardening factor, respectively (Hobbs and Ord [1988]; Karato et al. [2001]). We evaluate ${\left(\frac{\mathrm{\partial \sigma }}{\mathrm{\partial T}}\right)}_{\epsilon ,\dot{\epsilon }}$, and therefore Θ, by assuming quasi-steady state plastic flow laws for our slab model (Figure 4c; see the caption for the parameter values). The strain rate is assumed to be 5 ×10⁻¹⁵ s⁻¹ uniformly, but its effects are small. The strain hardening constant for the creep, $a\equiv \frac{1}{\sigma }{\left(\frac{\mathrm{\partial \sigma }}{\mathrm{\partial \epsilon }}\right)}_{\dot{\epsilon },T}$, is set to 0.5 following Hobbs and Ord ([1988]). The heat fraction κ is 1. According to the profiles of Θ across the slab along the line A-B in Figure 4b (bottom image of Figure 4c), for diffusion creep that is dominant when the grain size is 1 μ m, instability can only occur approximately 15 km around the coldest part of the slab (solid lines in Figure 4c). This is consistent with the observations (Figure 1b). When the grain size is larger, say, ≥5 μ m, dislocation creep dominates, but instability is still possible approximately 25 km around the coldest temperature (dashed lines in Figure 4c). Thus, thermal shear instability adequately explains the essential features of the the observation with the currently most acceptable rheological parameters.

In this study, we consider transformational faulting and thermal shear instability separately. If metastable olivine persists to the foci and the two mechanisms are concomitant, latent heat release due to the transformation increases the temperature, weakens the surrounding rocks, and enhances the instability. The stress regime around the deepest foci of this study, estimated by the focal mechanisms, is lateral compression parallel to WNW-ESE, nearly parallel to the direction of the maximum dip of the slab. The stress regime obtained by our simulation is consistent with the observation.