Thrusts control the thermal maturity of accreted sediments

Thermal maturity assessments of hydrocarbon-generation potential and thermal history rarely consider how upper-plate structures developing during subduction influence the trajectories of accreted sediments. Our thermomechanical models of subduction support that thrusts evolving under variable sedimentation rates and décollement strengths fundamentally influence the trajectory, temperature, and thermal maturity of accreting sediments. This is notably true for the frontal thrust, which pervasively partitions sediments along a low- and a high-maturity path. Our findings imply that interpretations of the distribution of thermal maturity cannot be detached from accounts of the length and frequency of thrusts and their controlling factors. Our approach takes these factors into consideration and provides a robust uncertainty estimate of maximum exposure temperatures as a function of vitrinite reflectance and burial depth. As a result, our models reduce former inconsistencies between predicted and factual thermal maturity distributions in accretionary wedges.

1 Introduction

Organic material transforms into coal, oil, and gas at rates primarily controlled by temperature. This transformation, critical for the hydrocarbon industry, is also useful to study the tectonic and sedimentary evolution of basins and orogens. The extent of this transformation in sediments, known as thermal maturity, can be measured as vitrinite reflectance, i.e. the percentage of incident light reflected from the surface of vitrinite particles in these sediments. Thermal maturity has been used to estimate the thermal evolution of igneous intrusions and seismic slip, the extent of low-grade metamorphism, porosity, and compaction in basin sediments, and the geothermal history of accreting material during subduction (e.g. Bostick and Pawlewicz, 1984; Rabinowitz et al., 2020; Fukuchi et al., 2017; Kamiya et al., 2017).

Inferences on the geothermal history of subduction margins based on thermal maturity depend on the trajectory followed by the accreting sediments (Miyakawa et al., 2019). Low-temperature, high-pressure metamorphic rocks in the subduction wedge are often attributed to the pressure maxima that typically predate the temperature maxima in accreted sediments undergoing diagenesis in the wedge (van Gool and Cawood, 1994). However, the existence of complicated patterns in sediment trajectories is supported by numerical models and field observations (Giunchi and Ricard, 1999). As the orogenic wedge evolves, sediments accreting along different paths reach different depths and velocities and are exposed to different regional peak temperatures. Miyakawa et al. (2019) proposed subdividing these trajectories based on their final characteristics, like thermal maturity. In this manner, the spatiotemporal evolution of sediments and their thermal maturity is regulated to the first order by the partition of incoming sediments along two endmember pathways: (i) a deeper path leading to elevated thermal maturities and constituted by under-thrusted material, the high thermal-maturity path, and (ii) a shallower path that typically lies closer to the surface or gets frequently exhumed to near-surface levels, the low thermal-maturity path.

Previous studies have used numerical and analogue approaches to study the trajectories of sedimentary particles, and their spatial and pressure–temperature evolution, as a function of changes in erosion, sedimentation, or décollement strength. The trajectory followed by under-thrusted sedimentary units is primarily determined by orogenic wedge dynamics and its controlling forces (Platt, 1986). Although these sediments may only be exhumed near the backstop of the wedge, the trajectories of other accreted sediments generally deflect toward the surface under the influence of erosion (Konstantinovskaia, 2005). In fact, sedimentary particle trajectories gradually shift from deflection toward the surface near the front of accretion to final exhumation near the wedge backstop (Wenk and Huhn, 2013). Still, even under-thrusted sediments, which would co-relate to high-maturity paths in our study, have variable pressure–temperature paths (Ruh, 2020). It is important to highlight that the majority of past studies have explored a snapshot of sediment trajectories, assuming that the general nature of trajectories remains relatively fixed with time or is stationary in nature. However, the intrinsic connection between thermal maturity and the comprehensive thermal exposure along the entire trajectory necessitates an in-depth investigation into the dynamic and transitory nature of sediment trajectories.

Although there is general consensus on the rate and extent of sediment trajectory transition from horizontal to vertical during accretion, the dynamic perturbations in sediment dynamics have yet to be adequately examined. For instance, while most studies show a great degree of correlation between the initial depth of incoming sediments and their final position in the wedge (e.g. Mulugeta and Koyi, 1992; Willett, 1992), a dynamic fluctuation in this correlation due to thrusting can result in non-stationary exhumation paths for accreting sediments in a wedge (e.g. Konstantinovskaia, 2005; Miyakawa et al., 2019). Much remains to be explored regarding the partition of high and low thermal-maturity paths and how sediments travel inside natural wedges, given the conventional assumption that accreting sediments remain at the same relative depth and translate along the adjacent “layers” without vertical mixing throughout the tectonic evolution of the wedge (Hori and Sakaguchi, 2011).

Our assessment identifies a primary gap in existing research: the prediction and mapping of the initial sediment influx to their location in the orogenic wedge. More specifically, the challenge lies in determining which portions of incoming sediment will predominantly constitute the core of the wedge and which will reside at comparatively shallower depths. Given that the maximum exposure temperature estimation from the thermal maturity is inherently reliant on the path of sediments inside the wedge, information on path diversity would inherently constrain the uncertainty in maximum exposure temperature used for the identification of paleo-thermal structures of subduction zones. Moreover, to better understand the time–depth paths of wedge sediments, their dependence on the initial state of undeformed sediments, and thus their thermal maturity, the factors that control the evolution of subduction–accretion systems, like sedimentation, erosion, and décollement strength, ought to be considered (Mannu et al., 2016; Simpson, 2010).

Here, we explore in detail the impact accretion in a subduction wedge has on the thermal maturity of its sediments. We simulate subduction–accretion using 2D finite-difference thermomechanical models incorporating empirical thermal conductivity values from the Nankai accretionary margin. We track the evolution of thermal maturity by computing vitrinite reflectance (%R_o) on each marker and throughout the model, using three well-established methods of %R_o computation, as accretion develops the wedge under different sedimentation rates and décollement strengths. These factors notably alter the trajectories and thermal maturities of incoming sediments. Particularly, thrusts define sharp thermal maturity boundaries leading to stark differences in the thermal maturity of sediments that accrete in different thrust blocks, even when they follow similar trajectories and lay nearby.

2 Geological settings and model generalisation

We use a generalised model for the subduction of an oceanic plate under a continental plate, with explicit integration of key parameters from the Nankai subduction margin off Kii Peninsula in south-west Japan. The Nankai subduction margin is a product of the ongoing, north-west-directed subduction of the Philippine Sea Plate beneath the Amurian Plate at a convergence rate of 4.1–6.5 cm yr⁻¹ (Seno et al., 1993; Heki et al., 1999; DeMets et al., 2010). Past studies posit the initiation of this subduction within the Nankai region at circa 6 Mya (Kimura et al., 2014). The accretionary wedge adjacent to the Nankai margin is marked by the accretion of thick sediment layers (>1 km), predominantly formed by overlying younger trench sediments atop Shikoku Basin sediments. Mean sedimentation rates of ∼0.4 mm yr⁻¹ for this area are calculated from on-land sediment data and may largely reach the trench through submarine channels (Korup et al., 2014).

Table 1 Properties of the different materials used for the model runs.

^∗ is strain-softened cohesion/coefficient of friction. ^# We have assumed flow law parameters such as A_D, E, V and n to be the same for dislocation and diffusion creep. T is temperature in Kelvin, and Z is the depth from the seafloor. The reference temperature for densities have been taken as the average temperature of the rock type. ^a Reference for densities (Turcotte and Schubert, 2002; Gerya and Meilick, 2011). ^b Reference cohesion values for sediments (Schumann et al., 2014). ^c Reference for angle of friction (Schumann et al., 2014; Ruh et al., 2014; Gerya and Meilick, 2011). ^d Reference for thermal conductivity (Clauser and Huenges, 1995; Sugihara et al., 2014). ^e Reference for flow laws and radiogenic heat production (Ranalli, 1995; Gerya and Meilick, 2011).

Another reason to select the Nankai subduction margin is that it is a particularly well-studied accretionary margin regarding its paleo-thermal history and thermal maturity distribution. For example, Underwood et al. (1992) and Sakaguchi (1999) used thermal maturity estimates from Shimanto accretionary wedge in the Nankai subduction margin to suggest that ridge subduction can explain the resulting paleo-heat flow. Following this, Ohmori et al. (1997) published a distribution of thermal maturity and maximum exposure temperature for the Shimanto accretionary wedge, identifying out-of-sequence thrusting in the region. The accretionary wedge adjacent to the Kumano forearc basin in the Nankai subduction margin has also been the subject of the NanTroSEIZE (Nankai Trough Seismogenic Zone) project, which drilled C0002 borehole during the 2012 Integrated Ocean Discovery Program Expedition 338. C0002 borehole is located south-west of Japan's Kii Peninsula in the Kumano Basin, within the Nankai accretionary margin, and extends 3348 m below the seafloor. Having data on both thermal maturity and thermal conductivity from the same borehole in subduction wedges is quite uncommon. To our knowledge, the C0002 borehole, located next to the Kumano forearc basin, is the only place where such data can be found in an accretionary wedge. Because of this unique characteristic, the C0002 borehole serves as an excellent dataset for validation purposes. We modify the thermal conductivity computation for sediments and décollement (see Table 1) to match the empirical relationship between depth and thermal conductivity, as measured on core samples in the borehole C0002 (Sugihara et al., 2014).

While these adjustments render our models somewhat specific to the Nankai accretionary wedge, we propose that the thermal conductivity values and trend are representative of patterns typically observed in forearc basins and accretionary wedges across the globe, making it broadly applicable to general subduction margins. For instance, in our simulations, the sediment thermal conductivity within our wedge steadily increases with depth from 0.96–4.0 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, which is within the range of thermal conductivity estimates of comparable depth in other subduction zones, such as the Hikurangi subduction margin, Japan Trench, and Taiwan subduction zone (Fig. S1 in the Supplement, Henrys et al., 2003; Lin et al., 2014; Chi and Reed, 2008). As a result, we compare our simulation results not only to thermal maturity values in the Nankai accretionary margin but also to those of the Miura–Boso plate subduction margin in central Japan and the fold-and-thrust belts of the Western Foothills complex in western Taiwan.

3 Methods

We employ I2VIS, a conservative finite-difference 2-D thermomechanical subduction–accretion model with visco-plastic/brittle rheology (Gerya and Yuen, 2003a, b). The code solves the governing equations for the conservation of mass, momentum, and heat as well as the advection equation with a non-diffusive marker-in-cell scheme constrained by thermal conductivity values inferred from Nankai accretionary wedge. Our numerical approach has several advantages over earlier attempts to simulate thermal maturity in an accretionary wedge, such as a more realistic geothermal profile, variable particle paths, and thermal evolution. In the following sections, we provide information regarding the governing equations, the modified thermal conductivity formulations based on the C0002 borehole, boundary conditions, the rheological model, model setup, surface processes, and the computation of thermal maturity.

3.1 Governing equations

The mass conservation is described by the continuity equation with the Boussinesq approximation of incompressibility.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial v_x}{\partial x}}+{\displaystyle \frac{\partial v_y}{\partial y}}=\mathrm{0},\end{array}$$

where v_x and v_y are horizontal and vertical components of velocity.

The equation for conservation of momentum with an incompressibility assumption is expressed in the 2D Stokes equation, for the x axis and y axis respectively,

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{xx}}{\partial x}}+{\displaystyle \frac{\partial {\mathit{\sigma }}_{xy}}{\partial y}}={\displaystyle \frac{\partial P}{\partial x}},\end{array}$$

where σ_xx, σ_xy, σ_yy are components of the deviatoric stress tensor; x and y denote the horizontal and vertical coordinates; and P is pressure.

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{yy}}{\partial y}}+{\displaystyle \frac{\partial {\mathit{\sigma }}_{xy}}{\partial x}}={\displaystyle \frac{\partial P}{\partial y}}-g\mathit{\rho }(T,P,C),\end{array}$$

where ρ is rock density and depends on rock type (C), temperature (T), and pressure as $\mathit{\rho }(T,P)={\mathit{\rho }}_{\mathrm{0}}(\mathrm{1}-\mathit{\xi }(T-T_{\mathrm{0}}\left)\right)(\mathrm{1}+\mathit{\varsigma }(P-P_{\mathrm{0}}\left)\right)$, where ξ is the coefficient of thermal expansion taken to be $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ K⁻¹ for all rock markers and 0 for air/water, ς is the coefficient of compressibility taken to be $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ MPa⁻¹ for all rock markers and 0 for air/water, ρ₀ is the reference density at reference temperature (T₀=298.15 K) and reference pressure (P₀=10⁵ K).

The thermal equation used in the model is as follows:

$$\begin{array}{}\text{(4)}& \mathit{\rho }{C}_P{\displaystyle \frac{\mathrm{D}T}{\mathrm{D}t}}={\displaystyle \frac{\partial q_x}{\partial x}}+{\displaystyle \frac{\partial q_y}{\partial y}}+H_{\mathrm{r}}+H_{\mathrm{a}}+H_{\mathrm{s}},\end{array}$$

where

$$\begin{array}{}\text{(5)}& {\displaystyle }{q}_x=-k(T,C,y){\displaystyle \frac{\partial T}{\partial x}},q_y=-k(T,C,y){\displaystyle \frac{\partial T}{\partial y}},\text{(6)}& {\displaystyle }\begin{array}{rl}& H_{\mathrm{a}}=T\mathit{\xi }{\displaystyle \frac{\mathrm{D}P}{\mathrm{D}t}},\\ & H_{\mathrm{s}}={\mathit{\sigma }}_{xx}{\dot{\mathit{\epsilon }}}_{xx}+{\mathit{\sigma }}_{yy}{\dot{\mathit{\epsilon }}}_{yy}+{\mathit{\sigma }}_{xy}{\dot{\mathit{\epsilon }}}_{xy}+{\mathit{\sigma }}_{yx}{\dot{\mathit{\epsilon }}}_{yx},\\ & H_{\mathrm{r}}=\text{const},\end{array}\end{array}$$

where $\frac{\mathrm{D}}{\mathrm{D}t}$ is the Lagrangian time derivative; ${\dot{\mathit{\epsilon }}}_{xx}$, ${\dot{\mathit{\epsilon }}}_{xy}$, ${\dot{\mathit{\epsilon }}}_{yy}$ are components of the strain rate tensor; q_x, q_y are the components of heat flux in the horizontal and vertical direction; g is the vertical gravitational acceleration; C_P is the isobaric heat capacity; and H_r, H_a, H_s denote the radioactive, adiabatic, and shear heat production respectively. $k(T,C,y)$ is the thermal conductivity, a function of composition, depth, and temperature (Table 1). The radioactive heat production H_r is constant for a rock type as mentioned in Table 1.

In order to accurately assess thermal maturity, it is crucial to consider the temperature distribution, which necessitates a realistic thermal conductivity profile when modelling thermal maturity. Many geodynamic models assume that thermal conductivity decreases as temperature increases, following a defined relationship (e.g. Clauser and Huenges, 1995). These models typically predict a decrease in thermal conductivity with depth within accretionary wedges, as geothermal profiles tend to increase in temperature with depth. However, empirical data reveal a different trend: thermal conductivity increases with depth, primarily due to sediment porosity influencing shallow thermal conductivity (Henrys et al., 2003; Lin et al., 2014). Additionally, the thermal conductivity values calculated using the Clauser and Huenges (1995) model are significantly higher than those observed at shallow depths (<3 km). To address these disparities, we incorporate the observed empirical relationship between depth and thermal conductivity from the IODP Site C0002 borehole in the Nankai accretionary wedge into our simulations. By adjusting the thermal conductivity formulation for sediments based on temperature and depth, we aim to replicate the empirical relationship observed in the core samples taken from the borehole at IODP Site C0002 (Sugihara et al., 2014) and account for the decrease in thermal conductivity near the surface caused by increased porosity. We modify the thermal conductivity formulation for sediments as a function of temperature and depth as follows.

$$\begin{array}{}\text{(7)}& k_{\text{sed}}=k_{\mathrm{0}}+{\displaystyle \frac{\mathrm{807}}{T+\mathrm{77}}}\left(\mathrm{1}-\exp \left({\displaystyle \frac{-Z^{\mathrm{2}}}{\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{7}}}}\right)\right)\end{array}$$

k₀=0.96 and 1.5 for the wedge sediment and décollement respectively. The larger thermal conductivity of the décollement emulates higher heat transfer in shear zones due to fluid advection (Fig. S1).

3.2 Rheological model

The expression for effective creep viscosities (η_eff) is computed as follows.

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\eta }}_{\text{disl}}=\mathrm{0.5}(\dot{{\mathit{\epsilon }}_{\text{II}}}{)}^{\frac{\mathrm{1}}{n}-\mathrm{1}}{A}_D^{\frac{\mathrm{1}}{n}}\exp \left(-{\displaystyle \frac{E_{\mathrm{a}}+V_{\mathrm{a}}P}{nRT}}\right),\text{(9)}& {\displaystyle }{\mathit{\eta }}_{\text{diff}}=\mathrm{0.5}{\displaystyle \frac{A_D}{S^{n-\mathrm{1}}}}\exp \left(-{\displaystyle \frac{E_{\mathrm{a}}+V_{\mathrm{a}}P}{RT}}\right),\text{(10)}& {\displaystyle }{\mathit{\eta }}_{\text{eff}}={\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\text{disl}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\text{diff}}}}\right)}^{-\mathrm{1}},\end{array}$$

where R is the gas constant (8.314 $\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{mol}}^{-\mathrm{1}}$) and A_D, n, m, E_a, and V_a are experimentally determined rheological parameters. A_D is the material constant (${\mathrm{Pa}}^{-n}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-m}$), n is the stress exponent, m is the grain size exponent, E_a is activation energy (J mol⁻¹), V_a is activation volume (J Pa⁻¹), and S is a stress factor for diffusion creep assumed to be 3×10⁴ Pa.

$$\begin{array}{}\text{(11)}& \dot{{\mathit{\epsilon }}_{\text{II}}}=\sqrt{{\displaystyle \frac{\dot{{\mathit{\epsilon }}_{ij}}\cdot \dot{{\mathit{\epsilon }}_{ij}}}{\mathrm{2}}}}\end{array}$$

The model uses visco-plastic rheology to account for both brittle rheology of the shallower and colder rigid lithosphere and deeper, hotter ductile lithosphere and asthenosphere. Using the plastic yield threshold as per the Drucker–Prager criterion we limit effective viscosity as

$$\begin{array}{}\text{(12)}& {\mathit{\eta }}_{\text{eff}}\le {\displaystyle \frac{P\cdot \sin \mathit{\phi }\cdot (\mathrm{1}-\mathit{\lambda })+C\cdot \cos \mathit{\phi }}{\mathrm{2}\dot{{\mathit{\epsilon }}_{\text{II}}}}},\end{array}$$

where C is cohesion and φ is an effective internal angle of friction or μ=tan φ where μ is the coefficient of internal friction and λ the fluid pressure ratio assumed to be 0 in all the simulations.

3.3 Boundary conditions

A free-slip boundary condition is implemented on all boundaries, except on the lower boundary, which is permeable in the vertical direction. On the lower boundary we implement an external free slip condition similar to where a free slip condition is satisfied at an external boundary such that

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial V_x}{\partial x}}=\mathrm{0}\text{ and }{\displaystyle \frac{\partial V_y}{\partial y}}={\displaystyle \frac{V_y}{\mathrm{\Delta }{Y}_{\text{external}}}},\end{array}$$

where V_x and V_y are the velocities in the horizontal and vertical directions at the boundary, ΔY_external is the depth that lies outside the modelling domain, and free slip condition is maintained. Similarly, we set thermally insulating boundary conditions on all sides except the lower one, where the external thermal boundary condition is implemented.

3.4 Surface processes

The rock–water/air boundary is simulated by an adaptive irregular grid that is advected horizontally and vertically and is coupled to the thermomechanical grid which controls the tectonic deformation of the surface. Apart from the tectonic changes, surface processes prescribed in the model can also change the topography. The surface process in the model is controlled by the conversion of rock markers to air/water and vice versa. All sedimentation in the model happens as a focused deposition of sediments from sea to land in morphological depressions (e.g. trench) and is modelled as follows (Fig. S2):

$$\begin{array}{}\text{(14)}& y_{\text{new}}=y_{\text{old}}+K\cdot y_{\text{fill}},\end{array}$$

where $K=min\left(\frac{V_{\text{budget}}}{V_{\text{basin}}},\mathrm{1}\right)$.

The shape of the basin and the resolution of the surface grid can lead to overfilling or underfilling when using Eq. (14) to fill the basin. To address this issue, we calculate the volume of deposited sediments and adjust for any deficit or overfill in the subsequent step. This ensures that, over time, the total value of sedimentation remains consistent with the prescribed value. However, it is challenging to ensure that all sediments added in a particular step are accommodated within the basins, especially in models with high sedimentation rates where significant runoff occurs. Therefore, the sedimentation rates mentioned in this study are computed as effective sedimentation rates after the model runs, rather than being predetermined. We perform multiple models runs (approximately 100) with sedimentation rates uniformly distributed in the range of 0.1–0.9 mm yr⁻¹. From these runs, we select models that exhibit appropriate sedimentation rates. This selection process ensures that the average sedimentation rates across all our models (ranging from 0.1 to 0.9 mm yr⁻¹) fall within the observed sedimentation rates in our chosen natural equivalent, the Nankai accretionary wedge in the south-western subduction margin of Japan (Korup et al., 2014). For more specific information about the model run and prescribed sedimentary conditions, please see Table 2.

Table 2 Model runs and their specific characteristic observations.

φ_b is décollement strength (internal angle of friction). φ is sediment strength. φ_ss is sediment strength (strain weakened)/(internal angle of friction). SR is average sediment rate (mm yr⁻¹). λ is pore–fluid pressure ratio. L is average length of the wedge (in kilometres) between ∼2.5 and 7.5 Myr. Length of the wedge is computed as the distance between trench and backstop (set at 1850 km from the right edge of the modelling domain). β is average basal dip angle β (in ^∘) between ∼2.5 and 7.5 Myr measured by fitting a line in the basal surface. α is average surface slope angle α (in ^∘) between ∼2.5 and 7.5 Myr measured computing the slope of fitting the best fitted line in the surface. D is average distance between the first and second frontal thrust between ∼2.5 and 7.5 Myr (in kilometres). The frontal thrust is always identified from the trench. The send thrust is identified by the high strain rate and deviation of the weak décollement material from the trend of oceanic plate. α predicted (φ_ss / φ) is the surface slope predicted using critical wedge theory using the β observed in the model and sediment strength (initial/strain weakened). T is average time a frontal thrust remains active between ∼3.5 and 7.5 Myr. 〈R_o%〉 is average vitrinite reflectance of the wedge between ∼3.5 and 7.5 Myr. %_top is proportion of >1 eventual R_o% (vitrinite reflectance at 7.5 Myr) in the shallow half of the incoming sediment at 2.5 Myr. %_bottom is proportion of >1 eventual R_o% (vitrinite reflectance at 7.5 Myr) in the deep half of the incoming sediments.
^∗ Please see Fig. S18 for details on the various measurement done on the wedge.

3.5 Thermal maturity calculation

The model computes the %R_o of each marker to estimate the thermal maturity of sediments during the model run using three widely used methods of thermal maturity modelling: Easy%R_o (Burnham and Sweeney, 1989; Sweeney and Burnham, 1990), Simple%R_o (Suzuki et al., 1993) and Basin%R_o (Nielsen et al., 2017). All the models presented here employ a simplified parallel Arrhenius reaction model, which accommodates an array of activation energies for every component of the kerogen, allowing it to estimate thermal maturity under varying temporal and thermal scales. The Easy%R_o model by Sweeney and Burnham (1990) can be described using the following equations:

$$\begin{array}{}\text{(15)}& {\displaystyle }{x}_i\left(t\right)=x_{\mathrm{0}i}\exp \left(-\int A\exp \left(-{\displaystyle \frac{E_{\mathrm{a}i}}{RT\left(t\right)}}\right)\mathrm{d}t\right),\text{(16)}& {\displaystyle }X\left(t\right)=\sum _{i=\mathrm{1}}^N{x}_i\left(t\right),\text{(17)}& {\displaystyle }F\left(t\right)=X(t=\mathrm{0})-X\left(t\right),\text{(18)}& {\displaystyle }\mathit{\%}{R}_{\mathrm{o}}=\mathit{\%}{R}_{\mathrm{o}\mathrm{0}}\exp \left(\mathrm{3.7}F\right),\end{array}$$

where x_oi are weights of reactions for ith component of the kerogen also described as the stoichiometric coefficient, A is the pre-exponential factor, E_ai is the activation energy of the ith component of the kerogen, R is the gas constant, T(t) is the temperature history, F is the percentage of fixed carbon, and %R_o0 is the vitrinite reflectance of the immature unaltered sediment. Sweeney and Burnham (1990) provided a set of 20 activation energies (E_ai) and the stoichiometric coefficient (x_oi) listed in Table 3. All thermal models used in this study use the same method of vitrinite reflectance computation albeit with different sets of activation energies, stoichiometric coefficient, pre-exponential factor, and %R_o0. Table 3 provides a comprehensive list of all these parameters.

Table 3 Parameters for Easy%R_o, Simple%R_o and Basin%R_o vitrinite reflectance model.

$A_{\text{Easy}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{13}}$ and %R_o0=0.2, $A_{\text{Simple}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{13}}$ and %R_o0=0.2, $A_{\text{Basin}}=\mathrm{9.7029}\times {\mathrm{10}}^{\mathrm{12}}$ and %R_o0=0.2104.

All these approaches for computing %R_o yield similar trends albeit with different absolute values. In the interest of clarity, we have mostly illustrated Easy%R_o, which is the most extensively used method for vitrinite reflectance computation and hereafter we refer Easy%R_o as %R_o, unless explicitly stated. %R_o is set to %R_oo in sediment markers at the start of the model until 2.5 Myr, while %R_o in markers for other rocks, air, and water is undefined at all times. After 2.5 Myr, the model computes %R_o on each marker as a function of temperature (T), time (t), and percentage of fixed carbon (F). The initial %R_o of newly deposited sediments is computed using an assumed water–sediment interaction temperature assumed to be the same as the thermocline. The thermocline used in the model has been estimated using the data obtained and made freely available by International Argo Program and the national programmes that contribute to it for the region near Nankai (Fig. S3; https://argo.ucsd.edu, last access: 2 August 2023, https://www.ocean-ops.org, last access: 2 August 2023).

3.6 Model setup

The modelling domain is 3500 km wide and 350 km deep and is divided into 3484×401 nodes populated with ∼125 million markers (Fig. 1). The high resolution of 220 m (horizontal)×130 m (vertical) that we assign to the site of accretionary wedge evolution decreases steadily toward the edges of the modelling domain to a minimum resolution of 3000 m×3200 m. The simulation consists of an oceanic plate converging with a velocity of ∼5 cm yr⁻¹ and subducting beneath a continental plate (Fig. 1). The convergence is prescribed internally using highly viscous nodes inside the oceanic and continental plates near the boundary of the models. The oceanic plate consists of a 1 km thick upper oceanic crust and a 7 km thick lower crust (Akuhara, 2018). The thickness of the oceanic lithosphere depends on its age, which is set to 20 Myr at the start of the simulation (Turcotte and Schubert, 2002). The initial age of the oceanic lithosphere corresponds to the age of the subducting lithosphere in the Nankai subduction margin (Zhao et al., 2021). Displacement along the megathrust, at the contact between subducting oceanic plate and the overriding continental plate, occurs in a relatively weak basal layer in accretionary wedges across the globe (Byrne and Fisher, 1990). We simulate this with a predefined configuration at the interplate, with a 350 m thick weak décollement below a sediment layer that is 1 km thick. The wedge forms above this interphase by the accretion of sediments against the continental plate. The continental plate consists of an upper and lower continental crust with thicknesses of ∼20 and ∼15 km respectively (Akuhara, 2018) and is underlain by a mantle lithosphere of ∼25 km. We use a thin (10 km) “sticky air” layer to overlay the top face of the rock strata inside the model which is a fluid with a low viscosity of 5×10¹⁷ Pa s, and a low density, similar to air (white in Fig. 1) or water (light blue in Fig. 1) (Crameri et al., 2012). The transition between the lithosphere and asthenosphere is prescribed to occur at 1300 ^∘C. A weak layer is emplaced at the junction of both plates, which fails mechanically and leads to subduction initiation. All sediments (light and dark brown in Fig. 1) are rheologically identical, but colours are alternated in time to allow tracking the development of different geological structures. Readers are referred to Table 1 for the rheological and thermal properties of all the materials used. Note that in our models, we measure all distances from the point where the continental and oceanic plates are initially located, which is situated 1850 km from the right boundary of the modelling area. The terms “landward” and “seaward” indicate the relative direction towards the continental plate or the oceanic plate respectively. The “backstop” refers to the edge of the continental plate that buttresses the wedge and acts akin to an indenter for the accretionary wedge. The “forearc high” represents the highest point in the forearc zone, which includes both the accretionary wedge and the forearc basin.

Figure 1 Initial model setup. (a) The lithological and geothermal map of the whole computational domain with boundary conditions. (b) The zoomed lithological and geothermal map of the inset illustrates the junction of continental and oceanic plates. The colours represent different lithology of the materials used in the models, with upper and lower crust represented by light and dark grey and upper and lower oceanic crust represented by dark and light green. The arrows around the computational domain represent the imposed boundary conditions, while the white contour lines (dashed in the zoomed panel) show the geothermal gradients used for the initial model. The numbers on the white contour lines represent the temperature values in ^∘C for the contour.

3.7 Experimental strategy

Here, we present a total of 10 simulations that vary in their effective basal friction or their effective sedimentation rate to discern patterns of thermal maturity evolution in wedge sediments. Models $M_{\mathrm{0}}^{\mathrm{4.5}}-M_{\mathrm{0}}^{\mathrm{14.5}}$ have no sedimentation and effective internal angle values for the décollement of φ_b=4.5, 7, 9.5, 12, and 14.5^∘ respectively. The chosen range of effective décollement strength is well within the range of values postulated by several studies for the Nankai accretionary wedge (Tesei et al., 2021). The rest of the models ($M_{\mathrm{0.1}}^{\mathrm{9.5}}-M_{\mathrm{0.9}}^{\mathrm{9.5}}$) have a medium-strength décollement and variable effective sedimentation rate ranging from 0.1 to 0.9 mm yr⁻¹. In all the models presented in this study, sedimentation is limited to the trench, extending from the sea to the land. Restricting sedimentation to the trench allows us to observe and analyse the length and frequency of thrust sheets, enabling comprehensive investigation of their role in determining sediment trajectories. With these models, we evaluate the particle trajectory and %R_o of accreting sediments as a function of décollement strength and sedimentation rate. To restrict the number of parameters influencing our observations, models have no erosion. Moreover, all models lack surface processes during the first ∼2.5 Myr and have sedimentation thereafter. Strain-softening has been modelled as a linear decrease in angle of friction (φ) and cohesion between a cumulative strain of 0.5 and 1.5. Sediments used in the model have an angle of friction (φ) of 30^∘ before a cumulative strain of 0.5 and a strain-softened value of 20^∘ after a threshold of 1.5 cumulative strain. Strain softening has been used in wedges to mimic the weakening of faults and shear zones due to lubrication, with values threshold taken from previous numerical studies (Hickman et al., 1995; Ruh et al., 2014).

4 Results

In our models, subduction begins at 0.1 Myr by failure of the weak material between continental and oceanic plate (Figs. 2 and S4–S13, also see supporting information videos). Continued and sustained accretion of sediments against the deforming continental crust forms the accretionary wedge from the interplate contact landwards. After ∼5 Myr, all models develop a distinct wedge in agreement with the critical wedge theory (Davis et al., 1983). Surface slopes, measured by fitting a line in the surface of the wedge for every time step between 2.5 and 7.5 Myr and reported as mean ± standard deviation, increase systematically, as effective basal friction increases from ∼4.5 to ∼14.5^∘ (Figs. 1 and S4–S13, Table 2, $M_{\mathrm{0}}^{\mathrm{4.5}}-M_{\mathrm{0}}^{\mathrm{14.5}}$). Whereas models with a relatively weaker décollement (as $M_{\mathrm{0}}^{\mathrm{4.5}}$, φ_b=4.5^∘) have surface slopes of 0.95 ± 0.3^∘, models with very strong décollement (as $M_{\mathrm{0}}^{\mathrm{14.5}}$, φ_b=14.5^∘) have slopes as steep as 5.9 ± 1^∘ (Table 2). Our estimations of surface slopes consistently exhibit an excess of approximately 1.5^∘ compared to the surface slopes predicted by the critical wedge theory (Table 2). This is probably due to the penetration of weaker décollement material into high shear zones, resulting in faults that are weaker than the strain-softened wedge material.

Figure 2 Typical thermomechanical evolution of the accretionary wedge for model. The illustrated figure is for the model $M_{\mathrm{0}}^{\mathrm{7}}$ at (a) 0.5, (b) 2.5, (c) 5.0, and (d) 7.5 Myr. Similar figures for other models have been illustrated in supplementary images. The colour map for the panels is same as Fig. 1.

Models without trench sedimentation grow solely by accretion of incoming seafloor sediments, with frequent nucleation of frontal thrusts. Models with weaker décollements develop thrust sheets that are lengthier but remain active for shorter periods. This is clear when comparing, for models with increasingly strong décollement ($M_{\mathrm{0}}^{\mathrm{4.5}},M_{\mathrm{0}}^{\mathrm{7}},M_{\mathrm{0}}^{\mathrm{9.5}},M_{\mathrm{0}}^{\mathrm{5}},M_{\mathrm{0}}^{\mathrm{14.5}}$), the average distance between first and second frontal thrusts, which are 15.5 ± 7.0, 12.1 ± 3.6, 8.8 ± 3.3, 8.7 ± 2.1, and 8.0 ± 1.8 km respectively. Increasing sedimentation rate also leads to an increase in thrust sheet length from 7.3 ± 1.1 km for model $M_{\mathrm{0.1}}^{\mathrm{9.5}}$ to 13.8 ± 7.8 km in model $M_{\mathrm{0.9}}^{\mathrm{9.5}}$.

In models with similar basal friction, models with higher sedimentation rates have lengthier thrust sheets that remain active for longer periods (Table 2). Steeper surface slopes with increased décollement strengths and change in thrush sheet length with sedimentation and décollement strength are well-known effects that have been confirmed by previous numerical (Ruh et al., 2012) and analogue (Malavieille and Trullenque, 2009; Storti and Mcclay, 1995) models. All the reported values are mean ± standard deviation values recorded between 2.5 and 7.5 Myr in individual models. All models exhibit a temperature gradient that corresponds well with the temperature profile observed in the boreholes at IODP Site C0002 in the Kumano forearc basin, on top of the Nankai accretionary wedge (Fig. S14).

4.1 Thermal maturity of the wedge

Sediments are more thermally mature in wedges that have a higher sedimentation rate or décollement strength. For example, the mean %R_o of simulations for wedges with the highest sedimentation is 12 % higher (0.75) than in those without sedimentation ($M_{\mathrm{0}}^{\mathrm{4.5}}$, Table 2, Fig. 3). Similarly, simulations of wedges with the strongest décollement have the highest mean %R_o (0.94) of all the simulations presented in this study.

Figure 3 Distribution of thermal maturity for different models at ∼6.0 Myr (3.5 Myr of thermal maturation). Panels A1–A5 show the thermal maturity distribution (computed using Easy%R_o) in subduction wedges of models as a function of décollement strength respectively. Panels A6–A10 show the thermal maturity distribution in subduction wedges of models function of sedimentation rate respectively. The grey colour of the markers indicate that no thermal maturity change in these sediments occurred. Panels B1–B10 and C1–C10 similarly show the thermal maturity distribution in subduction wedges computed using Simple%R_o and Basin%R_o, respectively.

Thermal maturity values increase with depth and landward distance from the trench to the forearc high irrespective of the décollement strength, sedimentation rates, and method of thermal maturity computation (Figs. 3 and 4). The absolute value of %R_o and the rate at which thermal maturity values increase landward from the trench are larger for wedges with high décollement strength (Fig. 4A). For wedges characterised by the same décollement strength but higher trench sedimentation, we observe that the rate of thermal maturity increases in a landward direction from the trench and remains consistent across these wedges (Fig. 4B). Comparing the values of %R_o along a horizontal marker at the depth of trench in several models emphasises this result; the model with the highest décollement strength reaches a maximum %R_o of 1.25 and has the highest rate of landward increase in thermal maturity (Fig. 4A). However, all models with similar décollement strength but different sedimentation do not visibly vary in their rate or magnitude of landward increase in thermal maturity. All models show a decrease in thermal maturity landward of the forearc high, commonly 0.2 %R_o. Other interesting observations that we explore below are the increased thermal maturity occurring in the vicinity of thrusts and the reversal in sediment maturity around out-of-sequence thrust active over longer times visible across several models (e.g. Fig. 3).

Figure 4 The variation in %R_o for a horizon as indicated by the orange band in the inset at 7.5 Myr. Panels A1 and A2 show all the models with different décollement strength. Panels B1 and B2 show all the models with different sedimentation rates. Horizons in panels A1 and B1 are located at 1 km depth from the surface, while in panels A2 and B2 the horizons are horizontal zones located at the trench depth. The horizontal distance from the backstop is normalised by the wedge length. Horizontal distance 0 represents the fixed backstop and 1 represents the trench.

The magnitude of %R_o varies consistently among Easy%R_o, Simple%R_o, and Basin%R_o. On average Easy%R_o had the smallest values, followed very closely by Basin%R_o (with an average difference of only 0.02). However, Simple%R_o had the highest average value of thermal maturity, being 0.16 and 0.13 higher than Easy%R_o and Basin%R_o (Fig. 3).

4.2 Sediment trajectory inside the wedge

In wedges with a higher décollement strength or sedimentation rate, sediments tend to follow high-maturity paths in larger proportions. We demonstrate this effect by creating a map of the thermal maturity of sediments at 7.5 Myr of the model run, mapped to their spatial position at 2.5 My of the model run to analyse the spatial correlation between sediment position (depth and distance) from the trench and thermal maturity (Fig. 5). We also show the mean thermal maturity attained by sediments at a given horizontal distance from the trench during this period by a dashed black line in Fig. 5. The scatter plot shows sharp changes in eventual thermal maturity with horizontal distance from the trench that relate to changes in sediment trajectory. The mean thermal maturity is also variable along the horizontal length of the wedge and has a periodicity (Λ) increasing in distance with higher sedimentation rate but relatively constant with changing basal friction (Fig. 5). The periodicity of mean %R_o was computed by finding the average wavelength of the auto-correlated mean %R_o. Whereas the mean thermal maturity has a short periodicity of ∼7.2 km for the model $M_{\mathrm{0}}^{\mathrm{9.5}}$ with no sedimentation rates, the model $M_{\mathrm{0.9}}^{\mathrm{9.5}}$ shows the longest periodicity of 21 km. However, for all models with no sedimentation ($M_{\mathrm{0}}^{\mathrm{4.5}}-M_{\mathrm{0}}^{\mathrm{14.5}}$), the periodicity remains relatively consistent in the range of 7–8 km.

Figure 5 Map of thermal maturity at 7.5 Myr mapped to sediments at 2.5 Myr. Panels A1–A5 and B1–B5 show the mapping for models $M_{\mathrm{0}}^{\mathrm{4.5}}$–$M_{\mathrm{0}}^{\mathrm{14.5}}$ and $M_{\mathrm{0.1}}^{\mathrm{9.5}}$–$M_{\mathrm{0.9}}^{\mathrm{9.5}}$. The vertical axis (distance from the oceanic plate) has been corrected for the bending of the plate. The horizontal axis represents the distance of sediments from the trench. The grey colour of the markers indicates that these sediments have been eroded/reworked due to slope failure. The dashed black line represents the mean %R_o attained sediment at a given distance from the trench. Λ represents the horizontal periodicity in mean %R_o for the given model.

Figure 3 also represents the distribution of trajectories that exist in an accretionary wedge and how these trajectories get impacted under trench sedimentation (a subset of these trajectories can be viewed in Fig. S15). Whereas in wedges with weak décollements ($M_{\mathrm{0}}^{\mathrm{4.5}}$), none of the shallowest half of incoming sediments reach %R_o>1 in 5 Myr, 2 % of sediments reach this value in wedges with strong décollement ($M_{\mathrm{0}}^{\mathrm{14.5}}$). The effects of décollement strength in the thermal maturity of sediments can be quantified as well at deeper levels, with one-eighth vs. more than one-half of the sediments surpassing values of %R_o=1 for the deepest half of incoming sediments (12 % and 54 % respectively) in weak vs. strong décollement wedges ($M_{\mathrm{0}}^{\mathrm{4.5}}$ vs. $M_{\mathrm{0}}^{\mathrm{14.5}}$) respectively. In wedges for the model without sedimentation ($M_{\mathrm{0}}^{\mathrm{9.5}}$), the top half of the incoming sediments fail to achieve %R_o>1, as opposed to ∼15 % of them reaching %R_o>1 in the models with a sedimentation rate of 0.9 mm yr⁻¹ ($M_{\mathrm{0.9}}^{\mathrm{9.5}}$). In sum, the proportion of sediments in the top half and bottom half of the wedge that reach high maturity steadily increases with both sedimentation rate and décollement strength (Table 2).

4.3 Patterns of trajectory and thermal maturity in incoming sediments

The diversity in the trajectory of sediments in the wedge leads to a plethora of pathways in which the sediments can become thermally mature and thus introduces epistemic uncertainty in the estimation of maximum exposure temperature. Figure 6 captures this uncertainty, where we plot the maximum exposure temperature as a function of %R_o for all the models simulated in this study. The colours of individual markers represent the depth of the markers normalised by the thickness of the wedge represented as Y_n (See Fig. S16 for more details). We find that almost all the models show a remarkable similarity in their relationship between maximum exposure temperature and %R_o (for individual models please see Fig. S16) and differ mostly in their proportion of sediments with extreme values of %R_o. We observe that the typical uncertainty in maximum exposure temperature increases with an increase in values of %R_o with the ∼15 ^∘C interval at around %R_o=0.2 compared to the ∼33 ^∘C interval at %R_o=3 (both for 95 % confidence interval, Fig. 6b). Moreover, we observe that incorporating information about the normalised depth of sediments (Y_n) significantly aids in constraining the maximum exposure temperature. For instance, although the overall uncertainty at %R_o=1, is ∼23 ^∘C, for sediments with a Y_n of 0.2–0.4 the uncertainty greatly reduces to only ∼10.5 ^∘C. Thus, the range of thermal maturity values for sediments clearly has a large correlation with their trajectories.

Figure 6 (a) Vitrinite reflectance (%R_o) vs. maximum exposure temperature in all models. The colours in panel a represent the depth of the sediments at 7.5 Myr normalised by the thickness of the wedge (Y_n). (b) Range of 95 % CI for Easy%R_o, Simple%R_o, and Basin%R_o. Y_n is the depth of the marker from the surface normalised by the thickness (vertical extent) of the wedge at the location of the marker. Please see Fig. S16b for computation of Y_n.

4.4 Comparison of Easy%R_o, Simple%R_o, and Basin%R_o

The usage of Easy%R_o, Simple%R_o, and Basin%R_o in our models provides us with a distinct perspective on the comparative (dis)advantages of each method in estimating thermal maturity values. The non-uniqueness of maximum exposure temperatures for the same values of %R_o arises from the variation in sediment trajectory and thermal exposure. This diversity among sediment markers results in multiple markers attaining the same level of thermal maturity. We refer to the range of maximum exposure temperatures corresponding to similar %R_o values as the uncertainty in maximum exposure temperatures. Uncertainty for all three models increases with increasing %R_o from ∼20 to ∼25 ^∘C at approx. %R_o=0.3 to ∼35 ^∘C at %R_o=3.5 (Fig. 6b). Easy%R_o, probably the best-recognised method of thermal maturity computation, yields the best constraint on uncertainty for very small changes nearing %R_o<1. For the values of %R_o between 1 and 3, all models yield very similar uncertainty, with Simple%R_o yielding the most constrained exposure temperatures (Fig. 6b). However, beyond %R_o=3, Simple%R_o becomes unreliable, with uncertainty in exposure temperatures as high as 55 ^∘C at %R_o=4. Easy%R_o yields an uncertainty range of ∼37 ^∘C until %R_o=4.4, and starts to be unreliable above this value. Basin%R_o remains consistent until a very high value of approx. %R_o=6, and thus provides the best constraint on the widest range of values of thermal maturity (Fig. 6b).

5 Discussion

The thermomechanical models presented in this study provide (a) an explanation for the trend in thermal maturity observed in accretionary wedges, (b) a new venue to explore the uncertainty in the estimation of maximum exposure temperature using vitrinite reflectance, and (c) an estimate of the minimum lateral distance between the trench and the location of a paleo-thermal anomaly on the subduction plate for it to identified after accretion.

5.1 Thermal maturity distribution and importance of thrusting in wedges

Collectively, our results support a general increase in thermal maturity with depth and landward in accretionary wedges. The thermal maturity increase with depth is primarily the result of progressively larger exposures to higher temperatures as depth of burial increases. Conversely, the landward increase in thermal maturity is caused by the long-term deformation of sediments accumulated at older times and the exhumation of sediments that were under-thrusted as they met the backstop. Our models demonstrate that the rate of landward thermal maturity increase is faster for thicker wedges, both for the case of sediment near the surface and deep inside the wedge (Fig. 4). This can be attributed to a larger proportion of sediments being exposed to higher temperatures over an extended duration within thicker wedges, but validating this result with natural observations remains challenging given the very limited availability of thermal maturity data across natural wedges. Accretionary wedges in our models can be simplified as a system where the subducting oceanic plate acts as the primary heat source, while the seafloor acts as a heat sink. The heat generated through other sources such as shear heating, radioactivity, and advection is relatively insignificant compared to the heat originating from the younger oceanic plate. In our simulations, we consider a relatively younger and hotter oceanic plate of approximately 20 Myr, which is consistent with the accretionary wedge in the Nankai region adjacent to the Kumano forearc basin (Zhao et al., 2021). Given that the convergence rate remains constant across all models, the heat received from the oceanic plate should remain relatively similar. However, as the wedge thickness increases, the temperature gradient between the boundaries of the wedge must become gentler, resulting in a larger portion of the wedge experiencing elevated temperatures. Moreover, frequent advection from the subduction channel also results in elevated temperatures in the core of the wedge. Finally, models with thicker wedges typically exhibit higher décollement strength, leading to increased shear heating at the base of the wedge. Observational studies conducted by Yamano et al. (1992) on the thermal structure of the Nankai accretionary prism have further highlighted that the landward increase in prism thickness is the most significant factor contributing to temperature variations within the wedge. Consequently, the sustained higher temperatures within thicker wedges over time would lead to a higher rate of landward thermal maturity.

Our models show two cases where the above-mentioned trend in thermal maturity is relevantly altered, which we name “on-fault increase” and “fault-block inversion”. For instance, Fig. 3 shows a steep rise in the thermal maturity of sediments at fault sites. Thermal maturity inversions by thrusting, which are commonplace in accretionary contexts, are the primary cause of thermal maturity differentiation among wedges with similar paleo-thermal structures. During fault-block inversions, the positive gradient of thermal maturity with depth is inverted as relatively mature sediments are thrusted over less mature sediments (Underwood et al., 1992). The strong differentiation in the trajectory of sediments led by thrusting has a larger influence over thermal maturity than their burial depth or their in-wedge location. This novel inference has probably remained concealed thus far due to the large number of parameters that condition thrust development, frequency, length, and thermal state and due to the lack of high-resolution thermal maturity data.

The thermal maturity that incoming sediments reach also varies periodically as a function of thrust frequency. By examining the lateral and vertical position of incoming sediments and their eventual thermal maturity, we can deduce that the overall movement of sediments in the wedge is predominantly layered but not stationary over time. Changes in the depth of the thermal maturity boundary are less frequent and have larger amplitudes with increased décollement strength, and especially, increased sedimentation rates (Fig. 5). The periodicity in the thermal maturity boundary marks the periodic oscillation of the predominant trajectory followed by incoming sediments, i.e. between accretion (low thermal-maturity path) and under-thrusting (high thermal-maturity path). As a result, it should also strongly correlate with the periodicity observed in the evolution of forearc topography (Menant et al., 2020) and the frequency of thrust formation in our models. This is expected given that thrusts are active over longer mean times and that they channel material toward the décollement more efficiently in wedges with stronger décollement or increased sedimentation. While sediments at internal and higher structural positions of the wedge are translated toward the surface and have a lower thermal maturity, sediments at external and lower structural positions are translated toward the décollement and have a relatively higher maturity. The entire cycle is repeated with the formation of new in-sequence thrust.

Figure 7 Mapping of eventual thermal maturity (vitrinite reflectance at 7.5 Myr) to the location of same markers at ∼4 Myr in model $M_{\mathrm{0.9}}^{\mathrm{12}}$. Panel a shows the values of thermal maturity for the markers, while the lithology of the wedge is shown in panel b. The half arrow represents the active frontal thrust. The sediments which were eroded by 7.5 Myr but exist at 4 Myr have been marked as eroded using grey dots.

This is a relevant observation for it typifies the causality of particular sediment grains following a high- or low-maturity path, a long-standing unanswered question (Miyakawa et al., 2019). We corroborate this observation by analysing the terminal thermal maturity of sediments across a frontal thrust active at a younger age. An example in Fig. 7 shows the thermal maturity of sediments at ∼7.5 Myr across a thrust active at ∼4 Myr. Whereas this occurs for all thrusts in the wedge, the frontal thrust is particularly pronounced in partitioning sediments into the high- and low-maturity paths. Thermal maturity correlates with sediment depth weakly near faults and more strongly away from them. The distance of sediment from the frontal thrust dictates the trajectory of sediment grains, and as a result, the pressure–temperature conditions to which they are exposed.

Our results show the need to consider all factors influencing fault frequency when inferring the geothermal history of contractional terrains by means of thermal maturity. In this study, we have considered solely how décollement strength and the rate of trench sedimentation vary the frequency, architecture, and overall behaviour of thrusts, and the frontal thrust, as the wedge evolves. Fortunately, this predictive exercise should be relatively straightforward, for the impact of these external factors on the fault structure of wedges has been established (Fillon and van der Beek, 2012; Mannu et al., 2016, 2017; Mugnier et al., 1997; Simpson, 2010; Storti and Mcclay, 1995), and the effect of each of these factors can be accounted for when assessing the trajectory of sediments and the distribution of thermal maturity in accretionary wedges. It is nevertheless important to note that the frequency of faults in a wedge can be impacted by many other factors, including hinterland sedimentation (Storti and Mcclay, 1995; Simpson, 2010; Fernández-Blanco et al., 2020), erosion (Konstantinovskaia, 2005; Willett, 1992), and seafloor topography (Dominguez et al., 2000).

5.2 Implications of thermal maturity evolution in a subduction wedge

The main implications of this contribution emerge from its predictive power. Our approach can predict to a precise degree the thermal maturity of sediments and the uncertainty associated with the maximum exposure temperature in accretionary contexts with known structuration. A more accurate quantification of the thermal evolution and thermal state of accreted sediments reduces the uncertainties attached to the location of temperature-led transformations of organic material into hydrocarbons in subduction margins and other accretionary contexts. Such increased accuracy in the distribution of thermally mature sediments may also be applied for improved assessments of the evolution in time of any other geothermal process, including seismic slip, magmatic and metamorphic extent, porosity, compaction, and diagenesis of sediments, as well as the reconstruction of convergent margins in general (Bostick and Pawlewicz, 1984; Ferreiro Mählmann and Le Bayon, 2016; Rabinowitz et al., 2020; Sakaguchi et al., 2011; Totten and Blatt, 1993; Underwood et al., 1992).

Figure 8 Position dependency of thermal maturity preservation. (A1) Distribution of %R_o at 2.5 Myr with a paleo-thermal anomaly emplaced at 130–145 km from the backstop. (A2) The evolution of the emplaced paleo-thermal anomaly from 2.5 to 6.5 Myr in case 1. (A3) Distribution of %R_o at 2.5 Myr. (B1) Distribution of %R_o at 2.5 Myr with a paleo-thermal anomaly emplaced at 145–160 km from the backstop. (B2) The evolution of the emplaced paleo-thermal anomaly from 2.5 to 6.5 Myr in case 2. (B3) Distribution of %R_o at 2.5 Myr with a paleo-thermal anomaly emplaced at 145–160 km from the backstop.

Our simulations also imply that the paleo-thermal information stored in the incoming sediments can only be retrieved if sediments are at appropriate locations with respect to emergent thrusts. We illustrate this using two runs of the same model and tracking an artificial thermal anomaly imposed on incoming sediments at two different locations (Fig. 8). This hypothetical thermal anomaly can be conceptualised as any alteration of the thermal maturity profile of incoming sediments, for example, elevated heat flows by an antecedent magmatic intrusion. While the change in %R_o associated with the short-lived thermal anomaly results in abnormally high values of thermal maturity in both sediment packages, it can only be retrieved for the end-model run of sediments located further from the trench (those in Fig. 8B3. Contrarily, the end-model run of sediments closer to the trench (those in Fig. 8A3) shows no signs of discontinuity in the thermal maturity distribution of the wedge. This is because we deliberately placed the thermal anomaly at sites that evolve at two structural locations during the model run, i.e. above and below a yet-undeveloped frontal thrust (Fig. 8). The sediment sector affected by the thermal anomaly closer to the trench is over-thrusted by the frontal thrust and remains in a footwall location thereafter (Fig. 8A1–A3). In contrast, the homologous sedimentary package further away from the trench is accreted by the frontal thrust and remains in a hanging-wall location (Fig. 8B3). Thus, the preservation of the record of an antecedent thermal anomaly is only possible in the former case. We further note that, in our simulations, the entire vertical column of sediments records the thermal anomaly, while in nature, the anomaly may affect only sediments at the deeper locations of the sedimentary pile, which are in turn the sediments that most likely to follow a high-maturity path. We thus regard the possibility of retrieving such antecedent geothermal information as minimal.

Finally, among the three methods of %R_o computation, Easy%R_o and Basin%R_o are more consistent and well-constrained on a wide range of thermal maturities in comparison to Simple%R_o, which seems to be particularly useful for a smaller range of thermal maturity values. This simply illustrates the fact that while Easy%R_o and Basin%R_o computation deals with several parallel reactions related to the maturity of kerogen (and hence multiple activation energies), Simple%R_o is based on best-fitted single activation energy, and hence yields large confidence intervals at the extreme %R_o values. Additionally, the inclusion of the higher activation energy reactions in Basin%R_o makes it the best-suited formulation for sediments at the deeper and shear zone sediments which usually get saturated using Easy%R_o.

5.3 Comparisons to previous numerical studies

The thermomechanical models presented in this study offer a dynamic representation of trajectories within the wedge. Although the averaged trends in thermal structure and sediment trajectories remain consistent, there are short-term dynamic fluctuations near the frontal thrust. These fluctuations contribute to a diverse range of sediment paths along the depth of the incoming sediments. Miyakawa et al. (2019) conducted a similar study, modelling vitrinite reflectance using Simple%R_o and a stationary thermal field, which also resulted in an increase in thermal maturity towards the continent and thermal maturity inversions due to thrusting. However, the use of Simple%R_o led to premature saturation and the disappearance of thermal maturity variations at a shallower depth in their model.

We can compare our findings with other geodynamic models that examine the thermal structure of the wedge, although there are only a limited number of numerical models of thermal maturity in wedges. Pajang et al. (2022) recently investigated the distribution of the brittle–ductile transition in wedges and proposed a region dominated by viscous shear near the backstop, with the wedge core reaching temperatures of 450 ^∘C and typically containing forearc basins. Although trench sedimentation in our model does not result in the formation of forearc basins, the overall flattening of the wedge slope and the high vitrinite reflectance in the core align with consistent structures. Moreover, the presence of highly mature sediments in the wedge core suggests compacted sediments with greater strength and higher P-wave velocity. Although empirical studies have shown a strong correlation between P-wave velocity and thermal maturity estimates for depths of up to 4 km (Baig et al., 2016; Mallick and Raju, 1995), the exact nature of this correlation may vary depending on the specific location. Nevertheless, the patterns of thermal maturity values in the wedge core in our models also correspond to the patterns of P-wave velocity observed in the Nankai and Hikurangi margins (Górszczyk et al., 2019; Nakanishi et al., 2018; Dewing and Sanei, 2009; Arai et al., 2020).

Two modes of sediment trajectory evolution, from incoming sediment to their position inside the wedge, are generally considered: depth dependence sediment trajectories, as observed in studies by Mulugeta and Koyi (1992) and Hori and Sakaguchi (2011), and crossover exhumation pathways, as illustrated by Konstantinovskaia (2005) and Miyakawa et al. (2019). We consider the latter as non-stationary sediment trajectories that vary with time and cut across sediment trajectories of sediments previously located at the same spatial position. Our models show that both modes of sediment trajectories are valid and that changes in trajectory patterns leading to path crossovers are controlled by the horizontal distance of sediments from the frontal thrust. Starting at a threshold distance from the trench, sediments at different depths follow laminar paths along different trajectories within the wedge. Laminar-type trajectories can be reproduced in a broad range of simulations and are particularly common in models with low sedimentation and décollement strengths. However, the depth dependence of sedimentary paths varies periodically as a function of distance from the trench of specific sedimentary packages (Fig. 5). This effect, which is particularly marked in the neighbourhood of the frontal thrust, explains the crossover paths for incoming sedimentary packages at similar depths but different horizontal locations (Konstantinovskaia et al., 2005). Therefore, thrust faults in the wedge act as the primary agent controlling whether sediments sustain depth-controlled laminar flow or sediment mixing.

5.4 Comparisons to natural wedges

Our models achieve thermal maturity distributions that are in good agreement with their natural analogues, despite several relevant assumptions. Our models are very simplified with regard to their natural analogues, with assumptions such as no elasticity, predefined décollement, no erosion, and simple and uniform rheology. Also, our models have an insufficient resolution for small-scale fault activity and lack empirical relations to simulate the compaction of sediments and multiscale fluid flow. Although these assumptions hinder a wholesale comparison between our simulations and natural examples of accretionary wedges, we still find an acceptable agreement between our model and natural observations, primarily due to simulations that have a temperature evolution assimilating empirical data and a fine spatiotemporal resolution. Our estimated %R_o values for the model are in very good agreement with those measured for the borehole C0002 Nankai accretionary wedge by Fukuchi et al. (2017) (Fig. 9). The maximum exposure temperature estimated from the observed thermal maturity for the C0002 borehole also strongly correlates with maximum temperatures recorded on markers in the model with similar thermal maturity with 95 % confidence (Fig. S17). However, our result is reliant on the empirical thermal conductivity profiles estimated for the C0002 borehole, which does not show any large thermal discontinuity between the forearc basin and inner wedge that has been observed in fossil accretionary wedges (e.g. Underwood et al., 1988).

Figure 9 Depth vs. thermal maturity (%R_o). The shaded (in violet) region shows the range of observed R_o% (mean ± 1 SD) from the C0002 borehole; coloured lines represent the values in models sampled from a 10 km wide hypothetical borehole 20 km seaward of the backstop as shown in the inset.

Landward increase in thermal maturity is well documented in studies of the Japan Trench, at the Miura–Boso plate subduction margin, the fold-and-thrust belts Western Foothills complex in western Taiwan, the Mesozoic accretionary prism in the Franciscan subduction complex in northern California, and the Cretaceous Shimanto accretionary complex in Nankai subduction margin (Kamiya et al., 2017; Sakaguchi et al., 2007; Underwood et al., 1989; Sakaguchi, 1999). The natural wedges mentioned above display vitrinite reflectance values with maximum %R_o values ranging from 0.2 to 4.0 near the surface, which is generally much higher than the near-surface %R_o values observed in our models. Underwood et al. (1988) suggested that this discrepancy is likely due to the ongoing process of progressive exhumation and erosion, leading to the exposure of deeper sections of the accretionary prism over time. As a result, younger wedges, such as those found in the Miura–Boso plate subduction margin, exhibit a much closer resemblance to the %R_o values near the surface of our models.

On-fault increases in vitrinite reflectance are well also documented in nature, as for boreholes C0004 and C0007, which sample the megasplay fault in Nankai accretionary margin (Sakaguchi et al., 2011). The vitrinite reflectance data from the megasplay and frontal thrusts in Nankai indicate the faults reach a temperature well in excess of 300 ^∘C during an earthquake, much larger than the background thermal field. Therefore, on-fault increases in thermal maturity are comparatively smaller in our simulations and lack the marked increase in %R_o observed at fault sites in nature. We consider this is due to a discrepancy in the rate of change in thermal diffusion occurring in simulated thrusts, given that our models develop much wider fault zones than their natural equivalents. For instance, the location of megasplay fault in C0007 borehole exhibits an unevenness within the high-reflectance zone with a maximum of approx. %R_o=1.9 (Sakaguchi et al., 2011). This is in line with the prediction by Fulton and Harris (2012) about the impact of fault thickness on change in vitrinite reflectance. Natural observations also exhibit a much higher incidence of on-fault increase in thermal maturity compared to our simulations, given that our models do not have sufficient spatial resolution to capture the large number of thin faults that develop inside the wedge.

Natural examples of fault-block inversion have been well-documented in natural settings, providing evidence of past thrust activity preserved in the shallower sections of the Nankai accretionary wedge. Sakaguchi (1999) reported the presence of step increments of thermal maturity, similar to increments in vitrinite reflectance in Figs. 3 and 4 across the faults. Other examples are the fault block inversion along the Fukase fault in the Shimanto accretionary wedge (Ohmori et al., 1997) and the inversion beneath the forearc basin in the Nankai accretionary wedge (Fukuchi et al., 2017).

Our study highlights that paleo-thermal anomalies that extend laterally beyond the average thrust spacing have a significantly higher likelihood of being retained in the final thermal maturity record of the wedge. This allows for several inferences. For example, the subduction of the Cretaceous ridge, as identified by Underwood et al. (1992) and Sakaguchi (1999), must have caused a substantial alteration in thermal maturity during the Kula–Pacific subduction in order to be discernible in vitrinite reflectance records. Likewise, we can anticipate the preservation of the paleo-thermal anomaly near Ashizuri in the southern Nankai wedge, which has high thrust frequency, in contrast to that at the Muroto transect, where thrust sheets are widely spaced. In the case of the accretionary wedge adjacent to the Boso Peninsula, Kamiya et al. (2017) proposed the emplacement of an ophiolite complex beneath the Miura group. Our findings indicate that the preservation of the thermal-advection heating event coincided with a decrease in trench sedimentation. This likely led to an increase in the thrust frequency, which facilitated the preservation of the thermal-advection heating event in the thermal maturity data.

6 Conclusion

This study demonstrates how contractional faults alter the paths of sediments as they accrete and how this fundamentally controls the distribution of the thermal maturity of sediments in accretionary wedges and emphasises the role that sedimentation rate and interplate contact strength have in such distribution. The increased resolution of our approach leads to findings that have relevant implications. For example, the geothermal history that can be retrieved from the thermal maturity of sediments in drills, i.e. at the shallow wedge, provides, at best, an incomplete record that is skewed towards the thermal evolution of sediments near the trench. Coevally, relevant sectors of sediments located further seaward, when not subducted, follow high-maturity paths that overwrite their antecedent thermal history. Finally, this study also provides a first-order uncertainty measure for the thermal maturity of sediments based on the diversity in their trajectory.