Relation between rheological properties and the stress state in subducting slabs

Abstract

The distribution of different stress states in the subducting slab indicated by centroid moment tensor solutions for intra-slab earthquakes can help constrain the rheological properties of the slabs. A comparison of slabs in the western and eastern Pacific realms shows contrasting patterns in the stress states down to depths of ~ 350 km. The majority of slabs in the western Pacific show a pair of down-dip compression (DC) and down-dip tension (DT) domains in the upper and lower parts of the slab reflecting the effects of the slab unbending during progressive subduction. In contrast, slabs in the eastern Pacific show predominantly in-plane DT stress irrespective of slab geometry. Two-dimensional numerical simulations assuming constant slab thickness and viscosity indicate that the development of slabs with in-plane DT stress at depths of 100–300 km requires the slabs to be thin and have a low viscosity (10²³ Pa s). Weak slabs bend easily and tend to fold when they encounter increased resistance to downward movement at the 660-km boundary. The associated DC stresses are not transmitted up the slab so negative buoyancy of the slab and DT stress dominates at intermediate depths for this type of slab. Most experimentally derived rheological parameters predict a high viscosity (> 10²⁴ Pa s) for such slabs. However, two-dimensional numerical simulations using temperature- and pressure-dependent viscosity show that a relatively low activation energy (~ 110 kJ/mol) for diffusion creep is a possible explanation for the observed distribution of stresses in the slabs. Such low activation energies are compatible with recent experimental work on diffusion creep of polyphase mantle materials in which a low effective activation energy for creep results from a slow grain growth due to pinning effect of the secondary phase. The simulations provide a mechanical explanation for the observed dominantly DT stress state at 100–300 km depths for young slabs and paired DT and DC stress states at the same depth range for old slabs.

Graphical Abstract

Introduction

The stress state in subducting slabs estimated from centroid moment tensor (CMT) solutions for intermediate-depth (< 350 km) earthquakes in the Pacific realm can be broadly divided into predominantly down-dip tension (DT) for east-directed subduction zones and predominantly down-dip compression (DC) for west-directed subduction zones. In contrast, most deep (> 350 km) earthquakes indicate DC (Alpert et al. 2010; Billen 2020; Chen et al. 2004; Isacks and Molnar 1969, 1971). The discovery of the double seismic zone and associated focal mechanisms indicating DC for the upper layer and DT for the lower layer (Hasegawa et al. 1978) suggests that the stress states of intermediate-depth slabs do not correspond to simple uniform in-plane DT or DC stresses. A detailed analysis of the global distribution of the stress states of intermediate-depth slabs requires a higher resolution of observations than those available for many areas with active subduction. For example, Alpert et al. (2010) analyzed average stress orientations for 50-km depth intervals and thus was unable to distinguish the stress orientations associated with double seismic zones. In addition, their analysis focused on earthquakes at z > 100 km in order to exclude earthquakes at the inter-plate megathrusts and consequently their analysis does not include intra-slab earthquakes at z < 100 km depths.

The intra-slab stress state is considered to be controlled by interaction between (i) negative buoyancy caused by the presence of a cold slab and/or the effect of an exothermic phase transition at the top of the transition zone; (ii) forces that act to resist slab penetration into the lower mantle as a result of increasing viscosity and/or endothermic phase transition at the bottom of the transition zone; and (iii) flexural stress associated with slab bending and unbending (Alpert et al. 2010; Chen et al. 2004; Hasegawa et al. 1978; Isacks and Molnar 1969, 1971; Ito and Sato 1992).

Sandiford et al. (2020) evaluate the relationship between slab geometry and stress state estimated from intra-slab earthquakes at depths < 300 km for six subduction segments (Fig. 1). Three western Pacific subduction segments (Tonga, Kuriles, Japan) show slab geometries that can be termed ‘normal subduction’ with one convex-up slab hinge (Fig. 1a–c). The three eastern Pacific subduction segments show different characteristics: the segments of Chile and Central America show slab geometries termed ‘low angle subduction’ with two convex-upward slab hinges, and the Peru segment shows a slab geometry termed ‘flat slab subduction’ with two convex-upward slab hinges and one concave upward slab hinge in between (Fig. 1d–f).

Fig. 1

Trench-perpendicular projection of CMT T-axes on vertical cross sections (modified from Sandiford et al. 2020). a Tonga; b Kuriles; c Japan; d Chile; e Peru; f Central America. The origin of the horizontal axis represents the location of the trench. The T-axes are plotted as projections of uniform-length vectors, with the length variation in the plotted T-axes reflecting the magnitude of the projected component of the vector. T-axes are colored according to the orientation relative to the slab midplane (shown with dashed black line): red represents earthquakes with a down-dip compressive (DC) sense, blue with a down-dip tensional (DT) sense

The stress state of slabs for the three western Pacific subduction segments is characterized by two pairs of DC and DT domains distributed on the upper and lower parts of the slab across the neutral plane: one pair occurs in the outer rise–trench region with DT in the upper part and DC in the lower part and the other occurs in an intermediate depth (50–300 km) region with DC in upper part and DT in the lower part. These stress states of the slab can be explained by bending (increase in the slab curvature) in the outer rise–trench region followed by unbending (decrease in the slab curvature) when the slab reaches deeper parts of the subduction zone (e.g., Sandiford et al 2020; Fig. 2b). However, the DC domain in the lower part of the slab is not well developed in the outer rise–trench region for either the Japan or Kuriles segments (Fig. 1a–c). The Chile and Peru segments also shows stress states associated with the bending–unbending processes near the trench. However, in more deeply subducted parts of the slab (> 100 km) DT dominates irrespective of the slab geometry (Fig. 1d, e). The plot for the Central America segment shows a mixed distribution of earthquakes indicating both DT and DC in the 50–100 km depth range. This apparent mixing is probably the result of projecting events onto one vertical plane ignoring along-strike lateral changes in slab geometry that likely result in diverse stress environments being shown in the same region. DT dominates in the more deeply subducted parts of the slab (Fig. 1f).

Fig. 2

Two processes affecting the stress state in the subducting slab: uniform in-plane stretching–shortening (a) and bending–unbending (b)

These features of intra-slab stress state suggest that the relative contributions of three forces due to negative buoyancy of the dense slab, resistance at the 660-km boundary and flexural bending of the slab varies among different subduction segments. The negative buoyancy causes an in-plane DT stress state that operates across the full width of subducting slabs and the resistance at 660-km boundary causes a DC slab stress state in the region immediately above the 660-km boundary. The balance between these two forces controls the depth of the transition from DT to DC stress states in the slab (Fig. 2a). A shallow transition is favored by greater resistance at the 660-km boundary, but is also influenced by the strength of the slab: weak slabs do not effectively transmit the compressive stress up dip along the slab from the 660-km boundary. A shallow transition may also be explained by a relatively small negative buoyancy due to the presence of a large volume of the metastable olivine (Kita et al. 2010). Kita et al. (2010) also suggest that such a change in buoyancy force may account for observed differences in the locations of the neutral planes between paired DT and DC domains associated with unbending processes. The above summary shows the stress state in slabs can be understood in terms of a combination of in-plane DT or DC stresses due to the effects of negative buoyancy and resistance at the 660-km boundary and paired DT and DC stresses due to bending–unbending processes (Fig. 2; Kita et al. 2010; Sandiford et al. 2020; Sippl et al. 2022). The relative importance of slab bending is strongly controlled by the rheological properties of the subducting slab, which determine the stresses associated with bending and hence the ability of the slab to transmit DC stress upwards along its length.

To examine the way in which the rheological properties of the subducting slab affect the distribution of DT and DC stresses, we calculated the stress state in a range of subducting slabs based on 2D numerical modeling. As oceanic lithosphere increases in age it cools, increases in thickness and density and the effective viscosity of any ductilely deforming regions increases. As a first step in our evaluation of the relative importance of such changes in thickness and associated rheological properties of subducting lithosphere, we compare models consisting of layers of constant thickness and constant viscosity (Model CV). Then, based on the results of Model CV, we conducted calculations using a model incorporating a temperature- and pressure-dependent viscosity (Model TPV). Model TPV uses a composite flow law including dislocation creep, diffusion creep and plastic yielding. Here, we focus on the effect of changes in the value of the effective activation energy for the diffusion creep.

Model

Model with constant viscosity (Model CV)

Our analysis uses the code I2VIS (Gerya and Yuen 2003), which is based on finite differences with a marker-in-cell technique and a fully staggered rectangular Eulerian grid. The numerical model domain (7000 km wide × 2900 km deep) includes non-uniform 751 × 318 rectangular grids with a resolution varying from 4 to 20 km. The upper central area (x = 3600–5400 km, z = 0–640 km) has the highest resolution of 4 × 4 km. The model includes 5.075 × 10⁶ markers.

In addition to upper mantle (UM), mantle transition zone (TM) and lower mantle (LM), the model includes a subducting oceanic plate and an overriding continental plate (Fig. 3a). The oceanic plate consists of three layers: one layer of oceanic crust (OC) and two layers of mantle (SP1, SP2). SP1 layer represents the shallow cold part of the oceanic mantle lithosphere and SP2 layer represents the deeper warmer part of oceanic mantle lithosphere. The continental plate consists of five layers: forearc crust (OPFC), crust (OPC), forearc mantle (OPFM), arc–back-arc mantle (OPM1), and cratonic mantle (OPM2). OPFC and OPFM layers represent the forearc cold nose. The initial geometry and physical properties are shown in Fig. 3 and Table 1. The oceanic plate (OC, SP1, SP2) has a positive density contrast of ∆ρ_slabU = 60 kg/m³ with respect to the reference density of 3300 kg/m³. The density of continental crust (OPFC and OPC) is 2900 kg/m³. The density of oceanic crust (OC) increases linearly from 3000 to 3360 kg/m³ for the depth range of 15–65 km (Fig. 3c). To incorporate the effect of the phase transitions at 410 and 660 km depth, the density contrasts (∆ρ₄₁₀ and − ∆ρ₆₆₀) are added to SP1 and SP2 for the depth ranges of 330–410 km and 660–690 km, respectively (Fig. 3c). We use the values of ∆ρ₄₁₀ = 124 kg/m³ and ∆ρ₆₆₀ = 193 kg/m³ to account for the olivine phase transitions (olivine to wadsleyite, and ringwoodite to bridgmanite and ferropericlase) and ~ 60% modal content of olivine in peridotite (Faccenda and Dal Zilio 2017).

Fig. 3

Model with constant thickness and viscosity (Model CV). a Model geometry and constituent materials. b Viscosity of materials. c Density of the oceanic crust and mantle materials as a function of depth. Density of the continental crust (OPC, OPFC) is 2900 kg/m³

Table 1 Geometrical and physical properties for model layers in Model CV

The governing equations for the flow are the conservation equations of mass and momentum. Conservation of mass is approximated by the incompressible continuity equation:

$$\frac{\partial {v}_{x}}{\partial x}+\frac{\partial {v}_{z}}{\partial z}=0,$$

(1)

where v_x and v_z are the horizontal and vertical components of the velocity vector, respectively.

The two-dimensional Stokes equations for creeping flow take the form:

$$\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{xz}}{\partial z}=\frac{\partial P}{\partial x},$$

(2a)

$$\frac{\partial {\sigma }_{zz}}{\partial z}+\frac{\partial {\sigma }_{xz}}{\partial x}=\frac{\partial P}{\partial z}-\rho g,$$

(2b)

where σ_xx, σ_zz, and σ_xz are the components of the deviatoric stress tensor, P is pressure, g is acceleration due to gravity, ρ is density.

The constitutive relationship is given by:

$${\dot{\varepsilon }}_{ij}=\frac{1}{2}\left(\frac{\partial {v}_{i}}{\partial {x}_{j}}+\frac{\partial {v}_{j}}{\partial {x}_{i}}\right),$$

(3a)

$${\sigma }_{ij}={2{\eta }_{\text{eff}}\dot{\varepsilon }}_{ij},$$

(3b)

where ${\dot{\varepsilon }}_{ij}$ is the deviatoric strain rate tensor. The effective viscosity, η_eff, is determined from the viscosity of each layer, η_vis, and the equivalent viscosity for the yield stress, η_yield, and corresponding strain rate, ${\dot{\varepsilon }}_{\text{vis}}$ and ${\dot{\varepsilon }}_{\text{yield}}$:

$${\eta }_{\text{eff}}={\left(\frac{1}{{\eta }_{\text{vis}}}+\frac{1}{{\eta }_{\text{yield}}}\right)}^{-1},$$

(4a)

$${\dot{\varepsilon }}_{\text{total}}={\dot{\varepsilon }}_{\text{vis}}+ {\dot{\varepsilon }}_{\text{yield}}.$$

(4b)

The yield stress and the equivalent viscosity are taken to be

$${\sigma }_{\text{yield}}={\text{min}}\left(c+\phi P,{\sigma }_{\text{max}}\right),$$

(5a)

$${\eta }_{\text{yield}}= \frac{{\sigma }_{\text{yield}}}{2{\dot{\varepsilon }}_{\text{yield}}},$$

(5b)

where c = 30 MPa and ϕ = 0.6 are the cohesion and effective coefficient of friction, respectively, and ${\sigma }_{\text{max}}=200 {\text{MPa}}$. The model uses a weak frictional strength (c = 1 MPa and ϕ = 0.005) for the oceanic crust at z < 120 km to facilitate decoupling between the subducting oceanic plate and the overriding plate.

To allow both the subducting and overriding plates to move freely, the plates taper toward both ends of the model. Material distribution in the regions x < 200 km, z < 200 km and x > 6500 km, z < 200 km (shown as black rectangles in Fig. 3a) is reset at each calculation step to simulate formation of a new oceanic plate. The velocity boundary conditions are free-slip everywhere.

To evaluate the effects of the thickness and viscosity of the subducting oceanic plate separately, we treat thickness and viscosity of the SP1 layer (h_SP1 = 30 or 50 km, η_SP1 = 10²³ or 10²⁴ Pa s) as variable parameters (Fig. 3a, b, Table 1). Geodynamic models suggest tectonic evolution of the subduction zones is strongly controlled by the motion of the overriding plate (Faccenna et al. 2017; Schepers et al. 2017; Yang et al. 2019) and absolute plate motion models indicate a fixed Eurasian plate and a westward motion for the South American plate (Becker et al. 2015; Wang et al. 2018). Therefore, we also compare models with an overriding plate either fixed to the sidewall or with a free trailing edge. To model a fixed overriding plate, we set internal boundary conditions of v = 0 at grid points of x = 300 km and z < 100 km. The conditions for each model cases are shown in Table 2.

Table 2 Parameters and conditions used for particular model cases

Low effective activation energy for diffusion creep accompanied by grain growth

Geodynamical modeling of the mantle beneath fracture zones and of flexure around seamounts suggests that the lithospheric mantle has a viscosity with an activation energy of 100–150 kJ/mol assuming a Newtonian rheology for the mantle (Cadio and Korenaga 2016; Courtney and Beaumont 1983; Watts and Zhong 2000; Watts et al. 2013). However, much higher activation energies (300–600 kJ/mol) for creep of mantle materials are predicted by many high-temperature deformation experiments (e.g., Hirth and Kohlstedt 2003) and numerical models with experimentally derived parameters estimate high slab viscosities of > 10²⁴ Pa s (Čížková and Bina 2013; Garel et al. 2014). Based on the results of creep and grain growth experiments for forsterite and enstatite aggregates (Nakakoji et al. 2018), Nakakoji and Hiraga (2018) proposed a rheological model for the lithosphere with a low effective activation energy due to diffusion creep accompanied by grain growth, and provide a possible solution to the apparent discrepancy between estimates based on observations of natural systems and those based on laboratory studies.

High-temperature deformation of polycrystalline materials follows the relationship:

$$\dot{\varepsilon }\propto \frac{{\sigma }^{n}}{{d}^{p}}{\text{exp}}\left(\frac{{Q}_{\text{creep}}}{RT}\right),$$

(6)

where d is the grain size, Q_creep is the activation energy for creep and the n and p are stress and grain size exponents, respectively. The grain size during static grain growth depends on time (t) and T:

$$d \propto \left[ {t{ } \cdot {\text{ exp}}\left( { - \frac{{Q_{{{\text{growth}}}} }}{RT}} \right)} \right]^{1/m} ,$$

(7)

where m is the grain growth exponent and Q_growth is the activation energy for grain growth. The grain growth rate for polymineralic rocks is significantly reduced (m = 3 or 4) compared to monomineralic rocks (m = 2). The secondary phase grains effectively pin grain boundaries of the primary phase. In order for the primary phase to grow, grain growth of the secondary phase due to Ostwald ripening and a corresponding decrease in the number of pinning grains are required.

Results for creep experiments of Nakakoji et al. (2018) showed three creep mechanisms operated at different stress levels: interface-controlled diffusion creep (n = 2 and p = 1) at low stress; grain boundary diffusion (Coble) creep (n = 1 and p = 3) at intermediate stress; and an intragranular dislocation-controlled process (n ≥ 3) at high stress. Nakakoji et al. (2018) obtained comparable activation energies of 430 kJ/mol for Coble creep and 468 ± 65 kJ/mol for grain growth, and the value of the grain growth exponent is estimated to be ~ 4 indicating grain growth proceeded by grain boundary diffusion in a polymineralic system. Nakakoji and Hiraga (2018) suggested that diffusion creep of forsterite, grain growth of forsterite and grain growth of enstatite are all controlled by Si diffusion along forsterite grain boundaries. Assuming Q_creep = Q_growth = Q_GB and p = 3 and m = 4, the viscosity during diffusion creep accompanied by grain growth can be expressed as:

$$\eta \propto t^{\frac{p}{m}} { } \cdot {\text{ exp}}\left( {\frac{{Q_{{{\text{creep}}}} - \frac{p}{m}Q_{{{\text{growth}}}} }}{RT}} \right) \propto t^{\frac{3}{4}} { } \cdot {\text{ exp}}\left( {\frac{{Q_{{{\text{GB}}}} }}{4RT}} \right),$$

(8)

where Q_GB is the activation energy for grain boundary diffusion (Nakakoji and Hiraga 2018). Figure 4 shows viscosity as a function of temperature and pressure using the relation:

$$\eta \propto t^{\frac{3}{4}} { } \cdot {\text{ exp}}\left( {\frac{{E_{{{\text{diff}}}} + V_{{{\text{diff}}}} P}}{RT}} \right),$$

(9)

with low effective values of activation energy and volume (E_diff = 110 kJ/mol, V_diff = 1 cm³/mol, Nakakoji et al. 2018), assuming constant t and a viscosity of 2 × 10²⁰ Pa s at P = 4 GPa and T = 1350 ºC. For a cold slab formed by fast subduction of an old plate, the slab viscosity is estimated to be > 10²⁴ Pa s. For a warm slab formed by slow subduction of a young plate, the slab viscosity is estimated to be < 10²⁴ Pa s and < 10²³ Pa s at the depth of the mantle transition zone.

Fig. 4

Viscosity as a function of pressure and temperature calculated based on Eq. 9 with the values of E_diff = 110 kJ/mol and V_diff = 1 cm³/mol, assuming constant t and the viscosity of 2 × 10²⁰ Pa s at P = 4 GPa and T = 1350 ºC. P–T conditions for cold and warm slabs are shown by black lines

Model with temperature- and pressure-dependent viscosity (Model TPV)

To estimate the stress state of the subducting slab using the rheological model of Nakakoji and Hiraga (2018), we performed simulations using models with temperature- and pressure-dependent viscosities. In this model, the rheology of the upper mantle is expressed by a combination of diffusion creep, dislocation creep and plastic yielding. The model uses the values of E_diff = 110 kJ/mol and V_diff = 1 cm³/mol (Nakakoji et al. 2018) for diffusion creep parameters. Here, we describe the model briefly and details of the model settings are shown in Additional file 1.

The model solves the conservation equations of mass, momentum, and energy, with the extended Boussinesq approximation. The energy equation includes latent heat for phase transitions, radiogenic heating, adiabatic heating and viscous dissipation. The model uses material- and depth-dependent thermal expansion and thermal conductivity. The Clapeyron slopes of the 410-km and 660-km boundaries are assumed to be γ₄₁₀ = 3 MPa/K and γ₆₆₀ = − 1 MPa/K, respectively (Faccenda and Dal Zilio 2017). The viscosity increases by ~ 25 at the 660-km boundary (Additional file 1: Fig. S1c).

The initial temperature distribution of the oceanic plate follows a half-space cooling model, and the initial temperature distribution of the continental plate follows a 1D steady-state conduction model incorporating radiogenic heat production. At the start of the calculation the age of the oceanic plate at the trench (A_tr) is set to be either 40 or 100 Ma. The subducting oceanic plate is free to move in response to the distribution of density and effective viscosity, except for the initial stage of the simulation (until the slab tip reaches 300 km depth) when a constant velocity is applied to the plate.

Results

Evolution of slab geometry and stress state of the overriding plate in Model CV

Figure 5a–d shows model results for the case of a thick (h_SP1 = 50 km) and viscous (η_SP1 = 10²⁴ Pa s) subducting plate with a fixed overriding plate (Fixed5024). As the overriding plate is fixed to the left side of the model, subduction occurs at a largely immobile trench and the overall trajectory of the slab is close to vertical (Fig. 5b–d). The resistance at the 660-km boundary causes folding of the slab to develop in the domain above the 660-km boundary and also an associated periodic change in velocity of the subducting plate (V_sp) (Fig. 5a). The horizontal stress in the overriding plate is extensional throughout the calculation period (Fig. 5i), and this is reflected in the very small but positive values of the trench velocity (V_tr) (Fig. 5a). The stress in the slab locally reaches the maximum yield stress (σ_max) in highly deformed regions, where the effective viscosity (η_eff in Eq. 4a) reduces to ~ 10²³ Pa s (Fig. 5b–d).

Fig. 5

Results for models with constant viscosity (Model CV). a–d Results for a model case with h_SP1 = 50 km, η_SP1 = 10²⁴ Pa s and a fixed overriding plate (Fixed5024). e–h Results for a model case with h_SP1 = 30 km, η_SP1 = 10²³ Pa s and a free overriding plate (Free3023). a, e Change in velocities of the subducting plate (V_sp), the overriding plate (V_op), and the trench (V_tr). b–d, f–h The time evolution of the slab geometry. The distribution of effective viscosity (η_eff) is shown by corresponding colors and arrows indicate the velocity field. i Distribution of horizontal stress (σ_xx) in the uppermost 100 km of the model

Figure 5e–h shows model results for a thin (h_SP1 = 30 km) and low viscosity (η_SP1 = 10²³ Pa s) subducting plate with a free-moving overriding plate (Free3023). In contrast to the case for the fixed overriding plate, the model with a free-moving plate results in trench retreat associated with trenchward motion of the overriding plate. The continuous trench retreat associated with slab rollback results in a gentle overall slab dip (Fig. 5f–h). After the slab first penetrates into the lower mantle, the average dip for the slab at depths < 100 km decreases continuously, and the slab geometry changes to flat-slab subduction after ~ 35 Myr subduction (Fig. 5h). The resistance at the 660-km boundary causes folding of the slab above the 660-km boundary and periodic change in the velocity of the subducting plate (V_sp) (Fig. 5e). The transition to flat-slab subduction is accompanied by a decrease in trench velocity (V_tr) (Fig. 5e) and a change in the horizontal stress in the overriding plate from extension to compression (Fig. 5i). This change in the stress state is reflected in the change in the relative velocities of the overriding plate and trench (V_op and V_tr) (Fig. 5e). The stress in the slab rarely reaches the maximum yield stress (σ_max), and the distribution of the effective viscosity (η_eff) in the subducting slab is almost uniform (Fig. 5f–h).

The calculated changes in velocities of the subducting plate, the overriding plate, and the trench for all eight cases are shown in Additional file 2: Fig. S2. All cases with a freely moving overriding plate exhibit a transition to flat-slab subduction and reduction of the trench velocity after ~ 35 Myr of subduction, irrespective of the thickness and viscosity of the subducting oceanic plate. The stress state of the overriding plate is also controlled by the movement of the overriding plate: tensional for a fixed overriding plate and exhibiting a switch from initial tensional to later compressional for a free-moving overriding plate. The period for the oscillation of V_sp corresponds to the period for slab buckling and is longer for a thick SP1 layer (h_SP1 = 50 km) than for a thin one (h_SP1 = 30 km) (Additional file 2: Fig. S2). The oscillation associated with a half period for the cases with a fixed overriding plate corresponds to the change between ocean-ward and continent-ward slab buckling.

Stress state in the subducting slab in Model CV

Figure 6 shows the stress state in the subducting plate after continuous subduction lasting ~ 45 Myr. For cases with a fixed overriding plate, the slab geometry changes periodically associated with slab folding, although the slab geometry during each phase does not change significantly after the slab first penetrates into the lower mantle. Therefore, for the cases with a fixed overriding plate we show the stress states for two different folding phases. For cases with a free overriding plate, the slab geometry at depths of z < 400 km is similar to that for cases with a fixed overriding plate until the slab penetrates into the lower mantle (Fig. 5f). The geometry subsequently changes to show flat-slab subduction (Fig. 5g, h). After ~ 40 Myr of subduction, the slab geometry at depths of z < 400 km does not show any major change. However, at depths of z > 400 km the slab geometry shows periodic changes due to folding.

Fig. 6

Distribution of stresses for the cases with a h_SP1 = 50 km, η_SP1 = 10²⁴ Pa s and a fixed overriding plate (Fixed5024); b h_SP1 = 30 km, η_SP1 = 10²⁴ Pa s and a fixed overriding plate (Fixed3024); c h_SP1 = 50 km, η_SP1 = 10²³ Pa s and a fixed overriding plate (Fixed5023); d h_SP1 = 30 km, η_SP1 = 10²³ Pa s and a fixed overriding plate (Fixed3023); e h_SP1 = 50 km, η_SP1 = 10²⁴ Pa s and a free overriding plate (Free5024); f h_SP1 = 30 km, η_SP1 = 10²⁴ Pa s and a free overriding plate (Free3024); g h_SP1 = 50 km, η_SP1 = 10²³ Pa s and a free overriding plate (Free5023); and h h_SP1 = 30 km, η_SP1 = 10²³ Pa s and a free overriding plate (Free3023). Color shows the square root of the second invariant of the stress tensor. Directions of the black lines indicate the orientation of the tension axes and the lengths of the lines are proportional to the stress magnitudes. For the cases with a fixed overriding plate, the stress states are shown at two different times with different slab geometries

For cases with a fixed overriding plate, the stress state of the slab changes periodically related to slab folding. During the period when the slab is straight in the depth range of 150–400 km, the stress state shows two pairs of DC and DT associated with bending at z < 150 km and unbending at z > 150 km (left panels of Fig. 6a–d). During the period when the oceanward convex slab hinge is descending to depths of z > 400 km, the stress magnitude in the intermediate depth range (z = 200–350 km) decreases, and for a case of h_SP1 = 30 km and η_SP1 = 10²³ Pa s, the stress state at z < 350 km is close to in-plane DT (right panels of Fig. 6d).

For cases with a free overriding plate, stress states in slabs at < 150 km depths show two pairs of DC and DT associated with the bending geometries of the flat-slab (Fig. 6e–h). The stress states for slabs at 150–350 km depths generally show a pair of DC and DT domains associated with unbending (Fig. 6e–g). One exception is for a case of h_SP1 = 30 km and η_SP1 = 10²³ Pa s, where the stress state at 150–350 km depths is almost entirely DT (Fig. 6h).

These results indicate that the stress state in the slab at z < 150 km is mainly controlled by the bending–unbending process related to the geometry of the slab, which in turn is controlled by the movement of the overriding plate. In contrast, the stress state in the slab at z = 150–350 km is mostly controlled by the thickness and viscosity of the slab: DT for a case of h_SP1 = 30 km and η_SP1 = 10²³ Pa s, and paired DC and DT domains for other model cases. In all eight cases, the stress state for z > 400 km is almost entirely DC.

Effect of the 410- and 660-km boundaries

To evaluate the effect of the phase transitions at 410- and 660-km depths and viscosity increase at the 660-km boundary, we conducted calculations where these conditions were varied. Additional file 2: Fig. S3a–c shows the results for a model that lacks the 410-km phase transition but is otherwise the same as Free3023 (Free3023n410). Except for slower subduction velocity, there is little difference in slab geometry and stress state in the slab compared with the case Free3023. A model with a larger buoyancy due to the 660-km phase transition (Free3023s660), shows that after the folded slab remained around the 660-km boundary, it gradually descends into the lower mantle and a flat-slab geometry is gradually established after ~ 60 Myr of subduction (Additional file 2: Figs. S3d, e). After the transition to the flat-slab subduction, the stress state in the slab in the upper mantle is almost same as for Free3023 (Additional file 2: Fig. S3f). The result for the case with a lower viscosity of the lower mantle (Free3023wLM) shows a larger subduction velocity and a weaker slab folding compared to the case Free3023 (Additional file 2: Fig. S3g, h). However, there is no major difference in stress state in the slab (Additional file 2: Fig. S3i).

Similar effects are found in the case of Fixed5024 (Additional file 2: Fig. S4). The absence of the 410-km phase transition (Fixed5024n410) reduced the subduction velocity (Additional file 2: Fig. S4a). A larger buoyancy due to the 660-km phase transition (Fixed5024s660) caused a tighter slab folding (Additional file 2: Fig. S4e). A lower viscosity of the lower mantle (Fixed5024wLM) resulted in a larger subduction velocity and a weaker slab folding (Additional file 2: Fig. S4g, h). However, these changes in the model did not result in any major difference in the overall stress state in the slab (Additional file 2: Fig. S4c, f, I).

Results of Model TPV

Figure 7 shows the results for Model TPV. For the case of a free overriding plate and a plate age at the trench (A_tr) of 40 Ma, the model shows that after the slab penetrates deep into the lower mantle, the subduction mode changes to a flat-slab subduction at ~ 50 Myr and the trench retreat rate slows down (Fig. 7a, b). The upper mantle slab deforms by diffusion creep (Fig. 7d) and is thin and less viscous—with particularly low values < 10²³ Pa s at z = 400–660 km (Fig. 7b)—due to small values of E_diff and V_diff. This low viscosity slab cannot transmit compressive stress caused by resistance at the 660-km boundary upward, leading to the dominance of DT stress state of the slab at z = 150–350 km (Fig. 7c).

Fig. 7

Results for models with temperature- and pressure-dependent viscosity (Model TPV) assuming E_diff = 110 kJ/mol and V_diff = 1 cm³/mol. a–d Results for the case with A_tr = 40 Ma and a free overriding plate. e–h Results for the case with A_tr = 100 Ma and a fixed overriding plate. a, e Changes in velocities of the subducting plate (V_sp), the overriding plate (V_op), and the trench (V_tr). b, f Slab geometry after deep slab penetration into the lower mantle. The distribution of the effective viscosity (η_eff) is shown in color and arrows indicate the velocity field. c, g Distribution of stress. Color indicates the square root of the second invariant of the stress tensor. Directions of the black lines indicate tension axes and the lengths of the lines are proportional to the stress magnitudes. d, h Distribution of the dominant deformation mechanisms in the upper mantle. Diff diffusion creep, Disl dislocation creep, Yield plastic yield. Dotted black lines indicate the 410-km and 660-km phase boundaries

For the case of a fixed overriding plate and A_tr = 100 Ma, the stress state for a thick and viscous slab in the upper mantle shows a pair of DC and DT domains at z = 100–350 km (Fig. 7f, g). The low-temperature core part of the slab is deformed by plastic yielding, and the adjacent domain is deformed by diffusion creep (Fig. 7h).

The temperature of the subducting slab is controlled by slab age and subduction velocity. For a free overriding plate and A_tr = 40 Ma, the maximum depth of the 600 ℃ isotherm is less than 200 km (Fig. 7c, d). This high slab temperature results from both the young slab age and the extremely low vertical component of the subduction velocity in the flat subduction domain. For a fixed overriding plate and A_tr = 40 Ma, the 600 ℃ isotherm extends to depths > 300 km and the slab has a stress state showing a pair of DC and DT domains at z = 100–350 km (Additional file 2: Fig. S5).

To assess the effect of changing the effective activation energy for creep on the results, we also examined the results for a model with E_diff = 375 kJ/mol and V_diff = 4 cm³/mol—figures that have been estimated from high-temperature experiments on olivine aggregates (Hirth and Kohlstedt 2003; Additional file 2: Fig. S6). Due to the large E_diff and V_diff, the upper mantle slab is deformed by plastic yielding and dislocation creep instead of diffusion creep (Additional file 2: Figs. S6d, h), resulting in a high slab viscosity (> 10²⁴ Pa s) regardless of slab age (Additional file 2: Figs. S6b, f). In addition, the stress state of the upper mantle slab shows two or four pairs of DC and DT domains at depths of < 350 km, the number of pairs depending on the slab geometry (Additional file 2: Figs. S6c, g).

Discussion

Evolution of slab geometry and stress state of the overriding plate in Model CV

The evolution of subduction processes with a freely moving overriding plate in the present model is consistent with previous studies with similar model settings (Schellart and Strak 2021; Strak and Schellart 2021). For the early stage of the subduction before the slab penetrates into the lower mantle, the return flow induced by subduction is mostly confined to the upper mantle near the trench and the trenchward velocity of the upper mantle decreases toward the left side of the model (Fig. 5f). This horizontal velocity gradient in the mantle flow exerts a drag on the overlying plate leading to an extensional stress state within the plate (Fig. 5f, i). In contrast, in the later stage deep penetration of the slab into the lower mantle is associated with the development of a large cell of return flow, which scales in size with the depth of the whole-mantle (Fig. 5h), and this flow exerts a drag over a correspondingly large portion of the overriding plate tending to drive it towards the trench. In addition, the penetration of the deep slab into the lower mantle produces an anchoring of the slab that resists lateral movement of the trench. The combination of the slab anchoring and large-scale convection that drags the overriding plate towards the trench results in compressive stress state in the overriding plate (Fig. 5i).

Several mechanisms have been suggested for the formation of flat slabs including buoyant features on the subducting plate (Arrial and Billen 2013; Axen et al. 2018; van Hunen et al. 2002), trenchward motion of the overriding plate (Arcay et al. 2008), and suction forces enhanced by various features (Rodríguez-González et al. 2012; Stevenson and Turner 1977; Taramón et al. 2015). Among these multiple mechanisms, trenchward motion of the overriding plate against the anchoring slab may be the most effective for the present model, because our model does not include any special buoyant features on the subducting plate and the flat-slab subduction occurs after the deep slab penetration into the lower mantle irrespective of slab thickness and viscosity. The development of flat-slab subduction associated with compressive stress state in the overriding plate after the slab has penetrated into the lower mantle are consistent with the geodynamic evolution of the South American subduction zone (Faccenna et al. 2017; Schellart 2020; Schepers et al. 2017).

For the cases of a fixed overriding plate, trenchward return flow leads to an extensional stress state in the overriding plate irrespective of the slab penetration depth and the size of the return flow cell. These features are consistent with the geodynamic evolution of the western Pacific margins, which is characterized by the formation of many Cenozoic back-arc basins (Sdrolias and Müller 2006) and the immobile nature of the Eurasian plate (Becker et al. 2015; Yang et al. 2019).

Stress state in the subducting slab in Model CV

The stress state in the subducting slab is considered to be a combination of two components, one associated with bending–unbending and the other with uniform in-plane compression or tension due to negative buoyancy and resistance at the 660-km boundary (Fig. 2; Kita et al. 2010; Sandiford et al. 2020; Sippl et al. 2022). The present results show a DC stress state in the slab at z > 400 km in all cases indicating that resistance at the 660-km boundary controls the stress state in the domain immediately above this boundary. At shallower levels of z = 150–350 km our model shows the stress state in the slab is controlled by the thickness and viscosity of the slab. The stress state is characterized by DT for h_SP1 = 30 km and η_SP1 = 10²³ Pa s, and by pairs of DC and DT domains for other cases. Because thin and less viscous slabs have relatively low resistances to bending, they respond to the compressive stress due to resistance at the 660-km boundary by folding and the stress is not transmitted up dip along the slabs. Therefore, at z = 150–350 km, the in-plane DT stress state due to negative buoyancy dominates over both the in-plane DC stress state due to the resistance at the 660-km boundary and the paired DT and DC stress state due to slab unbending. In addition, stress associated with the bending–unbending deformation is smaller for a thinner slab, because the strain rate associated with bending–unbending deformation increases with the distance from the neutral plane.

The change in the negative and positive buoyancy due to the phase transitions at the top and bottom of the transition zone can be estimated from the upward and downward displacement of the phase boundaries. For the results of Model TPV, the 2D triangular areas of upward and downward displaced phase boundaries are about 4600 km² and 1600 km², respectively, for the case of free overriding plate and A_tr = 40 Ma (Fig. 7d) and are about 8000 km² and 1800 km², respectively, for the case of fixed overriding plate and A_tr = 100 Ma (Fig. 7h). For Model CV with a constant viscosity, the corresponding areas are (h_SP1 + h_SP2) × 80 km and (h_SP1 + h_SP2) × 30 km, respectively. These values are 4800 km² and 1800 km², respectively, for the case of h_SP1 = 30 km and are 6400 km² and 2400 km², respectively, for the case of h_SP1 = 50 km. The displacement of the phase boundaries depends on the temperature distribution and thus the slab age and the subduction velocity in addition to the values for the Clapeyron slopes. Therefore, the conditions for Model CV roughly correspond to the values for the Clapeyron slopes of γ₄₁₀ = 3 MPa/K and γ₆₆₀ = − 1 MPa/K.

Rheological properties of the subducting slab

The results of Model CV (Figs. 5, 6) indicate that a thin and low viscosity (10²³ Pa s) slab is required to account for the dominantly DT stress state of the slab at depths of 100–300 km observed for the eastern Pacific. This is consistent with the young slab ages of the slabs of the eastern Pacific. However, relatively high slab viscosities of > 10²⁴ Pa s are expected from model calculations using experimentally derived rheological parameters for olivine aggregates even for young slab ages (Additional file 2: Fig. S6b; Garel et al. 2014). The results of Model TPV indicate that the first-order stress state of the slab estimated from intermediate-depth earthquakes in both the western and eastern Pacific subduction segments can be adequately explained by slab deformation controlled by diffusion creep of polyphase rocks (Fig. 7). The results of Model TPV are also consistent with the observation that there is a positive correlation between the thermal parameter (the product of the slab age and its vertical rate of descent) and the fraction of the intermediate-depth earthquake populations associated with the DC stress (Chen et al. 2004). The difference in stress states between the eastern Pacific and the western Pacific subduction segments is mainly caused by the age of subducting slab with a more limited contribution from the motion of the overriding plate.

The fine-grained mineral aggregates required for the polyphase model compatible with our model results may form in the oceanic lithosphere by several subduction processes including brittle fracturing and formation of hydrous minerals along faults developed in the outer rise–trench region (Boneh et al. 2019; Ranero et al. 2003) and the dehydration reactions that may be the cause of intermediate-depth earthquakes (Ferrand et al. 2017). If these fine-grained shear zones are formed at sufficiently narrow intervals, then on the scale of the whole slab deformation could be approximated by homogeneous diffusion creep with the conditions of slow grain growth for polyphase rocks.

For the dominant DT stress state in the domains of downward slab bending (z = 100–150 km) at the Chile and Peru segments, Sandiford et al. (2020) suggest that the seismic strain release would be restricted to the low temperature part of the slab above the neutral plane. For the deeper part of the slab in the Chile segment, which is associated with an unbending geometry (z > 150 km), in addition to the possibility of in-plane DT stress due to slab pull, Sandiford et al. (2020) suggest that progressive warming of the upper part of the slab and/or processes related to dehydration may leave it essentially aseismic at these depths, with the consequence that the seismically active zone lies beneath the midplane. The results for Model TPV with a free overriding plate and A_tr = 40 Ma using E_diff = 375 kJ/mol and V_diff = 4 cm³/mol show paired DT and DC stress states in the slab at intermediate depths (Additional file 2: Fig. S6c). The neutral plane in the slab at ~ 160 km depth is located below the coldest part of the slab and the temperature of the DT domain is lower than that of DC domain (Additional file 2: Fig. S7). This is consistent with the interpretation of Sandiford et al. (2020). However, the neutral plane in the slab at ~ 200–300 km depth is located along the coldest part of the slab and there is no major difference in temperature between DC and DT domains (Additional file 2: Fig. S7). This result suggests the simple explanation that dominant DT stress is caused by negative buoyancy of the slab is more likely.

Slabs beneath the western Pacific margins show variable geometries and some of them may reflect active or past back-arc spreading. Therefore, application of the model results that assume largely immobile trench to a specific subduction segment is difficult. However, the nearly vertical slab beneath the Mariana segment which shows dominant DT stress state at intermediate-depth (Alpert et al. 2010; Billen 2020; Chen et al. 2004) may be an example of the case shown in right panels of Fig. 6a–d.

Conclusions

To investigate the rheological properties of the subducting slab, we compared the stress state in the subducting slabs below the eastern and western Pacific margins with results of 2D numerical modeling. The simulation results of models with constant thickness and viscosity of the slab indicate that a thin and low viscosity (10²³ Pa s) slab is required to account for the dominantly DT (down-dip tension) stress state of the slab at depths of 100–300 km observed for the eastern Pacific. Such weak slabs only display limited resistances against bending and thus respond to increased resistance at the 660-km boundary by folding and cannot transmit the compressive stress upward.

High-temperature deformation experiments of olivine aggregates generally yield high activation energies for creep implying a high viscosity (> 10²⁴ Pa s) for even young slabs. However, the effective activation energy for diffusion creep of polyphase mantle materials has been estimated to be as low as 110 kJ/mol due to slow grain growth. Incorporating this low value for the activation energy for creep in our 2D numerical simulations with temperature- and pressure-dependent viscosity results in DT stress domains dominating at 100–300 km depths for young slabs and a pair of DT and DC stress domains at the same depths for old slabs corresponding to the unbending process. These results are compatible with the observed distributions of stress domains in the slabs of the Pacific realm.