3-D thermal regime and dehydration processes around the regions of slow earthquakes along the Ryukyu Trench

Several interplate seismic events, such as short-term slow slip events (S-SSEs) and low-frequency earthquakes (LFEs), have been identified in the Ryukyu Trench, southwestern Japan. As one of the specific characteristics of this seismicity, the depths at which S-SSEs occur at the plate interface beneath Okinawa Island are approximately 5–10 km shallower than those beneath the Yaeyama Islands. To elucidate the cause of this difference in depth, we constructed a three-dimensional, Cartesian thermomechanical subduction model and applied the subduction history of the Philippine Sea (PHS) plate in the model region. As a result, the interplate temperatures at which S-SSEs take place were estimated to range from 350 to 450 °C beneath Okinawa Island and from 500 to 600 °C beneath the Yaeyama Islands. The former temperature range is consistent with previous thermal modelling studies for the occurrence of slow earthquakes, but the latter temperature range is by approximately 150 °C higher than the former. Therefore, explaining how the depth difference in S-SSEs could be caused from the aspect of only the thermal regime is difficult. Using phase diagrams for hydrous minerals in the oceanic crust and mantle wedge, we also estimated the water content distribution on and above the plate interface of the PHS plate. Near the S-SSE fault planes, almost the same amount of dehydration associated with phase transformations of hydrous minerals from blueschist to amphibolite and from amphibolite to amphibole eclogite within the oceanic crust were inferred along Okinawa Island and the Yaeyama Islands, respectively. On the other hand, the phase transformations within the mantle wedge were inferred only beneath the Yaeyama Islands, whereas no specific phase transformation was inferred beneath Okinawa Island around the S-SSE occurrence region. Therefore, we conclude that dehydrated fluid derived from the oceanic crust at the plate interface would play a key role in the occurrence of S-SSEs.

. This mobility of the source of LFEs would indicate the existence of fluid near the plate boundary 6 . By incorporating a modified network inversion filter into GNSS time series data, the spatiotemporal source processes of five S-SSEs in the Yaeyama Islands region during the period from March 2010 to February 2013 were revealed, and LFE and VLFE activity was inferred to have initiated after the evolution of S-SSEs 7 . Considering the source mobility of LFEs and the synchronization among their occurrences, S-SSEs would activate the LFEs and VLFEs associated with fluid transport near the plate interface. The thermomechanical subduction modelling in the preceding study implied that the transition in the slip mode may be constrained by temperature-pressure conditions at the plate interface 8 . Previous studies discussing the relationship between permeability and fault strength [9][10] suggested that low permeability in the hanging wall results from serpentinization or deposition of silica due to dehydration from the slab near the mantle wedge corner, which would prompt the formation of a high pore fluid pressure region at the plate interface and reduce the effective normal stress there. Then, episodic tremor and slip (ETS) may occur at an isolated friction zone within this region. Low effective normal stress in the ETS zone resulting from high pore pressure and chemical dehydration reactions around the plate interface may induce semi-brittle deformation of serpentinite, which is governed by both brittle and ductile responses and has a low healing rate 11 . A variation in the slip behaviour from velocity strengthening to velocity weakening has been identified there, which is reasonable for the slip mechanism of slow earthquakes 12 . Another friction model 13 also supports the presence of the high pore pressure region to reproduce realistic properties of slow slip events. In addition to frictional properties, the plastic flow of antigorite-rich serpentinite influences the dynamics of fault slip 14 .
To investigate the similarity or difference between the Okinawa Island and Yaeyama Islands regions with respect to the thermal regimes and the dehydration processes near the S-SSEs and LFEs, we performed several numerical simulations of temperature fields associated with the subduction of the PHS plate using a 3-D Cartesian thermomechanical subduction model and estimated the water content distributions on and above the subducting PHS plate using phase diagrams for hydrous mid-ocean ridge basalt (MORB) of the oceanic crust and ultramafic rocks of the mantle wedge. We took certain parameters as unknowns and tested their sensitivity to the thermal regime. Then, we determined the most suitable model to fit the calculated thermal regime with the observed heat flow data. Following Ji and Yoshioka (2015) 15 , we constructed a 3-D Cartesian thermomechanical subduction model and performed numerical simulations during the period from 15 Ma to the present (0 Ma) to calculate the spatiotemporally changing temperature structure and mantle flow velocity field associated with the subduction of the PHS plate. The dimensions of the model were 500 km, 800 km, and 400 km in the x-, y-, and z-axis directions, respectively (Fig. 2). The numbers of nodes used were 50, 80, and 80 in the x, y, and z directions, respectively. The geometry of the subducting PHS plate was created by combining the Slab 2.0 model 16 17 , which is based on a seismic velocity structure survey around the Yaeyama Islands. The subduction of the PHS plate was modelled kinematically, and the PHS plate was constrained to subduct along a prescribed guide with an imposed convergence rate (Fig. 2). The model consisted of the subducting oceanic plate (within the prescribed guide), continental upper crust, continental lower crust, accretionary prism, and mantle (Fig. S1). The domains of the upper crust, lower crust, and accretionary prism were rigid; the other domains were composed of convective regions. For the tectonic setting associated with the subduction of the PHS plate along the Ryukyu Trench, we referred to a global plate rotation model 18 (Fig. S2) and a trench motion model 19 . From 15 to 3 Ma, the PHS plate moved in a northeasterly direction. Subsequently, the PHS plate changed its motion to a northwesterly direction at 3 Ma. In addition, the Ryukyu Trench retreated southeastward 19 ; the Okinawa Trough, which is an active backarc basin of the Ryukyu arc, initiated spreading associated with subduction of the PHS plate at approximately 2 Ma 20-21 . The spreading rate increased towards the southwest 22 , although its exact value from 2 Ma to the present is unclear. In this study, we evaluated the validity of the estimated thermal structures by comparing calculated heat flows with observed heat flows. For heat flow data, we used the Global Heat Flow Database (GHFD) 23 and land boreholes [24][25] . To compare the observed heat flows with the calculated heat flows, we used the following equation:

Model and tectonic setting
where hf obs i and hf cal i are the ith observed and calculated heat flow values, respectively, and N is the total number of data points from land boreholes and marine heat probes within the model region. We constructed a best-fit model to minimise S in Eq. (1).
We investigated several parameter combinations of the trench retreat rate and initiation age of the trench retreat, and compared the calculated heat flow values with the observed heat flow values to evaluate the parameter dependency on the trench retreat effect. For the parameter search range, the trench retreat rate was assumed to range from 15 to 20 mm/year 19 , and the initiation age of the trench retreat was set to range from 8 to 5 Ma, following a previous thermal modelling study in this region 26 . Furthermore, we considered the spreading rate of the Okinawa Trough, which ranged from 5 to 20 mm/year 22 uniformly in the + x direction (Fig. 2) from 2 Ma to the present, including the initiation of rifting in the Okinawa Trough. In this study, this effect was dealt with by simply increasing the trench retreat rate of the Ryukyu Trench only from 2 Ma to the present as an additional velocity constraint in the numerical simulation. All parameters for the physical properties and subduction history of the PHS plate are tabulated in Tables S1 and S2, respectively.   (Fig. S3), we obtained optimal values for the retreat rate of the Ryukyu Trench, spreading rate of the Okinawa Trough, and initiation age of trench retreat: 20 mm/year, 15 mm/year, and 7 Ma, respectively. This model is hereafter referred to as a best-fit model. In the best-fit model, the interplate temperature for the occurrence of S-SSEs and LFEs beneath Okinawa Island ranges from 350 to 450 °C (Fig. 3a). This temperature range for the occurrence of S-SSEs is considered to correspond to the temperature at which the slip mode changes from unstable sliding to stable sliding 8 . On the other hand, temperatures beneath the Yaeyama Islands range from 500 to 600 °C (Fig. 3a), and this range is not very consistent with the above-mentioned transition zone for the slip mode. However, it is not possible to discuss in detail how significant this temperature difference is, taking account of stress field which is not calculated in our kinematic subduction model. The temperature-depth gradients along Okinawa Island (Fig. 4b) are higher than those along the Yaeyama Islands because of the higher dip angle at depths greater than 40 km. Incidentally, the interplate temperature and the mantle wedge temperature, which is taken 5 km above the plate interface, decrease by at most 200 °C when incorporating the effect of trench retreat and the recent spreading rate of the Okinawa Trough into the model; these changes are due to the development of lateral flow in the mantle wedge. The influence of lateral flow in the mantle wedge tends to be stronger after the initiation of spreading of the Okinawa Trough at 2 Ma.
Additionally, the low-temperature region near Okinawa Island in Fig. 3a results from the older seafloor within the model region (Fig. S2(d)). On the other hand, the subduction velocity around the Yaeyama Islands is faster than that around Okinawa Island. Thus, the cold slab can easily reach greater depths around the Yaeyama Islands. The former and latter effects have a trade-off factor, although the former effect has a greater influence on the thermal regime, considering the lower temperature-depth gradient around Okinawa Island (Fig. 4a). However, lateral mantle flow in the along-arc direction associated with oblique subduction could also be a factor causing the difference in interplate temperature at the same depth between the two regions. However, such flow almost

Water content distributions.
We also calculated the water content distributions near the plate interface using the calculated thermal structure and phase diagrams for hydrous minerals in the oceanic crust and the mantle wedge, which were obtained from the Perple_X program 27 (Fig. 4a,b). We assumed a sharp dehydration transition when a certain hydrous mineral within the oceanic crust or the mantle wedge transformed into another higher P-T mineral. The results show that the shallower region where S-SSEs and LFEs occur beneath Okinawa Island is adjacent to the deeper region where the phase transformation from blueschist to amphibolite within the oceanic crust takes place (Figs. 3b and 4a). In contrast, the phase transformation from amphibolite to amphibole eclogite within the oceanic crust takes place close to the S-SSE fault plane around the Yaeyama Islands (Figs. 3b and 4a). We also investigated a temperature-depth path at a depth of 5 km above the slab surface (Fig. 4b). Consequently, phase transformations from brucite to antigorite and antigorite to chlorite in ultramafic rocks were identified beneath the Yaeyama Islands (Figs. 3c and 4b). The phase transformation from brucite to antigorite was identified at depths 5-10 km deeper than the S-SSE occurrence region along profile A, passing through Okinawa Island (Figs. 3c and 4b).
Comparing our results with the V p /V s ratio estimated from seismic tomography 28 (Fig. S4), we found that the region with a V p /V s ratio appropriate for the occurrence of slow earthquakes (1.75 < V p /V s < 1.85 29 ) corresponded well with the region where the phase transformation from brucite to antigorite within the slab mantle was identified in our results. In this study, we did not introduce a particle trace system for slab-derived fluid, but some preceding studies performed numerical simulations and showed that wet, slab-derived water or slab melt upwelled, passing through the mantle wedge [30][31] . In summary, almost the same amount of slab-derived water was supplied near the S-SSE occurrence regions beneath both Okinawa Island and the Yaeyama Islands and was associated with phase transformations in the oceanic crust despite a depth difference in the occurrence of S-SSEs. The phase transformations within the mantle wedge along Okinawa Island were shallower than those along the Yaeyama Islands (Fig. 4b) because the steeper dip angle along Okinawa Island may have led to mantle upwelling. However, phase transformations within the mantle wedge occurred in the region slightly deeper than the S-SSE fault planes.
In our future work, there are some challenges for improving our model to make it more realistic. For example, we need to trace the particle motion of the derived fluid, which would change the pore fluid pressure on the plate interface and the spatiotemporal frictional properties [9][10] . In addition, it is important to consider the changes in the physical properties of mantle rocks associated with their mineral metamorphism, such as serpentinization 9,14 .  27 . The red-and blue-solid lines denote the temperature-depth paths at the surface of the oceanic crust of the PHS plate along profiles A and B (Fig. 1), respectively. The thick red and blue lines on the temperature-depth paths denote the temperature-depth ranges of the occurrence of S-SSEs along profiles A and B, respectively. BS, LwsBS, AMP, AmpEC, ZoEC, LwsEc, DryEc, and GR represent the phases of blueschist, lawsonite blueschist, amphibolite, amphibole eclogite, zoisite eclogite, lawsonite eclogite, dry eclogite, and granulite, respectively. (b) Calculated temperature-depth paths on the phase diagram of hydrous ultramafic rock 27 . The red-and blue-solid lines denote the temperature-depth paths passing through a depth of 5 km above the surface of the PHS plate along profiles A and B (Fig. 1), respectively. The others are the same as in (a). www.nature.com/scientificreports/ We also need to introduce the effects of chemical reactions associated with equilibrium or non-equilibrium phase transformations considering the latent heat effect on solid-solid reactions in crustal and mantle rocks. As pressure and temperature conditions change, the physical properties of crustal and mantle rocks are affected by structural changes in the crystal lattices of minerals and the absorption or release of latent heat associated with their phase transitions 32 .

Conclusions
In this study, we performed three-dimensional thermomechanical subduction modelling of the Ryukyu subduction zone to investigate the cause of depth differences in S-SSEs and LFEs in terms of temperature field and dehydration processes. The results of this study can be summarised as follows: (1) The most suitable parameter values to explain the observed heat flow data for the retreat rate of the Ryukyu Trench, spreading rate of the Okinawa Trough, and initiation age of trench retreat are estimated to be 20 mm/year, 15 mm/year, and 7 Ma, respectively. (2) The interplate temperatures where S-SSEs take place range from 350 to 450 °C and 500 to 600 °C beneath Okinawa Island and the Yaeyama Islands, respectively. The latter range is approximately 150 °C higher than the former. Because the effect of older seafloor around Okinawa Island has a larger influence on the thermal regime than the increasing subduction velocity towards the southwest, the interplate temperature around Okinawa Island is more likely to become colder than that around the Yaeyama Islands.

Methods
In this study, we used an anelastic liquid approximation. For the governing equations, we used mass conservation, momentum, and energy equations to calculate the temperature and mantle flow velocity fields. The mass conservation equation can be written as where ρ s (z) and v are the fluid density and flow velocity vectors, respectively. The subscripts denote the adiabatic state. Following Ji et al. 33 , the momentum equation can be written as follows, where P is the pressure deviation from the hydrostatic pressure, τ ij is the stress tensor, R a is the Rayleigh number, α is the thermal expansivity, T is the temperature, and δ ij is Kronecker's delta. T can be written as where g is the acceleration of gravity, D is the thickness of the model, and C p is the specific heat at constant pressure. The energy equation is expressed as where t is the time, k is the thermal conductivity, v is a vertical component of the flow velocity, H r is the heat generation by a radioactive element, and Q f is the frictional heating on the plate interface. The density ρ is dependent only on temperature. Following Burkett and Billen 34 , we use viscosity, which is represented by a composite of diffusion creep and dislocation creep: with where η diff and η disl are the viscosities of diffusion and dislocation creeps, respectively, d is the grain size, p is the grain index, A is a coefficient, C OH is the water content, r is the water content index, n is the stress index, ε is the second invariant of the strain rate tensor, E is the activation energy, P lc is the lithostatic pressure defined by the depth gradient of compressibility, V is the activation volume, and R is the gas constant. where τ b and τ d are the shear stresses at the plate boundary in the brittle and ductile regimes, respectively, v t is the subducting velocity, and w is the thickness of the plate boundary deformation zone. The shear stress τ b in the brittle regime 36 can be written as where σ n is the normal stress and is the pore pressure ratio. In the best-fit model (Figs. 3 and 4), the pore pressure ratio is estimated to be 0.995. The shear stress in the ductile regime 37 can be written as where A is a coefficient, ζ is the shear strain rate and E ′ is the activation energy. By solving Eqs. (2), (3), and (5) as a coupled problem, we can obtain unknown parameters v , P , and T at each time step.
For the boundary condition for the mantle flow velocity, we imposed permeable conditions at the + x, −x, + y, −y, and + z planes (Fig. 2). When trench retreat was introduced into the model, the model surface was divided into two parts bounded by the time-dependent moving trench axis. On the oceanic-plate side, the oceanward retreat velocity was given on the model surface to realize trench retreat. On the other hand, on the continental-plate side, a flow velocity of (0, 0, 0) was given. As a result, a gap arises between the edge of the continental plate and the trench axis moving away with time. Therefore, we considered a growing accretionary prism with a triangular prism shape increasing with time during the trench-retreat period. The growth rate of the accretionary prism is set to be equal to the trench retreat rate to compensate for the accretionary-prism materials between the trench axis and the continental plate.
For the boundary condition for the temperature fields, we assigned adiabatic conditions at the − x, + y, − y, and + z planes. A fixed temperature of 0 °C was imposed at the − z plane, and the temperature T(z, t) at the + x plane was determined by the plate cooling model 38 : where T 0 and T 1 are the surface and basal temperatures of the lithosphere, respectively, z is the depth, d l is the depth at which the mantle temperature reaches T 1 , n is an index, κ is the thermal diffusivity, and t is the age of the PHS plate at the Ryukyu Trench. We used the RT1 model 39 , in which values of T 1 and d l are given as 1402 °C and 134.9 km, respectively.
The initial condition in the model domain is a steady state without mantle flows. Following Yoshioka and Sanshadokoro (2002) 40 , the temperature distribution for the half-space cooling model can be expressed as and where T 0 is a potential temperature, t ′ is the age of the continental plate, and z 0 is the depth below which the effect of adiabatic compression is considered. When z ≥ z 0 , we used Eq. (14) instead of Eq. (13), taking into account adiabatic compression.