Stress-driven fluid flow controls long-term megathrust strength and deep accretionary dynamics

The heterogeneity of frictional strength along the megathrust earthquake zone critically controls plate coupling and long-term subduction dynamics. However, the persistence and distribution of high-friction segments through space and time remain poorly constrained. Here, we show that accretion processes, such as tectonic underplating (i.e., basal accretion of material below the fore-arc region), can be used as a proxy to characterize the long-term frictional zonation of the subduction interface. We carry out numerical thermo-mechanical experiments, which predict a first-order control of tectonic-stress variations on fluid transport in deep fore-arc regions. Accordingly, positive feedback between fluid distribution and effective stress favours the stability of the interface frictional properties at Myr-scale which, in turn, controls the deep accretionary dynamics. We propose that the recognition of thick duplex structures resulting from successive underplating events over tens of Myr, allows for tracking subduction segments exhibiting an increasing frictional behaviour. Our numerical results help ascertain the long-term hydro-mechanical properties and distribution of coupling/decoupling segments of megathrust earthquake zones worldwide where active tectonic underplating is recognized.

Results from the cold subduction model. The reference experiment cold30-5-1 reproduces the subduction of an old and cold oceanic lithosphere below a 30-km-thick continental lithosphere at a convergence rate of 5 cm yr −1 . Accretion dynamics of this long-lived subduction zone is dominated by a succession of tectonic underplating events of sedimentary and basaltic material inserted at the base of the fore-arc crust as discrete tectonic slices ( Fig. 1 and Supplementary Movie S1). This leads to the formation of an ~80-km-wide dome-shaped structure, so-called duplex, after ~56 Myr of convergence. The most superficial terrigenous sedimentary layer of the subducting plate is underplated at ~16-17 km depth, whereas basaltic slices detach preferentially at ~24-25 km depth, forming two sub-duplex structures (i.e., underplating loci U Z1 and U Z2 , respectively; Fig. 1b). Interestingly, these underplating loci localize in the frontal part of persistent low-angle compressional shear zones rooting into the subduction channel (Fig. 1c). After each underplating event, the main active shear zone is successively stepping down and previously accreted slices are passively exhumed through the fore-arc crust. We stress that the along-dip positions of U Z1 and U Z2 remain remarkably stable over the entire duration of the experiment (see bar chart in Fig. 2a and Supplementary Fig. S2). At near-surface conditions, high-angle normal faults and erosion control the final exhumation and the unroofing of the duplex (Fig. 1c). Deeper down, the nappe pile gets narrower due to the stiff lower fore-arc crust. At mantle depth, restricted serpentinization atop the subduction channel and weakening of the dehydrating basaltic crust are responsible for the transient formation of a thick basaltic slice, which is not accreted to the overriding plate due to the persistence of the weak serpentinite layer (Fig. 1b,c).
As shown in Fig. 1d,e, fluid markers released from the subducting sediments and upper basaltic crust mostly migrate vertically upward throughout the duplex, creating fluid-oversaturated conditions. This transport is driven by fluid buoyancy, which prevails over low deviatoric stresses associated with passive exhumation of the duplex. In the shallow parts of the duplex, high-angle normal faults locally generate a low-stress region (i.e., tectonic underpressure) channelizing fluids until they reach the surface as localized vents (Fig. 1d). Tectonic underpressure also drives fluids downward throughout the gabbroic crust and underlying lithospheric mantle through tensile fractures related to slab bending 17 . In contrast, local tectonic overpressure patches controlled by low-angle shear zones hamper upward fluid transport atop the plate interface (Fig. 1d). Note that our model does not predict significant updip (or downdip) fluid flow along the plate interface. Near the tip of the mantle wedge corner, high compressional stresses in the fore-arc lower crust are responsible for a relative scarcity of free fluids and the transmission of shear stress to the subduction channel (Figs 1e and 2a). High stresses (up to ~300 MPa) extend deeper down into the subcontinental mantle, restricting upward fluid flow from the plate interface and, consequently, the serpentinization degree of the mantle wedge over the entire duration of the experiment (Fig. 1d). The crossing of the 300 °C isotherm activates thermodynamically-constrained metamorphic reactions 11,16 , resulting in fluid oversaturation and a significant decrease of the depth-integrated shear stress of the subduction channel at mantle depth (Figs 1e, 2a and Supplementary Fig. S2).
Fluid transport and the overall margin dynamics are consistently reproduced in additional experiments with a variable plate convergence rate (i.e., models cold30-4-1 and cold30-8-1; Supplementary Figs S3 and S4) and thicker overriding crust (i.e., model cold40-5-1; Supplementary Fig. S5), suggesting that these results are robust for a wide range of first-order parameters characterising active margins worldwide. Alternatively, variations in fluid transport properties (i.e., v perc ) result in significantly different fluid distribution and accretion dynamics (Fig. 3). Thus, a low v perc (i.e., model cold30-5-0.1; Supplementary Fig. S6) restricts fluid flow and promotes fluid-oversaturated conditions along the plate interface and immediately above, resulting in a low-friction channel and a poorly efficient fore-arc accretion. Alternatively, increasing v perc leads to a higher-friction channel and more efficient basal then frontal accretion processes (i.e., reference model cold30-5-1 and model cold30-5-10, respectively; Fig. 3 and Supplementary Fig. S7). For details on the results from these additional experiments, the reader is referred to Supplementary Text.

Modulation of fluid transport and accretion dynamics in the warm subduction model. Modelling
the subduction of a younger and warmer oceanic lithosphere (i.e., model warm30-5-1) results in a similar evolution as for the cold subduction model during the first ~20 Myr, with the formation of an early duplex by successive underplating events ( Fig. 4 and Supplementary Movie S2). Latter in the model evolution, the predominance of basal erosion over tectonic underplating along the plate interface at ~20-30 km depth leads the partial consumption of the accreted structure, which is only maintained by persistent underplating at ~15-16 km depth (underplating locus U Z1 ; Figs 2b, 4b and Supplementary Fig. S8). On top of this small duplex, the fore-arc margin is affected by a long-term subsidence and forward and backward thrusting (Fig. 4b,c). Despite such discrepancies in margin dynamics, the fluid circulation pattern in the warm subduction model is relatively consistent with the cold subduction setting (Fig. 4d,e), except that (i) thermodynamically-constrained release of slab-derived fluids occurs at shallower depth (i.e., ~25-30 km depth), resulting in a low-friction, fluid-oversaturated subduction channel near the base of the fore-arc crust (Fig. 2b) and (ii) thermal weakening of the mantle wedge allows for lower stress accumulation and vertical fluid migration up to the continental crust (Fig. 4d).

Discussion
High stress and fluid trapping. All our numerical experiments evidence the formation of a more or less developed duplex by accretion of sedimentary and basaltic slices at the base of the fore-arc crust, independently

Higher-friction segment
Higher-friction segment   www.nature.com/scientificreports www.nature.com/scientificreports/ (Central-Southern Chile) 21,22 . Geophysical imaging of duplex structures has also been reported from several actives subduction zones 23-26 (e.g., Alaska, central Japan, Cascadia, New Zealand), leading to suggest that tectonic underplating may be ubiquitous along most of active margins worldwide.
Predicted fluid distribution in the cold and warm subduction models (Figs 1e, 4e and 5) is globally consistent and supported by numerous field and seismological observations from ancient and active subduction zones, which provide robust evidences of pore fluid pressures near lithostatic values along the plate interface 27,28 . These include hydrofracture networks into fossilized subduction channels from 10 to >40 km depth 29   www.nature.com/scientificreports www.nature.com/scientificreports/ well as seismic imaging of a 2-8-km-thick, landward-dipping low-velocity zone (LVZ) with high V P /V S ratio near the mantle wedge corner along several active margins displaying contrasting thermal regimes 27 (e.g., SW Japan, Cascadia, Alaska). It is noticeable that fluid release from thermodynamically-constrained metamorphic reactions are activated at shallower depth (i.e., ~25-30 km depth) in the warm subduction model, resulting in a fluid-oversaturated subduction channel at the base of the fore-arc crust (Figs 4e and 5e), which is in line with the shallow LVZs monitored along warm active margins 31,32 (e.g., Cascadia, SW Japan). Independently from the thermal regime, these high fluid pressures must be achieved by undrained, low-permeability conditions at the top of the plate interface 31,[33][34][35] . Several hypotheses to explain this phenomenon have been proposed such as impermeable lithologies 36 , deformation-controlled grain-size reduction processes 37 or mineral precipitation 38 . According to our numerical experiments, we argue that tectonic-stress variations and associated dynamic pressure gradient on and above the plate interface exert a first-order control on fluid transport (Figs 1d, 4d and 5). Thus, the passively-exhuming duplex displays low deviatoric stresses, which allow for buoyancy-driven vertical fluid migration, in line with field observations worldwide showing abundant, nearly vertical, vein systems emplaced during the exhumation of paleo-accretionary duplexes [39][40][41] (Fig. 5a). Deeper down, the high differential stresses in the lower fore-arc crust of the cold subduction model partly close the channel at 30-35 km depth (in line with the conceptual vision of Cloos and Shreve 42 ) and avoid extensive fluid influx above the interface (Fig. 1d). The Dent Blanche massif (NW Alps) exposes a natural analogue corresponding to this situation, which has been subsequently exhumed during the Alpine orogeny 43 (Fig. 5b). Here, a strongly-hydrated 100-200-m-thick fossilized subduction interface is overlain by plastically-deformed orthogneiss almost devoid of evidence for hydrofracturing (see inset in Fig. 5b). Higher up in the Dent Blanche massif, granulite-facies rocks show numerous pseudotachylyte-bearing faults formed at ~30 km depth under differential stresses of several hundreds of MPa 44 , in agreement with estimates yield by our numerical results (Fig. 1d). At 35-55 km depth, the cold subduction model predicts a 2-4-km-thick fluid-oversaturated layer corresponding to the basaltic crust and the overlying channel (Fig. 1e). The formation of this fluid-rich zone is the consequence of persisting high stresses in the subcontinental mantle and massive fluid release via metamorphic reactions, which is supported by seismic imaging of the LVZs worldwide 27 and by field evidences of intense veining along the plate interface at similar depths 30,45 (Fig. 5c).   www.nature.com/scientificreports www.nature.com/scientificreports/ The lack of updip fluid flow along the subduction channel in our experiments apparently departs from numerical, experimental and field-based observations suggesting interface-parallel migration 43,46,47 . The discrepancy can be explained by the calculation of fluid transport in our models, which does not consider any porosity and permeability variations that are thought to be anisotropic and to occur at multiple spatial and temporal scales 46,48 . Indeed, a high-permeability subduction channel would favour along-dip fluid flow with respect to dynamic pressure gradients. However, it is noteworthy that massive upward fluid percolation from the deep subduction interface is suspected notably in warm subduction settings, based on seismically-imaged low-Poisson's ratio anomalies into the fore-arc crust 38,49 , as well as on the mantle-derived isotopic signature of mineral springs located at the surface 50,51 . This statement is supported by our modelling results, which predict an upward fluid infiltration through a persistent low-stress, thermally-weakened mantle wedge above a young and warm oceanic lithosphere (Fig. 4d). Alternatively, the maintaining of high pore fluid pressures in the deep fore-arc region is suspected from ancient and active subduction zones 27,28 , indicating important fluid accumulation over long timescales, responsible for partial serpentinization of the mantle wedge 52 . All these sensibly different observations make it difficult to decipher the respective contribution of interface-parallel versus across-interface fluid flow at Myr timescales. In view of our set of experiments, we suggest that the magnitude of the along-dip transport may be marginal with respect to the long-term fluid budget.
Tectonic underplating as a window on plate-interface frictional properties. The aforementioned heterogeneous distribution of fluid-oversaturated zones within the subduction channel in a cold regime results in along-dip shear-stress variations in the range of ~5-18 MPa with two spatially and temporally stable higher-friction regions ( Fig. 2a and Supplementary Fig. S2). A clear correlation is observed between these persistent higher-friction zones and the two U Z1 and U Z2 segments where tectonic underplating proceeds. Indeed, underplating appears to be triggered by the transition from low to higher subduction-channel shear stresses (and hence, from a wetter to a drier interface segment; Fig. 5d). This correlation is supported by similar results from models with a variable plate convergence rate ( Supplementary Figs S3 and S4), a thicker overriding crust ( Supplementary Fig. S5) or a younger and warmer subducting oceanic lithosphere (Figs 2b and 5e). Thus, it seems that the genetic link between subduction segments exhibiting an increasing frictional behaviour and tectonic slicing is robust, independently of the kinematics and the thermal regime of the active margin. However, it is noteworthy that, in the case of a warm subduction zone, the higher-friction zone near U Z1 is temporally stable while the patch near U Z2 progressively disappears ( Fig. 2b and Supplementary Fig. S8). Accordingly, margin dynamics is modified as underplating at U Z1 is overthrown by basal erosion at U Z2 , leading to the partial consumption of the former duplex structure (Fig. 4 and Supplementary Movie S2). These results support the earlier proposed conceptual model suggesting that high pore fluid pressures along the interface (and associated hydrofracturing) promote long-term basal erosion, as shown by seismic imaging at erosional margins 53,54 (e.g., Costa Rica, Ecuador). Furthermore, additional experiments show that variations in fluid transport properties strongly control the frictional state of the plate interface and, therefore, the efficiency and the depth of accretion processes (Fig. 3 and Supplementary Figs S6 and S7; see also Supplementary Text). Thus, the more drained (and hence strengthened) the subduction channel is, the more accretive the margin is, because of an increasing coupling between the subducting material and the top of the subduction channel.
The positive feedback existing between stresses and fluids 2 allows to maintain this correlation pattern between fluid distribution, high-friction subduction segments and deep accretion, both in space and time. Indeed, low-stress domains favour fluid influx which, in turn, enhances rheological weakening and hence stress lowering. Conversely, higher stresses cause fluid escape, which further increases stresses on the plate interface, ultimately triggering local underplating. The long-term consequence of this self-supported mechanism is the growth of accretionary wedges and duplex structures formed over tens of millions of years (Figs 1, 3 and 5d), strikingly recalling those identified in analogue modelling experiments 21,55 , previous numerical studies 56 and field-based conceptual models 13,22,57 . However, the natural record of very large duplexes such as the one formed after 56 Ma of convergence in our reference model (Fig. 1) are rare, as it is unlikely to maintain a stable strength regime of the interface over such a long period of time. Changes in the rheological properties of the subduction channel (e.g., subduction of topographic highs, variably-hydrated incoming plate) may lead to a switch from an accretionary to an erosional mode and hence to the basal erosion of the previously-formed duplex structure 11,53 .
A major implication from our study is that accretionary processes, and especially tectonic underplating, may be used as a proxy to ascertain the frictional properties of the plate interface over Myr timescales. Only a limited number of geophysical studies tentatively link active tectonic slicing, plate coupling and stress distribution and all of them tackle with a loss of resolution on plate interface processes deeper than 15 km 24,26,58,59 . In addition, GPS-based locking maps show a heterogeneous distribution of coupled subduction segments with a common decrease approaching the subcontinental mantle 8,10 but the link with long-term underplating is still challenging to establish. Our finding implies that regions where active duplexing is recognized could be viewed, on the short-term point of view, as the consequence of the transition from decoupled to coupled regions along the interface (Fig. 5). In active subduction settings, these patches would not be systematically visible as GPS-monitored frictional heterogeneities, which only reflect transient, short-term processes associated with the seismic cycle.
Last, variations in the frictional properties (and pore fluid pressure) along the plate interface have been invoked to explain slow slip events (SSE) updip and downdip the seismogenic zone 5,60 . Our results suggest that a rheological link between tectonic underplating regions and subduction segments exhibiting SSE must exist. Accordingly, stress must play a fundamental role for controlling fluid distribution and thus, the location of SSE phenomena along the interface. www.nature.com/scientificreports www.nature.com/scientificreports/ Methods Governing equations. The two-dimensional numerical experiments are carried out with the I2ELVIS code, which solved the continuity, momentum and heat conservation equations, based on a finite difference scheme and a non-diffusive marker-in-cell technique 15 . The continuity equation describes the conservation of mass, assuming a visco-elasto-plastic compressible fluid. It is solved on a staggered Eulerian grid and has the form: where ρ eff is the effective rock density calculated in Eq. (7), t the time, v i the viscous velocity and x i the spatial coordinates x and y. The momentum of the compressible fluid is then solved using the Stokes equation: where P is the pressure, σ ij the components of the deviatoric stress tensor and g i the gravitational acceleration (g x = 0 and g y = 9.81 m s −2 ). The heat conservation equation is formulated in a Lagrangian form to avoid numerical diffusion of temperature: where C p is the isobaric heat capacity, T the temperature, H r the radiogenic heat production, H a the adiabatic heat production, H s the shear heating and q i the heat flux solved as: where k is the thermal conductivity depending pressure, temperature and rock type (Supplementary Table S2).

Fluid implementation.
In our experiments, fluids are initially prescribed in the subducting oceanic lithosphere as (i) pore water in sediments and basaltic crust ( = X 1 w pore wt.%) and (ii) mineral bound water in sediments, basaltic crust and gabbroic crust. Pore water release by compaction is assumed constant from 0 to 75 km depth, also mimicking dehydration from low-temperature metamorphic reactions 29 (e.g., smectite-illite, opal-quartz transformations). This linearly-decreasing pore water content accounts for the kinetics of metamorphic reactions and the natural heterogeneity of rocks that result in distributed fluid release over a temperature range, rather than in discrete pulses 2,61 . Bound water release is calculated by free-energy minimization 16 as a function of pressure, temperature and rock type 11 . Resulting free water is then transported as newly-formed Lagrangian markers. Importantly, our models do not take into account the coupling between pore pressure and hydraulic properties (i.e., porosity and permeability) as considered in real two-phase flow systems 62 . Instead, we make the following assumptions: (i) fluid pressure approaches the solid dynamic pressure and (ii) hydraulic properties remain constant due to the balance between porosity enhancement and destruction 17 . Accordingly, fluid marker transport is controlled by the Stokes velocity field v i , the fluid buoyancy and the dynamic pressure gradients, such that: where v i water is the velocity of fluid markers, v perc the reference percolation velocity (varying from 0.1 to 10 cm yr −1 in our experiments; see Supplementary Text and Supplementary Table S1 for details on the parametric study) and k i the pressure-dependent coefficient calculated as: Fluid effects on rock density and rheology. In our experiments, bound water content in rock markers affects rock density as follow: rock 0 www.nature.com/scientificreports www.nature.com/scientificreports/ where ρ 0solid is the standard density of rocks, X fluid the mass fraction of fluid, α the thermal expansion and β the compressibility (see Supplementary Table S2 for a complete list of rock properties). Non-Newtonian visco-elasto-plastic rheologies are employed in these experiments, implying that the deviatoric strain rate tensor ε ij includes the three respective components: Details on the calculation of the rheological constitutive equations are available in ref. 15 . For our purpose, it is important to note that free fluids affect rock rheology by modifying the plastic strength σ yield , which limits the creeping (i.e., viscous) behaviour such that:  (Fig. S1b). The oceanic crust is initially underthrusted below the continental margin and a 10-km-thick weak zone is prescribed at the interface between the two plates to initiate subduction. The thermal structure of the oceanic lithosphere is initially defined by an oceanic geotherm with a cooling age evolving from 10 kyr (x = 0) to 53 Myr (x = 854 km) for the cold subduction setting (model sub53-1). To limit the size of the computational domain, the cooling of the oceanic lithosphere is prescribed as 10 times faster for 0 ≤ x ≤ 200 km. This high cooling zone is located at ~600 km away from the subduction zone, which allows to avoid any thermal or mechanical effect on modelled fore-arc dynamics. In the case of a warm subduction setting (model sub20-1), the maximum cooling age of the oceanic lithosphere is 20 Myr and its cooling is only 2 times faster for 0 ≤ x ≤ 200 km. A geothermal gradient of ~15° km −1 down to 90 km is defined for the continental lithosphere. Below, the asthenospheric geothermal gradient is 0.5° km −1 .
Velocity boundary conditions are free slip for the left, right and upper boundaries, while the lower boundary is open to ensure mass conservation in the computational domain. The top of the lithosphere is calculated as an internal free surface by using a low-viscosity layer simulating air (y < 10 km) or water. Resulting large viscosity contrast minimizes shear stresses at the top of the lithospheres, which efficiently approximates a free surface 64 . In our experiments, sedimentation and erosion processes are implemented independently by applying the following equation at the surface 65 : where y surf is the y coordinate of the surface, v x and v y the horizontal and vertical velocity components of the Stokes velocity field at the surface and v erosion and v sedim the global erosion and sedimentation rates, respectively, defined as (i) v erosion = 0.3 mm yr −1 and v sedim = 0 mm yr −1 for y < 10 km and (ii) v erosion = v sedim = 0 mm yr −1 for y > 10 km. A increased erosion/sedimentation rate of 1 mm yr −1 is applied to regions with steep surface slopes (i.e., > 17°) for smoothing the topographic surface. This is particularly relevant for the offshore fore-arc region where the increased sedimentation rate counterbalances the absence of global sedimentation rate prescribed in our experiments. Newly-formed sedimentary rocks are labelled as terrigenous sediments and display the same properties than the pelagic sediments (Supplementary Table S2).