the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The effect of temperaturedependent material properties on simple thermal models of subduction zones
Iris van Zelst
Cedric Thieulot
Timothy J. Craig
To a large extent, the thermal structure of a subduction zone determines where seismicity occurs through controls on the transition from brittle to ductile deformation and the depth of dehydration reactions. Thermal models of subduction zones can help understand the distribution of seismicity by accurately modelling the thermal structure of the subduction zone. Here, we assess a common simplification in thermal models of subduction zones, i.e. constant values for the thermal parameters. We use temperaturedependent parameterisations, constrained by lab data, for the thermal conductivity, heat capacity, and density to systematically test their effect on the resulting thermal structure of the slab. To isolate this effect, we use the welldefined, thoroughly studied, and highly simplified model setup of the subduction community benchmark by van Keken et al. (2008) in a 2D finiteelement code. To ensure a selfconsistent and realistic initial temperature profile for the slab, we implement a 1D plate model for cooling of the oceanic lithosphere with an age of 50 Myr instead of the previously used halfspace model. Our results show that using temperaturedependent thermal parameters in thermal models of subduction zones affects the thermal structure of the slab with changes on the order of tens of degrees and hence tens of kilometres. More specifically, using temperaturedependent thermal parameters results in a slightly cooler slab with e.g. the 600 ^{∘}C isotherm reaching almost 30 km deeper. From this, we infer that these models would predict a larger estimated seismogenic zone and a larger depth at which dehydration reactions responsible for intermediatedepth seismicity occur. We therefore recommend that thermo(mechanical) models of subduction zones take temperaturedependent thermal parameters into account, especially when inferences of seismicity are made.
 Article
(11484 KB)  Fulltext XML

Supplement
(60421 KB)  BibTeX
 EndNote
The thermal structure of subduction zones plays a vital role in controlling many geological and petrological processes, including the dehydration of the subducting plate (Peacock, 2001; Hacker et al., 2003a), the subsequent hydration of the mantle and overriding plate (Peacock, 1993a; Abers et al., 2017), and mineralogical variations, including serpentinisation (Hyndman and Peacock, 2003). Furthermore, seismicity can often be related to both the thermal structure and the various processes controlled by the pressure and temperature evolution of the slab (Scholz, 2019). For example, intermediatedepth earthquakes are associated with a process called dehydration embrittlement (e.g. Green and Houston, 1995; Peacock, 2001; Hacker et al., 2003b; Yamasaki and Seno, 2003; Jung et al., 2004; Wang et al., 2017), in which water is released during the compositional evolution of the slab as hydrous minerals progressively transform to less hydrous phases (e.g. from blueschist to eclogite; van Keken et al., 2011). The addition of free fluids to the system acts against the pressure of the surrounding rock, permitting earthquakes to occur at depths where the confining pressure is otherwise too large. These phase transitions are linked to specific temperature and pressure conditions, suggesting that a thorough grasp of those conditions at depth could indicate where intermediatedepth seismicity would be likely to occur (e.g. Hacker et al., 2003a). Similarly, megathrust earthquakes occur within the seismogenic zone, the downdip limit of which is thought to be the transition from brittle to ductile deformation (Peacock and Hyndman, 1999; Scholz, 2019) and is again controlled, directly or indirectly, by temperature, with isotherms of 350–450 ^{∘}C typically linked to this change (Hyndman and Wang, 1993; Hyndman et al., 1997; Gutscher and Peacock, 2003).
From these examples, it becomes clear that it is important to have a thorough understanding of the thermal structure of a slab in order to better understand the distribution of the full spectrum of seismicity associated with the subduction process. However, it is hard to obtain direct observational data on the thermal structure of the slab due to the inaccessibility of subduction zones and the difficulty of obtaining data at great depths (i.e. larger than 10 km).
The dependence of seismic wave speeds on temperature allows seismic tomography studies to give a broad overview of the largescale thermal structure of the subduction zone as a whole, but such studies typically lack the resolution to infer the thermal structure of the slab itself in great detail (e.g. Abers et al., 2006; Pozgay et al., 2009). In addition, the observed velocity anomalies in tomographic models are not exclusively due to temperature, and wave speed variations are also related to other factors, particularly composition, density, mineralogy, and the presence of fluids (e.g. Hacker et al., 2003a; Blom et al., 2017). Whilst borehole experiments and marine heat flow measurements can provide vital insights into the thermal state of the shallow seismogenic zone (e.g. Hyndman and Wang, 1993; Chang et al., 2010; Fulton et al., 2013; Harris et al., 2013; Yabe et al., 2019), such measurements are extremely local and fail to give a good overview of the conditions of the subduction zone as a whole, especially the finer details of the temperature structure within the slab.
In light of the limited available data on the thermal structure of subduction zones, geodynamic numerical modelling provides a way of investigating the complete temperature field of subduction zones in relation to the thermal and dynamic evolution of the slab (see Peacock, 2020, for an overview). The starting point for thermal models of subduction zones are onedimensional models of the cooling of the oceanic lithosphere that define the thermal structure of the slab for a certain plate age, including halfspace cooling models and more advanced plate models (McKenzie and Sclater, 1969; Parsons and Sclater, 1977; Stein and Stein, 1994; Hillier and Watts, 2005; McKenzie et al., 2005; Crosby et al., 2006; Emmerson and McKenzie, 2007; Richards et al., 2018). Extending this thermal modelling to two dimensions to study the thermal evolution of a subduction zone has provided insights into the predicted location of dehydration and melting processes linked to intermediatedepth seismicity (Ponko and Peacock, 1995; Peacock and Wang, 1999; van Keken et al., 2002; Abers et al., 2006; Syracuse et al., 2010; van Keken et al., 2012, 2019). Apart from pure thermal models, thermomechanical models with various complexities such as melting and dehydration reactions have also been employed (e.g. Gerya and Meilick, 2011; Gerya, 2011; Faccenda et al., 2012; Arcay, 2017; Beall et al., 2021), leading to insights into subduction dynamics and estimates of the depth of intermediatedepth seismicity and the geometry of the megathrust. When these types of models additionally account for an inertia term in socalled seismothermomechanical models (van Dinther et al., 2013a, b), megathrust slip events are resolved, allowing for estimates of the maximum size of the seismogenic zone and the distribution of seismicity in a given subduction geometry (van Dinther et al., 2014; Van Zelst et al., 2019; Petrini et al., 2020; Brizzi et al., 2020). These types of modelling have the advantage that the temperature can be calculated across the entire subduction zone with arbitrary resolution. However, their results depend on their initial and boundary conditions and the assumptions that enter the models at various stages (Van Zelst et al., 2022).
Indeed, numerical models of the temperature structure of subduction zones are subject to a range of simplifications. One, which we seek to address here, is that the thermal parameters in the model, i.e. the thermal conductivity, heat capacity, and density, are commonly assumed to be constant or merely material dependent. In contrast, laboratory experiments have shown that these parameters actually depend on temperature and can differ by as much as a factor of 2 depending on the temperature (e.g. Berman, 1988; Berman and Aranovich, 1996; Seipold, 1998; Hofmeister, 1999; Xu et al., 2004; Wen et al., 2015; Su et al., 2018). The inclusion of such parameters into models for the cooling of the oceanic lithosphere has made a significant difference to both the resulting thermal structure and its interpretation and implications (Denlinger, 1992; McKenzie et al., 2005; Richards et al., 2018). Previous onedimensional (Emmerson and McKenzie, 2007) and threecomponent (Chemia et al., 2015) slab studies have highlighted the potential for a similar impact on the more complex thermal structure of subduction zones.
Given the sensitivity of the various processes mentioned above to smallscale variations in the temperature evolution of the slab, we therefore seek to quantify the potential impact of the temperature dependence of thermal parameters on subduction zone thermal structure and to build towards their routine incorporation.
In order to assess the effect of temperaturedependent thermal parameters on the resulting thermal structure of the slab, we perform a systematic study by using the welldefined, highly simplified setup of the subduction community benchmark by van Keken et al. (2008) with the addition of temperaturedependent functions for the thermal conductivity, heat capacity, and density, as constrained by laboratory experiments (Sect. 2). We show that using temperaturedependent parameters in geodynamic models changes the resultant thermal structure of the slab relative to models with fixed values (Sect. 3). To relate this change in thermal structure to expected seismicity and mineralogical changes in the slab, we discuss the change in the expected depth of intermediatedepth seismicity when temperaturedependent thermal parameters are taken into account (Sect. 4). Going forward, we recommend the inclusion of temperaturedependent thermal parameters in future thermal models of subduction zones, especially when inferences on seismicity are made.
We base our models on the subduction zone community benchmark presented by van Keken et al. (2008). We use the tailormade twodimensional finiteelement Python code xFieldstone (https://github.com/irisvanzelst/xFieldstone, last access: 22 May 2023; Van Zelst, 2023) to solve the incompressible Stokes equations with Crouzeix–Raviart (Crouzeix and Raviart, 1973) elements (also used in the MILAMIN code; Dabrowski et al., 2008) and the conservation of energy using quadratic triangular elements. xFieldstone is based on Fieldstone #68, which is part of the opensource Fieldstone collection of educational finiteelement codes in computational geodynamics (https://cedrict.github.io/, last access: 22 May 2023). The exact version of xFieldstone used to produce the results presented in this work can be found in Van Zelst (2023).
In the following, we first discuss the governing equations (Sect. 2.1) and rheology (Sect. 2.2) of the physical model. We then present the model setup (Sect. 2.3), our formulation for the thermal structure of the oceanic plate at the trench on the left side of the model (Sect. 2.4), and the different functions we consider for the temperature dependence of the thermal parameters (Sect. 2.5). Based on these functions, we define the parameter space of this study (Sect. 2.6) and detail the model diagnostics used in this work (Sect. 2.7).
2.1 Governing equations
Following van Keken et al. (2008), we solve the incompressible formulation of the conservation of mass and momentum (i.e. the Stokes equations) for velocity v and pressure p:
where σ^{′} is the deviatoric stress tensor. In this formulation of the Stokes equations, we implicitly assume zero gravitational acceleration, such that we have purely kinematically driven subduction. This allows us to have the same (fixed) subduction geometry in all models.
We also solve for temperature T using the steadystate conservation of energy without external heat sources such as radiogenic or shear heating to simplify the model as follows:
where ρ is density, C_{p} is the heat capacity, and k is the thermal conductivity. Unlike van Keken et al. (2008), we make these thermal parameters temperature dependent instead of constant, as described in Sect. 2.5.
2.2 Rheology
We consider a purely viscous rheology, where we relate the deviatoric stress tensor σ^{′} to the viscosity η and the deviatoric strain rate tensor $\dot{\mathit{\epsilon}}$ according to
Here, the deviatoric strainrate tensor $\dot{\mathit{\epsilon}}$ is defined as
Initially, we run sets of models with different viscous rheologies to successfully reproduce the different benchmark cases presented in van Keken et al. (2008) (Sect. S1 in the Supplement; Figs. S1–S7 in the Supplement). In the following, we confine ourselves to a rheology that combines the diffusion and dislocation creep mechanisms used in van Keken et al. (2008). We implement this temperaturedependent rheology through an effective viscosity η_{eff}.
For the diffusion creep rheology, we use the simplified diffusion creep viscosity formulation η_{diff} for olivine, where we set the activation volume to zero as it multiplies with the full pressure, whereas the code only computes the excess pressure because we neglect buoyancy forces (Eq. 2). Additionally, we ignore any effects caused by hydration and grain size dependence, resulting in the following rheology:
where A_{diff} is a prefactor, E_{diff} is the activation energy, and R is the universal gas constant. Similarly, we use the following expression for a dislocation creep rheology:
where A_{disl} is a prefactor, E_{disl} is the activation energy, n is the powerlaw exponent, and ${\dot{\mathit{\epsilon}}}_{\mathrm{II}}=\sqrt{{\dot{\mathit{\epsilon}}}_{xx}^{\mathrm{2}}+{\dot{\mathit{\epsilon}}}_{xy}^{\mathrm{2}}}$ is the square root of the second invariant of the deviatoric strain rate tensor.
We combine these formulations for diffusion and dislocation creep into one rheology by assuming two viscous dampers in series (Schmeling et al., 2008):
To avoid unrealistically high stresses, we limit the maximum viscosity in the model to η_{max}=10^{26} Pa s for both the diffusion and dislocation creep rheology, such that the effective viscosity η_{eff} becomes
2.3 Model setup
We use the twodimensional model setup of the community benchmark for subduction zone modelling presented by van Keken et al. (2008) (Fig. 1). This setup is a highly simplified representation of a subduction zone. Its simplicity and the fact that this setup is well defined and thoroughly studied as a benchmark allows us to isolate the effect of the temperaturedependent thermal properties in this study. Indeed, the simplicity of the model geometry allows us to add complexities in other parts of the model, i.e. in this case, temperaturedependent thermal properties. The modelling philosophy of this study is best described as generic according to Van Zelst et al. (2022) as we do not aim to obtain realistic results directly comparable to any specific subduction zone.
We consider a domain that is L_{x}=660 km wide and L_{y}=600 km deep, with the origin of the coordinate system at the lower left corner and with the y axis being positive upwards. We discretise the domain by means of a structured triangular grid with a uniform resolution of 2.5 km, resulting in 528×480 triangular elements. We define a simple slab geometry with a 45^{∘} dip angle originating at the top left corner and a 50 km thick overriding plate at the top of the model. The remaining part of the model is the mantle wedge. Our chosen resolution ensures that the computational grid aligns with the bottom of the overriding plate and the wedge corner.
We fix the overriding plate by prescribing no slip (i.e. zero velocity in both the x and y direction) at its bottom boundary with the mantle wedge. We define the plate kinematics such that the downgoing slab subducts with a constant velocity of 5 cm yr^{−1} by prescribing this velocity at the top of the slab from the corner point at x=50 km and y=550 km to the bottom of the domain. At the corner point itself, we prescribe zero velocity.
For the conservation of energy, we apply a constant 0 ^{∘}C temperature boundary condition along the top of the model domain. At the righthand boundary, we apply a linear temperature gradient in the overriding plate from T=0 ^{∘}C at the top to 1300 ^{∘}C at the bottom of the overriding plate at y=550 km. Below that, incoming material (i.e. v_{x}<0) is assigned the maximum temperature in the model T_{max}=1300 ^{∘}C. At the left boundary, we apply either a halfspace cooling model with a slab age t_{s} of 50 Myr and constant thermal parameters (as used in the benchmark of van Keken et al., 2008) or a temperature profile extracted at 50 Myr from a onedimensional cooling plate model (following Richards et al., 2018, and discussed further in Sect. 2.4). The initial temperature field is constant, with T=0 ^{∘}C.
We first solve the Stokes equations across the entire domain. As we are only interested in the velocity field in the mantle wedge, we overwrite the resulting velocity solution in the slab and overriding plate by our boundary conditions, i.e. no slip in the overriding plate and a constant subduction velocity of 5 cm yr^{−1} in the slab. With the velocity solution determined, the heat equation is solved next. We then iteratively solve the Stokes and heat equation until convergence is reached, i.e. when the horizontal and vertical components of the velocity and the temperature compared to the previous iteration change less than a given tolerance. We choose a relative tolerance of 10^{−5} in our model runs for both velocity and temperature, although we also impose a maximum number of 50 iterations to limit the wall time of the model. Tests show that employing a tolerance of 10^{−3} (reached before 50 iterations) changes the model diagnostics from Sect. 2.7 by less than 1 ^{∘}C and has no effect on the reported isotherm depths. To prevent numerical oscillations in the solution that inhibit convergence of the temperature field, we limit the change of each new iteration solution via a relaxation parameter r after the first iteration according to
where ϕ is the solution for v_{x}, v_{y}, and T. This relaxation step is applied to the velocity components after solving the Stokes equations and to the temperature after solving for the heat equation. After trial and error, we choose r=0.8, which prevents numerical oscillations in the solution towards convergence.
2.4 Thermal structure of the oceanic plate at the trench
As the lefthand boundary condition for the conservation of energy, we prescribe the thermal structure of the incoming oceanic plate. In the original community benchmark, van Keken et al. (2008) used a simple halfspace cooling model (Turcotte and Schubert, 2002) with a plate age of 50 Myr with constant values for the thermal conductivity, heat capacity, and density:
However, the halfspace cooling model does not satisfy petrological constraints and fails to satisfy heat flow and bathymetric data for plate ages greater than ∼80 Myr (Richards et al., 2018). Hence, we follow Richards et al. (2018) by including a more complex and realistic plate model as input for the temperature profile of the incoming oceanic plate. This plate model also has the advantage that it easily incorporates temperaturedependent thermal parameters, which results in consistency between the thermal structure of the incoming plate and the thermal structure we solve for in the rest of the domain.
We calculate the structure of the incoming oceanic plate in a linked, separate Python script, with the coordinate convention that the y axis is positive downwards. The thermal structure of the oceanic plate is based on the following onedimensional heat equation (Turcotte and Schubert, 2002; McKenzie et al., 2005):
Following McKenzie et al. (2005) and Richards et al. (2018), we discretise this equation using a onedimensional time and spacecentred Crank–Nicolson finitedifference scheme that is stable in both space and time and solve it numerically with a predictor–corrector step (Press et al., 1992) according to
with
where m=n for the predictor step and $m=n+\frac{\mathrm{1}}{\mathrm{2}}$ for the corrector step. Additionally, we have
for the predictor step and
for the corrector step.
As input parameters, we choose a constant Δz of 1000 m and a constant time step Δt=1000 years. We have the same temperature boundary conditions as the 2D model domain for consistency, with a surface temperature of 0 ^{∘}C and a maximum temperature (mantle potential temperature) of 1300 ^{∘}C. We choose a plate thickness of 106 km in accordance with the optimum plate thickness found by Parsons and Sclater (1977), Sclater et al. (1980), and McKenzie et al. (2005) based on heat flow observations. We recognise that this plate thickness diverges from the results of Richards et al. (2018), but their results involved the inclusion of compositional variability in addition to the thermal dependence of material properties, which we do not include here. Hence, we use the older plate thickness value determination of McKenzie et al. (2005).
We solve for the temperature evolution of the incoming oceanic plate with the desired thermal parameters (Sect. 2.5) for 200 Myr, which we store in a lookup table (Figs. S8–S17 in the Supplement). The main part of the code then extracts the relevant temperature profile for a plate age of 50 Myr (van Keken et al., 2008) as input for the left boundary of the model domain, taking into account the different coordinate system conventions and the cubic interpolation between the 1D finitedifference coordinates and finiteelement nodes in case of differing resolutions. We then solve the entire system using tridiagonal elimination.
2.5 Temperaturedependent thermal parameters
We use temperaturedependent expressions for the thermal conductivity, heat capacity, and density using parameterisations based on observational experimental data for the way in which these values change with temperature.
Following McKenzie et al. (2005), we approximate the analytical expression for temperaturedependent thermal conductivity k (Fig. 2a) by Hofmeister (1999) with
where k has units of watts per metre kelvin (W m^{−1} K^{−1}), although T is in ^{∘}C in this expression, and where b=5.3 W m^{−1}, c=0.0015, ${d}_{\mathrm{0}}=\mathrm{1.753}\times {\mathrm{10}}^{\mathrm{2}}$ W m^{−1} K^{−1}, ${d}_{\mathrm{1}}=\mathrm{1.0365}\times {\mathrm{10}}^{\mathrm{4}}$ W m^{−1} K^{−2}, ${d}_{\mathrm{2}}=\mathrm{2.2451}\times {\mathrm{10}}^{\mathrm{7}}$ W m^{−1} K^{−3}, and ${d}_{\mathrm{3}}=\mathrm{3.4071}\times {\mathrm{10}}^{\mathrm{11}}$ W m^{−1} K^{−4} are constants. This expression considers both heat transport and the radiative heat transfer by phonons.
Like McKenzie et al. (2005), we also implement the temperaturedependent conductivity for olivine proposed by Xu et al. (2004) to account for the large uncertainties in the temperature dependence of the thermal conductivity:
where T is in ^{∘}C, k_{298}=4.08 W m^{−1} K^{−1} and n=0.406.
For the heat capacity C_{p} (Fig. 2b), we follow Berman (1988) to calculate the heat capacity of both fayalite and forsterite (McKenzie et al., 2005) such that we have
where C_{p} is the heat capacity in joules per kilogram kelvin (J kg^{−1} K^{−1}) and T is in K. We use updated values for the constants according to Berman and Aranovich (1996), resulting in ${a}_{\mathrm{0},\mathrm{fa}}=\mathrm{252}$ J mol^{−1} K^{−1}, ${a}_{\mathrm{1},\mathrm{fa}}=\mathrm{20.137}\times {\mathrm{10}}^{\mathrm{2}}$ J mol^{−1} K^{−½}, and ${a}_{\mathrm{3},\mathrm{fa}}=\mathrm{6.219}\times {\mathrm{10}}^{\mathrm{7}}$ J K^{2} mol^{−1} for fayalite and ${a}_{\mathrm{0},\mathrm{fo}}=\mathrm{233.18}$ J mol^{−1} K^{−1}, ${a}_{\mathrm{1},\mathrm{fo}}=\mathrm{18.016}\times {\mathrm{10}}^{\mathrm{2}}$ J mol^{−1} K^{−½}, and ${a}_{\mathrm{3},\mathrm{fo}}=\mathrm{26.794}\times {\mathrm{10}}^{\mathrm{7}}$ J K^{2} mol^{−1} for forsterite. To obtain the heat capacity in the correct unit of joules per kilogram kelvin, we multiply the equation from Berman (1988) where C_{p} is in joules per mole kelvin (J mol^{−1} K^{−1}) with $\frac{\mathrm{1000}}{{m}_{\mathrm{fa}\mathrm{fo}}}$, where m_{fafo} is the molecular mass of fayalite (fa) or forsterite (fo). We then obtain the effective heat capacity in the model by assuming a molar fraction f=0.11 of fayalite in the mantle according to McKenzie et al. (2005):
For the dependency of density on temperature (Fig. 2c), we follow the parameterisation of McKenzie et al. (2005) based on the integration of the parameterisation of the temperature dependence of the thermal expansivity according to Bouhifd et al. (1996):
where T is in K, ρ_{0}=3330 kg m^{−3}, T_{0}=273.15 K, ${\mathit{\alpha}}_{\mathrm{0}}=\mathrm{2.832}\times {\mathrm{10}}^{\mathrm{5}}$ K^{−1}, and ${\mathit{\alpha}}_{\mathrm{1}}=\mathrm{3.79}\times {\mathrm{10}}^{\mathrm{8}}$ K^{−2}.
Formulations for the temperature dependence of the thermal conductivity, heat capacity, and density other than the ones described here are also available (e.g. Berman and Brown, 1985; Seipold, 1998; Wen et al., 2015; Su et al., 2018), but here we limit ourselves to the formulations described in this section to test the effect of such variability.
In our preferred formulations for the temperature dependence of the thermal parameters, the thermal conductivity (Hofmeister, 1999; McKenzie et al., 2005) varies by a factor of 2.5 for the temperature range in our subduction zone models (i.e. 0–1300 ^{∘}C) with k=5.3 W m^{−1} K^{−1} for T=0 ^{∘}C and k=2.1 W m^{−1} K^{−1} for T=1300 ^{∘}C (Fig. 2). Similarly, the heat capacity (89 % forsterite; Berman and Aranovich, 1996) varies by a factor of 1.65 over the temperature range 0–1300 ^{∘}C. The temperature dependence of the density (Bouhifd et al., 1996) is less pronounced, with a variation of approximately 9 % in density for temperatures common to subduction zones. However, the formulations used here for the thermal conductivity, heat capacity, and density are mainly applicable to mantle rocks, which do not typically experience low temperatures of T<150 ^{∘}C. Instead, crustal rocks with different thermal properties are a better approximation at these temperatures, which is addressed further in the next section.
2.6 Parameter space
To systematically test the effect of temperaturedependent parameters on the thermal structure of subduction zones, we run the suite of simulations outlined in Table 1. We start with a reference model case2c_PvK based on the van Keken et al. (2008) benchmark models. Note that this model is not in the original suite of benchmark models of van Keken et al. (2008), the difference being that the rheology employed is a combination of diffusion and dislocation creep. Then, we first test the effect of adding the more complex temperature profile of the plate model instead of the halfspace cooling model in case2c_bc. We test the effect of the two different functions for the thermal conductivity with models case2c_k1 and case2c_k2, where the approximation of the function by Hofmeister (1999) is our preferred function, following McKenzie et al. (2005). Our preferred model for the heat capacity is the one where we use the function of Berman (1988) with the updated values of Berman and Aranovich (1996) for a composition of 89 % of forsterite and 11 % fayalite (case2c_Cp6). We also test the effect of using the older values of Berman (1988) (case2c_Cp3), and a composition of 100 % forsterite (case2c_Cp4) and 100 % fayalite (case2c_Cp5). Here, the numbers behind k and C_{p} in the model names refer to the flags used in the code to select different options for the temperaturedependent thermal parameters. We test the temperaturedependent density in case2c_rho. Finally, we combine our preferred functions of the thermal parameters in simulation case2c_all.
McKenzie et al. (2005)Hofmeister (1999)Xu et al. (2004)Berman and Aranovich (1996)Berman (1988)Berman and Aranovich (1996)Berman and Aranovich (1996)McKenzie et al. (2005)Bouhifd et al. (1996)McKenzie et al. (2005)Hofmeister (1999)Berman and Aranovich (1996)McKenzie et al. (2005)Bouhifd et al. (1996)McKenzie et al. (2005)Hofmeister (1999)Berman and Aranovich (1996)McKenzie et al. (2005)Bouhifd et al. (1996)McKenzie et al. (2005)Hofmeister (1999)Berman and Aranovich (1996)McKenzie et al. (2005)Bouhifd et al. (1996)^{a} None of the models listed in this table include the oceaniccrust parameterisation in the slab or an overriding plate. Instead, the first 10 models listed here also represent three additional batches of models. Models with the addition _mc (mimic crust) include the oceaniccrust parameterisation in the slab without an overriding plate. Models with the addition _op (oceanic plate) include the oceaniccrust parameterisation in the slab and a parameterisation for an oceanic upper plate. Models with the addition _cp (continental plate) also include the oceaniccrust parameterisation in the slab and a continental overriding plate parameterisation. So, for example, model case2c_PvK does not include any oceaniccrust parameterisation or overriding plate, while model case2c_PvK_cp does include the oceaniccrust parameterisation in the slab and a parameterisation for a continental upper plate. ^{b} The constant values used for the thermal conductivity k, heat capacity C_{p}, and density ρ are taken from van Keken et al. (2008) to be k=3 W m^{−}1 K^{−1}, C_{p}=1250 J kg^{−1} K^{−1}, and ρ=3300 kg m^{−3}.
To illustrate how the different simulations differ in terms of the temperature dependence of the thermal parameters, we show the thermal diffusivity κ in Fig. 3, calculated according to
Hence, when all thermal parameters k, C_{p}, and ρ are temperature dependent in model case2c_all, the overall temperature dependency of the model is largest. Compared to the constant thermal diffusivity used in the benchmark by van Keken et al. (2008), values are up to 319 % larger and up to 28 % smaller for the temperature range of our model. Large values for the diffusivity translate to rapid heat transfer, meaning that cold regions heat up faster and that hot regions cool down faster. In general, Fig. 3 shows that the thermal diffusivity is higher for low temperatures, meaning that the cold top of the slab will be heated faster.
The top of the oceanic plate, where temperatures are low and hence where thermal parameters differ most from the constant values used in van Keken et al. (2008), is the area of the plate where our assumption of a uniform mantle composition is most likely to be inappropriate. Crustal materials also have temperaturedependent thermal properties, but, with different mineralogical compositions to the olivinedominated mantle, they can have quite different values, with conductivities for crustal materials at temperatures <150 ^{∘}C typically being lower than those of mantle rocks (Fig. 2; e.g. Grose and Afonso, 2013; Richards et al., 2018). Given the effect that such variation would have on the rate of temperature change in the top few kilometres of the slab and the potential for this to insulate the deeper sections of the slab from the effects of the mantle wedge, we test the impact that incorporating a crustal layer in the subducting slab in our models might have using a simple parameterisation (see Sect. 4.2 and 4.3 for a discussion on this simplification and alternatives). In order to approximate an oceaniccrustal layer in the slab in our models, we define a crustal thickness of 7 km, for which we use half the value of the thermal conductivity for a given temperature (e.g. Grose and Afonso, 2013; Fig. 2). We then run our main set of 10 models with the addition of this oceaniccrust parameterisation in the slab to assess how the inclusion of explicitly modelled crustal heterogeneity might affect the results presented in this work. We indicate this set of models with the postfix _mc (mimic crust) in our simulation names.
Another important aspect in a subduction zone is the compositional structure of the overriding plate. In our standard set of models, there is no crustal layer within the overriding plate. To assess the effect of the crust of the overriding plate on the thermal structure of the slab, we run another two sets of models based on our main set of 10 models and including the oceaniccrust parameterisation in the slab of the _mc models but with an additional crustal layer in the overriding plate. The first set of models (denoted _op – oceanic plate) incorporates a crustal layer with a thickness of 7 km and half the thermal conductivity values for a given temperature, similarly to our parameterisation of oceanic crust in the slab, as if the overriding plate were oceanic in origin. The second set of models considers a continental upper plate (denoted _cp – continental plate) with a crustal thickness of 35 km, halved thermal conductivity values, and a density that is multiplied by 0.79 to obtain realistic values for crustal densities at low temperatures.
Note that for these three sets of models that include a parameterisation for oceanic crust in the slab, we also include this crustallayer approximation in the onedimensional plate model that calculates the left temperature boundary condition. In contrast, the halfspace cooling model used by van Keken et al. (2008) does not allow for the addition of layers. Therefore, the boundary condition in these case2c_PvK models is inconsistent due to the limitations of the halfspace cooling model deployed by van Keken et al. (2008).
To illustrate the applicability of our results to the variety of subduction zones observed on Earth, we also run two endmember models with constant and temperaturedependent thermal parameters for a model with a younger (t_{s}=20 Myr) and older (t_{s}=80 Myr) slab age compared to our reference slab age of t_{s}=50 Myr.
2.7 Model diagnostics
To assess our models and quantify their differences, we use the three diagnostics defined in the community benchmark by van Keken et al. (2008), as well as the maximum depth of certain isotherms and the surface heat flux. Following van Keken et al. (2008), we define a uniform rectangular grid of 111×110 points with 6 km spacing starting in the top left corner and stored rowwise. On this grid, we interpolate the discrete temperature field T_{ij} in ^{∘}C in a postprocessing step. Using this grid, we output (1) the temperature T_{x=60 km} at the top of the slab at y=540 km and x=60 km, just downdip of the mantle wedge corner; (2) the l^{2} norm of the temperature along the top of the slab T_{slab} between y=600 km and y=390 km, as defined by
and (3) the l^{2} norm of the temperature in the tip of the mantle wedge between 54 and 120 km depth, as defined by
In addition to the diagnostics previously used in van Keken et al. (2008), we further report additional diagnostics that relate more closely to changes in the thermal structure of the slab that impact other processes, particularly the main processes governing seismogenesis. We report the maximum depth of the 350 and 450 ^{∘}C isotherms within the slab, which are associated with the brittle–ductile transition and hence the downdip limit of the seismogenic zone of megathrust seismicity (Hyndman et al., 1997; Gutscher and Peacock, 2003). We also report the maximum depth of the 600 ^{∘}C isotherm in the slab, which, together with the 350 and 450 ^{∘}C isotherms, is associated with the main dehydration reaction fronts within the slab, and the associated intermediatedepth seismicity between 70 and 350 km depth (Peacock, 2001; Yamasaki and Seno, 2003; Kelemen and Hirth, 2007). In addition, we report the values of the surface heat flux in our models since the observed surface heat flux is a frequently used constraint for models. Besides that, we provide snapshots of relevant variables, such as the temperature, viscosity, and velocity.
Postprocessing and visualisation are primarily done using MATLAB scripts (available at https://doi.org/10.5281/zenodo.7957674; Van Zelst et al., 2023) with additional touchups in Adobe Illustrator. We use scientific colour maps by Crameri (2018a) and Crameri et al. (2020) to avoid visual distortion of the data and exclusion of readers with colourvision deficiencies (Crameri, 2018b). To compare the thermal parameters and initial temperature conditions of the different models, we colour the models according to the optimal qualitative colour palette by Tsitsulin (2019).
3.1 Models with constant thermal parameters
The results from the reference model case2c_PvK with constant thermal parameters are shown in Fig. 4. It shows a subducting plate with a relatively cold core and a cold overriding plate with the base of the overriding plate that spills into the mantle wedge. There is flow in the mantle wedge around the base of the overriding plate which reaches the tip of the mantle wedge at x=50 km and y=550 km and heats up the subducting plate from the top.
The reference model has a combined dislocation and diffusion creep rheology in contrast to the original cases presented in van Keken et al. (2008), which are either isoviscous (Figs. S1–S4 in the Supplement), purely diffusion creep (Fig. S5 in the Supplement), or purely dislocation creep (Fig. S6 in the Supplement). Despite the difference in rheology, the model diagnostics of our reference model do not change significantly with respect to the model with a pure dislocation or diffusion creep rheology presented in van Keken et al. (2008) (Fig. S7 in the Supplement). However, looking at the snapshots presented in Fig. 4 and comparing them to the benchmark models of van Keken et al. (2008) (Figs. S5, S6 in the Supplement), there are distinct differences between our reference model and the benchmark cases presented in van Keken et al. (2008) in terms of the viscosity and velocity field in the mantle wedge, as well as the temperature field within the slab. These differences are not evident from our quantitative model diagnostics as the differences manifest themselves at high temperatures in the mantle wedge. These high temperatures and the region of the mantle wedge are not included in our model diagnostics as they principally affect the area of the model domain outside the main focus of our study, i.e. the slab.
^{∗} See Table S1 in the Supplement for model diagnostics of the models including an oceaniccrust approximation in the slab without an upper plate and including an oceanic upper plate.
In model case2c_bc, we build upon our reference model and change the initial and boundary temperature condition of the subducting oceanic plate at the left of the model from a halfspace cooling model to the plate model. This does not incur major changes in the model diagnostics (Table 2), which is consistent with the similarity between the temperature profiles of the halfspace cooling model and the plate model with constant values for the thermal parameters (Fig. 1b).
3.2 Models with temperaturedependent thermal conductivity
Using the temperaturedependent thermal conductivity according to Hofmeister (1999) and McKenzie et al. (2005) in model case2c_k1 results in an overall colder model with the slab isotherms penetrating deeper into the mantle. This effect increases with temperature, with the 350 ^{∘}C isotherm reaching 20 km deeper into the mantle and the 600 ^{∘}C isotherm reaching almost 90 km deeper into the mantle compared to the reference model (Fig. 7). We observe a similar but less pronounced trend when we use the thermal conductivity by Xu et al. (2004).
3.3 Models with temperaturedependent heat capacity
When using a temperaturedependent heat capacity, the model diagnostics show larger temperatures in the mantle wedge compared to the reference model with a constant heat capacity value. Similarly, the subducting slab is warmer, and isotherms penetrate less deeply into the mantle. For our preferred heat capacity model with 89 % forsterite and values from Berman and Aranovich (1996), the temperature diagnostics in the mantle wedge are larger by up to 37.7 ^{∘}C, and the maximum depths reached by the isotherms in the slab are shallower by 13.7–50 km (Fig. 7).
Using the values of Berman (1988) instead of the updated values of Berman and Aranovich (1996) only incurs minor changes in the model diagnostics (a maximum of 0.9 ^{∘}C and 1.3 km). Changing the ratio of forsterite and fayalite to 100 % forsterite in model case2c_Cp4 results in a slightly warmer mantle wedge by up to 2.8 ^{∘}C and shallower slab isotherms identical to the isotherm depths obtained in model case2c_Cp3 with values from Berman (1988) (Fig. 7). In the purely fayalite model, case2c_Cp5, the heat capacity is lower, resulting in a model that is cooler than the model with 89 % forsterite (case2c_Cp6) but that is still warmer than the reference model with a constant value for the heat capacity. Disregarding the temperaturedependent aspect of heat capacity tested, the overall magnitudes of the heat capacity used in the four C_{p} models from Fig. 2b also differ. For example, the purefayalite heat capacity model has the lowest overall heat capacity. This trend in changing magnitude of the heat capacity is also consistently visible in the model results, with models with lower heat capacity exhibiting lower temperatures and models with higher heat capacity resulting in higher temperature diagnostics. However, it is not straightforward to include the model with constant heat capacity values in this trend as well. For example, model case2c_Cp5 with fayalite values consistently has a lower heat capacity than the reference model with constant values, but the overall model diagnostics still show larger temperatures like the models with both larger and smaller heat capacity magnitudes depending on the temperature. Hence, the temperature dependence of the heat capacity has nonlinear effects on the resulting temperature field.
3.4 Models with temperaturedependent density
When we use a temperaturedependent density in model case2c_rho, the model is cooler than the reference model, case2c_PvK, but the effect is less pronounced than for the thermal conductivity (Table 2; Fig. 7). This results in isotherms that reach deeper into the mantle.
3.5 Models including all three temperaturedependent thermal parameters
We show the results for the model case2c_all in Fig. 5. In this model, we include the temperaturedependent function for the thermal conductivity by Hofmeister (1999) and McKenzie et al. (2005), the function of the heat capacity for 89 % forsterite with values from Berman and Aranovich (1996), and the temperaturedependent density. We show the differences between this model and the reference model (Fig. 4) in Fig. 6a–b for easy comparison. Based on our model diagnostics (Table 2), the model is overall colder than the reference model, and the slab has a colder core with isotherms that reach deeper into the mantle when we use temperaturedependent expressions for all thermal parameters (Fig. 7). However, the effect is less pronounced than for the models where we only use a temperaturedependent expression for the thermal conductivity. This is likely because the warming effect of the temperaturedependent heat capacity cancels part of the cooling effect of using temperaturedependent thermal conductivity and density. The largest difference between the two models is the temperature structure of the overriding plate, which is colder in the temperaturedependent model. Although we specifically focus on the change in thermal structure in the slab in this work, the difference in temperature in the overriding plate due to the use of temperaturedependent thermal parameters instead of constant values affects the thermal structure of the slab. Since the overriding plate is colder in models including temperature dependence of the thermal parameters, the heating of the interface between the slab and the overriding plate is delayed, which likely plays an important role in timedependent models of the thermal evolution of slab dynamics.
Within the slab itself, Fig. 6a–b shows that the largest temperature differences are approximately −65 ^{∘}C in the shallow part. The top of the slab is colder in model case2c_all, allowing isotherms to reach deeper into the mantle. The difference in isotherm depth is 3.8 km for the 350 ^{∘}C isotherm, 10 km for the 450 ^{∘}C isotherm, and 28.8 km for the 600 ^{∘}C isotherm (Fig. 7). At the base of the lithosphere, the bottom of the slab is warmer by up to 35 ^{∘}C compared to the reference model. This difference in slab temperature between the two models is partly due to the difference in boundary condition as the plate cooling model in case2c_all is cooler than the halfspace cooling model of case2c_PvK at shallow depths due to increased conductivity at low temperatures and is warmer at larger depths due to the imposition of a lower thermal boundary condition in the plate model. This also depends on the age of the subducting slab (see Sect. 3.6).
To summarise the effect of using temperaturedependent thermal parameters for all our models with a plate age of 50 Myr, we plot the maximum depth of the 350, 450, and 600 ^{∘}C isotherms for each model in Fig. 7. With respect to the reference model with constant values, adding the temperaturedependent thermal conductivity by Hofmeister (1999) and McKenzie et al. (2005) results in the largest changes in isotherm depth, with the isotherms reaching greater depths. Using a temperaturedependent density also results in a colder top of the slab with deeper isotherms. In contrast, using a temperaturedependent heat capacity yields a warmer slab, with isotherms penetrating the mantle less deeply than the reference model. Combining the effect of temperaturedependent thermal conductivity, heat capacity, and density results in an overall effect of slab cooling, with the isotherms reaching greater depths.
The surface heat flux in the models varies across the surface, with values of ∼17–49 mW m^{−2} observed in our models (Fig. S28 in the Supplement). This is in line with typical surface heat flux values found by Davies and Davies (2010). At the boundaries of the models, unrealistically high surface heat fluxes are obtained as a result of artificial boundary effects. The surface heat fluxes of all models generally show the same trend and similar values. One exception to this is the heat flux observed in case2c_all, which is 10 mW m^{−2} larger near the trench than the heat flux observed in case2c_PvK and case2c_bc. Considering that typical surface heat flux values are in the range of tens of mW m^{−2} (Davies and Davies, 2010), this is a significant change in modelpredicted surface heat flux. So, our results show that including temperaturedependent thermal parameters increases the predicted surface heat flux.
3.6 Models with different slab ages
Similarly to the models with a plate age of 50 Myr, we see a cooling effect in the temperaturedependent thermal models for the models with different plate ages (Fig. 8), with the changes to the thermal structure of the slab being more pronounced with increasing slab age. Similarly, from Fig. 9, we see that there is a particularly strong trend when it comes to larger isotherms such as the 600 ^{∘}C isotherm, with the differences between the reference models including constant thermal parameters and the models with variable properties increasing when the plate gets older. Hence, including temperaturedependent thermal parameters has a larger effect for old, cold subducting plates. This is expected, as the functions used in this paper for temperaturedependent thermal properties (Fig. 2) have their most extreme values at lower temperatures, which are more prevalent in older, and hence colder, slabs.
3.7 Models with a parameterised oceanic crust in the slab
When we include a parameterisation for a crustal layer at the top of the oceanic plate in the model, the diagnostics show a warmer top of the slab and mantle wedge. Within the slab, temperatures are also warmer, resulting in shallower maximum depths of the isotherms within the slab (Table S1 in the Supplement; Fig. S29 in the Supplement; results are similar to the results of the overriding plate models shown in Fig. 7 – see next section).
The difference between isotherm depths for models with and without the oceaniccrust parameterisation becomes more pronounced with increasing temperature. So, for example, between models case2c_bc and case2c_bc_mc, the difference in depth of the 350 ^{∘}C isotherm is 16.3 km, whereas it is 47.5 km for the depth of the 600 ^{∘}C isotherm. Even though the isotherm depths are shallower and the wedge is warmer compared to the models without any oceaniccrust parameterisation in the slab, the observed trends concerning the effects of temperaturedependent thermal parameters are the same as for the models without an oceaniccrust parameterisation in the slab (Fig. S29 in the Supplement).
An exception to these trends is the reference model case2c_PvK_mc, which uses the halfspace cooling model as a lefthandside temperature boundary condition that does not account for the 7 km oceaniccrust layer in the slab. Using this simplified boundary condition actually results in a colder slab with deeper isotherms (Fig. S29 in the Supplement). The large differences between the model case2c_PvK_mc and the other models in the _mc model batch show that the choice of temperature boundary condition and the use of a consistent temperature boundary condition is crucial.
3.8 Models with a parameterised overriding plate
The models that include the parameterisation for an overriding plate also include the oceaniccrust parameterisation in the slab. The models including an oceanic or continental plate all show the same trends as the models that do not include an overriding plate parameterisation but do include the oceaniccrust parameterisation. Typically, the models including a continental overriding plate result in a warmer slab and shallower isotherms compared to the other models. Figure 7 shows the results of the models including a continental overriding plate, which are similar to the models with an oceanic overriding plate and the models without an overriding plate (Sect. 3.7; Fig. S29 in the Supplement). Indeed, the largest difference in isotherm depth is 8.75 km for the 600^{∘} isotherm between models case2c_k1_cp and case2c_k1_mc, but the results are predominantly similar enough such that they plot on top of each other (Fig. S29 in the Supplement). This indicates that the nature of the overriding plate, and indeed the inclusion of an overriding plate at all, does not significantly affect the temperature field in the subduction zone in our models and hence does not affect the conclusions of this work. Beyond the slab, the inclusion of a continental overriding plate results in a warmer base of the continental crust and a colder mantle wedge (Fig. 6c–d).
We find that using temperaturedependent thermal parameters does not significantly change the firstorder thermal structure of a subduction zone in our models, with the largescale features remaining the same. Hence, when considering largescale subduction dynamics, the use of temperaturedependent thermal parameters is likely not an important factor. However, our results clearly show that the temperaturedependent thermal parameters do affect the thermal structure of the slab on the order of tens of degrees and hence tens of kilometres in these simple models of subduction zones. On a similar scale, our results show that the temperature boundary condition on the lefthandside of the model influences the temperature field of the slab. On the other hand, the inclusion of an overriding plate does not significantly affect the temperature field of the modelled slab.
Using temperaturedependent thermal parameters also affects the predicted surface heat flux in our models. As surface heat flux is one of the principal observables used in constraining subduction zone thermal models, all regional subduction zone thermal models are therefore affected by the inclusion of temperaturedependent thermal parameters. Our models with different plate ages show that our conclusions are valid regardless of the slab age, but they still lack realism in terms of model geometry and the inclusion of many important processes relevant for the development of a realistic thermal structure of the slab.
The found effect of using temperaturedependent thermal parameters on the order of tens of degrees could be significant when assessing the detailed thermal structures of local slabs for seismological applications. In this section, we first discuss the implications of our results for modelling the thermal structure of subduction zones in light of megathrust, intermediatedepth, and deep seismicity while taking into account the simple and geometrically unrealistic nature of our models (Sect. 4.1). As our models are conceptual calculations for the impact of including temperaturedependent variables, these implications are generic rather than applicable directly to any given subduction zone (Van Zelst et al., 2022). We further discuss the potential implications of our thermal models for the geochemical and mineralogical evolution of the slab and the impact this may have on the flux of fluids through subduction zones (Sect. 4.2). We then discuss how realistic our models are, their limitations, and how future work could further improve both the models and their applicability (Sect. 4.3).
4.1 Implications for seismicity
The temperature structure of a slab determines to a large extent where seismicity is expected to occur through its effect on both the mode of failure and the onset of ductile behaviour and its control on geochemical transitions within the slab and along the megathrust interface, including dehydration reactions. Here we summarise those effects and highlight how the results presented in Sect. 3 translate to influences on the distribution and extent of intermediatedepth and deepfocus earthquakes and potentially on the extent of megathrust and related shallow seismicity.
4.1.1 Intermediatedepth seismicity
Although the shallow slab geometry in our model is clearly a simplification, at depths consistent with intermediatedepth seismicity, the slab dip of 45^{∘} in our models is realistic, with an average slab dip of 45.5^{∘} reported by Syracuse et al. (2010) in nature, although it remains highly variable between different subduction zones. Other studies also find that the slab dip is steeper away from the interface between the slab and the overriding plate (e.g. Jarrard, 1986; King, 2001; Hu and Gurnis, 2020). Therefore, we can make some inferences about the expected depth of intermediatedepth seismicity using our models. Intermediatedepth seismicity at depths of 75–300 km is commonly associated with a temperature range between 600 and 800 ^{∘}C, where dehydration embrittlement of the metamorphosed minerals in the slab occurs (e.g. Jung et al., 2004; Wang et al., 2017). Focusing on the 600 ^{∘}C isotherm in our models (Table 2; Fig. 7), we see that its depth changes significantly throughout our parameter space with a depth of 203.8 km for the reference model case2c_PvK, 232.5 km for model case2c_all, and endmembers of 291.3 km depth for model case2c_k1 and 152.5 km depth for models case2c_Cp3 and case2c_Cp4. Hence, the depth at which dehydration reactions are expected to occur varies by almost 140 km within our parameter space. When comparing the reference model, case2c_PvK, and case2c_bc with constant thermal properties to the most complex model, case2c_all, which includes temperature dependence of all thermal properties, the differences in depth are smaller (∼28.8 km). However, this is still a significant difference in predicted seismicity depth that could be important when comparing to data. Therefore, the predicted depth of intermediatedepth seismicity in thermal models of subduction should be viewed in light of the assumptions regarding the thermal parameters.
Additionally, previous thermal models (e.g. Syracuse et al., 2010; van Keken et al., 2012) that use constant values for the thermal parameters and reproduce a thermal structure that fits observed seismicity are expected to change when temperaturedependent thermal parameters are used, with implications for the thermochemical changes then ascribed to control intermediatedepth seismicity. Depending on the choices of the functions describing the thermal parameters and their interaction, the fit with observed seismicity can change. To accurately determine the depth of intermediatedepth seismicity and the relationship between the thermal structure of the slab and intermediatedepth seismicity, we recommend the use of temperaturedependent thermal parameters constrained by the insights into rock behaviour. Neglecting temperaturedependent thermal parameters could result in errors of tens of kilometres in the estimated depth of intermediatedepth seismicity or a misinterpretation of the relation between the thermal structure of a slab and observed intermediatedepth seismicity. When only one of the thermal parameters is temperature dependent (i.e. just the thermal conductivity or just the heat capacity) instead of all three, the differences in temperature profile are even larger as their separate effects cancel each other out to an extent when considered simultaneously. Hence, to accurately model the thermal structure of the subduction zone, especially at intermediate depth, the temperature dependence of all three thermal parameters should be taken into account simultaneously.
4.1.2 Deep seismicity
The cause of deep earthquakes (> 300 km) is hotly debated, with proposed mechanisms such as dehydration embrittlement, transformational faulting, and (grainsizeassisted) thermal runway as a result of shear heating (see Green and Houston, 1995; Frohlich, 2006; Zhan, 2020, for an overview). Regardless of the exact mechanism responsible for deep earthquakes, it is clear that the thermal structure of the slab plays a large role by providing the optimal conditions for each of these mechanisms to occur. In fact, recent studies by Jia et al. (2020) and Liu et al. (2021) have shown that local slab temperature likely controls the rupture of deep earthquakes. Our results show that the effect of using temperaturedependent thermal parameters instead of constant values grows more pronounced with depth. Therefore, we expect that adding temperaturedependent thermal parameters will significantly affect models of the thermal structure of slabs at depths between 300–600 km, relevant to deep earthquakes.
4.1.3 Megathrust seismicity
The spatial extent of megathrust seismicity and, in particular, the maximum potential earthquake magnitude for interface events in a subduction zone are determined by the size of the seismogenic zone (and hence the maximum rupture width) bounded by empirically derived updip and downdip limits (e.g. Heuret et al., 2011). The updip limit of the seismogenic zone is usually thought to be determined by the transition from velocitystrengthening behaviour in the shallowest part of the subduction zone to velocityweakening behaviour on the megathrust. However, it has been observed that large earthquakes such as the 2011 TohokuOki earthquake can rupture through velocitystrengthening materials, thereby increasing the seismogenic zone size (Fujiwara et al., 2011; Lay et al., 2011). Similarly, tsunami earthquakes are thought to rupture the shallowest part of the subduction zone (Kanamori, 1972; Satake and Tanioka, 1999; Lay et al., 2012; Satake, 2015).
The downdip limit of the seismogenic zone is typically associated with the transition from brittle behaviour at lower temperatures to the onset of ductile behaviour at higher temperatures. The exact isotherms corresponding to this change in deformation style are still debated, with estimates ranging from 250 to 550 ^{∘}C depending on the mineralogy (Tichelaar and Ruff, 1993; Peacock and Hyndman, 1999; Scholz, 2019). Most commonly though, the downdip limit of the seismogenic zone is associated with the 350 and 450 ^{∘}C isotherms (Hyndman and Wang, 1993; Hyndman et al., 1997; Gutscher and Peacock, 2003).
Using the constant slab dip of 45^{∘} in our model setup, we can then make a quick and rough calculation of the seismogenic zone width in our models if we assume that the upper limit of the seismogenic zone is the surface. For the reference model case2c_PvK, we observe depths of the 350 and 450 ^{∘}C isotherms of 77.5 and 110 km, respectively, which translates to a seismogenic zone size of 109.6–155.6 km. These rough estimates for the seismogenic zone size of our models are within the observed range (i.e. ∼75–250 km ; Heuret et al., 2011). In the model where we take all temperaturedependent parameters into account, model case2c_all, the isotherm depths change by 3.8 and 10.0 km for 350 and 450 ^{∘}C, respectively, resulting in a maximum change in seismogenic zone size of 14.1 km.
However, our models here are of limited use in assessing the sensitivity of the temperature along the shallow subduction interface to the inclusion of temperaturedependent thermal properties for several reasons. First of all, in our simplified model geometry, the shallow dip of our interface is significantly larger than that typically seen in the interface seismogenic zone of most subduction zones (typically 23±8^{∘}; e.g. Jarrard, 1986; Heuret et al., 2011; Schellart and Rawlinson, 2013). This translates to a typical depth of the downdip limit of the seismogenic zone of 40±5 km (Tichelaar and Ruff, 1993). Besides that, we refrain from including the compositional complexity necessary to appropriately model the thermal structure of the overriding plate and a sedimentary forearc. Lastly, we do not include the effects of shear heating and fluid circulation on the shallow interface (England, 2018).
However, noting the impact that the variation in thermal properties at low temperatures (e.g. Fig. 3) has on the rates at which cold material heats up near the top of the downgoing plate and in the wedge of the forearc, we recommend using temperaturedependent thermal parameters in thermal models of subduction zones, in addition to the other influences mentioned, for when physically realistic estimations of the seismogenic zone size are desired. Similarly, when observations are linked to the behaviour of the interface (e.g. limits on seismogenesis, on coupling, on slow slip), the inclusion of temperaturedependent thermal parameters may alter the inferred mineralogical and rheological controls on such transitions. Considering that we use an unrealistically steep subduction angle, our results on the changes in the maximum depth of the 350 and 450 ^{∘}C isotherms might underestimate the effect on the seismogenic zone size in realistic settings with lower megathrust angles.
4.2 Mineralogical evolution of the slab
As the subducting plate descends, it typically undergoes a range of mineralogical transitions relating to the increase in pressure and temperature. These mineralogical changes, particularly the location at which dehydration reactions release fluids into the slab system, play a controlling role in determining the location of intraslab seismicity and also in influencing a range of other geophysical phenomena, from the internal impedance and velocity contrasts within the slab (e.g. Abers, 2000; Rondenay et al., 2008) to the occurrence of slow slip events on the subduction interface (e.g. Peacock, 2009) to the development of a hydrated mantle forearc (e.g. Abers et al., 2017). The preservation of volatilehosting, lowertemperature material into the deeper mantle also plays a role in global geochemical cycles (e.g. Rüpke et al., 2004).
Whilst the kinematic constraints we impose on the slab in our models mean there is little variation in lithostatic pressure between models, we have shown that including the temperature dependence of thermal parameters in the modelling of slab thermal structures has an impact on the resultant temperature field. Whilst these changes are small relative to the total change in temperature experienced by the slab during subduction, they lead to a slightly different pressure–temperature evolution for the slab material. An additional crustal layer in our models, for which we currently only employ a simple parameterisation, further alters the temperature field. We note that the changes in the temperature evolution of the uppermost ∼7 km of the slab is particularly susceptible to the temperature dependence of thermal properties given the rapid variation of such values at low temperatures (Fig. 3). The mineralogical evolution of the shallowest part of the slab is therefore likely to be altered by the incorporation of temperaturedependent thermal properties, with initially more rapid heating at low pressures giving way to slower heating at higher pressures, in comparison to models using fixed, temperatureindependent thermal properties. Dehydration reactions in hydrated basaltic oceanic crust typically take place between 350–450 ^{∘}C, whilst those in the serpentinised oceanic mantle concentrate between 600 and 800 ^{∘}C (Hacker et al., 2003a). In linking geophysical observations to thermal models, we again note that the variation in depth of the 350 and 450 ^{∘}C isotherms in our models of up to 38.8 km with respect to the reference model case2c_PvK (Fig. 7) would translate for most subduction zones into a significant trenchperpendicular lateral shift. To a lesser extent, the same is true when merely considering the depth variation of ∼3.8–10 km for the 350 and 450 ^{∘}C isotherms between case2c_PvK and case2c_bc with constant thermal parameters and case2c_all with temperaturedependent thermal parameters, which represent our most simple and most complex models. This will have a significant impact on the source location of phenomena such as the migration of fluids from the slab to the forearc mantle and/or updip along the subduction interface.
Lastly, the model diagnostics we focus on here centre around the maximum depth of a given isotherm. However, the variation in thermal structure that we study will also impact on the thermal cross section of the slab at any given depth – with marginally colder slabs having a significantly greater volume of material below a given temperature at a given depth and hence potentially altering the volatile fluxes within slabs into the middle mantle.
4.3 Model limitations and future work
With the exception of a different rheology in the mantle wedge, where we combine both diffusion and dislocation creep, we use the same model setup as the subduction zone community benchmark presented by van Keken et al. (2008). We choose this model setup as it is well defined and documented and reproduced by many codes in the geodynamics community (see codes used in van Keken et al., 2008). Hence, we are able to study the effect of temperaturedependent thermal parameters on the thermal structure of subduction zones in an isolated, welldefined manner, although, as discussed, this does limit the direct applicably to observational data.
The model setup is greatly simplified, and many complexities that are known to influence the thermal structure of the slab are ignored. As illustrated in the benchmark of van Keken et al. (2008) itself, one of the largest influences on the thermal structure of the subducting slab is the employed rheology. The temperature model diagnostics in van Keken et al. (2008) change up to 189 ^{∘}C when changing from an isoviscous to a dislocation or diffusion creep rheology. To a lesser extent, the difference between a purely dislocation creep and diffusion creep rheology is noticeable with variations on the order of 10 ^{∘}C in model diagnostics. We find that employing a combined dislocation and diffusion creep rheology does not significantly change the model diagnostics compared to a purely dislocation or diffusion creep rheology. However, our approximation of combining a dislocation and diffusion creep rheology is simplistic. Using a composite rheology of diffusion and dislocation creep to properly account for the nonlinearity of the two rheologies would be more physically appropriate (Ranalli, 1995; Karato, 2008; Gerya, 2019). This would likely introduce changes to the temperature field of the slab on the same order as the differences observed between a pure diffusion and a pure dislocation model, as in van Keken et al. (2008). Similarly, it would be more appropriate to take into account the effects of the activation volume, which we currently set to zero. The activation volume is not well constrained (Karato and Wu, 1993; Dixon and Durham, 2018) but would affect the viscosity of the system and therefore indirectly change the flow patterns and thermal structure of the subduction zone.
Hence, the effect of using temperaturedependent thermal parameters in thermal models instead of constant values is a secondary effect to rheology when comparing isoviscous and nonlinear rheologies (i.e. compare Figs. S1–S4 to S5–S7 in the Supplement). However, when comparing nonlinear rheologies, using temperaturedependent thermal parameters instead of constant values will likely have a greater effect on the thermal structure of the slab than changing the details of the rheology formulation. Note that these conclusions relate to the thermal structure of the slab; the rheology plays a major role in the thermal structure of the mantle wedge and overriding plate, as is evident from the original benchmarks presented in van Keken et al. (2008). Although our parameterisation of the overriding plate captures some of these complexities, our results cannot address all the complexities introduced by a rheology tailored to continentalcrust rocks in the overriding plate. Similarly, our overriding plate parameterisation does not capture the effect of crustal radiogenic heat production on the temperature field. In the absence of this, based on our results, we predict that changes in the model with regards to the overriding plate will not significantly affect the temperature field of the slab, although the overriding plate thickness would change the temperature depth distribution.
Apart from a simplified rheology, we also employ a simplified geometry in our model setup. Although the model serves as a good benchmark and although we can infer some implications for seismicity from this simple setup, a strictly 45^{∘} dipping slab is not realistic. In nature, the slab dip changes with depth, with low dipping angles of 23±8^{∘} for the megathrust region (Heuret et al., 2011) and larger dip at depth (e.g. Isacks and Barazangi, 1977; King, 2001; Cruciani et al., 2005; Syracuse et al., 2010; Klemd et al., 2011; Hu and Gurnis, 2020). Therefore, more realistic models of the thermal structure of subduction zones include these complex geometries (e.g. Syracuse et al., 2010; van Keken et al., 2012). Our results indicate that, in these complex models of the thermal structure of the slab, it is important to take the temperature dependence of thermal parameters into account as well. Even though including them will likely not change the largescale subduction evolution, it is important to include the temperaturedependent thermal parameters for accurate comparison with (earthquake) data.
Although we focus here on the effect of using temperaturedependent thermal parameters, there are numerous other processes relevant to the developing thermal structure of a subduction zone (see van Keken et al., 2019, for an overview). For example, frictional (or shear) heating along the plate interface (e.g. Peacock, 1992; England and Molnar, 1993; Peacock, 1993b; Peacock and Wang, 1999; Burg and Gerya, 2005; Gao and Wang, 2014, 2017) and radiogenic heating in the overriding plate (e.g. Gao and Wang, 2014; England, 2018) introduce additional heat sources to the system and result in warmer slabs in line with petrological estimates of the temperatures of rocks in a subduction zone (PennistonDorland et al., 2015). Typically these processes are included in models where a temperaturedependent density formulation is used, although the conductivity and heat capacity are often still taken to be constants. We deliberately do not include these additional heat sources when including the temperaturedependent density to isolate its effect on the thermal structure of a subduction zone. However, we recognise that this may lead to thermodynamic inconsistencies similar to those introduced through inconsistent thermodynamic potentials calculated from the thermal parameters (Schubert et al., 2001; Van Zelst et al., 2022). Phase changes, such as serpentinisation in the mantle wedge corner (e.g. Hyndman and Peacock, 2003) and the transition from blueschist to hydrous eclogite (e.g. Hacker et al., 2003a), also play an important role in establishing the thermal structure of the slab as they are paired with the release of latent heat, density and subsequent volume changes, fluid production, and heat advection (see Peacock, 2020, for an overview of petrologic models). Fluid flow and hydrothermal circulation within the upper part of the slab efficiently transport heat updip towards the trench (e.g. Spinelli and Wang, 2008; England and Katz, 2010; Faccenda et al., 2012; Rotman and Spinelli, 2013; Harris et al., 2017). Depending on the subduction velocity, this can significantly reduce the temperature of the subduction interface (Rotman and Spinelli, 2013). In line with this, processes such as melting and melt transport at the top of the slab and in the mantle wedge corner (e.g. England and Katz, 2010; Bouilhol et al., 2015; Perrin et al., 2016), magmatism (e.g. Jones et al., 2018), erosion (e.g. Royden, 1993; England, 2018), sedimentation (e.g. England, 2018), anisotropy (e.g. Morishige and Tasaka, 2021), and threedimensional complexities (e.g. Gerya, 2011; Plunder et al., 2018; Wada, 2021; Qu et al., 2023) also play a role in the thermal structure of a subduction zone. In addition, subduction is an inherently timedependent process, with the temperature structure of the subducting slab likely changing throughout its evolution, which is not captured by the steadystate thermal models presented here (King, 2001; Holt and Condit, 2021). Here, we deliberately choose to ignore these complexities to isolate the effect of temperaturedependent thermal parameters on the thermal structure of the slab. Future studies could focus on these processes to explore their effect on models of the thermal structure of the slab.
Although we consider a set of models where we include a simple parameterisation as an approximation of a crustal layer in the slab, our models are still restricted to a single composition. Hence, we do not explicitly include a crustal layer, and we neglect the impact of the mineralogical evolution of the slab on the temperature structure, both through the variation in thermal parameters with evolving mineralogy and through the latent heat of mineralogical transformation. Our results suggest that including compositional heterogeneity, specifically oceanic crust, does not change the observed trends concerning the effect of temperaturedependent thermal parameters on the models. Hence, the main conclusions presented in this work are not affected by compositional heterogeneity. However, thermal parameters do vary greatly with composition (e.g. Whittington et al., 2009; Miao et al., 2014). Therefore, compositional variation plays a big role in controlling the thermal structure of the incoming plate (e.g. Richards et al., 2018) and therefore also the subduction zone dynamics (e.g. Gerya et al., 2004). Going forwards, the development of models that include both composition and temperature dependence in the thermal parameters is likely required to further progress the accurate modelling of the thermal structures of subduction zones.
Lastly, there are numerous functions describing the temperature and pressure dependence of the thermal parameters in the literature, and existing functions are continuously updated with improved values for constants to better fit laboratory data. It is outside of the scope of this work to test all different formulations, and here we follow McKenzie et al. (2005) and Richards et al. (2018), who used temperaturedependent thermal parameters for plate models of the cooling oceanic lithosphere. However, other possible functions of the temperature dependence of thermal parameters include formulations from studies like e.g. Berman and Brown (1985), Seipold (1998), Hofmeister (2007a), Wen et al. (2015), and Su et al. (2018). In addition, pressuredependent formulations have been proposed by e.g. Hofmeister (2007b), and studies have shown that the residual misfit between the model and the data is reduced upon including the pressure dependence of the thermal parameters (e.g. Grose and Afonso, 2013; Korenaga and Korenaga, 2016; Richards et al., 2018).
In this work, we look at the effect of using temperaturedependent thermal parameters in thermal models of subduction zones compared to using constant values.
We find that the use of temperaturedependent thermal parameters does not significantly affect the firstorder thermal structure of a subducting slab, like the choice of rheology or plate age would. However, the local thermal structure of the slab does change on the order of tens of degrees when using temperaturedependent thermal parameters. This secondorder effect could be significant for specific modelling applications that consider, for instance, seismicity and phase changes.
Using only a temperaturedependent conductivity decreases the temperature in the slab and results in a larger predicted seismogenic zone width and deeper intermediatedepth seismicity, with the maximum depth of the 600 ^{∘}C isotherm changing up to 87.5 km for a model using the thermal conductivity function of Hofmeister (1999) and McKenzie et al. (2005) compared to a reference model with constant values.
Solely employing a temperaturedependent heat capacity has the opposite effect and results in a warmer slab with a shallower downdip limit of the seismogenic zone and predicted depth of dehydration reactions responsible for intermediatedepth seismicity.
Using only a temperaturedependent density has the least effect on the thermal structure of the slab when compared to the reference model with constant values, although the slab is overall colder.
Combining the temperature dependence of the three thermal parameters negates the effect on the thermal structure of the slab slightly, but the strong cooling of the slab produced by both the temperaturedependent thermal conductivity and density dominates. Therefore, the modelled slab is colder than a slab modelled with constant thermal parameters with, e.g. the maximum depth of the 600 ^{∘}C isotherm changing by 28.8 km. In addition, the predicted surface heat flux increases when temperaturedependent thermal parameters are included. The importance of including temperaturedependent thermal parameters increases for increasing slab age as the functions of the thermal parameters used in this study have their most extreme values for lower temperatures.
Models including a parameterisation of oceanic crust in the slab and an oceanic or overriding plate show the same trends. We find that choosing consistent temperature boundary conditions is crucial and can otherwise lead to large differences in the temperature field of the slab. In contrast, the nature, or existence, of an overriding plate does not significantly affect the temperature field of the slab.
We caution that the conclusions drawn here stem from a highly simplified model of a subduction zone that is not suitable for direct comparisons to nature and does not aim to reproduce observations. Hence, the details of our findings will likely change for more complex and realistic model setups. However, even considering the simplifications in our model setup, our results indicate that the changes in the modelled thermal structure of the slab, as well as the predicted surface heat flux, will have important implications for the estimated size of the seismogenic zone in these kinds of thermal models and for predictions where intermediatedepth seismicity might occur. These implications for the inclusion of temperaturedependent thermal parameters are of a lower order than the established firstorder controls on the subduction zone thermal structure, such as the employed rheology and plate age. Nevertheless, we suggest that, for optimal comparison to data and to avoid misinterpretations, temperaturedependent thermal parameters are also an important modelling ingredient that should be taken into account when using thermal(mechanical) models of subduction zones, especially in studies that have seismological applications.
The models were run with the opensource code xFieldstone (https://github.com/irisvanzelst/xFieldstone; https://doi.org/10.5281/zenodo.7957374, Van Zelst, 2023). The data used to reproduce the van Keken et al. (2008) benchmark and the MATLAB scripts used for the postprocessing of the results and generation of the figures are shared in https://doi.org/10.5281/zenodo.7957674 (Van Zelst et al., 2023).
The supplement related to this article is available online at: https://doi.org/10.5194/se146832023supplement.
IvZ and TJC conceived the study. IvZ designed and ran the models, analysed the results, and wrote the article. CT and IvZ wrote the code xFieldstone. TJC supervised IvZ and contributed to the analysis of the models. All authors discussed the results and contributed to the model design and the final article. The author order of the second and third authors was decided by a coin flip.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We thank the editor Taras Gerya and three anonymous reviewers at Solid Earth as well as three previous anonymous reviewers at JGR: Solid Earth for their constructive feedback that helped to improve this article.
We thank Peter van Keken for providing the original data from van Keken et al. (2008) to benchmark the code presented here and for providing additional insights on the model setup.
Elements of this work were undertaken on ARC4, part of the HighPerformance Computing facilities at the University of Leeds, UK. We also ran simulations on the Plejades work stations of the German Aerospace Center (DLR), Germany. Iris van Zelst and Timothy J. Craig were funded by the Royal Society (UK) through University Research Fellowship no. URF\R1\180088 and Research Fellows Enhancement Award no. RGF\EA\181084 and CEC19\100101. Timothy J. Craig was also supported through COMET, the UK Natural Environment Research Council's Centre for the Observation and Modelling of Earthquakes, Volcanoes, and Tectonics. Iris van Zelst additionally received financial support and endorsement from the DLR Management Board Young Research Group Leader Program and from the Executive Board Member for Space Research and Technology.
This research has been supported by the Royal Society (grant nos. URF\R1\180088, RGF\EA\181084, and CEC19\100101), COMET (the UK Natural Environment Research Council's Centre for the Observation and Modelling of Earthquakes, Volcanoes, and Tectonics), the DLR Management Board Young Research Group Leader Program and the Executive Board Member for Space Research and Technology.
The article processing charges for this openaccess publication were covered by the German Aerospace Center (DLR).
This paper was edited by Taras Gerya and reviewed by three anonymous referees.
Abers, G., van Keken, P., and Hacker, B.: The cold and relatively dry nature of mantle forearcs in subduction zones, Nat. Geosci., 10, 333–337, 2017. a, b
Abers, G. A.: Hydrated subducted crust at 100–250 km depth, Earth Planet. Sc. Lett., 176, 323–330, 2000. a
Abers, G. A., van Keken, P. E., Kneller, E. A., Ferris, A., and Stachnik, J. C.: The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow, Earth Planet. Sc. Lett., 241, 387–397, 2006. a, b
Arcay, D.: Modelling the interplate domain in thermomechanical simulations of subduction: Critical effects of resolution and rheology, and consequences on wet mantle melting, Phys. Earth Planet. In., 269, 112–132, 2017. a
Beall, A., Fagereng, Å., Davies, J. H., Garel, F., and Davies, D. R.: Influence of subduction zone dynamics on interface shear stress and potential relationship with seismogenic behavior, Geochem. Geophy., Geosy., 22, e09267, https://doi.org/10.1029/2020GC009267, 2021. a
Berman, R. G.: Internallyconsistent thermodynamic data for minerals in the system Na_{2}OK_{2}OCaOMgOFeOFe_{2}O_{3}Al_{2}O_{3}SiO_{2}TiO_{2}H_{2}OCO_{2}, J. Petrol., 29, 445–522, 1988. a, b, c, d, e, f, g, h, i
Berman, R. G. and Aranovich, L. Y.: Optimized standard state and solution properties of minerals: I. Model calibration for olivine, orthopyroxene, cordierite, garnet, and ilmenite in the system FeOMgOCaOAl_{2}O_{3}TiO_{2}SiO_{2}, Contrib. Mineral. Petr., 126, 1–24, 1996. a, b, c, d, e, f, g, h, i, j, k, l, m, n
Berman, R. G. and Brown, T. H.: Heat capacity of minerals in the system Na_{2}OK_{2}OCaOMgOFeOFe_{2}O_{3}Al_{2}O_{3}SiO_{2}TiO_{2}H_{2}OCO_{2}: representation, estimation, and high temperature extrapolation, Contrib. Mineral. Petr., 89, 168–183, 1985. a, b
Blom, N., Boehm, C., and Fichtner, A.: Synthetic inversions for density using seismic and gravity data, Geophys. J. Int., 209, 1204–1220, 2017. a
Bouhifd, M. A., Andrault, D., Fiquet, G., and Richet, P.: Thermal expansion of forsterite up to the melting point, Geophys. Res. Lett., 23, 1143–1146, 1996. a, b, c, d, e, f, g
Bouilhol, P., Magni, V., van Hunen, J., and Kaislaniemi, L.: A numerical approach to melting in warm subduction zones, Earth Planet. Sc. Lett., 411, 37–44, 2015. a
Brizzi, S., Van Zelst, I., Funiciello, F., Corbi, F., and van Dinther, Y.: How sediment thickness influences subduction dynamics and seismicity, J. Geophys. Res.Sol. Ea., 125, e2019JB018964, https://doi.org/10.1029/2019JB018964, 2020. a
Burg, J.P. and Gerya, T.: The role of viscous heating in Barrovian metamorphism of collisional orogens: thermomechanical models and application to the Lepontine Dome in the Central Alps, J. Metamorph. Geol., 23, 75–95, 2005. a
Chang, C., McNeill, L. C., Moore, J. C., Lin, W., Conin, M., and Yamada, Y.: In situ stress state in the Nankai accretionary wedge estimated from borehole wall failures, Geochem. Geophy. Geosy., 11, Q0AD04, https://doi.org/10.1029/2010GC003261, 2010. a
Chemia, Z., Dolejš, D., and SteinleNeumann, G.: Thermal effects of variable material properties and metamorphic reactions in a threecomponent subducting slab, J. Geophys. Res.Sol. Ea., 120, 6823–6845, 2015. a
Crameri, F.: Scientific colourmaps, Zenodo [code], https://doi.org/10.5281/zenodo.1243862, 2018a. a
Crameri, F.: Geodynamic diagnostics, scientific visualisation and StagLab 3.0, Geosci. Model Dev., 11, 2541–2562, https://doi.org/10.5194/gmd1125412018, 2018b. a
Crameri, F., Shephard, G. E., and Heron, P. J.: The misuse of colour in science communication, Nat. Commun., 11, 1–10, 2020. a
Crosby, A., McKenzie, D., and Sclater, J.: The relationship between depth, age and gravity in the oceans, Geophys. J. Int., 166, 553–573, 2006. a
Crouzeix, M. and Raviart, P.A.: Conforming and nonconforming finite element methods for solving the stationary Stokes equations I, Rev. Fr. Automat. Infor. Mathématique, 7, 33–75, 1973. a
Cruciani, C., Carminati, E., and Doglioni, C.: Slab dip vs. lithosphere age: no direct function, Earth Planet. Sc. Lett., 238, 298–310, 2005. a
Dabrowski, M., Krotkiewski, M., and Schmid, D.: MILAMIN: MATLABbased finite element method solver for large problems, Geochem. Geophy. Geosy., 9, Q04030, https://doi.org/10.1029/2007GC001719, 2008. a
Davies, J. H. and Davies, D. R.: Earth's surface heat flux, Solid Earth, 1, 5–24, https://doi.org/10.5194/se152010, 2010. a, b
Denlinger, R. P.: A revised estimate for the temperature structure of the oceanic lithosphere, J. Geophys. Res.Sol. Ea., 97, 7219–7222, 1992. a
Dixon, N. A. and Durham, W. B.: Measurement of activation volume for creep of dry olivine at uppermantle conditions, J. Geophys. Res.Sol. Ea., 123, 8459–8473, 2018. a
Emmerson, B. and McKenzie, D.: Thermal structure and seismicity of subducting lithosphere, Phys. Earth Planet. In., 163, 191–208, 2007. a, b
England, P.: On shear stresses, temperatures, and the maximum magnitudes of earthquakes at convergent plate boundaries, J. Geophys. Res.Sol. Ea., 123, 7165–7202, 2018. a, b, c, d
England, P. and Molnar, P.: The interpretation of inverted metamorphic isograds using simple physical calculations, Tectonics, 12, 145–157, 1993. a
England, P. C. and Katz, R. F.: Melting above the anhydrous solidus controls the location of volcanic arcs, Nature, 467, 700–703, 2010. a, b
Faccenda, M., Gerya, T. V., Mancktelow, N. S., and Moresi, L.: Fluid flow during slab unbending and dehydration: Implications for intermediatedepth seismicity, slab weakening and deep water recycling, Geochem. Geophy. Geosy., 13, Q01010, https://doi.org/10.1029/2011GC003860, 2012. a, b
Frohlich, C.: Deep earthquakes, Cambridge university press, ISBN 9781107297562, https://doi.org/10.1017/CBO9781107297562, 2006. a
Fujiwara, T., Kodaira, S., No, T., Kaiho, Y., Takahashi, N., and Kaneda, Y.: The 2011 TohokuOki earthquake: Displacement reaching the trench axis, Science, 334, 1240–1240, 2011. a
Fulton, P., Brodsky, E. E., Kano, Y., Mori, J., Chester, F., Ishikawa, T., Harris, R., Lin, W., Eguchi, N., Toczko, S., and Expedition 343, 343T, and KR1308 Scientists: Low coseismic friction on the TohokuOki fault determined from temperature measurements, Science, 342, 1214–1217, 2013. a
Gao, X. and Wang, K.: Strength of stickslip and creeping subduction megathrusts from heat flow observations, Science, 345, 1038–1041, 2014. a, b
Gao, X. and Wang, K.: Rheological separation of the megathrust seismogenic zone and episodic tremor and slip, Nature, 543, 416–419, 2017. a
Gerya, T.: Introduction to numerical geodynamic modelling, Cambridge University Press, ISBN 9781316534243, https://doi.org/10.1017/9781316534243, 2019. a
Gerya, T. V.: Future directions in subduction modeling, J. Geodynam., 52, 344–378, 2011. a, b
Gerya, T. V. and Meilick, F.: Geodynamic regimes of subduction under an active margin: effects of rheological weakening by fluids and melts, J. Metamorph. Geol., 29, 7–31, 2011. a
Gerya, T. V., Yuen, D. A., and Maresch, W. V.: Thermomechanical modelling of slab detachment, Earth Planet. Sc. Lett., 226, 101–116, 2004. a
Green, H. W. and Houston, H.: The mechanics of deep earthquakes, Annu. Rev. Earth Pl. Sc., 23, 169–213, 1995. a, b
Grose, C. J. and Afonso, J. C.: Comprehensive plate models for the thermal evolution of oceanic lithosphere, Geochem. Geophy. Geosy., 14, 3751–3778, 2013. a, b, c
Gutscher, M.A. and Peacock, S. M.: Thermal models of flat subduction and the rupture zone of great subduction earthquakes, J. Geophys. Res.Sol. Ea., 108, 2009, https://doi.org/10.1029/2001JB000787, 2003. a, b, c
Hacker, B. R., Abers, G. A., and Peacock, S. M.: Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H_{2}O contents, J. Geophys. Res.Sol. Ea., 108, 2029, https://doi.org/10.1029/2001JB001127, 2003a. a, b, c, d, e
Hacker, B. R., Peacock, S. M., Abers, G. A., and Holloway, S. D.: Subduction factory 2. Are intermediatedepth earthquakes in subducting slabs linked to metamorphic dehydration reactions?, J. Geophys. Res.Sol. Ea., 108, 2030, https://doi.org/10.1029/2001JB001129, 2003b. a
Harris, R., Yamano, M., Kinoshita, M., Spinelli, G., Hamamoto, H., and Ashi, J.: A synthesis of heat flow determinations and thermal modeling along the Nankai Trough, Japan, J. Geophys. Res.Sol. Ea., 118, 2687–2702, 2013. a
Harris, R. N., Spinelli, G. A., and Fisher, A. T.: Hydrothermal circulation and the thermal structure of shallow subduction zones, Geosphere, 13, 1425–1444, 2017. a
Heuret, A., Lallemand, S., Funiciello, F., Piromallo, C., and Faccenna, C.: Physical characteristics of subduction interface type seismogenic zones revisited, Geochem. Geophy. Geosy., 12, Q01004, https://doi.org/10.1029/2010GC003230, 2011. a, b, c, d
Hillier, J. and Watts, A.: Relationship between depth and age in the North Pacific Ocean, J. Geophys. Res.Sol. Ea., 110, B02405, https://doi.org/10.1029/2004JB003406, 2005. a
Hofmeister, A.: Mantle values of thermal conductivity and the geotherm from phonon lifetimes, Science, 283, 1699–1706, 1999. a, b, c, d, e, f, g, h, i, j, k, l, m
Hofmeister, A. M.: Thermal conductivity of the Earth's deepest mantle, in: Superplumes: Beyond Plate Tectonics, Springer, 269–292, https://doi.org/10.1007/9781402057502, 2007a. a
Hofmeister, A. M.: Pressure dependence of thermal transport properties, P. Natl. Acad. Sci. USA, 104, 9192–9197, 2007b. a
Holt, A. F. and Condit, C. B.: Slab temperature evolution over the lifetime of a subduction zone, Geochem. Geophy. Geosy., 22, e2020GC009476, https://doi.org/10.1029/2020GC009476, 2021. a
Hu, J. and Gurnis, M.: Subduction duration and slab dip, Geochem. Geophy. Geosy., 21, e2019GC008862, https://doi.org/10.1029/2019GC008862, 2020. a, b
Hyndman, R. D. and Peacock, S. M.: Serpentinization of the forearc mantle, Earth Planet. Sc. Lett., 212, 417–432, 2003. a, b
Hyndman, R. D. and Wang, K.: Thermal constraints on the zone of major thrust earthquake failure: The Cascadia subduction zone, J. Geophys. Res.Sol. Ea., 98, 2039–2060, 1993. a, b, c
Hyndman, R. D., Yamano, M., and Oleskevich, D. A.: The seismogenic zone of subduction thrust faults, Isl. Arc, 6, 244–260, 1997. a, b, c
Isacks, B. L. and Barazangi, M.: Geometry of Benioff zones: Lateral segmentation and downwards bending of the subducted lithosphere, in: Island Arcs, Deep Sea Trenches and BackArc Basins, edited by: Talwani, M. and Pitman III, W. C., American Geophysical Union, 1, 99–114, ISBN 9781118665756, https://doi.org/10.1029/ME001, 1977. a
Jarrard, R. D.: Relations among subduction parameters, Rev. Geophys., 24, 217–284, 1986. a, b
Jia, Z., Shen, Z., Zhan, Z., Li, C., Peng, Z., and Gurnis, M.: The 2018 Fiji Mw 8.2 and 7.9 deep earthquakes: One doublet in two slabs, Earth Planet. Sc. Lett., 531, 115997, https://doi.org/10.1016/j.epsl.2019.115997, 2020. a
Jones, D. W. R., Katz, R. F., Tian, M., and Rudge, J. F.: Thermal impact of magmatism in subduction zones, Earth Planet. Sc. Lett., 481, 73–79, 2018. a
Jung, H., Green Ii, H. W., and Dobrzhinetskaya, L. F.: Intermediatedepth earthquake faulting by dehydration embrittlement with negative volume change, Nature, 428, 545–549, 2004. a, b
Kanamori, H.: Mechanism of tsunami earthquakes, Phys. Earth Planet. In., 6, 346–359, https://doi.org/10.1016/00319201(72)900581, 1972. a
Karato, S.I.: Deformation of earth materials, vol. 463, Cambridge Press Cambridge, https://doi.org/10.1017/CBO9780511804892, 2008. a
Karato, S.I. and Wu, P.: Rheology of the upper mantle: A synthesis, Science, 260, 771–778, 1993. a
Kelemen, P. B. and Hirth, G.: A periodic shearheating mechanism for intermediatedepth earthquakes in the mantle, Nature, 446, 787–790, 2007. a
King, S. D.: Subduction zones: observations and geodynamic models, Phys. Earth Planet. In., 127, 9–24, 2001. a, b, c
Klemd, R., John, T., Scherer, E., Rondenay, S., and Gao, J.: Changes in dip of subducted slabs at depth: petrological and geochronological evidence from HP–UHP rocks (Tianshan, NWChina), Earth Planet. Sc. Lett., 310, 9–20, 2011. a
Korenaga, T. and Korenaga, J.: Evolution of young oceanic lithosphere and the meaning of seafloor subsidence rate, J. Geophys. Res.Sol. Ea., 121, 6315–6332, 2016. a
Lay, T., Ammon, C. J., Kanamori, H., Xue, L., and Kim, M. J.: Possible large neartrench slip during the 2011 Mw 9.0 off the Pacific coast of Tohoku Earthquake, Earth Planets Space, 63, 32, https://doi.org/10.5047/eps.2011.05.033, 2011. a
Lay, T., Kanamori, H., Ammon, C. J., Koper, K. D., Hutko, A. R., Ye, L., Yue, H., and Rushing, T. M.: Depthvarying rupture properties of subduction zone megathrust faults, J. Geophys. Res.Sol. Ea., 117, B04311, https://doi.org/10.1029/2011JB009133, 2012. a
Liu, H., Gurnis, M., Leng, W., Jia, Z., and Zhan, Z.: Tonga Slab Morphology and Stress Variations Controlled by a Relic Slab: Implications for Deep Earthquakes in the TongaFiji Region, Geophys. Res. Lett., 48, e2020GL091331, https://doi.org/10.1029/2020GL091331, 2021. a
McKenzie, D. and Sclater, J.: Heat flow in the eastern Pacific and sea floor spreading, Bulletin Volcanologique, 33, 101–117, 1969. a
McKenzie, D., Jackson, J., and Priestley, K.: Thermal structure of oceanic and continental lithosphere, Earth Planet. Sc. Lett., 233, 337–349, 2005. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab
Miao, S., Li, H., and Chen, G.: Temperature dependence of thermal diffusivity, specific heat capacity, and thermal conductivity for several types of rocks, J. Therm. Anal. Calorim., 115, 1057–1063, 2014. a
Morishige, M. and Tasaka, M.: Limited impact of anisotropic thermal conductivity in the mantle wedge on the slab temperature in the Tohoku subduction zone, Northeast Japan, Tectonophysics, 820, 229110, https://doi.org/10.1016/j.tecto.2021.229110, 2021. a
Parsons, B. and Sclater, J. G.: An analysis of the variation of ocean floor bathymetry and heat flow with age, J. Geophys. Res., 82, 803–827, 1977. a, b
Peacock, S. M.: Blueschistfacies metamorphism, shear heating, and PTt paths in subduction shear zones, J. Geophys. Res.Sol. Ea., 97, 17693–17707, 1992. a
Peacock, S. M.: Largescale hydration of the lithosphere above subducting slabs, Chem. Geol., 108, 49–59, 1993a. a
Peacock, S. M.: The importance of blueschist > eclogite dehydration reactions in subducting oceanic crust, Geol. Soc. Am. Bull., 105, 684–694, 1993b. a
Peacock, S. M.: Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle?, Geology, 29, 299–302, 2001. a, b, c
Peacock, S. M.: Thermal and metamorphic environment of subduction zone episodic tremor and slip, J. Geophys. Res.Sol. Ea., 114, B00A07, https://doi.org/10.1029/2008JB005978, 2009. a
Peacock, S. M.: Advances in the thermal and petrologic modeling of subduction zones, Geosphere, 16, 936–952, 2020. a, b
Peacock, S. M. and Hyndman, R. D.: Hydrous minerals in the mantle wedge and the maximum depth of subduction thrust earthquakes, Geophys. Res. Lett., 26, 2517–2520, 1999. a, b
Peacock, S. M. and Wang, K.: Seismic consequences of warm versus cool subduction metamorphism: Examples from southwest and northeast Japan, Science, 286, 937–939, 1999. a, b
PennistonDorland, S. C., Kohn, M. J., and Manning, C. E.: The global range of subduction zone thermal structures from exhumed blueschists and eclogites: Rocks are hotter than models, Earth Planet. Sc. Lett., 428, 243–254, 2015. a
Perrin, A., Goes, S., Prytulak, J., Davies, D. R., Wilson, C., and Kramer, S.: Reconciling mantle wedge thermal structure with arc lava thermobarometric determinations in oceanic subduction zones, Geochem. Geophy. Geosy., 17, 4105–4127, 2016. a
Petrini, C., Gerya, T., Yarushina, V., van Dinther, Y., Connolly, J., and Madonna, C.: Seismohydromechanical modelling of the seismic cycle: methodology and implications for subduction zone seismicity, Tectonophysics, 791, 228504, https://doi.org/10.1016/j.tecto.2020.228504, 2020. a
Plunder, A., Thieulot, C., and van Hinsbergen, D. J. J.: The effect of obliquity on temperature in subduction zones: insights from 3D numerical modeling, Solid Earth, 9, 759–776, https://doi.org/10.5194/se97592018, 2018. a
Ponko, S. C. and Peacock, S. M.: Thermal modeling of the southern Alaska subduction zone: insight into the petrology of the subducting slab and overlying mantle wedge, J. Geophys. Res.Sol. Ea., 100, 22117–22128, 1995. a
Pozgay, S. H., Wiens, D. A., Conder, J. A., Shiobara, H., and Sugioka, H.: Seismic attenuation tomography of the Mariana subduction system: Implications for thermal structure, volatile distribution, and slow spreading dynamics, Geochem. Geophy. Geosy., 10, Q04X05, https://doi.org/10.1029/2008GC002313, 2009. a
Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: Numerical recipes in C++, The art of scientific computing, Cambridge University Press, 2, 1002, ISBN 9780521880688, 1992. a
Qu, R., Zhu, W., Ji, Y., Xie, C., Zeng, D., and Zhang, F.: Subduction thermal regime, petrological metamorphism and seismicity under the Mariana arc, Sci. Rep., 13, 1948, https://doi.org/10.1038/s41598023290041, 2023. a
Ranalli, G.: Rheology of the Earth, Springer Science & Business Media, ISBN 9780412546709, 1995. a
Richards, F., Hoggard, M., Cowton, L., and White, N.: Reassessing the thermal structure of oceanic lithosphere with revised global inventories of basement depths and heat flow measurements, J. Geophys. Res.Sol. Ea., 123, 9136–9161, 2018. a, b, c, d, e, f, g, h, i, j, k
Rondenay, S., Abers, G. A., and van Keken, P. E.: Seismic imaging of subduction zone metamorphism, Geology, 36, 275–278, 2008. a
Rotman, H. M. and Spinelli, G. A.: Global analysis of the effect of fluid flow on subduction zone temperatures, Geochem. Geophy. Geosy., 14, 3268–3281, 2013. a, b
Royden, L. H.: The steady state thermal structure of eroding orogenic belts and accretionary prisms, J. Geophys. Res.Sol. Ea., 98, 4487–4507, 1993. a
Rüpke, L. H., Morgan, J. P., Hort, M., and Connolly, J. A.: Serpentine and the subduction zone water cycle, Earth Planet. Sc. Lett., 223, 17–34, 2004. a
Satake, K.: 4.19 – Tsunamis, in: Treatise on Geophysics, edited by: Schubert, G., 2nd edn., 477–504, Elsevier, Oxford, https://doi.org/10.1016/B9780444538024.000865, 2015. a
Satake, K. and Tanioka, Y.: Sources of tsunami and tsunamigenic earthquakes in subduction zones, Pure Appl. Geophys., 154, 467–483, https://doi.org/10.1007/s000240050240, 1999. a
Schellart, W. P. and Rawlinson, N.: Global correlations between maximum magnitudes of subduction zone interface thrust earthquakes and physical parameters of subduction zones, Phys. Earth Planet. In., 225, 41–67, 2013. a
Schmeling, H., Babeyko, A., Enns, A., Faccenna, C., Funiciello, F., Gerya, T., Golabek, G., Grigull, S., Kaus, B., Morra, G., Schmalholz, S. M., and van Hunen, J.: A benchmark comparison of spontaneous subduction models – Towards a free surface, Phys. Earth Planet. In., 171, 198–223, 2008. a
Scholz, C. H.: The mechanics of earthquakes and faulting, Cambridge university press, https://doi.org/10.1017/9781316681473, 2019. a, b, c
Schubert, G., Turcotte, D. L., and Olson, P.: Mantle Convection in the Earth and Planets 2 Volume Set, Cambridge University Press, https://doi.org/10.1017/CBO9780511612879, 2001. a
Sclater, J., Jaupart, C., and Galson, D.: The heat flow through oceanic and continental crust and the heat loss of the Earth, Rev. Geophys., 18, 269–311, 1980. a
Seipold, U.: Temperature dependence of thermal transport properties of crystalline rocks – a general law, Tectonophysics, 291, 161–171, 1998. a, b, c
Spinelli, G. A. and Wang, K.: Effects of fluid circulation in subducting crust on Nankai margin seismogenic zone temperatures, Geology, 36, 887–890, 2008. a
Stein, C. A. and Stein, S.: Constraints on hydrothermal heat flux through the oceanic lithosphere from global heat flow, J. Geophys. Res.Sol. Ea., 99, 3081–3095, 1994. a
Su, C., Liu, Y., Song, W., Fan, D., Wang, Z., and Tang, H.: Thermodynamic properties of San Carlos olivine at high temperature and high pressure, Acta Geochimica, 37, 171–179, 2018. a, b, c
Syracuse, E. M., van Keken, P. E., and Abers, G. A.: The global range of subduction zone thermal models, Phys. Earth Planet. In., 183, 73–90, 2010. a, b, c, d, e
Tichelaar, B. W. and Ruff, L. J.: Depth of seismic coupling along subduction zones, J. Geophys. Res.Sol. Ea., 98, 2017–2037, 1993. a, b
Tsitsulin, A.: Optimal qualitative colour palettes, http://tsitsul.in/blog/coloropt/ (last access: 10 May 2021), 2019. a
Turcotte, D. and Schubert, G.: Geodynamics, Cambridge University Press, https://doi.org/10.1017/CBO9780511843877, 2002. a, b
van Dinther, Y., Gerya, T. V., Dalguer, L. A., Corbi, F., Funiciello, F., and Mai, P. M.: The seismic cycle at subduction thrusts: 2. Dynamic implications of geodynamic simulations validated with laboratory models, J. Geophys. Res.Sol. Ea., 118, 1502–1525, 2013a. a
van Dinther, Y., Gerya, T. V., Dalguer, L. A., Mai, P. M., Morra, G., and Giardini, D.: The seismic cycle at subduction thrusts: Insights from seismothermomechanical models, J. Geophys. Res.Sol. Ea., 118, 6183–6202, 2013b. a
van Dinther, Y., Mai, P. M., Dalguer, L. A., and Gerya, T. V.: Modeling the seismic cycle in subduction zones: The role and spatiotemporal occurrence of offmegathrust earthquakes, Geophys. Res. Lett., 41, 1194–1201, 2014. a
van Keken, P. E., Kiefer, B., and Peacock, S. M.: Highresolution models of subduction zones: Implications for mineral dehydration reactions and the transport of water into the deep mantle, Geochem. Geophy. Geosy., 3, 1056, https://doi.org/10.1029/2001GC000256, 2002. a
van Keken, P. E., Currie, C., King, S. D., Behn, M. D., Cagnioncle, A., He, J., Katz, R. F., Lin, S.C., Parmentier, E. M., Spiegelman, M., and Wang, K.: A community benchmark for subduction zone modeling, Phys. Earth Planet. In., 171, 187–197, 2008. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag, ah, ai, aj
van Keken, P. E., Hacker, B. R., Syracuse, E. M., and Abers, G. A.: Subduction factory: 4. Depthdependent flux of H_{2}O from subducting slabs worldwide, J. Geophys. Res.Sol. Ea., 116, B01401, https://doi.org/10.1029/2010JB007922, 2011. a
van Keken, P. E., Kita, S., and Nakajima, J.: Thermal structure and intermediatedepth seismicity in the TohokuHokkaido subduction zones, Solid Earth, 3, 355–364, https://doi.org/10.5194/se33552012, 2012. a, b, c
van Keken, P. E., Wada, I., Sime, N., and Abers, G. A.: Thermal structure of the forearc in subduction zones: A comparison of methodologies, Geochem. Geophy. Geosy., 20, 3268–3288, 2019. a, b
Van Zelst, I.: irisvanzelst/xFieldstone: Final version xFieldstone used in Van Zelst et al. (2023, Solid Earth), Zenodo [code], https://doi.org/10.5281/zenodo.7957374, 2023. a, b, c
Van Zelst, I., Wollherr, S., Gabriel, A.A., Madden, E. H., and van Dinther, Y.: Modeling megathrust earthquakes across scales: Oneway coupling from geodynamics and seismic cycles to dynamic rupture, J. Geophys. Res.Sol. Ea., 124, 11414–11446, 2019. a
Van Zelst, I., Crameri, F., Pusok, A. E., Glerum, A., Dannberg, J., and Thieulot, C.: 101 geodynamic modelling: how to design, interpret, and communicate numerical studies of the solid Earth, Solid Earth, 13, 583–637, https://doi.org/10.5194/se135832022, 2022. a, b, c, d
Van Zelst, I., Thieulot, C., and Craig, T. J.: Data & scripts – The effect of temperaturedependent material properties on simple thermal models of subduction zones, Zenodo [data set], https://doi.org/10.5281/zenodo.7957674, 2023. a, b
Wada, I.: A simple picture of mantle wedge flow patterns and temperature variation, J. Geodynam., 146, 10848, https://doi.org/10.1016/j.jog.2021.101848, 2021. a
Wang, J., Zhao, D., and Yao, Z.: Seismic anisotropy evidence for dehydration embrittlement triggering intermediatedepth earthquakes, Sci. Rep., 7, 1–9, 2017. a, b
Wen, H., Lu, J.H., Xiao, Y., and Deng, J.: Temperature dependence of thermal conductivity, diffusion and specific heat capacity for coal and rocks from coalfield, Thermochim. Acta, 619, 41–47, 2015. a, b, c
Whittington, A. G., Hofmeister, A. M., and Nabelek, P. I.: Temperaturedependent thermal diffusivity of the Earth’s crust and implications for magmatism, Nature, 458, 319–321, 2009. a
Xu, Y., Shankland, T. J., Linhardt, S., Rubie, D. C., Langenhorst, F., and Klasinski, K.: Thermal diffusivity and conductivity of olivine, wadsleyite and ringwoodite to 20 GPa and 1373 K, Phys. Earth Planet. In., 143, 321–336, 2004. a, b, c, d, e
Yabe, S., Fukuchi, R., Hamada, Y., and Kimura, G.: Simultaneous estimation of in situ porosity and thermal structure from core sample measurements and resistivity log data at Nankai accretionary prism, Earth Planets Space, 71, 1–15, 2019. a
Yamasaki, T. and Seno, T.: Double seismic zone and dehydration embrittlement of the subducting slab, J. Geophys. Res.Sol. Ea., 108, https://doi.org/10.1029/2002JB001918, 2003. a, b
Zhan, Z.: Mechanisms and implications of deep earthquakes, Annu. Rev. Earth Pl. Sc., 48, 147–174, 2020. a