Thermal structure of continental subduction zone: high temperature caused by the removal of the preceding oceanic slab

The thermal structure of the continental subduction zone can be deduced from high‐pressure and ultra‐high‐pressure rock samples or numerical simulation. However, petrological data indicate that the temperature of subducted continental plates is generally higher than that derived from numerical simulation. In this paper, a two‐dimensional kinematic model is used to study the thermal structure of continental subduction zones, with or without a preceding oceanic slab. The results show that the removal of the preceding oceanic slab can effectively increase the slab surface temperature of the continental subduction zone in the early stage of subduction. This can sufficiently explain the difference between the cold thermal structure obtained from previous modeling results and the hot thermal structure obtained from rock sample data.


Introduction
The thermal structure of subducted plates plays a dominant role in negative driving buoyancy force, phase transition, dehydration and partial melting during subduction (Peacock, 1996(Peacock, , 2003van Keken et al., 2011). These processes can cause large earthquakes and volcanic activities in subduction zones (Hacker et al., 2003;Grove et al., 2012). Therefore, understanding the thermal structure of subducted plates has always been a focus of plate tectonics. Analytical and numerical methods have been used to study the thermal structure of slabs (McKenzie, 1969;Molnar and England, 1990;Peacock and Wang KL, 1999;van Keken et al., 2002;Syracuse et al., 2010). The results show that the thermal structure of slabs is mainly affected by conduction rate (McKenzie, 1969;Molnar and England, 1990). Higher vertical subduction velocities and older plate ages result in colder slab thermal structures, and vice versa (Peacock and Wang KL, 1999;Syracuse et al., 2010). Other factors such as phase transition, partial melting, and geometry, may also lead to the variations of slab thermal structure in differ-ent subduction zones (Leng W and Mao W, 2015 and references therein;Leng W and Huang LZ, 2018).
Previous studies have focused on oceanic subduction zones (van Keken et al., 2002;Syracuse et al., 2010), partly because most of the subduction zones at present are oceanic subduction zones. However, along the Tethys orogeny, from the Alps to the Iranian Plateau to the Qinling-Dabie, continental subduction zones have been widely observed since the Mesozoic (Zheng YF, 2012;Zheng YF and Chen YX, 2016). Therefore, it is crucial to establish a model of the continental subduction zone and analyze the thermal structure of the subducted plate. In particular, high-pressure (HP) and ultra-high-pressure (UHP) metamorphic rock samples indicate that plate thermal structures in continental subduction zones are generally hotter than predicted by numerical models at depths less than 120 km (Penniston-Dorland et al., 2015). HP/UHP rocks provide an ultra-hot slab thermal gradient of up to 20°C/km, while the cold slab thermal gradient predicted by the numerical model is less than 10°C/km (Gerya et al., 2002;Syracuse et al., 2010). This difference between the experimental and numerical results needs to be resolved. (Kohn et al., 2018). However, other studies challenge this conclusion (van Keken et al., 2019) and propose that these HP/UHP rocks are preferentially exhumed from warm subduction environments (van Keken et al., 2018). On the other hand, during the transitional period from oceanic subduction to continental subduction, it is generally proposed that oceanic slab breakoff occurs before subduction of continental plate (Davies and von Blanckenburg, 1995). Based on seismic observations (Kounoudis et al., 2020), laboratory experiments (Fernández-García et al., 2019) and numerical simulations (Pusok et al., 2018), slab breakoff could effectively remove the preceding oceanic slab at the transition from oceanic to continental subduction. Furthermore, the numerical model results show that when the oceanic slab tip first pierces the specific depth horizon, the slab surface temperature reaches a maximum value, then gradually decreases with time and approaches to a steady value (Kincaid and Sacks, 1997). Therefore, the removal of the preceding oceanic slab can effectively increase the slab surface temperature of a continental subducted plate before it reaches a steady state. This dynamic effect may explain the high temperature observed in HP/UHP rock samples, which is worthy to be further explored.
Here we use a numerical model to study the thermal structure of the continental subduction zone. Specifically, we focus on the thermal structure of the continental subduction zone without a preceding oceanic slab and compare it with the results from HP/UHP metamorphic rock samples.

Method
We use the two-dimensional finite element code CITCOM to solve the conservation equations of mass, momentum and energy (Moresi and Solomatov, 1995). The geometry of the model is shown in Figure 1a. In the horizontal and vertical directions, the box size is 1143 × 600 km, and the grid resolution is 705 × 641.
Ours is a kinematic model driven by an imposed subduction velocity which does not include the thermal buoyancy. Compared with the fully dynamic model, the kinematic model allows us to better control the subduction process, such as the subduction ve-locity and dip angle. Besides the imposed subduction velocity, the top boundary is no-slip and other boundaries are stress-free which allows material to flow in and out ( Figure 1a). Since the slab thermal structure is mainly controlled by the vertical velocity of the subduction (Peacock and Wang KL, 1999;Syracuse et al., 2010), we chose to fix the dip angle of θ = 30° in all models, and only vary the subduction velocity. Following Syracuse et al. (2010), the decoupling point between the subducted plate and the overriding plate is set at 80 km depth, where the flow velocity becomes continuous between the slab and the mantle wedge. We benchmark this subduction model with the results of van Keken et al. (2008).
We have determined the initial temperature profiles of oceanic and continental subduction zones. For the oceanic slab preceding the continental plate, the initial temperature is calculated using a half-space cooling model (Turcotte and Schubert, 2002). For the subducted continental plate, the calculation method of the initial temperature before subduction is similar to that of Syracuse et al. (2010): the crust is divided into upper crust and lower crust with different thickness (h 1 = 7 km, h 2 = 30 km), density (ρ 1 = 2700 kg/m 3 , ρ 2 = 2900 kg/m 3 ) and radiogenic heating rate (H 1 = 1.3 μW/m 3 , H 2 = 0.27 μW/m 3 ). Then, given the surface heat flux q s , the temperature T at depth y can be calculated as: where k = 3 W/mK is the thermal conductivity. Note that the temperature of the lithospheric mantle is calculated from the heat flux at the Moho surface and linearly increases to the mantle temperature T 1 = 1300°C. Figure 1b shows the obtained initial temperature profiles of oceanic plates with different plate ages and contin- (3) Fixed subduction velocity  3) and (4) represent the preceding oceanic slab, subducting continental plate, overriding plate and mantle wedge, respectively. For areas below the slab surface, the velocity is fixed as the subduction velocity. The top boundary is fixed with zero velocity. Other boundaries are stress-free which allows material to flow in and out. θ = 30° is the dip angle. (b) The initial thermal structure of continental subducting plates with different surface heat fluxes (solid lines) and oceanic subducting slabs with different slab ages (dashed lines). ental plates with different surface heat fluxes. The temperature profile for the overriding plate is calculated in the same way, with the averaged continental heat flux q s = 65 mW/m 2 .
Model viscosity is temperature-dependent: η 0 where = 1.8×10 21 Pa·s is the reference viscosity; T 0 = 1573 K is the reference temperature; E = 335 kJ/mol is the diffusion activation energy, where the dislocation creep is ignored since our model is kinematic; R = 8.31 J/mol K is the gas constant.

Results
We are most interested in the effect of oceanic slab breakoff on the thermal structure of the continental subduction zone. Therefore, during continental subduction, we run two types of models, with or without the preceding oceanic slab.
Cases A1 and A2 are the two cases with the preceding oceanic slab. The age of the oceanic slab is 40 Ma and the subduction velocity is 5 cm/yr. After 11.9 million years (Myr) of subduction, the front end of the oceanic slab reached a depth of 300 km (Figure 2a), and the slab surface temperature at the top surface (SST) reached a steady state (Figure 3a, black curve). It can be observed that SST of oceanic subduction approaches the thermal gradient of 5°C/km in the shallow part, and then rises after decoupling from the overriding plate at the depth of 80 km. This thermal structure of the cold oceanic slab is consistent with previous numerical results (Syracuse et al., 2010;Penniston-Dorland et al., 2015). At 12 Myr, oceanic subduction transited to continental subduction. The surface heat flux q s of the subducted continental plate is set as 30 mW/m 2 , which represents the cold craton plate. Based on geological observations, due to the continental plate's buoyancy, the transition from oceanic subduction to continental subduction may significantly reduce the subduction velocity, so we adopt two cases: the continental subduction velocity is set at 1 cm/yr (case A1) and 5 cm/yr (case A2). For case A1, after another 40 Myr of subduction, the front end of the continental plate reaches a depth of 200 km (Figure 2b). During the transition to continental subduction, the reduced subduction velocity of 1 cm/yr resulted in a slight increase in SST compared with the preceding rapid oceanic subduction (Figure 3a, red curve). Then, after the transition, the SST gradually decreases to a steady state (Figure 3a, green and blue curves). As a comparison, for case A2, we maintain the same continental subduction velocity as the preceding oceanic subduction velocity of 5 cm/yr. During the transition from oceanic subduction to continental subduction, SST gradually decreases to a steady state (Figure 3b). These two cases show that the thermal structure of the continental subduction zone is similar to that of the oceanic subduction zone in the presence of the preceding oceanic slab.
Cases B1 and B2 are the same as cases A1 and A2, except that we removed the preceding oceanic slab and let the continental plate subduct at the beginning. It can be regarded as a slab breakoff event at the transition stage from oceanic subduction to continental subduction (Davies and von Blanckenburg, 1995). Figure 2c and 2d show the temperature and velocity fields of case B1 at 19.2 and 40.2 Myrs. The key difference between cases B1 and A1 is that, without the preceding oceanic slab in case B1, the SST is much higher when the continental plate first subducted to a shallow depth (Figure 3c, black, red and green curves). The continental plate directly penetrates the hot mantle and forms a thermal structure greater than 10°C/km. After a long period of evolution, the SST gradually cools down, and finally reaches a steady state (Figure 3c, blue curve), which is very close to the steady state SST of case A1 with the preceding oceanic slab. When the continental subduction velocity is large (5 cm/yr in case B2), the SST results show a similar trend (Figure 3d). That is, in the absence of the preceding oceanic plate, the continental subduction zone exhibits a hot thermal structure when the plate first subducts to a shallow depth. Only when the continental plate subducted to deep depths (>200 km), can its thermal structure approach a cold steady state. The final cold steady states are the same for cases with or without the preceding oceanic slab (comparing blue curves in Figure 3a and 3c, or in Figure 3b and 3d).
All of these above cases use a cold initial temperature profile for the subducted continental plate. We increased surface heat flux q s from 30 to 65 mW/m 2 and ran a new case C1. This case reflects the subduction of the normal continental lithosphere at a velocity of 1 cm/yr without the preceding oceanic slab. Figure 4 shows the evolution of SST with time in case C1. We also plot the pressuretemperature data of the HP/UHP rock samples from different continental subduction zones (Penniston-Dorland et al., 2015). It can be observed that the modeled continental subduction plate shows a hot thermal structure before approaching the final cold steady state, especially at shallow depths where the front end of the continental plate first reaches (Figure 4, black, red and green curves). Most rock data can be covered by our modeled temperature profiles at different evolutionary times. Our model predicts that the hot thermal structure (>10°C/km) of the continental subduction zone can only occur at shallow depths; this is consistent with the rock sample data. When the continental plate subducted to a deeper depth, its thermal structure became cold and finally   reached a steady state. Therefore, a subducted continental plate without the preceding oceanic slab may well explain the discrepancy between the cold thermal structure from the previous modeling results and the hot thermal structure from the HP/UHP rock sample data: the removal of the preceding oceanic slab due to slab breakoff can effectively increase the slab surface temperature during the early stage of continental subduction.

Discussion and Conclusion
We did not consider shear heating effects, which are believed to significantly increase the slab temperature at shallow depths (Kohn et al., 2018). Compared with the shear heating model, our model proposes that the hot continental slab temperature is a transient feature of the initial stage of continental subduction. For a specific continental subduction zone, if more accurate ages of rock samples can be determined which show the transition from hot to cold thermal structure, our model will become more credible.
Our model assumes that the preceding oceanic slab completely breaks at the beginning of continental subduction. If the preceding oceanic slab did not break until the continental plate subducted to a deep depth, then the thermal structure of the continental plate would be similar to previous studies showing a cold temperature profile (Figure 3a, 3b).
In conclusion, we use a two-dimensional kinematic numerical model to study the thermal structure of the continental subduction zone. We find that the shallow thermal structure of the sub-ducted continental plate can be hot in the early stage of subduction without a preceding oceanic slab. This explains the observed hot thermal structure from HP/UHP rock samples. Our results provide a new view for the thermal tectonic evolution of the continental subduction zone with a slab breakoff event.