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

Armel Menant, Samuel Angiboust, Taras Gerya

► To cite this version:

Armel Menant, Samuel Angiboust, Taras Gerya. Stress-driven fluid flow controls long-term megathrust strength and deep accretionary dynamics. Scientific Reports, 2019, 9, 10.1038/s41598-019-46191-y. hal-03508181

HAL Id: hal-03508181

https://hal.science/hal-03508181

Submitted on 3 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_2.png" /></figure>

OPEN

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

Armel Menant<fnref n="1" />,², Samuel Angiboust<fnref n="1" /> & Taras Gerya<fnref n="1" />,<fnref n="3" />

Received: 22 January 2019

Accepted: 17 June 2019

Published online: 04 July 2019

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.

Subduction fluids are intricately linked to tectonic processes<fnref n="1" />,² and have paramount effects on the frictional strength of subduction interfaces and, therefore, on the seismogenic behaviour of megathrust faults<fnref n="3" />⁻⁵. Based on seismological and geodetic observations, interfaces appear highly segmented with a heterogeneous distribution of coupled and decoupled patches due to pore fluid pressure variations⁶⁻⁸, but also geological and structural complexities of the lower and upper plates⁹. However, the persistence of this short-term frictional pattern is difficult to confirm despite crucial implications at various timescales, such as predicting potential ruptures on the interface¹⁰ and assessing margin dynamics over thousands to millions of years¹¹,¹².

The accretion of downgoing material at the base of the fore-arc region¹³ (i.e., tectonic underplating) has been proposed to be triggered by local variations of plate coupling¹⁴ and may, to this end, shed light on hydro-mechanical processes and stress distribution along the plate interface. However, the exact connection between interface properties, deformation mechanisms and detachments of tectonic slices remains elusive and requires to deal with pore fluid pressure and effective stress interrelations. In this context, we assess for the first time the relations between long-term fluid flow, stress regime and tectonic underplating events in deep fore-arc regions, using high-resolution, two-dimensional, thermo-mechanical experiments.

## Modelling Tectonic Underplating and Fluid Transport

In this study, we reproduce a self-consistent ocean-continent subduction zone governed by conservation laws and visco-elasto-plastic rheologies¹⁵, solved in a high-resolution domain (i.e., 0.5 and 0.4 km between each node in x and y directions, respectively). The model takes into account hydration and dehydration processes by implementing compaction and thermodynamically constrained metamorphic reactions¹⁶, as well as fluid transport driven by fluid buoyancy and dynamic pressure gradients¹¹,¹⁷. As a result, fluids defined as Lagrangian markers move from high- to low-stress regions, depending on calculated deviatoric stresses (see Methods and Supplementary Fig. S1 for details on the modelling procedure and the initial setup).

**Modelling strategy.** To assess the variability of subduction zones worldwide, we performed a set of numerical experiments by varying the cooling age of the subducting oceanic lithosphere (i.e., 53 Myr and 20 Myr), the

<footnote n="1">¹Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005, Paris, France. ²Present address: GFZ Helmholtz Centre Potsdam, German Research Centre for Geosciences, Telegrafenberg, 14473, Potsdam, Germany.</footnote>

<footnote n="3">³Institute of Geophysics, Swiss Federal Institute of Technology (ETH), Zürich, Switzerland. Correspondence and requests for materials should be addressed to A.M. (email: armel.menant@gmail.com)</footnote>
plate convergence rate (i.e., from 4 to 8 cm yr⁻¹) and the overriding continental crust thickness (i.e., 30 and 40 km thick; Supplementary Table S1). Aside from these parameters, the main uncertainty on fluid transport calculation is the reference percolation velocity (v<sub>perc</sub>; see Methods). In our reference model, v<sub>perc</sub> = 1 cm yr⁻¹, which is of same order of magnitude than estimations for long-term fluid flow in accretionary wedges<sup>18,19</sup>. However, these estimations have been obtained from *in-situ* measurements and numerical modelling on the most frontal part of the accretionary wedge and the uncertainty on the fluid velocity remains high in the deep fore-arc region. Accordingly, additional experiments have been carried out by varying v<sub>perc</sub> to explore the effect of fluid drainage on plate interface strength and accretion processes (Supplementary Table S1).

**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⁻¹. 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*<sub>Z1</sub> and *U*<sub>Z2</sub>, 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*<sub>Z1</sub> and *U*<sub>Z2</sub> 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¹⁷. 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¹,¹⁶, 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*<sub>perc</sub>) result in significantly different fluid distribution and accretion dynamics (Fig. 3). Thus, a low *v*<sub>perc</sub> (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*<sub>perc</sub> 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*<sub>Z1</sub>; 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).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Cold subduction zone and fluid flow pattern (reference model *cold30-5-1*). (a) Compositional map. Overview of the subduction zone. Rock types with initially-prescribed pore water content are italicized. (b) Compositional map. Zoom on the deep fore-arc region. (c) Strain-rate map. (d) Differential-stress map. The black streamlines indicate fluid flow pattern calculated according to dynamic pressure gradients. (e) Fluid-distribution map. Fluid markers (in blue) indicate local fluid oversaturation. See Fig. 2a for details on the identification of the preferential loci for underplating *U*<sub>Z1</sub> and *U*<sub>Z2</sub>. Thick dashed lines (in white and black on panels *d* and *e*, respectively) depict the subduction channel.</figcaption></figure>

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
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Relations between shear-stress variations along the subduction channel and discrete underplating events in (a) cold and (b) warm subduction models. Line charts show the shear-stress evolution within the subduction channel integrated over a 1-Myr-long period (i.e., 55-56 Myr and 43-44 Myr for the cold and warm subduction models, respectively). Bar charts show the horizontal distribution of the maximum upward vertical velocity component at the top of the subduction channel for each time step, which is thought to reflect the distribution of underplating events along the plate interface. The main underplating loci $U_{Z1}$ and $U_{Z2}$ are highlighted by grey bands. In the warm subduction model, the underplating locus $U_{Z2}$ is only active until ~35 Myr (see Supplementary Fig. S8 for details).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Correlation between fluid distribution, frictional properties of the plate interface and fore-arc accretion dynamics. (a) Model-based sketches show the first-order geometry of the fore-arc margin for three simulations displaying different reference percolation velocity ($v_{perc}$). Accreted material is in yellow. Fluid markers (in blue) indicate local fluid oversaturation. See Fig. 1 and Supplementary Figs S6 and S7 for details on models *cold30-5-1*, *cold30-5-0.1* and *cold30-5-10*, respectively. (b) Line chart shows shear-stress variations within the subduction channel integrated over 1-Myr-long. Green, orange and red lines correspond to models with $v_{perc}$ equal to 0.1, 1 and 10 cm yr$^{-1}$, respectively.</figcaption></figure>

from the age of the oceanic lithosphere, the plate convergence rate and the thickness of the overriding crust (Figs 1, 4 and Supplementary Figs S3, S4 and S5). Similar nappe stacking has been largely recognized from paleo-accretionary complexes, including the Franciscan complex (Western US)²⁰ and the Western Series
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Warm subduction zone and fluid flow pattern (model *warm30-5-1*). (a) Compositional map. Overview of the subduction zone. Rock types with initially-prescribed pore water content are italicized. (b) Compositional map. Zoom on the deep fore-arc region. (c) Strain-rate map. (d) Differential-stress map. The black streamlines indicate fluid flow pattern calculated according to dynamic pressure gradients. (e) Fluid-distribution map. Fluid markers (in blue) indicate local fluid oversaturation. See Fig. 2b for details on the identification of the preferential locus for underplating *U*<sub>ZI</sub>. Thick dashed lines (in white and black on panels **d** and **e**, respectively) depict the subduction channel.</figcaption></figure>

(Central-Southern Chile)²¹,²². Geophysical imaging of duplex structures has also been reported from several actives subduction zones²³⁻²⁶ (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²⁷,²⁸. These include hydrofracture networks into fossilized subduction channels from 10 to >40 km depth²⁹,³⁰, as
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Fluid distribution, hydro-mechanical properties of the plate interface and preferential location of tectonic underplating in cold and warm subduction settings. Field evidences for fluid distributions at various depths from ancient subduction zones allow for supporting model results (see location of panels a,b,c on panels d and e). (a) Foliation-parallel and nearly vertical quartz veins crosscutting greenschist-facies metasediments (Western Series, Central Chile, high thermal gradient). (b) Interpreted panoramic view of the Dent Blanche basal tectonic contact (Western Alps, intermediate thermal gradient). Foliated orthogneiss located several hundreds of meters above the interface displays no evidence of fluid percolation (see inset). (c) Hydrofractured blueschist-facies metavolcanics (Zagros, Iran, low thermal gradient). (d,e) Model-based sketches linking hydro-mechanical properties of the plate interface with tectonic underplating events in cold and warm subduction zones, respectively. Line charts show shear-stress variations along the interface. Bar charts show the horizontal distribution of underplating events along the plate interface for the entire model duration (see Fig. 2 for details on these charts). Blue bands depict the relative proportion of fluid markers (i.e., fluid-oversaturated conditions) into the subduction channel. Blue arrows depict the dominant fluid flow pattern. LVZ: Low-velocity zone.</figcaption></figure>

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²⁷ (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³¹,³² (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³¹,³³⁻³⁵. Several hypotheses to explain this phenomenon have been proposed such as impermeable lithologies³⁶, deformation-controlled grain-size reduction processes³⁷ or mineral precipitation³⁸. 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³⁹⁻⁴¹ (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⁴²) 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⁴³ (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⁴⁴, 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²⁷ and by field evidences of intense veining along the plate interface at similar depths³⁰,⁴⁵ (Fig. 5c).
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<sup>13,46,47</sup>. 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<sup>46,48</sup>. 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<sup>28,49</sup>, as well as on the mantle-derived isotopic signature of mineral springs located at the surface<sup>50,51</sup>. 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<sup>27,28</sup>, indicating important fluid accumulation over long timescales, responsible for partial serpentinization of the mantle wedge<sup>52</sup>. 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<sup>53,54</sup> (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² 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<sup>21,55</sup>, previous numerical studies<sup>56</sup> and field-based conceptual models<sup>13,22,57</sup>. 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<sup>11,53</sup>.

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<sup>24,26,58,59</sup>. In addition, GPS-based locking maps show a heterogeneous distribution of coupled subduction segments with a common decrease approaching the subcontinental mantle<sup>8,10</sup> 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<sup>5,60</sup>. 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.
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¹⁵. 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:

$$ \frac{D \ln \rho_{eff}}{Dt} + \frac{\partial v_i}{\partial x_i} = 0 \qquad (1) $$

where $\rho_{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:

$$ -\frac{\partial P}{\partial x_i} + \frac{\partial \sigma_{ij}}{\partial x_i} = - \rho_{eff} g_i, \qquad (2) $$

where P is the pressure, $\sigma_{ij}$ the components of the deviatoric stress tensor and $g_i$ the gravitational acceleration ($g_x=0$ and $g_y=9.81 \text{ m s}^{-2}$). The heat conservation equation is formulated in a Lagrangian form to avoid numerical diffusion of temperature:

$$ \rho_{eff} C_p \frac{DT}{Dt} = - \frac{\partial q_i}{\partial x_i} + H_r + H_a + H_s, \qquad (3) $$

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:

$$ q_i = -k \frac{\partial T}{\partial x_i}, \qquad (4) $$

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_{w_{pore}} = 1$ 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²⁹ (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²,⁶¹. Bound water release is calculated by free-energy minimization¹⁶ as a function of pressure, temperature and rock type¹¹. 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⁶². 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¹⁷. Accordingly, fluid marker transport is controlled by the Stokes velocity field $v_p$, the fluid buoyancy and the dynamic pressure gradients, such that:

$$ v_{i,water} = v_i + v_{perc} k_i, \qquad (5) $$

where $v_{i,water}$ is the velocity of fluid markers, $v_{perc}$ the reference percolation velocity (varying from 0.1 to 10 cm yr⁻¹ 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:

$$ k_i = \frac{\rho_{crust} g_i - \frac{\partial P}{\partial x_i}}{(\rho_{crust} - \rho_{fluid}) g_y}, \qquad (6) $$

where $\rho_{crust}$ and $\rho_{fluid}$ are respectively the reference crustal and fluid densities ($\rho_{crust} = 2300 \text{ kg m}^{-3}$ and $\rho_{fluid} = 1000 \text{ kg m}^{-3}$) and $|k_i| = 2$ for $|k_i| > 2$. Upward transport is defined by positive values. Once moving, fluid markers may be then consumed by rock markers (either as pore or mineral bound water), depending on their stable water content. Further details on the fluid implementation are provided by refs¹¹,¹⁷.

**Fluid effects on rock density and rheology.** In our experiments, bound water content in rock markers affects rock density as follows:

$$ \rho_{eff} = \rho_{rock}(1 - X_{fluid}) + \rho_{fluid} X_{fluid}, \qquad (7) $$

with

$$ \rho_{rock} = \rho_0 (1 - \alpha(T - 298)) (1 + \beta (P - 0.1)), \qquad (8) $$
where $\rho_{0solid}$ is the standard density of rocks, $X_{fluid}$ the mass fraction of fluid, $\alpha$ the thermal expansion and $\beta$ 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 $\varepsilon_{ij}$ includes the three respective components:

$$ \dot{\varepsilon}_{ij} = \dot{\varepsilon}_{ij_{viscous}} + \dot{\varepsilon}_{ij_{elastic}} + \dot{\varepsilon}_{ij_{plastic}} \quad (9) $$

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 $\sigma_{yield}$, which limits the creeping (i.e., viscous) behaviour such that:

$$ \eta \leq \frac{\sigma_{yield}}{2 \dot{\varepsilon}_{II}} \quad (10) $$

with

$$ \eta = \dot{\varepsilon}_{II} \frac{1-n}{n} A_D \frac{1}{n} \exp \left( \frac{E + P V}{n R T} \right) \quad (11) $$

and

$$ \sigma_{yield} = C + P \sin(\varphi_{dry}) (1 - \lambda_{fluid}) \quad (12) $$

where $\eta$ is the effective creep viscosity, $\dot{\varepsilon}_{II}$ the second invariant of the strain rate tensor, $n$ the creep exponent, $A_D$ the pre-exponential factor, $E$ the activation energy, $V$ the activation volume, $R$ the gas constant, $C$ the cohesion, $\varphi_{dry}$ the internal friction angle for dry rocks and $\lambda_{fluid}$ the pore fluid pressure factor. The latter is defined as $\lambda_{fluid} = 0$ for dry rocks and $\lambda_{fluid} = 0.99$ for fluid-oversaturated rocks (indicated locally by the presence of fluid markers) except at the surface where hydrostatic conditions are assumed (i.e., $\lambda_{fluid} = 0.4$). Because no fluid infiltration or alteration is evidenced a few hundreds of meters above fossilized subduction channels$^{43,63}$, fluid weakening effects have been deactivated for the upper and lower continental crusts.

**Numerical setup.** In all numerical experiments, the computational domain measures 1500 × 200 km in the $x$ and $y$ direction, respectively (Fig. S1a). Eulerian grid resolution is 2.0 × 1.5 km except at the centre of the domain where it is 0.5 × 0.4 km. Additionally, ~8 millions of randomly distributed markers are initially prescribed for advecting material properties and computing water release, transport and consumption. Convergence rate (defined within the convergence condition region) is 5 cm yr$^{-1}$. The initial setup is designed with a 30-km-thick overriding continental crust composed of 15 km of felsic upper crust and 15 km of mafic lower crust. The 7.5-km-thick subducting oceanic crust is made up of 0.5 km of pelagic sediments, 2 km of hydrated basaltic crust and 5 km of gabbroic crust (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 \le x \le 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 \le x \le 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}$:

$$ \frac{\partial y_{surf}}{\partial t} = v_y - v_x \frac{\partial y_{surf}}{\partial x} - v_{sedim} + v_{erosion} \quad (13) $$

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).

## References

1. Peacock, S. A. Fluid Processes in Subduction Zones. *Science* **248**, 329–337 (1990).

2. Saffer, D. M. & Tobin, H. J. Hydrogeology and Mechanics of Subduction Zone Forearcs: Fluid Flow and Pore Pressure. *Annu Rev Earth Planet Sci* **39**, 157–186 (2011).

3. Scholz, C. H. Earthquakes and friction laws. *Nature* **391**, 37–42 (1998).
4. van Dinther, Y. et al. The seismic cycle at subduction thrusts: Insights from seismo-thermo-mechanical models. *J Geophys Res Solid Earth* **118**, 6183–6202 (2013).

5. Gao, X. & Wang, K. Rheological separation of the megathrust seismogenic zone and episodic tremor and slip. *Nature* **543**, 416–419 (2017).

6. Audet, P. & Schwartz, S. Y. Hydrologic control of forearc strength and seismicity in the Costa Rican subduction zone. *Nat Geosci* **6**, 852–855 (2013).

7. Saffer, D. M. Mapping fluids to subduction megathrust locking and slip behavior: Fluids and Subduction Megathrust Locking. *Geophys Res Lett* **44**, 9337–9340 (2017).

8. Moreno, M. et al. Chilean megathrust earthquake recurrence linked to frictional contrast at depth. *Nat. Geosci* **11**, 285–290 (2018).

9. Kopp, H. Invited review paper: The control of subduction zone structural complexity and geometry on margin segmentation and seismicity. *Tectonophysics* **589**, 1–16 (2013).

10. Chlieh, M. et al. Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust. *Earth Planet Sci Lett* **400**, 292–301 (2014).

11. Gerya, T. V. & Meilick, F. I. Geodynamic regimes of subduction under an active margin: effects of rheological weakening by fluids and melts. *J Metamorph Geol* **29**, 7–31 (2011).

12. Tan, E., Lavier, L. L., Van Avendonk, H. J. A. & Heuret, A. The role of frictional strength on plate coupling at the subduction interface. *Geochem Geophys Geosyst* **13** (2012).

13. Platt, J. P. Dynamics of orogenic wedges and the uplift of high-pressure metamorphic rocks. *Geol Soc Am Bull* **97**, 1037 (1986).

14. Agard, P., Plunder, A., Angiboust, S., Bonnet, G. & Ruh, J. The subduction plate interface: Rock record and mechanical coupling (from long to short time scales). *Lithos* (2018).

15. Gerya, T. V. & Yuen, D. A. Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems. *Phys Earth Planet Inter* **163**, 83–105 (2007).

16. Connolly, J. A. D. Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation. *Earth Planet Sci Lett* **236**, 524–541 (2005).

17. Faccenda, M., Gerya, T. V. & Burlini, L. Deep slab hydration induced by bending-related variations in tectonic pressure. *Nat Geosci* **2**, 790–793 (2009).

18. Solomon, E. A. et al. Long-term hydrogeochemical records in the oceanic basement and forearc prism at the Costa Rica subduction zone. *Earth Planet Sci Lett* **282**, 240–251 (2009).

19. Screaton, E. J., Wuthrich, D. R. & Dreiss, S. J. Permeabilities, fluid pressures, and flow rates in the Barbados Ridge Complex. *J Geophys Res Solid Earth* **95**, 8997 (1990).

20. Kimura, G., Maruyama, S., Isozaki, Y. & Terabayashi, M. Well-preserved underplating structure of the jadeitized Franciscan complex, Pacheco Pass, California. *Geology* **24**, 75 (1996).

21. Glodny, J. et al. Internal dynamics of a paleoaccretionary wedge: insights from combined isotope tectonochronology and sandbox modelling of the South-Central Chilean forearc. *Earth Planet Sci Lett* **231**, 23–39 (2005).

22. Willner, A. P. Pressure–Temperature Evolution of a Late Palaeozoic Paired Metamorphic Belt in North–Central Chile (34°–35°30′S). *J Pet.* **46**, 1805–1833 (2005).

23. Moore, J. C. et al. EDGE deep seismic reflection transect of the eastern Aleutian arc-trench layered lower crust reveals underplating and continental growth. *Geology* **19**, 420 (1991).

24. Kimura, H., Takeda, T., Obara, K. & Kasahara, K. Seismic Evidence for Active Underplating Below the Megathrust Earthquake Zone in Japan. *Science* **329**, 210–212 (2010).

25. Calvert, A. J., Preston, L. A. & Farahbod, A. M. Sedimentary underplating at the Cascadia mantle-wedge corner revealed by seismic imaging. *Nat Geosci* **4**, 545–548 (2011).

26. Henrys, S. et al. SAHKE geophysical transect reveals crustal and subduction zone structure at the southern Hikurangi margin, New Zealand. *Geochem Geophys Geosyst* **14**, 2063–2083 (2013).

27. Bostock, M. G. The Moho in subduction zones. *Tectonophysics* **609**, 547–557 (2013).

28. Bürgmann, R. The geophysics, geology and mechanics of slow fault slip. *Earth Planet Sci Lett* **495**, 112–134 (2018).

29. Moore, J. C. & Vrolijk, P. Fluids in accretionary prisms. *Rev Geophys* **30**, 113 (1992).

30. Bebout, G. E. & Penniston-Dorland, S. C. Fluid and mass transfer at subduction interfaces—The field metamorphic record. *Lithos* **240–243**, 228–258 (2016).

31. Audet, P., Bostock, M. G., Christensen, N. I. & Peacock, S. M. Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing. *Nature* **457**, 76–78 (2009).

32. Shelly, D. R., Beroza, G. C., Ide, S. & Nakamula, S. Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip. *Nature* **442**, 188–191 (2006).

33. Hyndman, R. D., McCrory, P. A., Wech, A., Kao, H. & Ague, J. Cascadia subducting plate fluids channelled to fore-arc mantle corner: ETS and silica deposition. *J Geophys Res Solid Earth* **120**, 4344–4358 (2015).

34. Nakajima, J. & Hasegawa, A. Tremor activity inhibited by well-drained conditions above a megathrust. *Nat Commun* **7**, 13863 (2016).

35. Delph, J. R., Levander, A. & Niu, F. Fluid Controls on the Heterogeneous Seismic Characteristics of the Cascadia Margin. *Geophys Res Lett* (2018).

36. Kato, A. et al. Variations of fluid pressure within the subducting oceanic crust and slow earthquakes. *Geophys Res Lett* **37** (2010).

37. Menegon, L., Fusseis, F., Stütinitz, H. & Xiao, X. Creep cavitation bands control porosity and fluid flow in lower crustal shear zones. *Geology* **43**, 227–230 (2015).

38. Audet, P. & Bürgmann, R. Possible control of subduction zone slow-earthquake periodicity by silica enrichment. *Nature* **510**, 389–392 (2014).

39. Breeding, C. M. & Ague, J. J. Slab-derived fluids and quartz-vein formation in an accretionary prism, Otago Schist, New Zealand. *Geology* **30**, 499–502 (2002).

40. Okamoto, A., Kikuchi, T. & Tsuchiya, N. Mineral distribution within polymineralic veins in the Sanbagawa belt, Japan: implications for mass transfer during vein formation. *Contrib Miner. Pet.* **156**, 323–336 (2008).

41. Rowe, C. D., Meneghini, F. & Moore, J. C. Fluid-rich damage zone of an ancient out-of-sequence thrust, Kodiak Islands, Alaska. *Tectonics* **28** (2009).

42. Cloos, M. & Shreve, R. L. Subduction-Channel Model of Prism Accretion, Melange Formation, Sediment Subduction, and Subduction Erosion at Convergent Plate Margins: 2. Implications and Discussion. *Pure Appl Geophys* **128**, 501–545 (1988).

43. Angiboust, S. et al. Probing the transition between seismically coupled and decoupled segments along an ancient subduction interface. *Geochem Geophys Geosyst* **16**, 1905–1922 (2015).

44. Menant, A., Angiboust, S., Monié, P., Oncken, O. & Guigner, J.-M. Brittle deformation during Alpine basal accretion and the origin of seismicity nests above the subduction interface. *Earth Planet Sci Lett* **487**, 84–93 (2018).

45. Taetz, S., John, T., Bröcker, M., Spandler, C. & Stracke, A. Fast intraslab fluid-flow events linked to pulses of high pore fluid pressure at the subducted plate interface. *Earth Planet Sci Lett* **482**, 33–43 (2018).

46. Kawano, S., Katayama, I. & Okazaki, K. Permeability anisotropy of serpentinite and fluid pathways in a subduction zone. *Geology* **39**, 939–942 (2011).
47. Morishige, M. & van Keken, P. E. Along-arc variation in short-term slow slip events caused by 3-D fluid migration in subduction zones. *J Geophys Res Solid Earth* **122**, 1434–1448 (2017).

48. Ingebritsen, S. E. & Manning, C. E. Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism. *Geofluids* (2010).

49. Savard, G., Bostock, M. G. & Christensen, N. I. Seismicity, Metamorphism, and Fluid Evolution Across the Northern Cascadia Fore Arc. *Geochem Geophys Geosyst* **19**, 1881–1897 (2018).

50. Doğan, T., Sumino, H., Nagao, K. & Notsu, K. Release of mantle helium from forearc region of the Southwest Japan arc. *Chem Geol* **233**, 235–248 (2006).

51. McCrory, P. A., Constantz, J. E., Hunt, A. G. & Blair, J. L. Helium as a tracer for fluids released from Juan de Fuca lithosphere beneath the Cascadia forearc. *Geochem Geophys Geosyst* **17**, 2434–2449 (2016).

52. Hyndman, R. D. & Peacock, S. M. Serpentinization of the forearc mantle. *Earth Planet Sci Lett* **212**, 417–432 (2003).

53. Sage, F., Collot, J.-Y. & Ranero, C. R. Interplate patchiness and subduction-erosion mechanisms: Evidence from depth-migrated seismic images at the central Ecuador convergent margin. *Geology* **34**, 997 (2006).

54. Ranero, C. R. et al. Hydrogeological system of erosional convergent margins and its influence on tectonics and interplate seismogenesis. *Geochem Geophys Geosyst* **9** (2008).

55. Gutscher, M.-A., Kukowski, N., Malavieille, J. & Lallemand, S. Cyclical behavior of thrust wedges: Insights from high basal friction sandbox experiments. *Geology* **24**, 135–138 (1996).

56. Ellis, S., Beaumont, C. & Pfiffner, O. A. Geodynamic models of crustal-scale episodic tectonic accretion and underplating in subduction zones. *J Geophys Res* **104**, 15169–15190 (1999).

57. Angiboust, S. et al. A 100-m.y.-long window onto mass-flow processes in the Patagonian Mesozoic subduction zone (Diego de Almagro Island, Chile). *Geol Soc Am Bull* **130**, 1439–1456 (2018).

58. Bassett, D. & Watts, A. B. Gravity anomalies, crustal structure, and seismicity at subduction zones: 2. Interrelationships between fore-arc structure and seismogenic behavior. *Geochem Geophys Geosyst* **16**, 1541–1576 (2015).

59. Tsuji, T., Minato, S., Kamei, R., Tsuru, T. & Kimura, G. 3D geometry of a plate boundary fault related to the 2016 Off-Mie earthquake in the Nankai subduction zone, Japan. *Earth Planet Sci Lett* **478**, 234–244 (2017).

60. Saffer, D. M. & Wallace, L. M. The frictional, hydrologic, metamorphic and thermal habitat of shallow slow earthquakes. *Nat. Geosci* **8**, 594–600 (2015).

61. Schmidt, M. W. & Poli, S. Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation. *Earth Planet. Sci. Lett.* **163**, 361–379 (1998).

62. Spiegelman, M. Flow in deformable porous media. Part 1 Simple analysis. *J Fluid Mech* **247**, 17 (1993).

63. Bachmann, R. et al. Exposed plate interface in the European Alps reveals fabric styles and gradients related to an ancient seismogenic coupling zone. *J Geophys Res Solid Earth* **114** (2009).

64. Schmeling, H. et al. A benchmark comparison of spontaneous subduction models—Towards a free surface. *Phys Earth Planet Inter* **171**, 198–223 (2008).

65. Gorczyk, W., Willner, A. P., Gerya, T. V., Connolly, J. A. D. & Burg, J.-P. Physical controls of magmatic productivity at Pacific-type convergent margins: Numerical modelling. *Phys Earth Planet Inter* **163**, 209–232 (2007).

## Acknowledgements

This study has been funded with an IDEX-USPC research chair grant 16C538 and contributes to the IdEx Université de Paris ANR-18-IDEX-0001. Robin Lacassin is acknowledged for insightful discussions and early revision of the manuscript. The manuscript benefited from suitable comments by four anonymous reviewers. This is IPGP contribution #4050.

## Author Contributions

S.A. conceived the original idea, T.G. designed the 2D-thermomechanical code, A.M. elaborated the numerical study together with S.A. and T.G., A.M. performed all the numerical experiments and interpreted the results. All authors discussed the results and interpretations and participated in writing the paper.

## Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-019-46191-y.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2019