Gaspard Farge
email: farge@ipgp.fr

Claude Jaupart

Nikolai M Shapiro

Episodicity and Migration of Low Frequency Earthquakes modeled with Fast Fluid Pressure Transients in the Permeable Subduction Interface

We model fluid pressure diffusion in a channel along the subduction interface in which the permeability varies through a valving mechanism • Tremor is generated as valves break open and interactions between valves through fluid pressure transients shape realistic event patterns • The input fluid flux in the fault controls seismic activity, opening possibilities to diagnose hydraulic conditions using LFE patterns

Introduction

The extensive records of weak low-frequency seismicity detected in many fault zones over the last two decades provide powerful information on their dynamics. At intermediate depths (30-50 km) in several subduction zones, seismic activity takes the form of tremor (e.g., [START_REF] Obara | Nonvolcanic deep tremor associated with subduction in southwest japan[END_REF][START_REF] Kao | A wide depth distribution of seismic tremors along the northern Cascadia margin[END_REF][START_REF] Payero | Nonvolcanic tremor observed in the Mexican subduction zone[END_REF][START_REF] Husker | Temporal variations of non-volcanic tremor (NVT) locations in the Mexican subduction zone: Finding the NVT sweet spot[END_REF]. It now seems clear that tremor is made of swarms of low-frequency earthquakes (LFEs) and many authors use the term "tremor" as a synonym for LFE activity (e.g., [START_REF] Shelly | Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip[END_REF]Song et al., 2009;[START_REF] Peacock | High pore pressures and porosity at 35 km depth in the Cascadia subduction zone[END_REF][START_REF] Audet | Teleseismic constraints on the geological environment of deep episodic slow earthquakes in subduction zone forearcs: A review[END_REF][START_REF] Frank | Low-frequency earthquakes in the Mexican Sweet Spot: LFES IN THE MEXICAN SWEET SPOT[END_REF][START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF][START_REF] Rogers | Episodic Tremor and Slip on the Cascadia Subduction Zone: The Chatter of Silent Slip[END_REF][START_REF] Bostock | Low frequency earthquakes below southern Vancouver Island[END_REF]. These LFEs appear to be due to repetitive, impulsive and localized sources [START_REF] Beroza | Slow Earthquakes and Nonvolcanic Tremor[END_REF][START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF]. High rates of tremor activity are associated with geodetically observed slow deformation episodes called slowslip events (SSEs) (e.g., [START_REF] Kostoglodov | The 2006 slow slip event and nonvolcanic tremor in the Mexican subduction zone[END_REF][START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes: UNCOVERING SILENT SLIP[END_REF][START_REF] Rousset | Slow slip events in the roots of the San Andreas fault[END_REF][START_REF] Delbridge | Geodetic Measurements of Slow-Slip Events Southeast of Parkfield, CA[END_REF]. The correlation between LFE and SSE activity gave rise to the term Episodic Tremor and Slip (ETS) [START_REF] Rogers | Episodic Tremor and Slip on the Cascadia Subduction Zone: The Chatter of Silent Slip[END_REF]. In a more general acceptation, transient seismic and geodetic signals that are both generated during inter-seismic phases are often lumped together as slow earthquake phenomena. This new class of seismic and geodetic observations holds great promise for improving our understanding of the seismo-tectonic cycle.

According to the framework that prevails today, slow earthquakes are the "slow" analogs of "regular" earthquakes. Geodetic and seismic events are manifestations of the same process of slow slip at transitional depths in fault zones, in between unstable stickslipping and stable sliding areas. The transitional part is made of slow sliding patches that are responsible for geodetically detected SSEs and that contain many small embedded frictionally unstable asperities generating LFEs. This framework explains why SSEs and tremor activity are correlated in time and accounts for the S-wave dominated LFE radiation pattern (e.g., [START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF]. Numerical models of these phenomena (e.g., [START_REF] Shibazaki | On the physical mechanism of silent slip events along the deeper part of the seismogenic zone[END_REF][START_REF] Liu | Spontaneous and triggered aseismic deformation transients in a subduction fault model[END_REF][START_REF] Ben-Zion | Episodic tremor and slip on a frictional interface with critical zero weakening in elastic solid[END_REF][START_REF] Ide | Modeling fast and slow earthquakes at various scales[END_REF][START_REF] Luo | Rate-and-State Model Casts New Insight into Episodic Tremor and Slow-slip Variability in Cascadia[END_REF][START_REF] Barbot | Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault[END_REF] are often based on the rate-and-state friction (RSF) formalism developed for "regular" earthquakes [START_REF] Dieterich | Earthquake nucleation on faults with rate-and statedependent strength[END_REF]. This efficient modeling approach has two significant shortcomings, however (van den [START_REF] Van Den Ende | A comparison between rate-and-state friction and microphysical models, based on numerical simulations of fault slip[END_REF]. First, all the dynamic vari-ability of the interface and its surroundings is collapsed onto a planar boundary. Second, it depends on empirical RSF parameters that are not inferred directly from processes occurring within the fault zone. Thus, it is difficult to specify RSF parameters and account for internal fault zone dynamics in a self-consistent manner. These shortcomings may be limiting in the brittle-to-ductile transition zone of subduction faults because the deformation proceeds over an increasingly wide off-fault region as depth increases [START_REF] Platt | Rheology and stress in subduction zones around the aseismic/seismic transition[END_REF]. Source dimensions for tremor and LFEs are typically inferred to be a few hundred meters (e.g., [START_REF] Chestler | Evidence for a scale-limited low-frequency earthquake source process[END_REF][START_REF] Supino | Self-similarity of low-frequency earthquakes[END_REF][START_REF] Farge | Moment-Duration Scaling of Low-Frequency Earthquakes in Guerrero, Mexico[END_REF][START_REF] Sammis | A Granular Jamming Model for Low-Frequency Earthquakes[END_REF], close to the thickness of the highly deformed permeable fault zone (e.g., [START_REF] Angiboust | JGR: Solid Earth June). Probing the transition between seismically coupled and decoupled segments along an ancient subduction interface[END_REF]. This calls for more complete models dealing with processes in the interior of these zones.

There can be little doubt that pore fluid pressures and their variations play a key role in fault zone dynamics [START_REF] Saffer | Hydrogeology and Mechanics of Subduction Zone Forearcs: Fluid Flow and Pore Pressure[END_REF]. From a frictional perspective, high pore pressures reduce the fault's strength, bringing it closer to unstable sliding conditions. Numerous geophysical and geological observations provide evidence for fluid saturation in the subduction interface [START_REF] Saffer | Hydrogeology and Mechanics of Subduction Zone Forearcs: Fluid Flow and Pore Pressure[END_REF]. It is widely believed that fluids are supplied by metamorphic dehydration reactions in the downgoing plate [START_REF] Manning | Coupled reaction and flow in subduction zones: Silica metasomatism in the mantle wedge[END_REF][START_REF] Van Keken | Subduction factory: 4. Depth-dependent flux of H 2 O from subducting slabs worldwide[END_REF][START_REF] Frezzotti | The chemical behavior of fluids released during deep subduction based on fluid inclusions[END_REF][START_REF] Hyndman | Cascadia subducting plate fluids channelled to fore-arc mantle corner: ETS and silica deposition[END_REF]. High values of the V p /V s ratio and electrical conductivity indicate that pore fluids are at nearlithostatic pressures [START_REF] Peacock | High pore pressures and porosity at 35 km depth in the Cascadia subduction zone[END_REF][START_REF] Wannamaker | Segmentation of plate coupling, fate of subduction fluids, and modes of arc magmatism in Cascadia, inferred from magnetotelluric resistivity[END_REF]Audet & Kim, 2016, among others), implying that the fault should be very weak. This explains why seismic activity is highly sensitive to very small dynamic stress changes generated by passing teleseismic waves [START_REF] Rubinstein | Seismic wave triggering of nonvolcanic tremor, episodic tremor and slip, and earthquakes on Vancouver Island[END_REF] or tidal deformation [START_REF] Rubinstein | Tidal Modulation of Nonvolcanic Tremor[END_REF][START_REF] Thomas | Tremor-tide correlations and near-lithostatic pore pressure on the deep San Andreas fault[END_REF][START_REF] Beeler | Inferring fault rheology from low-frequency earthquakes on the San Andreas[END_REF][START_REF] Royer | Tidal modulation and triggering of low-frequency earthquakes in northern Cascadia[END_REF].

There is mounting evidence for large variations of pore pressure on the time scale of the slow-earthquake cycle. During slow slip events, seismic images, gravimetric studies as well as earthquake source mechanisms indicate that pore fluid pressure around the subduction interface drops markedly [START_REF] Tanaka | Temporal gravity anomalies observed in the Tokai area and a possible relationship with slow slips[END_REF][START_REF] Warren-Smith | Episodic stress and fluid pressure cycling in subducting oceanic crust during slow slip[END_REF][START_REF] Gosselin | Seismic evidence for megathrust fault-valve behavior during episodic tremor and slip[END_REF]. The migration of LFE activity is highly suggestive of pore-pressure diffusion (Frank, Shapiro, et al., 2015;[START_REF] Cruz-Atienza | Rapid tremor migration and pore-pressure waves in subduction zones[END_REF]. In exhumed subduction interfaces, the brittle-to-ductile transition seems to be intrinsically coupled to transient pore-pressure cycling and fluid transport [START_REF] Angiboust | JGR: Solid Earth June). Probing the transition between seismically coupled and decoupled segments along an ancient subduction interface[END_REF][START_REF] Platt | Rheology and stress in subduction zones around the aseismic/seismic transition[END_REF][START_REF] Behr | What's down there? The structures, materials and environment of deep-seated slow slip and tremor[END_REF]. Episodes of fluid accumulation and drainage have been doc-manuscript submitted to JGR: Solid Earth umented in shallow thrust faults and are recorded by metasomatic reactions in host rocks [START_REF] Angiboust | Channelized Fluid Flow and Eclogite-facies Metasomatism along the Subduction Shear Zone[END_REF][START_REF] Saffer | The permeability of active subduction plate boundary faults[END_REF][START_REF] Taetz | Fast intraslab fluid-flow events linked to pulses of high pore fluid pressure at the subducted plate interface[END_REF]. They have been attributed to changes of permeable pathways due to tectonic deformation, precipitation/dissolution or deposition/erosion processes [START_REF] Giger | Permeability evolution in quartz fault gouges under hydrothermal conditions[END_REF]Jäger et al., 2017a;Williams et al., 2019). In recent efforts to model ETS phenomena, the emphasis has shifted from a static view of the impact of fluid pressure on RSF parameters to dynamic coupling between time-dependent fluid pressure and slip (e.g. [START_REF] Bernaudin | Episodic Tremor and Slip Explained by Fluid-Enhanced Microfracturing and Sealing[END_REF][START_REF] Cruz-Atienza | Rapid tremor migration and pore-pressure waves in subduction zones[END_REF][START_REF] Zhu | Fault valving and pore pressure evolution in simulations of earthquake sequences and aseismic slip[END_REF].

With the recognition that pore fluid pressures are high and variable in the subduction interface, it may not be appropriate to attribute seismic events to seismic slip only, without considering other processes that are active in the fault zone. Indeed, there are indications that the high fluid pressure transients that seem fundamental in shaping tremor patterns may well generate seismic events on their own. Observations do not demonstrate unambiguously that slow slip events and low-frequency seismic radiation are generated by the same process. It is true that these two types of activity are frequently correlated in time, but the sources of seismic and geodetic events are not systematically co-located (e.g. [START_REF] Hirose | Slow Earthquakes Linked Along Dip in the Nankai Subduction Zone: Fig. 1[END_REF][START_REF] Kostoglodov | The 2006 slow slip event and nonvolcanic tremor in the Mexican subduction zone[END_REF]. In addition, in some subduction zones, LFEs do not conform to the moment-duration scaling of classical shear rupture [START_REF] Bostock | Magnitudes and moment-duration scaling of low-frequency earthquakes beneath southern Vancouver Island[END_REF][START_REF] Farge | Moment-Duration Scaling of Low-Frequency Earthquakes in Guerrero, Mexico[END_REF] and do not follow the Gutenberg-Richter size distribution characteristic of regular earthquake activity [START_REF] Sammis | A Granular Jamming Model for Low-Frequency Earthquakes[END_REF]. Finally, it should be noted that the predominance of S-waves in the seismic signals does not necessarily establish a shear rupture mechanism. S-wave dominated signals can arise from volumetric sources with non-spherical geometries (e.g. [START_REF] Melnik | Deep long period volcanic earthquakes generated by degassing of volatile-rich basaltic magmas[END_REF] or from singleforce sources associated with fast local pressure variations [START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF].

Strong arguments in favor of a fluid-related origin for LFEs comes from laboratory studies and volcano seismology. In volcanic systems, shear deformation and slip do not play a predominant role, and yet tremors and LFEs (labelled as long-period earthquakes) are generated in abundance. These events have been attributed to fluid related processes, such as resonance in fluid-filled cracks [START_REF] Chouet | Long-period volcano seismicity: Its source and use in eruption forecasting[END_REF], strong flow in tortuous channels [START_REF] Julian | Volcanic tremor: Nonlinear excitation by fluid flow[END_REF][START_REF] Benson | Laboratory simulation of volcano seismicity[END_REF] and rapid pressurization of volatiles [START_REF] Melnik | Deep long period volcanic earthquakes generated by degassing of volatile-rich basaltic magmas[END_REF][START_REF] Wech | Deep long-period earthquakes generated by second boiling beneath Mauna Kea volcano[END_REF], among others. In some cases, their sources are distributed over large depth intervals of a few tens of kilometers and appear to migrate on time scales of a few months or days [START_REF] Shapiro | Deep and shallow long-period volcanic seismicity linked by fluidpressure transfer[END_REF]. These migrations have been attributed manuscript submitted to JGR: Solid Earth to pore-pressure transients associated with hydraulic connections between different parts of a volcanic plumbing system.

A key piece of information on fault zone dynamics is the distribution of LFEs in space and time. Here, we explore the potential of pressure transients in a permeable subduction interface to generate this type of activity. The structure and hydraulics of such a channel have rarely been studied in relation to seismicity [START_REF] Segall | Dilatancy, compaction, and slip instability of a fluid-infiltrated fault[END_REF][START_REF] Aochi | Self-induced seismicity due to fluid circulation along faults[END_REF][START_REF] Cruz-Atienza | Rapid tremor migration and pore-pressure waves in subduction zones[END_REF][START_REF] Zhu | Fault valving and pore pressure evolution in simulations of earthquake sequences and aseismic slip[END_REF]. We investigate the behavior of plugs of low permeability in a heterogeneous channel. Pore pressure variations can arise in two different ways. They may be generated by a time-dependent fluid supply, but there are no known mechanisms for strongly fluctuating rates of metamorphic dehydration reactions. Alternatively, one can invoke changes of permeability, as in the classical faultvalve model [START_REF] Sibson | Subducting Slab Ultra-Slow Velocity Layer Coincident with Silent Earthquakes in Southern Mexico[END_REF]. In this model, an increase of pore fluid pressure can generate an earthquake rupture, which opens up or reactivates fluid pathways leading to an increase of permeability and a decrease of pore pressure. Clogging of pores due to precipitation from circulating fluids or to the trapping of particles then acts to decrease permeability, implying an increase of pore pressure and a new cycle of pressure change. The permeability increase of the first phase of the cycle can be achieved in two different ways.

One is shear displacement in association with fault dilation [START_REF] Im | Cyclic permeability evolution during repose then reactivation of fractures and faults[END_REF] -although it appears that, in some cases, this may in fact act to decrease permeability in the slip direction [START_REF] Auradou | Permeability anisotropy induced by the shear displacement of rough fracture walls[END_REF]. In the other mechanism, high pore pressures can induce fracturing and/or unclogging of small pores and pore throats. In the laboratory, impulsive flushing events similar to hydrofracturation have been observed in densely packed permeable beds due to sudden particle detachment and called erosive bursts (Jäger et al., 2017b;Bianchi, Thielmann, et al., 2018). These bursts are associated with sharp pore pressure drops and recur as particles are redeposited. Coupled variations of permeability and pore fluid pressure due to particle clogging and un-clogging have been observed in many laboratory experiments [START_REF] Candela | Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing[END_REF]Bianchi, Thielmann, et al., 2018) and have been surmised for fault zones [START_REF] Manga | Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms[END_REF]. These phenomena have also been documented in large-scale systems. In some oil fields, for example, long-term records show that the bulk permeability of the reservoir oscillates in relation to cycles of deposition and removal of asphalt aggregates [START_REF] Sahimi | The Effect of Asphalt Precipitation on Flow Behavior and Production of a Fractured Carbonate Oil Reservoir During Gas Injection[END_REF]. We therefore argue that these mechanisms are bound to be active in faults, where fluid circulates under high pressure gradients, carrying precipitates and cataclastic fines that can be deposited and dis-lodged. In this framework, fast local pore pressure variations [START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF] and/or rupture of brittle heterogeneities due to increased fluid pressure [START_REF] Fagereng | Mélange rheology and seismic style[END_REF][START_REF] Kotowski | Length scales and types of heterogeneities along the deep subduction interface: Insights from exhumed rocks on Syros Island, Greece[END_REF]) may be responsible for LFEs. We expect that the unclogging or rupture of an individual plug induces pore pressure variations that propagate to neighboring plugs and generates swarms of similar events.

Through the rapid pressure drop that occurs, the unclogging of an individual plug is responsible for a local impulsive force on the solid matrix which generates elastic waves with both compressional and shear polarisation components. [START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF] have shown that the radiation pattern of such a single-force source is consistent with seismic observations in subduction zones due to the constrained source-receiver geometry, and that it can account for the weak dependency of LFEs source duration on their magnitude, as observed in Cascadia and Guerrero [START_REF] Bostock | Magnitudes and moment-duration scaling of low-frequency earthquakes beneath southern Vancouver Island[END_REF][START_REF] Farge | Moment-Duration Scaling of Low-Frequency Earthquakes in Guerrero, Mexico[END_REF].

Here, we build on this study to evaluate the seismic patterns that can be generated by a series of interacting plugs in a permeable channel. The clogging/unclogging mechanism does not involve shear rupture explicitly and we will discuss a posteriori how it should reflect on slip and its sensitivity at high fluid pressures.

The paper is organized as follows. We first describe the main characteristics of tremor/LFE activity in subduction zones that must be accounted for by any physical model. We then develop a quantitative model, starting with the basic processes that are involved in the clogging and unclogging of a plug and the corresponding equations. Results are described for a range of fluid input rates and a regime diagram is obtained and justified by simple theoretical arguments. The activity patterns that are generated are compared to observed ones, and conditions that must be met for a successful model are specified.

Tremor patterns observed in subduction zones

We summarize the most salient features of tremor activity in subduction zones. The Guerrero, Mexico, area is particularly appropriate to study LFE activity because it stretches over a distance of more than 100 km along-dip, owing to a low dipping angle. A very large catalog of events is available there [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF].

manuscript submitted to JGR: Solid Earth

Spatial Segmentation

At a local scale, the fact that LFEs occur in families of recurring events with similar waveforms indicates that they are due to the repeated activation of individual subkilometric scale sources [START_REF] Chestler | Evidence for a scale-limited low-frequency earthquake source process[END_REF]. In subduction zones, activity is often segmented in the along-strike direction and concentrated in patches [START_REF] Brudzinski | Segmentation in episodic tremor and slip all along Cascadia[END_REF][START_REF] Poiata | Segmentation and interaction of low-frequency earthquake activity along the strike of subducting slab in shikoku, japan[END_REF]. In Shikoku, Japan, and elsewhere, events delineate plurikilometric striations that are aligned with the subduction direction and that have been linked to the subduction of seamount chains [START_REF] Ide | Striations, duration, migration and tidal response in deep tremor[END_REF]. Shorter striations have also been observed in the Olympic Peninsula, Cascadia [START_REF] Ghosh | Rapid, continuous streaking of tremor in Cascadia[END_REF]. In Guerrero, Mexico, activity is split into two well-separated patches along the dip of the interface [START_REF] Payero | Nonvolcanic tremor observed in the Mexican subduction zone[END_REF][START_REF] Kostoglodov | The 2006 slow slip event and nonvolcanic tremor in the Mexican subduction zone[END_REF][START_REF] Husker | Temporal variations of non-volcanic tremor (NVT) locations in the Mexican subduction zone: Finding the NVT sweet spot[END_REF][START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF] (Figure 1).

Tremor Migration

Bursts of tremor activity migrate along the subduction interface on time scales of seconds to weeks, and spatial scales from a few to a few hundreds of kilometers. Shortlived, short-distance migrations are often restricted to an individual tremor patch and proceed at rates in the 10-100 km/hr range (on time scales T of minutes to hours over length scales L of tens of kilometers) [START_REF] Ghosh | Rapid, continuous streaking of tremor in Cascadia[END_REF][START_REF] Cruz-Atienza | Rapid tremor migration and pore-pressure waves in subduction zones[END_REF]. In marked contrast, long-lived, long-distance migrations that connect spatially separated patches are slower, with velocities in the ∼ 1-10 km/day range (T ∼ weeks to months, L ∼ 100 km) (Frank, Shapiro, et al., 2015;[START_REF] Kao | Northern Cascadia episodic tremor and slip: A decade of tremor observations from 1997 to 2007: NORTHERN CASCADIA ETS[END_REF]. Such systematics between time and space scales are reminiscent of a diffusion process, such that L 2 ∝ T and therefore V ∼ L/T ∝ 1/L.

In Japan and Cascadia, the largest transients of activity (T ∼weeks, L ∼100 km) propagate along-strike in a belt lying in a 30 to 60 km depth range. It has been suggested that this depth interval delineates a zone rheologically conducive to slow-slip [START_REF] Ide | Variety and spatial heterogeneity of tectonic tremor worldwide[END_REF][START_REF] Gomberg | Reconsidering earthquake scaling[END_REF]. In Guerrero, Mexico, the active segment of the subducted slab is dipping at a low angle and stretches over as much as 200 km in a 40-50 km depth range [START_REF] Pardo | Shape of the subducted Rivera and Cocos plates in southern Mexico: Seismic and tectonic implications[END_REF][START_REF] Kim | Geometry and seismic properties of the subducting Cocos plate in central Mexico[END_REF]. Due to this particularly favorable configuration, it is clear that the largest migrations systematically originate from the deep end of the fault and slowly propagate updip (Frank, Shapiro, et al., 2015). Short-lived migrations (T ∼ minutes to hours) often proceed along striations [START_REF] Ide | Variety and spatial heterogeneity of tectonic tremor worldwide[END_REF], sometimes in rapid sweeps [START_REF] Ghosh | Rapid, continuous streaking of tremor in Cascadia[END_REF]. Remarkably, these rapid migrations sometimes occur within larger scale migrations, propagating in the reverse direction with respect to the large scale propagation. This has been observed in the alongdip direction in Guerrero, Mexico [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF][START_REF] Cruz-Atienza | Rapid tremor migration and pore-pressure waves in subduction zones[END_REF] and in the along-strike direction in Cascadia [START_REF] Houston | Rapid tremor reversals in Cascadia generated by a weakened plate interface[END_REF].

Temporal Patterns

Tremor activity occurs during episodic bursts of events, which are composed of shorter bursts. Temporal clustering of activity covers time scales of years for the largest bursts -often correlated with slow-slip events and long-distance migrations -to days or even minutes for smaller bursts (Wech & Creager, 2011;[START_REF] Idehara | Regional and global variations in the temporal clustering of tectonic tremor activity[END_REF] -called cascades. The imbrication of small bursts within larger ones makes it difficult to determine a characteristic duration. The shortest bursts proceed over a few minutes, whereas the largest ones can persist for days or even months (e.g. [START_REF] Ghosh | Rapid, continuous streaking of tremor in Cascadia[END_REF][START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF][START_REF] Shelly | Complexity of the deep San Andreas Fault zone defined by cascading tremor[END_REF] As a temporal point process, this activity is clustered and does not conform to simple memoryless realizations. It has been analyzed using non-Poissonian statistics [START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF][START_REF] Beaucé | Systematic Detection of Clustered Seismicity Beneath the Southwestern Alps[END_REF].

The recurrence interval between tremor bursts decreases with increasing depth along the subduction interface (Wech & Creager, 2011). In addition, the time-clustering of activity can change from one tremor patch to its neighbors [START_REF] Idehara | Regional and global variations in the temporal clustering of tectonic tremor activity[END_REF][START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF], and can evolve as a function of time [START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF]. This has been interpreted as evidence for variations of fault zone conditions in both space and time [START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF].

Tremor activity Patterns in Guerrero, Mexico

Figure 1 shows the spatio-temporal patterns of LFE activity in the Guerrero tremor subduction segment, which is concentrated in two separate patches. One, called the sweet spot, is located at the downdip end of the segment and has a high rate of activity with a high level of time clustering [START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF]. The other patch, called transient zone [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF], is further updip and on the whole more quiescent. It hosts episodic bursts of activity with a recurrence interval of about 3 months [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF], which coincide with short-term SSEs recurrence [START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes: UNCOVERING SILENT SLIP[END_REF]. One of these episodes is displayed in Figure 1c. Activity migrates away from the sweet spot over a hundred kilometers updip. These episodes are built from many shorter bursts that may drift in both updip and downdip directions (Figure 1d).

The Guerrero subduction zone is also known for its very large SSEs that recur at intervals of about 4 years (e.g., [START_REF] Radiguet | Slow slip events and strain accumulation in the Guerrero gap, Mexico[END_REF] and are accompanied by intense tremor activity in both the transient zone and sweet spot. These events persist for several months with slow slip occurring mostly updip of the tremor-generating segment. The unstable slipping segment of the fault encroaches upon the transient zone over a small distance, where short-term SSEs are concentrated. During both short and long-term SSEs, tremor activity in the transient zone switches from episodic to highly clustered [START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF].

The most striking feature of the Guerrero subduction zone is perhaps that the patch with maximum tremor activity (the sweet spot) is separated from the patch with largest amount of slow slip (the transient zone) by more than 20 km.

Model design

In this section, we develop a model for time changes of permeability in a heterogeneous porous channel that runs along the dip of the subduction interface. Based on physical arguments, we associate sudden changes of fluid pressure with the generation of LFE-like, elementary tremor events.

Diffusion equation for pore pressure variations

Pore connectivity is modeled at a macroscopic level with a permeability that varies along the channel k(x). The fluid mass flux per unit cross-sectional area, q(x, t), depends on pore pressure P (x, t) according to Darcy's law:

q(x, t) = - k(x)ρ η ∂P ∂x (x, t) + ρg sin α (1)
where ρ is the fluid's density, η is the fluid viscosity, α is the dip angle of the subduction interface and g the acceleration of gravity. Mass balance dictates that :

manuscript submitted to JGR: Solid Earth - ∂q ∂x (x, t) = ∂(ρφ) ∂t . ( 2
)
where φ is porosity. We consider that pores are saturated with fluid and neglect temperature changes, which is justified by the small thickness of the channel and the small flow rates that are expected. Changes of the local mass of fluid can be due to the fluid compressibility as well as expanding/contracting pores. This is accounted for by an effective rock-fluid compressibility β, such that:

β = 1 ρφ ∂(ρφ) ∂P =⇒ - ∂q ∂x (x, t) = ρφβ ∂P ∂t (x, t). (3)
For simplicity, we neglect density variations in the body force term in equation ( 1), because they are due to changes of flow-induced dynamic pressure, which account for a fraction of the total pressure, and because they are not expected to play a significant role.

We split pressure into a hydrostatic component and a time-dependent dynamic component due to flow noted p. Neglecting changes of fluid viscosity, we obtain the following equation for the dynamic pore pressure component:

∂p ∂t = 1 φβη ∂ ∂x k(x) ∂p ∂x . (4)
We assume that permeability is piecewise constant. Thus, in each segment, the dynamic pressure obeys a diffusion equation:

∂p ∂t = k φβη ∂ 2 p ∂x 2 , (5)
with diffusivity D = k/(φβη). Along the modeled channel, D is therefore piecewise constant.

Figure 2 shows how we simplify the subduction fluid transport system conceptually and how it is represented as a permeable channel with plugs of low permeability.

A mechanism for large changes of permeability

Our model has two main ingredients, which deal with the nature of the fluid phase and its interaction with the network of pores in the permeable medium. Subduction zone metamorphic fluids are complex silica-saturated mixtures [START_REF] Manning | Coupled reaction and flow in subduction zones: Silica metasomatism in the mantle wedge[END_REF][START_REF] Frezzotti | The chemical behavior of fluids released during deep subduction based on fluid inclusions[END_REF]. As they rise, they undergo decompression and, due to the pressure dependence of solubility, are expected to precipitate either amorphous silica (silica gel) or quartz nanoparticles as well as other mineral phases. Depending on the availability of nucleation sites, precipitation may proceed along the pore walls or in the bulk fluid. High temperature decompression experiments on silica-rich fluids demonstrate that the silica nanoparticles that are generated are readily mobilized by fluid flow [START_REF] Amagai | Silica nanoparticles produced by explosive flash vaporization during earthquakes[END_REF][START_REF] Okamoto | Formation of silica particles from supercritical fluids and its impacts on the hydrological properties in the crust[END_REF]. Small clasts generated by grinding motion along the fault are also likely involved (e.g. [START_REF] Han | Clay-clast aggregates in fault gouge: An unequivocal indicator of seismic faulting at shallow depths[END_REF]. As in other deep permeable systems, the interstitial fluid phase is thus expected to carry fines and colloids [START_REF] Manga | Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms[END_REF].

Natural permeable media are made of pores and fractures connected by thinner throats which have a strong effect on permeability because they determine the effectiveness of fluid pathways through a volume. Capture and/or precipitation of even small amounts of solids within throats can decrease permeability by orders of magnitude [START_REF] Mcdowell-Boyer | Particle transport through porous media[END_REF]. Calculations by [START_REF] Beckingham | Evaluation of macroscopic porosity-permeability relationships in heterogeneous mineral dissolution and precipitation scenarios[END_REF] show that, depending on whether pore throats are open or clogged, permeability can change by more than a factor of ten.

These changes are associated with negligible porosity variations because pore throats occupy a very small fraction of the total interstitial volume. Thus, the transition between high and low permeability states is very abrupt, as plugging thin throats between larger pores involves very small amounts of sediment. Even if they are much smaller than a pore throat, fines and colloids can get stuck through a jamming process at rates that depend on their concentration [START_REF] Agbangla | Collective dynamics of flowing colloids during pore clogging[END_REF][START_REF] Civan | Modified formulations of particle deposition and removal kinetics in saturated porous media[END_REF][START_REF] Delouche | Dynamics of progressive pore clogging by colloidal aggregates[END_REF]. Clogging can be achieved in a matter of minutes [START_REF] Manga | Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms[END_REF][START_REF] Candela | Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing[END_REF][START_REF] Delouche | Dynamics of progressive pore clogging by colloidal aggregates[END_REF].

Material that gets stuck in pore throats may be remobilized by the flow [START_REF] Civan | Modified formulations of particle deposition and removal kinetics in saturated porous media[END_REF][START_REF] Kudrolli | Evolution of Porosity and Channelization of an Erosive Medium Driven by Fluid Flow[END_REF]Jäger et al., 2017b;[START_REF] Bianchi | Tomographic Study of Internal Erosion of Particle Flows in Porous Media[END_REF][START_REF] Amagai | Silica nanoparticles produced by explosive flash vaporization during earthquakes[END_REF]. Flow exerts a shear stress on the encasing solid, noted σ, that scales with the average pathway diameter d and the pressure gradent ∂p/∂x. If this stress exceeds a certain threshold value noted σ c , infillings are eroded and transported by the fluid [START_REF] Civan | Modified formulations of particle deposition and removal kinetics in saturated porous media[END_REF]. This may occur in catastrophic fashion in events called erosive bursts (Bianchi, Thielmann, et al., 2018). Such bursts may be composites of several smaller unclogging events (Bianchi, Thielmann, et al., 2018). The process has been studied in detail by [START_REF] Kudrolli | Evolution of Porosity and Channelization of an Erosive Medium Driven by Fluid Flow[END_REF]. In these experiments, a pure fluid is fed into a polydisperse particle bed. When the flow rate exceeds a threshold value, such that the flow-induced stress exceeds the threshold value σ c , the flow dislodges and mobilizes the smaller particles. Interestingly, the authors observe small precursor events when the flow rate is just below critical and a large avalanche when the flow rate reaches the critical value. This is due to the collective behavior of particles, as some of them are not able to move if they are in contact with stronger ones that resist detachment. Thus, these particles are only dislodged when their neighbors are, which triggers an erosion avalanche. As long as σ > σ c , particles cannot remain stuck and incipient coatings of pore walls are eroded, ensuring that fluid pathways stay fully open. Clogging thus requires that σ ≤ σ c . The alternation between erosion and deposition phases leaves little opportunity for particles and colloids to become permanently cemented to one another and to pore walls. Clogging and unclogging events therefore proceed rapidly, and can be considered as instantaneous at the time-scales of LFE activity.

As shown above, these processes have been observed in laboratory experiments on different types of porous media. Owing to their local character, they have not been documented in the field but, on a large scale, simultaneous variations of permeability and fluid pressure in faults [START_REF] Manga | Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms[END_REF] and oil reservoirs [START_REF] Sahimi | The Effect of Asphalt Precipitation on Flow Behavior and Production of a Fractured Carbonate Oil Reservoir During Gas Injection[END_REF] have been attributed to erosion/deposition processes.

A valve mechanism

We model clogging and unclogging events as follows. Within a long permeable channel, let us consider a clogged plug where fluid pathways are narrow and contorted because they wrap around captured particles (closed state). In this case, the average path-

σ -= λ 1 d lo δp break c w = σ c (6)
where λ 1 is a coefficient and w is the plug length along the channel. As plug particles are set in motion by the hydraulic stress and can no longer adhere to throat walls, the effective pore aperture d increases from d lo to d hi . In the permeable medium, variables are defined as spatial averages over a representative volume, and this is true for porosity and pressure in particular. As the solid particles get flushed out of the thin throats that connect larger pores, porosity is barely affected and, by the same token, so is the fluid pressure, because it is the average pressure for a group of pores and pore throats, with the larger pores acting to buffer local pressure variations at the pore throat scale.

Thus, the pressure differential across the pore is maintained at δp break c as the plug opens up [START_REF] Jäger | Clogging at pore scale and pressure-induced erosion[END_REF]. Therefore, σ increases to σ + :

σ + = λ 1 d hi δp break c w (7)
which is larger than σ c . σ then decreases gradually as the pressure distribution adapts to the new permeability structure. The plug remains unobstructed as long as σ is larger than σ c . Clogging occurs when σ = σ c such that δp reaches a second threshold value of pressure differential noted δp clog c :

λ 1 d hi δp clog c w = σ c . (8)
Clogging causes the permeability to drop back to k lo , driving an increase in δp until unclogging occurs again, thus starting a new cycle. This is a valve behavior akin to the toggleswitch of [START_REF] Miller | Permeability as a toggle switch in fluid-controlled crustal processes[END_REF].

This framework allows an internally consistent mechanism for a succession of unclogging and clogging events. One requirement is that δp clog c < δp break c . From equations ( 6) and ( 8), we deduce that:

δp break c δp clog c = d hi d lo . (9)
Thus, the fact that d hi > d lo automatically ensures that δp clog c < δp break c

. One can go one step further using relationships between the average width of fluid pathways d and permeability k. A common parameterization is the Kozeny-Carman equation [START_REF] Civan | Modified formulations of particle deposition and removal kinetics in saturated porous media[END_REF]:

k = λ 2 φd 2 τ (10)
where λ 2 is a coefficient that depends on the nature of the fluid pathways, e.g. tubes or slits, φ is porosity and τ is tortuosity. We may neglect variations of porosity, as discussed above, and obtain:

k hi k lo = τ closed τ open
δp break c δp clog c = d hi d lo = τ open τ closed k hi k lo 1/2 . ( 12
)
For the small porosity media that are involved here, tortuosity ratio τ open /τ closed is likely to be in a 1/2-1/4 range [START_REF] Ghanbarian | Tortuosity in porous media: A critical review[END_REF] Figure 3 summarizes how these processes are implemented in our model. In valve segments of the channel, permeability responds to the local fluid pressure gradient according to the hysteretical cycle shown in Figure 3b.

Valves as elementary seismic sources

At each valve-opening, the fluid pressure gradient that has built up is suddenly dissipated in a rapid, localized diffusion transient. Variations of fluid pressure on each sides of the valve are responsible for temporal changes of the force applied to the solid walls of the permeable medium, which are able to generate seismic waves, as shown by [START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF]. Other seismic processes may also be triggered. For example, high fluid pressures may induce the seismic rupture of critically-stressed brittle heterogeneities (e.g. [START_REF] Taetz | Fast intraslab fluid-flow events linked to pulses of high pore fluid pressure at the subducted plate interface[END_REF], and an impulsive turbulent flow transient may generate resonance in tortuous pathways [START_REF] Benson | Laboratory simulation of volcano seismicity[END_REF]. We thus consider that as they open, valves act as elementary seismic sources and generate LFE-like events.

Permeability values, boundary and initial conditions

We will consider heterogeneous permeable channels with a series of randomly distributed valves. This corresponds to one type of random permeability distribution in-volving segments where fluid pathways are more contorted than elsewhere. The behavior of any valve depends on what happens in neighboring ones, so that pressures and flow rates must be calculated everywhere in the channel.

Metamorphic reactions occur at different depths depending on the dehydrating material, which includes sediments, serpentinized mantle and hydrated oceanic crust. Focussing on the oceanic crust, dehydration proceeds at depths exceeding the tremor source region [START_REF] Van Keken | Subduction factory: 4. Depth-dependent flux of H 2 O from subducting slabs worldwide[END_REF] and most of this fluid gets trapped in the high permeability subduction interface, sealed from above by the low permeability mantle wedge [START_REF] Hyndman | Cascadia subducting plate fluids channelled to fore-arc mantle corner: ETS and silica deposition[END_REF]. Thus, the channel is fed at its base by an input fluid flux q in .

For a first instalment of the model, it is instructive to maintain this flux at a constant value. At its updip end, the channel is connected to a high permeability region which may be either continental crust [START_REF] Hyndman | Cascadia subducting plate fluids channelled to fore-arc mantle corner: ETS and silica deposition[END_REF] or an accretionary prism. In this case, pore pressures are close to hydrostatic values, such that the dynamic pressure is zero (p out = 0).

For simplicity, we assume that in the channel outside of valves, permeability remains constant at the value k hi and that valves have the same length w in the x direction. Physical variables are made dimensionless using scales described in Table 1. Values adopted for the various valve parameters are listed in Table 2. Boundary and initial conditions are illustrated in Figure 2b,c. We start all simulations with an equilibrium profile corresponding to a low value of flux q in (t < 0) = 0.09, such that all valves are closed with permeability k lo . At t = 0, q in is changed to its value for the run.

Numerical implementation

Solutions for the pressure field in time were obtained using a Crank-Nicholson finite difference scheme, with a first order, centered discretization of the second order space derivatives of pressure [START_REF] Press | Numerical Recipes -The Art of Numerical Computing[END_REF], section 20.2, pp 1045-1048). The scheme is implicit, centered and second order in time, and is unconditionally stable. Valves are ei-

ther in an open or a closed state, as described in section 3.3 above, and when the pressure differential at their boundary reaches a threshold value, their open/closed state is changed accordingly. Each opening of a valve is counted as a seismic event.

In all simulations presented here, distance along the fault is discretized in dimensionless increments δx = 0.002 and the dimensionless time step is set to δt = h 2 /2D = manuscript submitted to JGR: Solid Earth 2.10 -6 . In order to test the accuracy of the numerical implementation, we calculated analytically the pressure variations that follow an unclogging event in and around an isolated valve, and compared them to numerical solutions. Convergence was then tested by lowering both space and time steps, and we verified that accuracy improved with increasingly fine discretization. Section S1 in the supplementary information file describes the methodology and results of the convergence and accuracy tests.

Time-dependent valve behavior

We first use a single valve to show how the activity regime is controlled by the input flux and how an intermittent valving behavior requires certain conditions (Figures 4 and5). We then illustrate how neighboring valves interact using 3 valves (Figure 6).

A single valve

As mentioned above, the valve mechanism depends on five control parameters, four that describe valve behavior, k lo , k hi , δp clog c , δp break c , and the input flux q in . We identify the conditions under which a steady state solution exists, with a valve that is either open or closed. In steady state, the flux through each valve is equal to the input flux:

q in = ρ η k v δp w ( 13
)
where k v , the valve permeability, may be equal to k lo or k hi , and where δp is again the be expressed as follows, using q in , the flux a valve is subjected to:

q in ≤ ρ η k lo δp break c w = q break c
for a closed valve, ( 14)

q in ≥ ρ η k hi δp clog c w = q clog c for an open valve. ( 15
)
It introduces two threshold values for the fluid flux, noted q break c and q clog c , respectively.

The two conditions in ( 14) and (15) define domains in (δp, k v ) space corresponding to a closed valve or an open one in steady-state (Figure 4). In this case, q break c < q clog c and there are no steady-state solutions in the q break c -q clog c interval (Figure 4a). In this interval, the input flux being higher than q break interval. In this interval, steady-state solutions are possible for the two valve states (Figure 4b). Which one is selected depends on the initial state of the valve. In this case, therefore, there is no valve activity.

This analysis leads to the following instability criterion for the valve:

q break c < q clog c ⇐⇒ δp break c δp clog c < k hi k lo . (16)
In a geological context, the tortuosity ratio is small and lies in a restricted range (see section 3.3, last paragraph) and in contrast, permeability is expected to change by one order of magnitude [START_REF] Mcdowell-Boyer | Particle transport through porous media[END_REF][START_REF] Miller | Permeability as a toggle switch in fluid-controlled crustal processes[END_REF][START_REF] Beckingham | Evaluation of macroscopic porosity-permeability relationships in heterogeneous mineral dissolution and precipitation scenarios[END_REF]. Using the argument developed for equation ( 12), where we derive the threshold δp ratio as a function of the permeability ratio, we conclude that the instability criterion should be met in geological porous media where erosion/deposition processes are active. Such permeable systems have been reported to self-organize around these critical erosion/deposition conditions [START_REF] Kudrolli | Evolution of Porosity and Channelization of an Erosive Medium Driven by Fluid Flow[END_REF].

Figure 5 shows results for three values of input flux, corresponding to steady-state with a closed valve (q red < q break c ), a permanent intermittent regime with an unstable valve that switches between closed and open states (q break c < q yellow < q clog c ), and finally steady-state with an open valve (q green > q clog c to allow the valve to open at all, no activity is recorded. In the last case, q green is so high that after an initial transient during which the valve opens and closes to accomodate the incoming flux, the valve ends up open, submitted to a flux that does not allow it to close.

In the intermediate case, q yellow is both high and low enough to allow the valve to open

Valve-valve interaction via pore pressure diffusion

A simulation with three valves shows how the opening of a valve may trigger the opening of adjacent valves (Figure 6). In this simulation, the middle valve v2 opens up, driving a sudden permeability increase. As a consequence, the δp that has built up across the valve re-equilibrates through diffusion, decreasing pressure updip of the valve and increasing it downdip of it. This induces an increase of pore pressure gradient through the two adjacent valves, eventually bringing them closer to unclogging conditions. We note that such constructive interaction causes migrations of activity in cascades, in both up-and downdip directions (see section 5.2).

Emergence of collective patterns for a large number of interacting valves

We now explore systematically a model involving 29 randomly distributed valves (Figure 2b). With such a valve density, there are patches with valves that are closer together than on average, thus strengthening interactions that may generate spatial patterns of activity. All the valves are identical. We take the same permeability and threshold δp as specified in Table 2, implying that the threshold flux q clog c and q break c are the same as in previous sections (equations ( 14) and ( 15)). In fault-zones, we expect that conditions for enhanced permeability variations are met (see section 4. 1, equation (16) in particular) and hence focus on the case in which q break c < q clog c .

Characteristics of valve activity in the permanent regime

Model outputs for runs with low (q in = 0.16) and high (q in = 0.81) input fluxes are shown in Figure 7 and Figure 8. At any given time, a simple diagnosis of the state of the system is given by the bulk equivalent permeability k eq . Considering that channel segments and valves constitute hydraulic resistors in series, k eq can be written:

k eq = L (L -N cl w)/k hi + N cl w/k lo (17)
where N cl is the number of closed valves, and L is the channel length. For any of the three regimes, activity rate, k eq and the pressure difference across the channel ∆p reach a dy-namic equilibrium around a stable value after an initial transient, thus defining a permanent regime (Figures 7 and8).

The analysis of section 4.1 for a single valve does not depend on the presence of neighboring valves, and holds for a channel with multiple valves. The three permanent regimes for the 29-valve system are illustrated in Figure 9: one with intermittent valve activity and two quiescent ones. Valve activity occurs when the flux is within the q break c q clog c range. Within that range, we find that the activity style changes as a function of the input flux (Figure 10). Near the two ends of that range, activity is episodic or timeclustered (see also Figures 7 and8), and almost periodic for the highest values of q in (as shown later on, Figures 11 and12). For intermediate flux values in the middle of the flux range, activity is random and Poisson-like, with a seemingly constant overall rate. This result holds for the valve distribution displayed in Figure 2b, and for 30 other random distributions of 29-valves with the same characteristics (parameters in Table 2). This will be discussed in more detail below.

Here, the catalog of synthetic events in space and time is analyzed in the form of a two-dimensional point process (t i , x i ) illustrated in time-location graphs (Figures 7d,8d,11 and 12c). To characterize the level of time clustering (episodicity) of synthetic events, we follow [START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF] and use statistical tools for the analysis of temporal point processes [START_REF] Lowen | Fractal-Based Point Processes[END_REF]. The main idea is to compute the spectrum of the event count signal autocorrelations, and to estimate the power-law exponent γ of its high fre- ), activity appears to be close to scale-invariant clustering with a burst recurrence rate that depends on the time scale of the measurement. For higher values of the input flux close to the upper bound q clog c , activity proceeds in quasi periodic bursts, with periods that seem to increase with increasing q in .

As the input flux is increased from the lower threshold value q break c to the upper one q clog c , the channel is increasingly open on average. This may be assessed from the average proportion of open valves which increases as a function of q in (Figure 10c). The proportion of open valves dictates the evolution of k eq , the bulk permeability of the channel as a whole. For a low input flux (q in ≈ q break c ), valves are closed most of the time and activity bursts are linked to phases of opening. In the other extreme, for q in ≈ q clog c , permeability is open most of the time. Swarms of events then occur mainly when the channel shuts down temporarily and valve pressures increase until unclogging conditions are met.

Cascades and migrations of synthetic activity

In the random distribution of valves that was adopted here, there are three patches where valves are closer to one another than elsewhere, around x = 0.1, x = 0.5 and x = 0.85 (Figure 11). In patches of closely packed valves, the constructive interaction we describe in Figure 6 is responsible for the occurrence of rapid cascades of events.These migrations may proceed in both downdip and updip directions, apparently with similar propagation rates. Obviously, the interaction is weaker and slower for valves that are separated by larger distances. Thus, the rates of cascading and event migration depend on the valve distribution. , q in = 0.93. Figure 11c shows nearly deterministic patterns of collective behavior at both short time and spatial scales. Cascades of triggering/opening events are associated with migrations of activity that show up as oblique dot alignments in the time-location chart. For this run as well as others, the migration rate of cascades takes similar dimensionless values v ≈ 4 -6 × 10 2 , with faster velocities (up to v ≈ 2 × 10 3 in dimensionless units) in the densest valve patches.

The fastest migrations appear to start and stop at the edges of densely populated patches. In such a patch, valves are close to one another and act in tandem. This shows up as an activity cluster which hosts its own internal migration patterns. Larger scale valve patches are activated sequentially in the updip direction. These large scale migrations proceed at slower rates than those within a single valve patch (Figure 11b). These large-scale migrations propagate at V ≈ 22 × 10 2 in dimensionless values from one end of the channel to the other.

Large-scale valving behavior

Segments of the channel containing a number of closely-spaced valves act like macroscopic valves, such that their equivalent permeability varies in hysteretical cycles with the total pressure difference across them. This behavior is achieved over a large range of scales, from small sets of a few valves to large valve patches and eventually for the full channel. Figure 12 depicts the cycles of equivalent permeability k eq and total pressure difference across the whole channel, noted ∆p, for three cycles of channel activity. Porepressure buildup and release varies in tandem with open/closed permeability phases, thus defining large-scale valving cycles coupled with the occurrence of seismic activity.

Random valve distributions

In order to evaluate the impact of the valve distribution on the previous results, and12). We conclude that characteristic patterns of synthetic activity do not depend on the specific valve distribution, but should rather depend on valve clustering in space and the average inter-valve distance in the simulated channel.

Discussion

The key feature of the model is that pore-pressure diffusion induced by the opening of a valve triggers neighboring valve openings. In spite of the model simplicity, the generated activity exhibits complex patterns of migration and time clustering, comparable to seismic observations. An important result is that, for a given valve distribution,

i.e. a given permeability structure of the channel, the style of activity depends on the input fluid flux.

In this section, we evaluate the relevance of our model in several ways. The most basic information that can be extracted from the seismic observations consists of activity patterns and the associated time scales. Therefore, we first compare simulated and observed activity patterns from Guerrero, Mexico. Model results have been derived in dimensionless form and we then show that various observations are consistent with the same values for the variables in physical units. Using all these results, we draw general conclusions about the relations between the fluid flux and LFE activity in relation to the characteristics of the permeable subduction interface. We then discuss the impact of assumptions on valve and channel characteristics and directions for future work. In a last part, we evaluate possible extensions of the model.

Comparison with LFE activity patterns at Guerrero, Mexico

We first observe that activity migrates in both the updip and downdip directions within the sweet spot patch at rates that are faster than that from the sweet spot to the transient zone (Figure 1), as in our simulations. On the longest timescale, LFE activity in Guerrero appears time-clustered [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF][START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF]. Activity bursts recur in episodic fashion and are mostly concentrated in the downdip sweet spot patch. At times, during the most vigorous episodes, activity crosses over into the less active updip transient zone. Our model generates a similar pattern when the input flux q in is near threshold values. At the lower end of the dimensionless operating flux range (i.e. q break c -q clog c ), the channel is mostly closed and activity emerges when it opens up in episodic fashion.

In contrast, at the upper end of the flux range, the channel is mostly open and activity bursts are concentrated at the times of partial channel closure. Therefore, two different sets of values for valve opening/closing δp thresholds, valve permeability and q in could correspond to the Guerrero activity. The first situation is reminiscent of a faultvalve mechanism, as high fault zone permeability correlates in time with seismicity. Several studies interpret tremor activity in Guerrero using this correlation (Frank, Shapiro, et al., 2015;[START_REF] Cruz-Atienza | Rapid tremor migration and pore-pressure waves in subduction zones[END_REF]. In the framework of our model, however, the second situation seems as likely as the first one. We thus argue that seismicity bursts may not be systematically linked to an increase of fault permeability and subsequent fluid pres-sure surge, but may also be due to episodic phases of clogging in the fault, leading to transient fluid pressure build-up and release.

In our calculations, the lowermost patch is always active and is at the origin of larger scale migrations towards the less active upper patch, which trigger activity there. This is reminiscent of what happens in Guerrero (Figure 1). This behavior is a consequence of the boundary conditions that are imposed in the model. At the base of the channel, the fixed input flux constantly stimulates activity whereas the fixed pressure that is imposed at the top acts to buffer pressure variations. Thus, activation of the upper valve cluster only occurs once sufficient fluid pressures and volumes have built up in the lower part of the channel.

Scaling numerical results to LFE activity at Guerrero

The main parameter to constrain is the fault-zone diffusivity D = k/φβη. Using values from the literature, we end up with a wide range which is not useful. The channel porosity is φ = 0.01-0.05 [START_REF] Peacock | High pore pressures and porosity at 35 km depth in the Cascadia subduction zone[END_REF], the fluid-rock effective compressibility is β = 10 -10 -10 -9 Pa -1 [START_REF] Wibberley | Hydraulic diffusivity of fault gouge zones and implications for thermal pressurization during seismic slip[END_REF][START_REF] Shvab | Thermophysical properties of supercritical water and bond flexibility[END_REF], and the dynamic viscosity of supercritical water is η = 10 -4 -10 -3 Pa.s -1 [START_REF] Shvab | Thermophysical properties of supercritical water and bond flexibility[END_REF].

Values for the background permeability k are very sensitive to the measurement method and, above all, to the spatial scale of the measurement. Here we do not consider laboratoryscale determinations on drill-core samples and focus on measurements that capture the large-scale channel permeability. This leads to permeability values that can be as small as 10 -18 m 2 and as large as 10 -11 m 2 [START_REF] Doan | In situ measurement of the hydraulic diffusivity of the active Chelungpu Fault, Taiwan[END_REF]Frank, Shapiro, et al., 2015;[START_REF] Saffer | The permeability of active subduction plate boundary faults[END_REF][START_REF] Hendriyana | Influence of structure and pore pressure of plate interface on tectonic tremor in the Nankai subduction zone, Japan[END_REF]. With all these values, we find that diffusivity D lies in a very large 10 -4 -10 5 m 2 .s -1 range.

In order to scale model results to the Guerrero, Mexico, conditions, we need values for the characteristic spatial scale X , fault-zone diffusivity D, and the characteristic time scale T . For this diffusion process, T = X 2 /D. For the Guerrero tremor zone, X = 100 km. Using D = 10 2 m 2 .s -1 , consistent with other estimates based on tremor migrations [START_REF] Hendriyana | Influence of structure and pore pressure of plate interface on tectonic tremor in the Nankai subduction zone, Japan[END_REF], we obtain a characteristic time scale T = 10 8 s ≈ 3 years.

Using values of the hydraulic parameters that are quoted above, a diffusivity of D = 10 2 m 2 .s -1 implies a permeability of k = 10 -14 -10 -12 m 2 .

We next use the simulated activity in Figure 11 as a generic case. Such a quasiperiodic behavior allows an unambiguous determination of the characteristic recurrence time between large activity bursts, T r ≈ 0.05, and their characteristic duration, τ ≈ 0.015, in dimensionless units. Both values are representative of highly time-clustered activity, which occur for both low and high fluid fluxes, and in a large number of 29-valve distributions. Scaling variables back to physical units, we find that T r ≈ 3 months and τ ≈ 2 weeks, which are close to the estimates of [START_REF] Frank | Daily measurement of slow slip from low-frequency earthquakes is consistent with ordinary earthquake scaling[END_REF] for Guerrero: T G r = 2-3 months, τ G = 5-20 days. In physical units, the migration velocities described in section 5.2 scale up to v = 15-20 km.hr -1 for short, rapid migrations, and to V = 7 km.hr -1 for the largest migrations that proceed from bottom to top of the channel. In Guerrero, short migrations propagate indeed at velocities of a few tens of kilometers per hour [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF][START_REF] Cruz-Atienza | Rapid tremor migration and pore-pressure waves in subduction zones[END_REF], but the larger scale migrations are slower than our model estimate, with a velocity of about 1 km.day -1 . This may be attributed to the distance that separates valve patches which has a large impact on the propagation time in a diffusion process.

Using the above physical scales, synthetic activity rates range from 1 event per day for low activity periods to about 40 events per day during bursts in near-critical fluxes, and are uniformly extremely high for intermediate fluxes (Figure 10a andb). In Guerrero, average values of the activity rate are typically 600 events per day in a background regime and up to a few thousands of events per day during bursts. The activity rate obviously depends on the number of valves that are accounted for and a channel with only 29 valves does not accurately represent the larger-scale, more complex Guerrero system.

Last, but not least, we evaluate pore pressure and fluid flux values. For the fluid flux scale, Q, we use the estimate of van [START_REF] Van Keken | Subduction factory: 4. Depth-dependent flux of H 2 O from subducting slabs worldwide[END_REF] for the Cocos subduction.

The fluid flux per unit subduction zone width is Q T = 10 Tg.Myr -1 .m -1 . Using the thickness of the subduction fault zone inferred from geological observations [START_REF] Angiboust | JGR: Solid Earth June). Probing the transition between seismically coupled and decoupled segments along an ancient subduction interface[END_REF][START_REF] Tarling | Fluid overpressure from chemical reactions in serpentinite within the source region of deep episodic tremor[END_REF], h = 200 m, we obtain a mass flux per unit area equal to q T = Q T /h = 3 10 -7 kg.s -1 .m -2 . This estimate is an average value for the whole subduction zone and the actual fluid flow should be focused into channels at several points along the strike of the fault [START_REF] Ague | 4.6 -Fluid Flow in the Deep Crust[END_REF][START_REF] Piccoli | Field-based evidence for intra-slab high-permeability channel formation at eclogite-facies conditions during subduction[END_REF]. Such a channeling effect may increase the fluid flux by up to one order of magnitude. With this caveat in mind, we assume that the lower fluid flux threshold q break c , which is equal to 0.1 in dimensionless units in our simulations, is of the same order of magnitude as the dehydration meta-manuscript submitted to JGR: Solid Earth morphic flux. We deduce that q break c ≈ 3 10 -7 --3 10 -6 kg.s -1 .m -2 . We then find that the characteristic pressure scale P = η/ρk Q × X = 30 MPa -30 kPa, using the parameter ranges previously estimated.

The pressure difference between the two ends of the channel ∆p is close to 1 dimensionless flux unit most of the time. Therefore, in physical units, ∆p ≈ 30 MPa-30 kPa.

Each closed valve gets unclogged for δp = 0.04P = 1 MPa-1kPa. The magnitude of the force that is applied to the plug walls depends on area. For a valve in a channel with a 10 3 m 2 cross-section, [START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF] were able to match the amplitude of seismic waves with a pressure drop of 5 MPa. Our channel might have a larger cross-section, and hence allow for large forces with comparatively lower pressure drops. The fact that these values are of the same order of magnitude demonstrates a certain self-consistency of our analysis, backed by independent observations of seismic patterns, geological flux measurements, and seismic characteristics of LFEs.

Perspectives for diagnosing hydraulic conditions in real fault zones

According to our model, the rate and style of seismic activity depends on the hydraulic regime of the permeable channel. Provided that tremor and LFEs are seismic expressions of valve-controlled pore pressure variations, tremor patterns can be used to investigate the hydraulics of fault zones, and eventually other natural systems such as hydrothermal fields and volcanoes where similar seismogenic mechanisms have been invoked [START_REF] Honda | Periodic magma movement in the conduit with a barrier: A model for the volcanic tremor[END_REF]).

For given valve parameters that meet the condition for instability (δp break The rate of fluid release by metamorphic reactions depends on the subduction rate and on the type of material that dehydrates. It cannot be the same everywhere as, for example, subduction rates vary by almost one order of magnitude. Therefore, comparing seismicity patterns in different subduction zones is likely to prove very useful. In particular, specific attention should be paid to subduction zone segments that do not seem to generate LFE activity (e.g. [START_REF] Bocchini | Does Deep Tectonic Tremor Occur in the Central-Eastern Mediterranean Basin?[END_REF]. As shown in this paper, this may be due to a fluid flux that is either very high or very small.

Perspectives for exploring model configurations

The present hydraulic model was designed to evaluate the potential of the clogging/unclogging valve mechanism and to illustrate the various types of activity that can be generated.

It can be made more complex in many ways. As discussed above, the migration and spatial patterns of activity depend on the distribution of valves in the domain. The density of valves can be related to the roughness of the subduction interface, which may be inherited from seafloor sedimentary cover or topography prior to subduction. Valve density impacts how strongly and quickly valves can communicate through pressure transients, and might be constrained by comparing magnitude-duration or magnitude-frequency scalings of synthetic and observed LFE activity. Those scalings at the LFE source level could also help constrain valve width and breaking/clogging criteria. We do not expect, however, to discover new activity regimes as solutions are determined by the four valve control parameters introduced above.

In all the cases reported here, all valves are identical, with the same width w, the . In natural conditions, the hydraulic and mechanical properties of the permeable fault zone should depend on the ambient pressure and hence on depth. The depth dependence of valve properties may explain -within the framework of our model -why tremor bursts recur more frequently at the deeper end of LFE areas (Wech & Creager, 2011;[START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF]. As pressure and temperature increase with increasing depth, both the background and valve permeabilities should decrease under pressure. Depending on how permeability changes with depth, a flux that is near-critical for updip valves may well be in the intermediate range for the downdip ones. This may account for the behavior of the Guerrero tremor area, where activity is much more frequent and continuous downdip than updip. With lower permeability, mi-manuscript submitted to JGR: Solid Earth grations should also be slower. With that in mind, the fast migration rates observed in the sweet spot in Guerrero would imply that valve density might increase in this depth range.

Model parameters that must be explored include the valve characteristics (length, permeability in a state of clogging), the distribution of valves and the total number of valves. However, the valve instability criterion (equation 16) constrains how permeability and pressure threshold values can vary. As discussed above, we expect that this condition is met in the pore-pore throat configuration that is relevant to many natural permeable systems [START_REF] Mcdowell-Boyer | Particle transport through porous media[END_REF][START_REF] Beckingham | Evaluation of macroscopic porosity-permeability relationships in heterogeneous mineral dissolution and precipitation scenarios[END_REF]. Ultimately, a more realistic description of the process should include time-dependent opening and closure processes, but very few constraints on the associated time-scales are available.

We have shown in section 5.3 that neighboring valves may act as a single macrovalve in some conditions (Figure 12), and this is worthy of a dedicated investigation. Determining which valve spacing and hydraulic conditions enable this behavior is likely to

provide key insight into the interaction mechanisms that generate large bursts of activity. One could also allow for some random distribution of valve properties. However, we may draw from our current understanding and predict that the end result would be that only a subset of valves get activated for a given flux. Other boundary conditions than those used here are also worth investigating. For example, the subduction channel may not be connected to a very large permeability medium and allowing for an output flux that depends on pressure may lead to large-scale cyclic activity. In addition, dehydration surges may occur due for example to kinetic constraints on metamorphic reactions or to highly hydrated portions of the downgoing plate crossing the phase-change boundary. Ongoing and future work with similar models will investigate these mechanisms and further explore parametric controls on valve behavior.

Coupling between hydraulic transients and fault slip

An important direction to explore in the future is the coupling between hydraulic transients and fault slip. The latter can be described using a RSF friction parameterization with coefficients that depend on pore pressure (e.g. [START_REF] Luo | Rate-and-State Model Casts New Insight into Episodic Tremor and Slow-slip Variability in Cascadia[END_REF]. In such a framework, the main challenge is to constrain the causal relationship between fluid pressure transients, in-channel damage (valve opening), and slip. A key factor is a high pore pressure because it lowers the fault strength. As discussed in the introduction, most models of RSF fault slip consider that slip comes first and induces a permeability increase through fracturing and dilation. Another perspective to consider is the effect of changes of permeability on frictional stability, which would allow to explore how slip could be influenced by valving in the fault. Unclogging acts to lower pore pressures and hence to increase the fault strength. At the same time, it lowers the effective friction coefficient because it reduces the solid to solid contact area, i.e. the roughness of the fault interface. In this manner, unclogging may induce slip. The fluid pressure component of fault strength may therefore be modulated by the openness of valves because closed valves act as contact points between fracture walls and hence may act as high effective normal stress barriers to slip. The interactions between fault slip, damage and permeability evolution should be explored in depth in coupled physical models.

Conclusion

In

  way aperture d and permeability k are both small, with values d lo and k lo , respectively. If and when δp, the pressure difference across the plug, reaches a critical value noted δp break c such that the shear stress on the plug particles σ is equal to critical value σ c , unclogging occurs. As a consequence, pore throats are no longer obstructed by particles and permeability increases to k hi k lo , in association with an increase of d from d lo to d hi (open state). This drives a lowering of pore pressures upstream of the plug, and an increase downstream. This pressure transient is very rapid initially, and then decelerates, in a diffusion-induced flow transient. A key point is that σ initially increases. Just be-fore unclogging, it is equal to:

  τ open and τ closed are the tortuosities of the open and closed plugs. Fluid pathways are more irregular and contorted in a closed plug than in an open one, implying that τ closed > τ open . Using equation (9, one deduces that:

  pressure difference across the valve. If δp ≤ δp break c in a closed valve, stuck particles cannot be removed and the valve remains closed. Conversely, if δp ≥ δp clog c in an open valve, particles cannot be captured and the valve remains open. These conditions can

c

  , and in turn, the input flux being lower than q clog c , the now open valve will eventually close. Therefore, the valve must constantly switch between open and closed states in order to handle the imposed input flux. This oscillatory regime is illustrated schematically in Figure 3. In the other case, called enhanced pressure variations, k hi /k lo < δp break c

  manuscript submitted to JGR: Solid Earth when closed, and to close when open. After an initial transient, the valve activates in a very regular fashion.

  quency decay. A flat spectrum with γ ≈ 0 corresponds to a nearly Poissonian process, whereas γ values that are significantly larger than 0 are indicative of time clustering. As shown in Figure 10a, activity is Poisson-like for intermediate values of the input flux, such that γ ≈ 0. For values of the input flux at both ends of the q break c -q clog c range, activity is characterized by high values of γ, indicating time clustering. Time-series of the activity rate reveal episodic bursts (Figure 10b). When the input flux is low (close to q break c

Figure 11

 11 Figure 11 depicts typical synthetic activity patterns. The input flux is close to the upper threshold q clog c

  Figure 10a and c show the averaged values of the clustering exponent γ, the mean activity rate and the time-averaged proportion of open valves at each value of the input flux. We observe similar spatio-temporal patterns within individual patches and at the scale of the inter-patch distance as shown in section 5.1 for the example valve distribution (Figures11 and 12). We conclude that characteristic patterns of synthetic activ-

  /k lo , equation 16), the activity regime depends on a dimensionless input rate of fluid into the channel. This dimensionless input rate is calculated with respect to a Darcy flux for some reference permeability value (chosen to be k hi here). For application to natural conditions, the same dimensionless input rate can be obtained for different pairs of fluid flux and permeability values. Conversely, the same fluid flux may be responsible for different activity styles in different parts of the subduction interface depending on the local permeability value. In addition, time changes of the fluid flux, due for example to the subduction of more or less hydrated parts of an oceanic plate, may induce changes of activity.

  same low permeability value k lo in a state of clogging and the same values of opening/closing pressure thresholds δp break c and δp clog c

1 .

 1 this work, we have developed a model to explain the tremor patterns in subduction zones as symptoms of rapid, strong, localized pressure transients emerging from intermittent fluid transport in a permeable fault interface. We model fluid transport in a 1D permeable channel, fed by a constant input flux. Based on field observations in faults, and laboratory experiments in porous media, we have argued that clogging and un-clogging events occur in the permeable channel, causing large changes of permeability. Based on simple physical and geological arguments, we have shown that such events can alternate in a simple valve mechanism. The strong pressure transient that occurs when a valve opens is able to generate a seismic event, and we have therefore considered valve activity as synthetic seismic activity. We have found that the behavior of a permeable channel containing valves depends on two valve parameters (pressure ratio δp break c /δp clog c and permeability ratio k hi /k lo ). In one case (δp break c /δp clog c > k hi /k lo ), the channel is in a stable steady-state regime with no valve activity for all values of the input flux. In the other case (δp break c /δp clog c < k hi /k lo ), the channel is unstable with valves that open and close repeatedly if the input flux is within a well-defined operating range. We have shown that only the latter is relevant to geological systems. Therefore, in this framework, the value of the fluid flux going through a fault zone is a major control on the observed seismicity. For values of the flux outside the operating range, the fault zone is permanently closed or open. For val-ues within the operating range, permeability valves opening and closing in cycles generate permanent activity and intermittent fluid transport. In this permanently unstable regime, spatiotemporal patterns of modeled tremor emerge from valve-to-valve interactions, and are very similar to those observed in many subduction zones. Within the operating range for permanent activity, the input flux controls the dynamics of the permeable system. The system is mostly closed for low flux values and mostly open for high values, with time-clustered bursts of activity. For intermediate flux values within the range, the system is mostly open, and closes in quasi-periodic bursts of activity. Furthermore, the activity rate is highest and almost constant through time. Using this model and more elaborate versions of it, one may hope to use microseismicity patterns to probe the hydraulics of permeable fault zones and how they change with time. manuscript submitted to JGR: Solid Earth Table Characteristic dimensions used to scale physical variables Pore pressure P = (ρ r -ρ)gZ Dynamic pressure difference across channel q Massic flow rate b Q = k hi ρ/η × P/X Flux resulting from characteristic pore pressure difference a k Permeability K = k hi Background channel permeability a along a fully open channel, D = k hi /(ρηβ) b per unit channel cross section area

Figure 1 .

 1 figs/map_patterns.pdf

Figure 4 .

 4 Figure 4. Requirements for steady-state valve regimes (equations (14) and (15)) in k-δp space (in log-log scale). The two values for permeability k and threshold values for the valve pressure difference δp are shown with vertical and horizontal dashed lines. In this space, iso-flux contours appear in gray lines of slope -1 in log-log scale (-45 o ). Blue and red segments indicate values of δp that are allowed for a closed valve and an open one in steady-state, respectively. Steady-state conditions are achieved in the green and red shaded areas. See explanations in the text. (a) Valve with enhanced permeability variations (k hi /k lo > δp break c

:

  Figure 1.

  Table 2 lists values for the pres-Depending on the four valve parameters, there are two different types of solutions. Let us first consider enhanced permeability variations, such that k hi /k lo > δp break

sure thresholds and permeabilities, as well as for the flux thresholds used for the valve set up of this study. manuscript submitted to JGR: Solid Earth c /δp clog c .

Table 2 .

 2 Valve parameters for all simulations presented in this study

		Parameter	Value (scaled)
	w	Valve width	0.02
	k hi	Open valve permeability	1
	k lo	Closed valve permeability	0.05

manuscript submitted to JGR: Solid Earth Acknowledgments This study was supported by the European Research Council under the European Union Horizon 2020 research and innovation program (Grant Agreement 787399-SEISMAZE)

and by the Russian Ministry of Education and Science (Grant 14.W03.31.0033). Numerical computations were performed on the S-CAPAD platform, IPGP, France. All figures were produced using Matplotlib (Hunter, 2007). The LFE catalog in Guerrero, Mexico is made available on the Slow-Earthquake database, under the name Frank2014 (http://wwwsolid.eps.s.u-tokyo.ac.jp/ sloweq/).

 2.