Gaspard Farge
email: gafarge@ucsc.edu

Claude Jaupart

William B Frank

Nikolai M Shapiro

Along-strike

come

Introduction 1.Source processes of tectonic tremor and low-frequency earthquakes

In the deeper parts of subduction faults (at about 40 km depth), seismic signals with a lower frequency than earthquakes are detected in a wide frequency band ranging from 0.01 Hz to 10 Hz. Impulsive events called low-frequency earthquakes (LFEs) and very-low frequency earthquakes (VLFEs) are detected in narrow high (1)(2)(3)(4)(5)(6)(7)(8)(9)(10) and low (0.02-0.05 Hz) frequency bands, respectively [START_REF] Obara | Nonvolcanic Deep Tremor Associated with Subduction in Southwest Japan[END_REF][START_REF] Katsumata | Low-frequency continuous tremor around the Moho discontinuity away from volcanoes in the southwest Japan: DEEP TREMOR IN THE SOUTHWEST JAPAN[END_REF][START_REF] Ito | Slow Earthquakes Coincident with Episodic Tremors and Slow Slip Events[END_REF]. Tectonic tremor is a long (10s to days), emergent seismic signal in which LFEs and VLFEs are almost systematically detected, from the same localized sources, which persist through time [START_REF] Rubin | Imaging slow slip fronts in Cascadia with high precision cross-station tremor locations: IMAGING SLOW SLIP FRONTS IN CASCADIA[END_REF][START_REF] Chestler | A Model for Low-Frequency Earthquake Slip[END_REF]. Tremor, LFEs and VLFEs appear to be systematically correlated with one another in space and time, suggesting that they are manifestations of the same broadband phenomenon in different frequency intervals (e.g. [START_REF] Masuda | Bridging the gap between low-frequency and very-low-frequency earthquakes[END_REF]. In this study, we will study the collective, long-term patterns of LFE activity to understand their underlying physical processes. We will therefore use the terms "tremor" or "tremor activity" as a synonym for LFE activity, considering that tremor is composed of many LFEs [START_REF] Shelly | Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip[END_REF][START_REF] Ide | Empirical Low-Frequency Earthquakes Synthesized From Tectonic Tremor Records[END_REF]) -although we recognize that this approach might gloss over potentially different source processes. Tremor activity is intermittent and proceeds in periods of intense activity formed of clusters of LFEs, separated by periods of quiescence [START_REF] Obara | Nonvolcanic Deep Tremor Associated with Subduction in Southwest Japan[END_REF][START_REF] Rogers | Episodic Tremor and Slip on the Cascadia Subduction Zone: The Chatter of Silent Slip[END_REF][START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF][START_REF] Brudzinski | Nonvolcanic tremor along the Oaxaca segment of the Middle America subduction zone[END_REF][START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF][START_REF] Idehara | Regional and global variations in the temporal clustering of tectonic tremor activity[END_REF]. The temporal clustering of activity, the duration, sizes and recurrence timescales of tremor bursts vary between subduction zones. Bursts typically last from less than an hour to a week and are seperated by hours to months (e.g. [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF]. Thus, the timescales of tectonic tremor activity can be several orders of magnitude longer than the average recurrence time and duration on an individual LFE.

The geological processes that generate tectonic tremor remain elusive. The strongest episodes of tremor occur during slow-slip events on the fault interface [START_REF] Obara | Episodic slow slip events accompanied by non-volcanic tremors in southwest Japan subduction zone: EPISODIC SLOW SLIP AND TREMOR IN JAPAN[END_REF][START_REF] Rogers | Episodic Tremor and Slip on the Cascadia Subduction Zone: The Chatter of Silent Slip[END_REF][START_REF] Shelly | Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip[END_REF][START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes: UNCOVERING SILENT SLIP[END_REF] and episodes of strong fluid pressure variations in their source region (Frank, Shapiro, et al., 2015;[START_REF] Nakajima | Repeated drainage from megathrusts during episodic slow slip[END_REF][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]. This hints at recurring cycles of accumulation and release of hydro-mechanical stress in the fault often called fault-valving R. [START_REF] Sibson | Implications of fault-valve behaviour for rupture nucleation and recurrence[END_REF]. As fluid is released through metamorphic dehydration reactions of the oceanic crust [START_REF] Anderson | Geophysical and Geochemical Constraints at Converging Plate Boundaries-Part I: Dehydration in the Downgoing Slab[END_REF], regions of low permeability above the slab can seal the interface, leading to high-pressure anomalies [START_REF] Shelly | Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip[END_REF][START_REF] Calvert | Sedimentary underplating at the Cascadia mantle-wedge corner revealed by seismic imaging[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]. These lead to decreased effective normal stress on the fault, and slip can therefore nucleate with minor changes of the stress field. As the slip front propagates, fueled and stabilized by the high fluid pressures [START_REF] Segall | Dilatant strengthening as a mechanism for slow slip events[END_REF], the induced fracturing is likely to open fluid pathways, thereby lowering the local fluid pressure (Frank, Shapiro, et al., 2015). Strong tremor activity seems to occur during this phase. When the permeability around the fault heals, fluid pressure rises again, and a new cycle begins (R. [START_REF] Sibson | Implications of fault-valve behaviour for rupture nucleation and recurrence[END_REF].

The geological study of outcrops of paleo-subduction zones indicates that fluid pressures are indeed highly heterogeneous and variable, with evidence of supra-lithostatic pressures triggering vein-opening [START_REF] Angiboust | Probing the transition between seismically coupled and decoupled segments along an ancient subduction interface[END_REF][START_REF] Muñoz-Montecinos | Episodic hydrofracturing and large-scale flushing along deep subduction interfaces: Implications for fluid transfer and carbon recycling (Zagros Orogen, southeastern Iran)[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][START_REF] Tarling | A common type of mineralogical banding in serpentine crack-seal veins[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]. The radiation patterns of LFEs, VLFEs and tremor is at least in part consistent with shear slip along the fault interface [START_REF] Ide | Mechanism of deep low frequency earthquakes: Further evidence that deep non-volcanic tremor is generated by shear slip on the plate interface[END_REF][START_REF] Royer | A comparative study of low frequency earthquake templates in northern Cascadia[END_REF]Ide & Yabe, 2014;[START_REF] Imanishi | Determination of focal mechanisms of nonvolcanic tremor using S wave polarization data corrected for the effects of anisotropy[END_REF], triggered by aseismic slow slip [START_REF] Ando | A slip pulse model with fault heterogeneity for low-frequency earthquakes and tremor along plate interfaces[END_REF][START_REF] Ben-Zion | Episodic tremor and slip on a frictional interface with critical zero weakening in elastic solid: NVT-ETS and criticality[END_REF][START_REF] Sammis | A Granular Jamming Model for Low-Frequency Earthquakes[END_REF]. The accumulation of geologic evidence of coupled deformation and fluid flow at the source scale suggests that LFEs are also possibly associated with fluid movement in the fault [START_REF] Muñoz-Montecinos | Episodic hydrofracturing and large-scale flushing along deep subduction interfaces: Implications for fluid transfer and carbon recycling (Zagros Orogen, southeastern Iran)[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]. The sudden fluid mass advection could generate a single-force source on the matrix [START_REF] Takei | Why have the single force and torque been excluded from seismic source models[END_REF], producing an LFE wavefield with a radiation pattern hardly distinguishable from a double-couple owing to the small aperture of seismic networks in subduction zones [START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF][START_REF] Ohmi | Deep low-frequency earthquakes beneath the focal region of the Mw 6.7 2000 Western Tottori earthquake: DLF BE-NEATH THE FOCAL REGION OF THE[END_REF]. A similar model has been suggested by [START_REF] Ukawa | A monochromatic earthquake suggesting deepseated magmatic activity beneath the Izu-Ooshima Volcano, Japan[END_REF] to explain low frequency volcanic earthquakes that may be generated by a sudden acceleration of fluid that would result in a changing viscous shear force acting on the conduit walls and oriented parallel to the flow. In this hydraulic process, the LFE source duration is controlled solely by the fluid diffusivity of the permeable medium, and does not depend on event magnitude, which is consistent with observations in Cascadia and Guerrero, Mexico [START_REF] Bostock | Magnitudes and moment-duration scaling of low-frequency earthquakes beneath southern Vancouver Island[END_REF][START_REF] Farge | Moment-Duration Scal-ing of Low-Frequency Earthquakes in Guerrero, Mexico[END_REF]. This key property does not hold in Japanese subduction zones, however, which led some authors to attribute it to an observational bias [START_REF] Ide | Detection of Low-Frequency Earthquakes in Broadband Random Time Sequences: Are They Independent Events?[END_REF][START_REF] Supino | Self-similarity of low-frequency earthquakes[END_REF].

The very small magnitude of LFEs and restrictions on the number of seismic stations that can hamper a clear-cut diagnosis of LFE source mechanisms. Thus, we take a different approach and focus on the role fluid flow processes. In the R. [START_REF] Sibson | Implications of fault-valve behaviour for rupture nucleation and recurrence[END_REF] valve mechanism, fault slip and changes of fluid pressure are coupled. The common assumption is that abrupt pressure changes are triggered by seismic rupture but we shall argue that they could well be due to processes internal to the permeable fault zone. One key feature of LFE activity is that tremor bursts involve seismic sources of limited extent (hundreds of meters at most) and take place over large distances (tens of kilometers), showing that some long-range interaction mechanism is at work. This has usually been attributed to slow-slip propagation but simple models of tremor and seismicity show that fluid-assisted interaction between neighbouring sources can also lead to long-range synchronization [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF][START_REF] Fukuda | Model for tectonic tremors: Enduring events, moment rate spectrum, and moment-duration scaling[END_REF]. The smallest patterns of LFE clustering last less than an hour, which is consistent with source-scale processes.

Segmentation of tremor intermittence in subduction zones

Tremor activity is segmented in subduction zones: the recurrence intervals and temporal clustering are variable along both the dip and strike of the fault. In many subduction zones, tremor activity is not evenly distributed and several segments with different time clustering characteristics can be delimited along the strike of the fault (Poiata, Vilotte, Shapiro, Supino, & Obara, 2021;[START_REF] Wang | Identifying the Recurrence Patterns of Nonvolcanic Tremors Using a 2-D Hidden Markov Model With Extra Zeros[END_REF][START_REF] Brudzinski | Segmentation in episodic tremor and slip all along Cascadia[END_REF][START_REF] Husker | Characteristic Tectonic Tremor Activity Observed Over Multiple Slow Slip Cycles in the Mexican Subduction Zone[END_REF]. The geologic scale (tens to hundreds of kilometers) and permanence of such segmentation indicates that it is likely due to large-scale, structural heterogeneity in the fault zone [START_REF] Brudzinski | Segmentation in episodic tremor and slip all along Cascadia[END_REF][START_REF] Maury | Spatiotemporal Variations in Slow Earthquakes Along the Mexican Subduction Zone[END_REF][START_REF] Ide | Striations, duration, migration and tidal response in deep tremor[END_REF]. Spatial variations of frictional stress and weakness are likely to shape both the source-scale distribution of tremor sources, and the large-scale segmentation of tremor episodicity [START_REF] Nakajima | Tremor activity inhibited by welldrained conditions above a megathrust[END_REF][START_REF] Kano | Strength of tremor patches along deep transition zone of a megathrust[END_REF][START_REF] Cattania | Precursory Slow Slip and Foreshocks on Rough Faults[END_REF]. The local plate dehydration rate, magnitude of fluid escape away from the fault zone, and thus the local fluid pressure and flow conditions are also bound to control the source-and regional-scale distribution of tremor [START_REF] Audet | Possible control of subduction zone slowearthquake periodicity by silica enrichment[END_REF][START_REF] Halpaap | Earthquakes track subduction fluids from slab source to mantle wedge sink[END_REF][START_REF] Gosselin | Seismic evidence for megathrust fault-valve behavior during episodic tremor and slip[END_REF][START_REF] Mclellan | Margin-wide variations in slab dehydration in Cascadia and their relationship to slow slip[END_REF]. One can anticipate that the characteristics of segment-scale activity are rooted in the spatial distribution of seismic sources (shorter-scale heterogeneity), which has been shown to be a stationary structural feature of the plate interface [START_REF] Rubin | Imaging slow slip fronts in Cascadia with high precision cross-station tremor locations: IMAGING SLOW SLIP FRONTS IN CASCADIA[END_REF]. If fluid circulation is a major control mechanism, however, the flow rate is also likely to be a key variable.

The tremor zone in Shikoku in the Nankai subduction zone (Japan) is particularly well suited for studying the structure of the fault interface at tremor depths. Long-term, high-resolution catalogs of tremor and LFE in this region reveal large-scale spatial structures on the subduction interface [START_REF] Ide | Striations, duration, migration and tidal response in deep tremor[END_REF]. They allow to resolve a temporal clustering and tremor burst recurrence at subdaily timescales, and their segmentation along the fault (T. [START_REF] Wang | Identifying the Recurrence Patterns of Nonvolcanic Tremors Using a 2-D Hidden Markov Model With Extra Zeros[END_REF]Poiata, Vilotte, Shapiro, Supino, & Obara, 2021;[START_REF] Ide | Forecasting tectonic tremor activity using a renewal process model[END_REF]. Studies have attributed the segmentation of temporal activity to regional scale variations in frictional strength of the megathrust, due to locally drained or undrained conditions in the interface [START_REF] Kano | Strength of tremor patches along deep transition zone of a megathrust[END_REF][START_REF] Nakajima | Tremor activity inhibited by welldrained conditions above a megathrust[END_REF]. In a previous study [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF], we have designed a model suitable to investigate in more details how dynamic fluid circulation processes could shape tremor activity in this region. By modeling source distribution and interactions in the fault, it allows us to link temporal characteristics of seismicity to the dynamic of fluid circulation in the fault, and its structural controls.

In this study, our aim is therefore to understand how the temporal patterns of tremor activity in Shikoku can provide an information about the local structure of hydraulic conditions in the fault zone. To do so, we will investigate how tremor sources can synchronize, and how they conspire to shape the observed temporal clustering of activity, long recurrence timescales, and coherence of activity over wide spatial scales. We focus our analysis on the role of the dynamics of permeability, fluid flux and fluid pressure in the fault zone. We build upon an earlier study of along-dip fluid circulation in a dynamically permeable fault zone, where low-permeability plugs open and close in a valve-like manner [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF]. The valve mechanism simulates local changes of permeability and fluid pressure that may be associated with LFE source processes. We have previously shown that in such a system, realistic, tremor-like patterns of activity emerge due to the dynamics of fluid pressure. In the present work, we systematically explore how the hydraulic forcing (the input fluid flux) and the spatial distribution of low-permeability valves control the intermittence of seismicity, and how in turn, the temporal characteristics of activity can be used to infer the hydraulic conditions in the fault interface.

The article is organized as follows. (1) We describe a new method to characterize the temporal clustering and recurrence timescales of LFE activity. We apply this method to characterize the segmentation of tremor intermittence in Shikoku, Japan. (2) We summarize key aspects of the valve model and the hydraulic and seismic behavior it produces.

(3) We explore the interaction between two permeability valves. We show that its strength is controlled by the input fluid rate as well as the valve separation. (4) In a fourth part, we systematically explore how the spatial distribution of valves -their number and spatial clustering -and the input flow rate set the intermittence of tremor activity. Using simulations ot explore the flux/valve distribution phase space, we explore how they determine the strength of long-range interactions, and hence control how synchronized the activity is. (5) In a last part, we elaborate on the physical origin of tremor activity segmentation in the Shikoku (Japan) subduction zone. We end the paper with a discussion of the scope and limitations of our model.

Characterizing tremor intermittence

Clustering and characteristic timescales of LFE activity

Tremor activity can be described in two complementary ways: either through a continuous prism by measuring how long or how often a tremor source is active in time, or through a more discrete point process description, looking at the occurrences of individual events, LFEs. In a point process description of the activity, LFEs occur in a "bursty" manner, with very short and very long inter-event delays much more probable than in a constant-rate, Poisson process [START_REF] Goh | Burstiness and Memory in Complex Systems[END_REF]. They form clusters of rapid LFE occurrences (sub-hourly delays), separated by long, quiet periods (one event per day or less) [START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF]. LFE clusters -i.e. tremor bursts -last between an hour and a week. Similarly to LFEs, their occurrences cluster in time, more likely accurring soon before or after another burst. The degree of temporal clustering of activity varies between subduction zones, and between segments of a subduction zone. Some zones can exhibit an extremely clustered activity in time, with almost no background seismicity, and larger, less frequent bursts dominating the activity, while other regions can produce less clustered, more continuous activity, formed of smaller, shorter, more frequent bursts. Single bursts of activity can be spaced by a few hours in the most active periods, which recur on the scale of months to years. This superimposition of timescales of recurrence -individual events, clusters, clusters of clusters, etc -is characterized as a "scale-free" phenomenon by studies on tremor [START_REF] Idehara | Regional and global variations in the temporal clustering of tectonic tremor activity[END_REF] and LFE activity [START_REF] Frank | The evolving interaction of low-frequency earthquakes during transient slip[END_REF]Poiata, Vilotte, Shapiro, Supino, & Obara, 2021). However, when looking at the largest bursts of activity in a given region, one or several characteristic timescales of recurrence can be estimated. In the tremor zone of Guerrero (Mexico) for instance, the largest tremor bursts, associated with large (geodetic magnitude M7.5) slow-slip events recur every 4 to 5 years [START_REF] Radiguet | Slow slip events and strain accumulation in the Guerrero gap, Mexico[END_REF]. Smaller bursts occur every three months or so, activating smaller parts of the fault, associated with M6.4 slow-slip events [START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes: UNCOVERING SILENT SLIP[END_REF][START_REF] Frank | Daily measurement of slow slip from low-frequency earthquakes is consistent with ordinary earthquake scaling[END_REF]. [START_REF] Ide | Forecasting tectonic tremor activity using a renewal process model[END_REF] proposed a systematic characterisation of long-term and short-term recurrence of tremor episodes in the Nankai subduction zone, highlighting periods of 10-100 days and 1-10 hours. It is likely that the period of recurrence of tremor bursts reflects characteristic spatial scales of the dynamics at play: while the shorter period might reflect source-to-source local interactions, the longest detectable periods could be caused by largescale synchronization of whole segments of the subduction interface.

The segmentation of LFE activity in Shikoku

In Figure 1e, we represent the LFE activity in time along the strike of the Shikoku tremor zone, in Japan (catalog from Poiata, Vilotte, Shapiro, Supino, & Obara, 2021).

Three zones with different clustering and recurrence behavior can be visually identified, from least to most temporally clustered: zone 1, the westernmost part between 50-140 km along strike, zone 2, the easternmost part between 190-250 km along strike, and zone 3 in between, 140-190 km. In zone 1, activity is characterized by a relatively low level of temporal clustering, and a high background rate of events. Bursts are frequent, and recur on timescales of 2-3 months, although with a strong variability. In the easternmost zone, activity is more clustered, and bursts are slightly less frequent, recurring every 3 months in the eastern and western halves of the zone. Finally, in the boundary region between the two previous zones, the activity is least intense, with a very low background rate. Most of the events occur during large transients of activity starting in neighboring segments which propagate into the middle one, every 6 months to slightly less than In all examples presented here (Shikoku in Figure 1, Guerrero in Figure S1 and Cascadia in Figure S2), the variations of temporal clustering, regularity of the recurrence and timescale of recurrence seem to be correlated: the stronger the clustering, the longer and more regular the recurrence, and therefore the more periodic the activity. In the light of this, large clusters, quasi-periodic behavior and long recurrence intervals can be interpreted as effects of the synchronization of the activity of tremor sources through physical interactions between them. Given a constant rate of activity for many sources, the stronger the synchronization between them, the more likely events are to occur within clusters, the larger the clusters, and therefore the longer their recurrence times. Measuring synchronization through temporal clustering and the timescale and regularity of the recurrence of activity therefore provides important insight into the physical mechanisms that allow tremor sources to interact, and shape their collective behavior.

Clustering and recurrence analysis

In order to study the synchronization of LFEs in models and observations, we design simple measures of the level of clustering of activity, and of timescales of recurrence of clusters of activity, and how variable they are. Because we want to capture the largerscale dynamics of tremor generating processes, those measures are designed to capture the longer periods of tremor episodicity. They are based on an explicit detection of clusters, to avoid the ambiguity of spectral measurements -based on autocorrelation or Fourier analysis of activity time series for instance (e.g. [START_REF] Beaucé | Systematic Detection of Clustered Seismicity Beneath the Southwestern Alps[END_REF] -regarding which specific patterns the measured timescales are related to. The measures and algorithm are designed to be used automatically and with minimal tuning on both real and simulated catalogs of events, with a focus on being able to characterize both relatively constant and very clustered activity. This will allow us to gain a direct understanding of the characteristics of synthetic seismicity for large batches of simulations (several tens of thousands), with no changes of parameters. in Shikoku, Japan (2013[START_REF] Sibson | Tensile overpressure compartments on low-angle thrust faults[END_REF]. (a) 2D histogram of LFE activity in 2.7 km by 2.7 km horizontal bins, colored according to the number of sources projected onto the surface (catalog from Poiata, Vilotte, Shapiro, Supino, and Obara (2021)). Depth contours of the slab are taken from a model by [START_REF] Iwasaki | Fundamental structure model of island arcs and subducted plates in and around Japan[END_REF]. The dotted white line indicates the along-strike direction. The first measure is a temporal clustering index c, which quantifies the extent to which events occur within clusters, rather than as a more homogeneous process in time.

q) q) q) q) q2 q2 q2 q2 q2 O Q O Q O Q IZIRXWTIVFMRWUVXRSVQ (MWXERGIEPSRKWXVMOIOQ )ZIRXGSYRXIZFMR
To do so, we measure the proportion p peaks ev of events occurring during "peaks" in the time series of event counts -the number of events per time bin δt. For later reference, we chose δt = 3e -4 to count events produced in our simulations, and δt = 2 days to count events in observed tremor in Shikoku (Figure 1). Peaks of activity are defined as all consecutive bins in which the event count is higher than 66% of the time -following a similar methodology as Frank, Shapiro, et al. (2015) for LFE burst detection. Peaks therefore cover about 34% of the time period. For a constant activity rate (for instance a Poisson process), events are homogeneously distributed in time, thus about 34% of the events occur during peaks. If the activity is clustered, events occur mostly during peaks, and p peaks ev is higher than 34%. The clustering index c is constructed by mapping linearly the proportion of events in peaks p peaks ev to 0-1:

c = p peaks ev -0.34 1 -0.34 (1)
If 100% of events are in peaks, events occur only in clusters and c = 1. If only 34% of events occur in peaks, events are more or less as likely to occur in peaks as outside of peaks, the activity is not clustered and c = 0. Figure S4 and Section 2 of the supplementary materials illustrate and describe this step of the analysis in more details.

In a second step, we measure the characteristic timescales of recurrence and their regularity. We perform this analysis only if the activity is temporally clustered, when c > 0.25 -i.e. more than 50% of events are in peaks. To do so, we design an automatic identification of clusters of events in time, based on the proximity of their occurrence times t i . Clusters of events are detected using a density-based clustering algorithm (DBSCAN, [START_REF] Ester | A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise[END_REF], scikit-learn implementation [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF], choice of parameters is detailed in Section S2.3 of the Supporting Information). Clusters of events are then classified according to the number of events that make them up, using DBSCAN on cluster population size. In some cases, more than one characteristic size of clusters can be detected, and recurrence timescales should be measured only for a given size of cluster. Measuring the average of inter-cluster time intervals can give an estimation of the characteristic timescale of recurrence for a given class of clusters, and the ratio of the standard deviation to the average of inter-cluster times gives an estimate of how variable or regular the recurrence is.

By explicitly detecting clusters, and classifying them by size, we implicitly assume that the clustering of activity exhibits characteristic scales of cluster size, duration or recurrence delay, which is not always the case for LFE activity [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]Poiata, Vilotte, Shapiro, Supino, & Obara, 2021), or intermittent microseismic activity in general [START_REF] Beaucé | Systematic Detection of Clustered Seismicity Beneath the Southwestern Alps[END_REF][START_REF] Beaucé | Microseismic Constraints on the Mechanical State of the North Anatolian Fault Zone 13 Years After the 1999 M7.4 Izmit Earthquake[END_REF], where scale-free behavior is often observed. The algorithm is tailored to focus on the longest timescales of recurrence, that is on the largest clusters of activity detected in a time series. This timescale is in itself a characteristic scale of the system, and it is likely controlled of by large-scale structure of the interface.

We apply this method to the activity in Shikoku, to compute the clustering index c (Figure 1c), the recurrence timescale and its variability (Figure 1d). The automatic measurements are performed for the whole period (2013)(2014)(2015)(2016)(2017), in bins of 20 km along strike, with 10 km overlap. The segmentation in three zones that we established visually using the time-dip representation of activity (Figure 1e) is confirmed by the automatic measurements: zone 1 produces the least clustered activity (c ≈ 0.6), with the most frequent bursts (about every 2 months), zone 2 produces slightly less frequent bursts (about every 3 months), and an activity that is overall more clustered (c ≈ 0.8), and zone 3 (140-190 km) produces the most clustered activity (c ≈ 0.9), with the longest timescales of recurrence (around 6 months). This novel method therefore produces accurate and simple measures of intermittency of real patterns of seismicity, and its efficiency will be most useful to characterize synthetic activity for large batches of simulations, later in this study.

3 Model description

A valve model to describe fault zone permeability

Many geophysical measurements (e.g. [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] Wannamaker | Segmentation of plate coupling, fate of subduction fluids, and modes of arc magmatism in Cascadia, inferred from magnetotelluric resistivity[END_REF][START_REF] Rubinstein | Seismic wave triggering of nonvolcanic tremor, episodic tremor and slip, and earthquakes on Vancouver Island[END_REF] reveal that the regions of the fault zone where tremor originates has a high porosity and is saturated with fluid at near-lithostatic pressures. The fluid is mainly composed of supercritical water, freed from local dehydration reactions of the oceanic plate minerals [START_REF] Tarling | Fluid overpressure from chemical reactions in serpentinite within the source region of deep episodic tremor[END_REF], or originating from deeper sources and then channeled along the fault zone to the tremor source region [START_REF] Hyndman | Cascadia subducting plate fluids channelled to fore-arc mantle corner: ETS and silica deposition[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]. The permeability of the fault zone is characterized by a higher permeability than the overriding plate, and a strong planar anisotropy, because of the development of a fabric in the shear zone [START_REF] Evans | Permeability of fault-related rocks, and implications for hydraulic structure of fault zones[END_REF][START_REF] Audet | Teleseismic constraints on the geological environment of deep episodic slow earthquakes in subduction zone forearcs: A review[END_REF].

This permeability structure traps fluid and channels it along the plate interface, preferentially through localized channels of high permeability [START_REF] Piccoli | Field-based evidence for intra-slab high-permeability channel formation at eclogite-facies conditions during subduction[END_REF][START_REF] Eymold | Constraining fault architecture and fluid flow using crustal noble gases[END_REF][START_REF] Angiboust | Channelized Fluid Flow and Eclogite-facies Metasomatism along the Subduction Shear Zone[END_REF][START_REF] Ague | 4.6 -Fluid Flow in the Deep Crust[END_REF]. Updip of the tremor source region, fluid can partly escape the plate interface through fracture networks in the overriding plate, or at the root of the accretionary wedge [START_REF] Hyndman | Cascadia subducting plate fluids channelled to fore-arc mantle corner: ETS and silica deposition[END_REF], allowing the pressure to drop slightly from lithostatic levels. As thoroughly detailed in [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF], Section 3, we describe fluid circulation processes in a one-dimensional, high-permeability channel along the dip of the fault zone, saturated with high-pressure fluid. This channel is fed at its deep end by a constant metamorphic fluid flux q in , and is bound by a constant fluid pressure at its outer end, where it outputs to a highly connected fracture network towards the surface. We represent this system in Figure 2a andb. In a fluidsaturated porous medium, Darcy's law dictates that the mass flux of fluid per unit area of channel cross-section q(x, t) is proportional to the local gradient of fluid pressure p(x, t)

along dip, and to permeability k(x, t):

q(x, t) = -k(x, t) × ρϕ η × ∂p ∂x (x, t) (2)
where η is the fluid viscosity, ρ the fluid's density, ϕ the fault zone porosity. The fluid pressure p(x, t) evolves along the dip of the channel x in time t according to a diffusion equation:

∂p ∂t = ∂ ∂x D(x, t) ∂p ∂x (3)
with a diffusivity

D(x, t) = k(x, t) ηβϕ , ( 4
)
where β is a composite, fluid-matrix compressibility.

In the tremor source region, the combination of a constant fluid input and low-permeability seal over the fault zone creates the near-lithostatic fluid pressure. Along the fault-zone channel, permeable segments filled with high pressure fluid are separated by impermeable barriers, caused by the fault-zone heterogeneity -for instance in grain/block size, fracture density and/or thickness of the permeable zone. At equilibrium, the heterogeneous distribution of permeability results in a decrease of fluid pressure in a stepwise fashion towards the surface along the dip of the fault : hydrostatic fluid pockets separated by fluid pressure discontinuities across barriers [START_REF] Gold | Fluid ascent through the solid lithosphere and its relation to earthquakes[END_REF][START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF].

Because the input fluid flux in the fault zone is a constant stressing rate, the fluid pressure can progressively overcome the lithostatic burden just downdip of permeability barriers, forcing the dynamic increase of permeability [START_REF] Hubbert | Mechanics of Hydraulic Fracturing[END_REF][START_REF] Etheridge | High fluid pressures during regional metamorphism and deformation: Implications for mass transport and deformation mechanisms[END_REF][START_REF] Sibson | Tensile overpressure compartments on low-angle thrust faults[END_REF], and the fluid influx and rapid fluid pressure adjustment observed with both geophysical and geological means, at the source [START_REF] Muñoz-Montecinos | Episodic hydrofracturing and large-scale flushing along deep subduction interfaces: Implications for fluid transfer and carbon recycling (Zagros Orogen, southeastern Iran)[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] and subduction scale [START_REF] Angiboust | Probing the transition between seismically coupled and decoupled segments along an ancient subduction interface[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][START_REF] Nakajima | Repeated drainage from megathrusts during episodic slow slip[END_REF][START_REF] Tanaka | Temporal gravity anomalies observed in the Tokai area and a possible relationship with slow slips[END_REF].

The dynamic opening of permeability could occur because of various mechanisms: hydraulic fracturing caused by the pressure difference across a barrier generating a force overcoming the rock's cohesive strength, opening of tensile cracks as fluid pressure overcomes the lithostatic burden [START_REF] Hubbert | Mechanics of Hydraulic Fracturing[END_REF], or dilatant slip triggered by the reduced effective normal stress on friction surfaces [START_REF] Mitchell | Experimental measurements of permeability evolution during triaxial compression of initially intact crystalline rocks and implications for fluid flow in fault zones[END_REF]. After a rapid opening of permeability and relaxation of the fluid pressure, the permeability heals shut through a combination of clogging by cataclastic grains carried in the fluid [START_REF] Candela | Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing[END_REF], and mechanical and chemical cementation of fractures (e.g. [START_REF] Yasuhara | Compaction of a Rock Fracture Moderated by Competing Roles of Stress Corrosion and Pressure Solution[END_REF]. At the P-T conditions of the tremor source region, hot-pressing experiments suggest that even chemical healing could occur on the scale of days or less (e.g. [START_REF] Giger | Permeability evolution in quartz fault gouges under hydrothermal conditions[END_REF].

In our model, the small-scale heterogeneity of permeability in the channel is represented in a discrete, binary manner: short segments of the channel have a lower permeability k lo than the surrounding k hi = λk lo (λ >> 1), as represented in Figure 2b.

The low-permeability segments behave as valves. We represent the opening of those segments in a simple manner: when the pressure difference δp across them exceeds a thresh- dient and pore sizes when a particle-laden fluid circulates in it [START_REF] Jäger | The Mechanism behind Erosive Bursts in Porous Media[END_REF][START_REF] Bianchi | Critical Bursts in Filtration[END_REF][START_REF] Souzy | Transition from clogging to continuous flow in constricted particle suspensions[END_REF][START_REF] Candela | Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing[END_REF]. More details on the modeling and its consistency can be found in [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF].

A key aspect of the model is that the along-dip width of the valve w v is chosen arbitrarily, as the smallest scale of permeability heterogeneity we describe. Indeed, valves are elementary, dynamic segments of permeability, and we have shown in a previous study that their interactions seem to be creating macro-valves: segments of the channel of larger width, with a similar valving behavior [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF]. Here, our goal is to understand what controls this phenomenon, and how it affects the characteristics of seismic activity in the fault zone. The valves we use should therefore be thin enough that the effects on the seismicity are dominated by their collective behavior. Their width is set at a hundredth of the channel length along-dip, to be convenient for our numerical approach. In essence, each elementary valve can however be composed of many smaller valving segments, resulting in a macro-valve of scale w v .

In the model, permeability k changes occur without changes in porosity ϕ, which can occur if changes occur mostly to pore throats, changing the connectivity between pores [START_REF] Steinwinder | Role of Pore and Pore-Throat Distributions in Controlling Permeability in Heterogeneous Mineral Dissolution and Precipitation Scenarios[END_REF]. Permeability is therefore the only control on fault zone transport properties in our model. All physical variables are scaled to characteristic values for our system: x = 1 is the length of the channel, k = 1 is the open permeability of the channel, t = 1 is the diffusive timescale across the channel when

it is fully open, a unit of fluid pressure p corresponds to the lithostatic burden over the depth of the channel, a unit of fluid q corresponds to the flux obtained when a fluid pressure difference ∆p = p(x = 0) -p(x = 1) = 1 is applied across the full length of the channel, when its permeability is k = 1 everywhere.

We assume that each valve opening is associated with an elementary tremor event,

i.e. an LFE, emitted from the position of the valve, at the time of opening. It simulates the opening of permeability that is presumably associated with an LFE triggered by locallyhigh fluid pressure, through a rupture, unclogging, or hydrofracture-like event (e.g. [START_REF] Kotowski | Length scales and types of heterogeneities along the deep subduction interface: Insights from exhumed rocks on Syros Island, Greece[END_REF][START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF][START_REF] Muñoz-Montecinos | Episodic hydrofracturing and large-scale flushing along deep subduction interfaces: Implications for fluid transfer and carbon recycling (Zagros Orogen, southeastern Iran)[END_REF]. The cumulated activity of those synthetic LFEs could build up into signals that could be detected as large magnitude LFEs, VLFEs, and tremor. This description is simplistic, but our fo-cus is not on an accurate physical description of the source of LFEs, VLFEs or tremor, but rather on the sensitivity to and effect on the fluid pressure field of elementary tremor sources, that allows source-to-source interaction.

In practice, we built and optimized a stable, accurate algorithm to solve the diffusion equation, compute fluid flux, valve state and therefore permeability throughout the channel [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF]. The main control parameters of this system are the value of the input fluid flux q in , and the valve distribution in the channel. We retrieve a catalog of valve openings, which we consider to be an analog of LFE catalogs. Figure 2c displays the activity in time and space emitted by valves in a channel, as fluid is forced through it. To our knowledge, this model is the first to explicitly and causally describe the link between the hydraulic dynamics in the fault zone and the activity of tremor in time and space, opening a rich scope of exploration.

Valve behavior for different values of the input flux

Closed valves act as barriers to fluid flow in the channel, and as they impede the flow in the channel, the pressure difference across them can increase to the point they break. However, the fluid input rate q in into the system controls if a valve can open when closed, and close when open. Two threshold flux values can be derived from the δp thresholds, using Darcy's law for the flux through a valve in a permanent regime. The flux needed to reach the breaking threshold δp break c when a valve is closed is:

q break c = ρ η k lo δp break c w v , (5)
where w v is the valve's width. Therefore, when q in > q break , and therefore q in has to be lower than: . In this case, values of flux determine three distinct valve behaviors:

q clog c = ρ η k hi δp clog c w v . ( 6
• When q in < q break c , the flux is lower than the breaking threshold, and valves will stay closed, or close if they are open.

• When q in > q clog c , the flux is higher than the closing threshold, and valves will stay open, or open if they are closed.

• When q break c < q in < q clog c , a closed valve will open because q in > q break c , and as q in < q clog c , the now-open valve will eventually close. In these conditions, valves are permanently unsteady, opening and closing in cycles, generating sustained activity.

In experiments with model porous media, the permeability state (clogged, unclogged or variable) undergoes similar phase transitions depending on the fluid input rate in the system [START_REF] Jäger | The Mechanism behind Erosive Bursts in Porous Media[END_REF][START_REF] Bianchi | Critical Bursts in Filtration[END_REF].

Figure 3 upper panels (a-f) display the cycle of pressure and permeability for a valve, for three values of q in in the range allowing sustained opening and closing

q break c < q in < q clog c
. The closer the input flux q in is to either threshold values, the closer the valve will be to a permanently closed or open state, as seen in the time series of valve permeability (Figure 3b,d,f). For instance in Figure 3a it. We will show that as two valves of width w v are close enough, they start behaving as a macro-valve of larger width. Therefore, the width w v should be considered the elementary width at which the heterogeneity of permeability is defined, but other scales of heterogeneity will emerge in systems of many valves where spatial cluster of valves behave as macro-valves, and we fix w v = 0.01 for all valves in the rest of the study. As permeability is the only transport property that is variable in space, the value of λ describes how quickly fluid pressure diffuses inside a valve relative to outside of it in the high permeability channel. When this ratio is very high, diffusion within the valve is much slower than outside the valve, and the effects of the pressure variations outside the valve will thus dominate the dynamics. In this study, we take λ = k hi /k lo = 20, for which closed valves virtually behave as barriers to the fluid diffusion, and the exact value of this ratio essentially does not influence valve dynamics. Our final goal being to understand how large-scale, complex patterns of valve activation emerge from interactions between elementary valves, all valves in the rest of the study have fixed parameters (Table 1). The main control parameters on valve activity will therefore be the value of the input flux q in relative to the opening and closing thresholds, and the valve distribution -their number N v and position relative to one another in the channel.

4 Characterizing source interaction in a two-valve system 4.1 Valves interact through pressure transients

In order to understand how valves interact through pressure transients and which parameters control how strong and fast this interaction is, we first design systems with two valves, and observe how they evolve. Valves are at equal distance from the center of the channel x = 0.5 in normalized length units. The valve width is w v = 0.01, and the distance between them is d v , normalized to w v (Figure 4a).

In Figure 4, we describe a simulation in which a system with two valves at a distance of d v = 0.8w v is subjected to a low flux q in = 0.13 -close to q break c , the critical flux above which closed valves will open. In the permanent regime and independently of the initial conditions, we observe that the cycles of pressurization and release of both The triggering interaction between valves is therefore carried by the fluid pressure field in the permeable system. Because it is a diffusive system, the further the triggering pressure transient has to travel to the neighboring valve, the more the interaction between the valves is damped and slow. In other words, the distance between valves d v should be an essential control on the interaction strength and synchronization.

The sensitivity of valves to small variations of pressure around them is also directly related to their criticality: their being close or far to the threshold of δp that will make them change state -open or close. As shown in the previous section and Figure 3, the amount of time at the valve spends close to threshold is controlled by the value of the input flux q in . The closer q in is to critical values q break c or q clog c , the closer to thresholds the valves are on average. In Figure 4b andc, the input flux is q in = 0.13 ≈ q break c = 0.1 and it is visible that both v1 and v2 approach δp break c tangentially. Thus, the closer q in is to critical values, the stronger the interactions between two valves should be, as the smallest change of pressure is likely to trigger opening of the valve when it is closed, and closing of the valve when it is open.

Valve spacing controls activity synchronization

In order to characterize the influence of the intervalve distance d v on the synchronization of valve activity, we design two-valve systems with various d v , spanning d v = 0 to 20w v . All systems are subjected to a constant flux, and we characterize the event recurrence timescales when the permanent regime is reached. Various values of flux are tested for all systems, and in Figure 5, we display the results of this experiment for three domains of q in , and for the full range of intervalve distances. As the two valves are identical, during the permanent regime of all simulations, the events occur in repeated sequences of alternating events: a two-event burst (Figure 5).

Two timescales therefore characterize the activity:

• δt intra the time interval between the first and second opening of the burst

• δt inter the time interval in between the beginnings of two neighboring bursts, which can be thought of as the burst recurrence time.

In Figure 5b andc, we show the evolution of these two timescales, normalized by δt 0 the period of the cycle of an isolated, identical valve subjected to the same flux.

The first observation is that at large distances (d v > 5w v , Figure 5b), burststherefore events for each valves -occur with a recurrence corresponding to an isolated valve: the cycle of either valve is not affected by the other valve.When the distance between the valves gets shorter, the interaction between valves becomes stronger, producing two effects on the valves' cycle. (1) In Figure 5c, we observe that for all fluxes, the closer the valves are (the lower d v ), the shorter δt intra is. As valves get closer, the events are closer in time, and valves synchronize.

(2) In Figure 5b, we observe that as valves get very close (d v < 0.5w v for q in = 0.11), δt inter increases for lower and lower d vfor all three domains of flux. In other words, the closer the valves, the longer the recurrence delay between bursts, and the less frequently valves activate. For such close proximity between valves, the inter-event time is so short that the two valves together can be considered as a twice-wider macro-valve. Interestingly, in Figure 5, we see that for all fluxes, the macro-valve consistently activates with about 4 times the delay of an isolated, elementary valve. We reported the emergence of similar macro-valving behavior Valve spacing effect on valve-to-valve interaction a b in systems of many strongly-interacting valves in our previous study (Farge et al., 2021, 574 Section 5.3).

575

A comparison of the curves corresponding to the low, medium and high input flux 576 ranges in Figure 5b and c allows to assess the effect of the q in on the interaction strength.

577

In Figure 5b, we see that for low and high q in ranges (red and green curves), which are be made for the delay between events in bursts δt intra : shorter valve spacings are nec-essary in the medium flux range to reach a similarly short delay between events in bursts, compared to the high and low, closer to criticality, flux ranges.

Here, we have demonstrated that the interaction between two identical valves are stronger when (1) they are closer together, (2) the flux q in into the channel is close to the thresholds of opening or closing, putting valves in a near-critical state. As interactions get stronger, the valves synchronize in time, and their behavior starts resembling that of a larger, less active macro-valve, with different effective permeability and triggering thresholds.

5 Emergence and variability of synchronization in complex valve systems

Control parameters and simulation setup

In reality, the heterogeneous permeability of a fault-zone channel should be composed of a complex spatial distribution of valves. In such systems, the synchronization of sources is more complex, and could rely on intricate hydraulic interactions. Based on the results for the two-valve experiments, we investigate how the input flux in the system and the valve distribution shape interactions in complex N v -valve systems, and if they are conducive to intermittent, clustered activity, and quasi-periodic bursts. To do so, we run simulations in complex channels with many identical valves, systematically testing how the valve distribution along-dip and values of q in relative to the opening and closing thresholds (q break c , q clog c ) affect the modeled tremor activity.

In the previous section, we have shown that the distance between two valves d v controls the strength of interactions for a given input flux. In a system with N v >> 1 valves, valves can be close together either because there are many of them in the system, or because of their spatial clustering which results in dense patches separated by less populated regions. In order to test the effect of density and spatial clustering, we design systems where valves are distributed using a Weibull distribution for the inter-valve distance

d v .
The probability density function of the distribution is:

p(d v ) = u d 0 d v d 0 u-1 × exp ((-d v /d 0 ) u ) , (7)
and its mean

d v = d 0 Γ 1 + 1 u . ( 8
)
where d 0 is the scale parameter, u the shape parameter of the distribution, and Γ is the gamma function.

The Weibull distribution provides a straightforward parameterization of valve density and spatial clustering in the system, allowing to choose between perfectly regular and highly clustered distributions solely with the choice of the shape parameter u in equa- A valve system is built by first specifying a buffer zone on each sides of the system 630 to reduce edge effects, common for all distributions, then by randomly drawing inter-valve 631 distances using the chosen parameters (u and d v ), and finally distributing them on the 632 regular space grid, by simply rounding down the inter-valve distance to the closest dis-633 crete distance possible. We also ensure that the spacing of the last valve to the one be-634 fore that is not too far from the target inter-valve distance, so that all valves do not end 635 up near the input for the most clustered valve distributions. Figure 6 We run simulations for distributions described by

v = 9.5w v N v = 9 u = 0.3GPYWXIVIH u = 14SMWWSR u = 9VIKYPEV d v = 5.0w v N v = 15 d v = 2.0w v N v = 29 d v = 1.0w v N v = 43 x d v = 0.2w v N v = 71 x x 1SVIGPYWXIVIH (IRWIV

Parameterizing valve distributions

N v = 9-71 (d v = 9.5w v -0.2w v , 640 16
values), u = 0.2-9 (8 values, 128 theoretical distributions in total). We randomly 641 generate 30 distributions for each (N v , u), in order to average out the effects of specific 642 valve arrangements on our results, and obtain statistical significance. In all systems, valves 643 have identical width, closed permeability and opening/closing δp thresholds (Table 1). sure is kept at p out = 0.

648

In addition to pressure, flux and permeability in time and space, we record valve 649 states in time and space, and build a catalog of valve activations to simulated a seismic-650 ity for each simulation. This catalog is then automatically characterized by computing:

• the activity rate, measured by the average event count per valve in a given time increment δt,

• the extent to which activity occurs in clusters with the clustering index c (c = 0 events are not clustered at all, c = 1 all events occur within clusters),

• the recurrence timescales of clusters in cases when the activity is clustered in time (c > 0.25), measured as the average delay mean k (∆T k ) between clusters of a same size,

• and how stable in time the recurrence is, using as an indicator the inverse of the coefficient of variation of the measured delays mean k (∆T k )/std k (∆T k ). When mean k (∆T k )/std k (∆T k ) > 2, a simulation is labelled "periodic".

The details of how the latter three measures are computed can be found in Section 2.3.

Control of the input flux on activity synchronization

In the previous sections we showed that when the input flux q in is close to the thresholds of opening (q break c ) and closing (q clog c ), valves are more sensitive to small changes of pressure. In Section 4, we have shown that this results in strong interactions between valves in two-valve systems subjected to near-threshold q in . With stronger interactions, valve activity tends to synchronize, producing less frequent, more synchronous bursts of events. In [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF], we demonstrated that near-threshold fluxes produce a more clustered seismicity in valve systems with N v = 29 valves, with u = 1 (Poissonian valve distribution). Here, we will generalize this analysis to the full spectrum of densities and spatial clustering of valves.

We start by focusing on the activity produced in a clustered valve system, u = 0.3,

with N v = 29 valves (d v = 2w v )
, for three values of the flux: (a) q in = 0.13, near the opening threshold of valves q break c , valves are mostly closed in average [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF], (b) q in = 0.57, far from thresholds, (c) q in = 1.02 near the closing threshold of valves

q clog c
. Figure 7 displays the results of the experiment. The first observation is that the results from the two-valve system generalize to this system with N v > 2. In (a) and (c), q in is close to a closing or opening threshold, and activity is intermittent, organized in clusters of events in time, which recur quasi-periodically, whereas when q in is far from both thresholds (in (b)), activity is more continuous, and less clustered. In this system at least, the input flux controls the synchronization of valve activity: the closer to an opening or closing threshold it is, the more clustered and periodic the activity seems to be.

A second observation is that coherent patterns of activity build over larger distances when the flux is closer to threshold. Indeed, when q in is far from both thresholds, activity is intermittent, but it synchronizes on a much smaller spatial scale than when q in is nearthreshold. When valves are more critical, triggering interactions can cascade farther in the system, as most valves are in the same state: close to threshold, and therefore sensitive to the opening or closing of their neighbor. When the flux is farther from threshold, valves cannot synchronize as well in time and space, as there is less chance that the pressure transient when a valve opens will reach a near-threshold valve. As seen in the two-valve system, the long-range interaction allowed by near-critical fluxes allows largescale, long-period cycles of activity for large clusters of valves (0.5 < x < 1 in Figure 7a.ii).

As valves synchronize in space and time, a large-scale valving behavior emerges, favored by valve criticality.

Using the full range of our simulations, we extend the analysis to all valve systems, to observe how they respond to values of the flux in between the closing and opening thresholds. Figure 8 synthesizes our results. For all valves systems, dense or sparse, clustered or regular, the closer q in is to critical values of the flux, the closer valve systems are to staying closed (q in ≈ q break c ) or open (q in ≈ q clog c ). In Figure 8a, we see that valves therefore produce less events just as a result of flux conditions, but also because of valveto-valve interactions. Interactions between neighboring valves force them into a less fre- quent activation (Figure 5). The closer q in is to thresholds, the more synchronized the activity gets. Valve activity synchronizes in clusters (c > 0.5), and long-range interactions produce macro-valving behavior, which tends to make the activity more periodic, and with longer periods (Figure 8b andc). As in real systems, periods of burst recurrence (T ∼ 10 -2 -1) are orders of magnitude longer than the period of individual valve cycles (T ∼ 10 -4 -10 -3 ). When valves in the system experience a flux that is further from critical conditions, they open and close rapidly, without time to interact and build synchronicity: activity is high, less clustered, less periodic, with shorter periods.

In subduction zones (e.g. in Shikoku, Figure 1), most regions exhibit strong temporal clustering (c > 0.7), timescales of burst recurrence several orders of magnitude longer than individual source recurrence, and coherent activity in space, over scales larger than the individual source scale. Such conditions are only found in our system for input fluxes into the fault zone that are near-threshold, both low and high. However, it seems unlikely that the permeable system in the fault zone could stay mostly open. Because of the high temperature and pressure, ductile deformation of pores and crystallization processes in the pores would rapidly shut most of the pathways [START_REF] Tarling | A common type of mineralogical banding in serpentine crack-seal veins[END_REF], and form pockets of fluid, separated by low-permeability barriers [START_REF] Gold | Fluid ascent through the solid lithosphere and its relation to earthquakes[END_REF]. The near-threshold, low-flux regime (q break c < q in < 0.3) therefore seems to be the more realistic part of the flux domain. According to our model, in real subduction zones, a small decrease in flux (between spatial regions or in time) could therefore result in a monotonic decrease in seismicity, and favor interactions between permeability valves, thus further lowering seismicity rates, and producing a more clustered, more periodic activity, synchronizing over larger distances and longer periods. For the rest of this study, we will thus focus our analysis on low, near-opening-threshold values of the input flux, q in < 0.3.

Control of the valve distribution on activity synchronization

In Section 4.2, we showed using a simplified system that the interaction between two valves gets stronger and faster as the distance between them decreases. Those effects should persist in a more complex system: the closer the valves get in the system, the more synchronized in time and long-period the activity should be. In order to independently demonstrate the effect of an increase in valve density (d v and N v ) and spatial clustering (u), we use the activity in three valve systems, (a) a rather homogeneous system, with tered valve distribution u = 0.3, and finally (c) a denser system (N

N v = 29 valves (d v = 2w v )
v = 43, d v = 1w v ),
with the same spatial clustering as in the previous system (u = 0.3). All three valve systems are subjected to a low input flux q in = 0.13 ≈ q break c . Figure 9 . This shows that as valves get closer in the system, their interactions generate stronger interactions, which build clusters of activations. As elementary valves form spatial clusters in the channel, their collective behavior builds macro-valves, which produce the bursts of activity seen in system (b)

for instance (x ≈ 0.1 and x > 0.5). The most striking result here is that even though valve clusters in the channel are far apart, the collective effect of a macro-valve on the pressure field is such that it builds long-range interaction between valve clusters, which synchronize valve clusters throughout the channel. This effect occurs in system (c), where a large-scale valving behavior emerges at the scale of the whole channel, as the four valve clusters synchronize their activity. In our previous study, we have shown that such cycles of activity are associated with cycles of permeability opening and closing, fluid pressure accumulation and release, similar to the cycle of a single valve [START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF].

The different scales of macro-valving allowed by interactions -cluster of valves (system (b), x < 0.2), cluster of valve clusters (system (c), full channel) -build increasingly long periodicity of activity, as fluid has to diffuse from the input throughout the closed system to load all valves to the point of breaking before a collective reopening of the system.

This test is therefore consistent with our hypothesis. The closer the valves get, either due to spatial clustering (lower u) or more numerous valves (higher N v , lower d v ), the more synchronized their activity is: they all activate at the same time during bursts, and are inactive outside of those episodes. In addition, the activity seems to also display longer recurrence timescales as the valves are closer locally in the system.

A simple measure of how many valves are interacting in the system N int v should capture how u and N v both affect the intermittence of activity, and conveniently reduce the dimension of the problem. We define N int v as the number of valves that have a neighbor closer than 0.5w v , on either side. The value of d v < 0.5w v as an interacting distance is chosen on the basis of the the two-valve experiments described in Section 4: in Figure 5c, it is visible that valves interact very strongly at distances lower than d v = 1w v and that this effect is even clearer for d v > 0.5w v . We choose the lower bound, d v < 0.5w v for the interacting distance. When the activity characteristics are plotted along N int v , we capture the effects of the local proximity of valves in the system, either due to spatial clustering or overall density. Figure 10 shows systems with different N int v , show-ing that a similar number of valves can be at interaction distance when the system is dense and not very clustered, and when it is less dense but with stronger spatially clustered valves.

The results outlined in the previous paragraphs stand for all valve distributions.

In Figure 11, the activity rate, temporal clustering and periodicity are plotted as a function of the number of interacting valves N int v , for all valve distributions, subjected to a low value of flux, close to the closing threshold q in = 0.13. The three simulations from 9 are shown as colored dots in panels a, b, and c, and represented again in d, e and f. The signature of increasing interactions in the system is visible as the number of valves at interacting distance rises in the system, due to both more clustered or denser distribution. The closer valves get in the system, the stronger their interactions, and the more clustered their activity gets (c > 0.5, Figure 11b). As they interact strongly, they form macro-valves, their activity rate drops (Figure 11a), and they activate in clusters with a more and more periodic behavior (Figure 11c). Finally, as more and larger clusters form with increasing density and spatial clustering of valves, the activity synchronizes on larger distances, and becomes periodic on longer timescales (Figure 11c). Strikingly, the various space scales of synchronization (macro-valves) seem to be reflected in the detected periodicities of activity: both short-(T ≈ 3 × 10 -2 ) and long-period (T ≈ 3 × 10 -1 ) are detected, corresponding roughly to periodicities of small valve clusters (of width ∼ 0.2), and of the full channel -e.g. system (c) in Figure 9. Once more, it is noteworthy that the periods of activity bursts (5×10 -3 to 1 scaled time units for most detected periods) that emerge are orders of magnitude longer than the period of an isolated valve cycle (8.1 × 10 -4 scaled time units for an isolated valve with q in = 0.13).

The noticeable drop in number of periodic simulations for N int v = 45-50 is an artifact due to the way we design valve systems. In order to test perfectly regular valve systems, the tested values for N v jump from 43 to 57. As 45-50 are systems which degree of spatial valve clustering is low (u > 1), putting valves far apart. In this case, valve clusters in space are few and small, which results and a less periodic activity, and lower periodicities -when they emerge.

N int v ≤ N v , systems with N int v = x (IRWIV QSVIVIKYPEV N v = 29ZEPZIWd v = 2.0w v u = 1.0 x N v = 43ZEPZIWd v = 1.0w v u = 1.0 x N v = 57ZEPZIWd v = 0.5w v u = 2.1 x 0IWWHIRWI QSVIGPYWXIVIH N int v = 12 N v = 13ZEPZIWd v = 6.0w v u = 0.3 x N int v = 24 N int v RYQFIVSJZEPZIWEXMRXIVEGXMRKHMWXERGId v < 0.5w v N v = 29ZEPZIWd v = 2.0w v u = 0.3 x N int v = 41 N v = 43ZEPZIWd v = 1.0w v u = 0.3

Parameterizing valve proximity

Discussion

Linking model parameters and subduction zone hydraulic properties

The spatio-temporal patterns of intermittence of tremor in subduction zones gives the most direct access to minute deformation and hydraulic processes occurring within or around the subduction interface (e.g. [START_REF] Bernaudin | Episodic Tremor and Slip Explained by Fluid-Enhanced Microfracturing and Sealing[END_REF]Cruz-Atienza et al., 2018;[START_REF] Luo | Rate-and-State Model Casts New Insight into Episodic Tremor and Slow-slip Variability in Cascadia[END_REF]. In this work, we propose a simple representation of the hydraulic processes in the fault based on elementary permeability-valve processes, in order to investigate the role that hydraulic dynamics may play in shaping the intermittence of tremor. The main characteristics of tremor we seek to reproduce and understand are the emergence of the strong temporal clustering of events [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]Poiata, Vilotte, Shapiro, Supino, & Obara, 2021), the long, quasi-periodic recurrence of tremor bursts [START_REF] Brudzinski | Segmentation in episodic tremor and slip all along Cascadia[END_REF][START_REF] Frank | Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico[END_REF][START_REF] Husker | Characteristic Tectonic Tremor Activity Observed Over Multiple Slow Slip Cycles in the Mexican Subduction Zone[END_REF]Poiata, Vilotte, Shapiro, Supino, & Obara, 2021), and how these two phenomena vary along-strike in tremor source regions of subduction zones, at scales of tens to hundreds of kilometers [START_REF] Brudzinski | Segmentation in episodic tremor and slip all along Cascadia[END_REF]Poiata, Vilotte, Shapiro, Supino, & Obara, 2021). The model we developed reproduces such patterns of seismicity because of interactions between elementary seismic sources through fluid pressure transients in the permeable fault interface. We showed that these interactions are the basis of the emergence of time and space scales of coherent activity much wider than the time and space scale of activation of a single, isolated source. This feature of our results is valid for a wide range of the control parameters in our system, and it therefore seems that large-scale, spatio-temporal patterns of LFE activity in reality could be influenced by slow hydraulic interaction between sources. In the system we present here, interactions are controlled by the criticality of valves (q in vs valve thresholds) and valve distribution.

We argue that in reality valves should spend most of their times closed, and the input flux be above, but near their opening threshold q break c , producing highly intermittent, on-off activity (cf Section 5.2). In our system, the distance between q in and q break c controls the intermittence of activity. It can change due to changes of (1) q in , (2) valve breaking threshold δp break c , and finally of (3) the permeability of the valve when closed k lo . More generally, those parameters represent (1) a hydraulic stressing rate, (2) a hy-draulic weakness (or strength), and (3) how much the valve translates flow into stress, a kind of stress-leakage term. The segmentation of tremor intermittence in a real subduction zone could therefore rely on spatial variations of those parameters. The variations of fluid input is due in part to spatial variations in focusing of the flow in largescale permeable channels [START_REF] Piccoli | Field-based evidence for intra-slab high-permeability channel formation at eclogite-facies conditions during subduction[END_REF][START_REF] Eymold | Constraining fault architecture and fluid flow using crustal noble gases[END_REF][START_REF] Angiboust | Channelized Fluid Flow and Eclogite-facies Metasomatism along the Subduction Shear Zone[END_REF][START_REF] Ague | 4.6 -Fluid Flow in the Deep Crust[END_REF], perhaps dug by past subduction of seamounts [START_REF] Ide | Striations, duration, migration and tidal response in deep tremor[END_REF]. The thickness of the subducting crust, speed and angle of the subduction also controls dehydration rate [START_REF] Manea | A geodynamical perspective on the subduction of Cocos and Rivera plates beneath Mexico and Central America[END_REF]. Variations in temperature can affect how well the fault interface is sealed from above by silica deposits [START_REF] Audet | Possible control of subduction zone slowearthquake periodicity by silica enrichment[END_REF], and therefore the importance of sinks and leaking along the fault-zone, which competes with loading rate, inhibiting valve activity (e.g. [START_REF] Halpaap | Earthquakes track subduction fluids from slab source to mantle wedge sink[END_REF]. Finally, the intrinsic properties of the fault zone govern the strength and permeability of valves, and can strongly be affected by the geology of the subducting and overriding plates, and local level of damage, controlling both the availability of pathways for the fluid to circulate and the plasticity of the fault zone material. Other studies have suggested that local increase in flux and availability of high fluid pressures might increase the frequency of tremor bursts [START_REF] Mclellan | Margin-wide variations in slab dehydration in Cascadia and their relationship to slow slip[END_REF][START_REF] Wech | A continuum of stress, strength and slip in the Cascadia subduction zone[END_REF], but our work stresses that the intensity of q in has to be compared to the relative contribution of strength and sinks in fluid pressure accumulation in the fault zone.

Indeed, the emergence of temporal clustering and large time and space scales of tremor activity also relies on the number of interacting sources, which in our description is directly based on the permeability structure of the fault zone. Our results tend to show that long-range synchronization and longer timescales are built through the collective behavior of dense patches of permeability valves in the channel (Section 5.3). Macro-valves built out of these patches have a strong effect on the pressure field. They are wide, mostly closed segments of the permeable channel in the fault zone, and by opening all at once, they generate strong, long-period, long-wavelength variations of fluid pressure, which can diffuse farther than the transients generated by elementary valve openings, and therefore can synchronize the valving and seismic activity of wide parts of the fault zone. In Section 5.3, we show that for more than N int v = 20 valves at interacting distance, the produced seismicity is highly temporally clustered, quasi-periodic, on timescales orders of magnitude longer than single-valve cycles. For the valve widths used here w v = 0.01, this corresponds to at least 20 % of the system behaving in a valve-like manner, with significant spatial clustering of those elementary segments. If it were the case, the observed patterns of tremor activity could be driven by large-scale, long-period fluid pressure and permeability transients in the fault zone, built on elementary valving processes. In other words, N int v is a rough description of the heterogeneity of permeability of the fault zone.

It is linked with the amount of larger scale irregularities of permeability, which end up behaving as macro-valves. Ridges, seamounts or fracture zones dominate the kilometric roughness of the subduction plane, and should define macro structures of permeability in the fault. Such structures generate kilometric-scale heterogeneities of stress, damage, permeability, and mechanical properties of the fault zone. As tremor sources seem to also be spatially clustered in patches [START_REF] Rubin | Imaging slow slip fronts in Cascadia with high precision cross-station tremor locations: IMAGING SLOW SLIP FRONTS IN CASCADIA[END_REF][START_REF] Ide | Striations, duration, migration and tidal response in deep tremor[END_REF][START_REF] Chestler | A Model for Low-Frequency Earthquake Slip[END_REF], the spatially-clustered Weibull distribution (u < 1) seems to reasonably approximate the physical and observational characteristics of the tremor zone.

Plate topography and tremor segmentation in Shikoku

Long-wavelength geophysical observations in subduction zones, of magnetic or gravimetric anomalies for instance, reveal that subduction zones are structurally segmented on a scale of tens to a hundred of kilometers. This segmentation seems to emerge mostly from the topographic and internal structure of the subducting oceanic plate [START_REF] Blakely | Subduction-zone magnetic anomalies and implications for hydrated forearc mantle[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][START_REF] Wang | Invited review paper: Fault creep caused by subduction of rough seafloor relief[END_REF][START_REF] Shillington | Link between plate fabric, hydration and subduction zone seismicity in Alaska[END_REF][START_REF] Bassett | Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief: CRUSTAL STRUCTURE AND SEISMICITY: 1[END_REF], which displays heterogeneous structures -ridges, seamounts, fracture zones -on such a scale. As rough terrain enters the subduction, it modifies stress, damages both the subducting and overriding plate, carries more sediments and water into the subduction. Those factors can directly affect how seismicity manifests in the rough segment by modifying the mechanical properties of the fault zone. And indeed, the structural heterogeneity of the subducting plate often correlates with the spatial variations of seismicity and tremor along strike in subduction zones, on a very similar scale as the observed topographic features of the incoming seafloor [START_REF] Blakely | Subduction-zone magnetic anomalies and implications for hydrated forearc mantle[END_REF][START_REF] Wang | Invited review paper: Fault creep caused by subduction of rough seafloor relief[END_REF][START_REF] Shillington | Link between plate fabric, hydration and subduction zone seismicity in Alaska[END_REF][START_REF] Bassett | Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief: CRUSTAL STRUCTURE AND SEISMICITY: 1[END_REF][START_REF] Ide | Striations, duration, migration and tidal response in deep tremor[END_REF]. The topography of the incoming plate could directly affect the permeability structure, valving properties and channeling of the fluid in the fault zone. In the next paragraph, we use the results of our model to try to link the segmentation of activity in Shikoku with geologicscale properties of the plate interface, and interpret how the topography of the subducting Philippine Sea plate could be the underlying cause of these spatial variations.

In Section 2 (Figure 1), we used our novel characterization techniques to propose a kilometric-scale segmentation of tremor intermittence in Shikoku. Three segments can be identified: the first two segments (segment 1 and 2) at each along-strike extremity (d < 140 km, and d > 190 km, Figure 1, Figure 12), are characterized by a high level of clustering and relatively short timescales of recurrence (T ≈ 3 months). Segment 3 in between is characterized by low, almost exclusively clustered activity, with a longer, clearer period (T ≈ 6 months). Interestingly, [START_REF] Ide | Universality of slow earthquakes in the very low frequency band[END_REF] show that the dominant focal mechanism of very-low-frequency earthquakes (VLFEs) are different between each of those region -a segmentation that our model cannot account for. In Figure 12, we see that the along-strike extent of those segments and spots of high tremor activity within them is quite similar to the dimension of seamounts on the seafloor of the Philippine Sea plate, or the ones observed in the accretionary wedge [START_REF] Yamazaki | Subducting seamounts and deformation of overriding forearc wedges around Japan[END_REF]: about a few tens of kilometers. The alignment of the Kinan seamount chain with the segmented tremor zone along the convergence direction could indicate that such structures on the subducting slab, or their lasting effects on the medium as they plow at depth, are the main factor shaping the segmentation of activity in the Shikoku tremor zone.

The segmentation of tremor intermittence in Shikoku can be interpreted as contrasts of activity synchronization in time and space, and therefore of interaction strength between sources. Stronger or weaker source-to-source interactions produce more (segment 3) or less (segments 1 and 2) clustered, periodic activity, with more or less coherence on large scales of time and space. Our work goes beyond this general interpretation by showing that such interactions can occur in a dynamically permeable channel in the fault zone, and that the criticality of sources and their spatial distribution directly control the interaction strength, and therefore shape the intermittence of activity.

The input flux q in partly controls how close valves are to threshold, and therefore how triggerable they are. If valves are relatively similar for the three segments in Shikoku, the lower activity, higher temporal clustering and clearer-defined, longer periodicity in the buffer segment (segment 3, 140 km < d < 190 km along-strike) can be due to a relatively lower metamorphic flux in the region compared to the neighboring regions 1 and 2. The locally lower q in would impose a slower build up of pressure behind the permeability valves, which would spend more time in a closed state, close to their opening threshold. Segment 3 would therefore produce a lower seismicity, and a more clustered, longperiod activity than the end segments 1 and 2, as valves would be closer to threshold.

The same criticality contrast between segments can also arise from a difference in valve strength, assuming the input fluid flow -and transport properties of the channelsdoes not vary significantly along-strike. The valves in the buffer, segment 3, could indeed be stronger -higher breaking threshold and/or lower permeability -but undergoing the same fluid input rate. For a constant input flux along-strike, this results in the same variation of criticality along strike: valves in the buffer region are closer to threshold, producing a more clustered and long-period activity, and valves in the segments 1 and 2 framing it are further from threshold, producing a less clustered, more intense and shorter-period activity. Both these effects can stem from the presence of, or lasting effects of one or several subducting seamounts [START_REF] Ide | Striations, duration, migration and tidal response in deep tremor[END_REF], that strongly fracture the medium in the end segments of the tremor zone. The amount of fluid that comes through the interface around the tremor source region would be enhanced by a channeling effect, and simply by a higher volume of dehydrating crust at depth. In Figure 12, the alignment of the Kinan seamount chain with the subducting seamount detected by [START_REF] Yamazaki | Subducting seamounts and deformation of overriding forearc wedges around Japan[END_REF] and the tremor zone could indicate that the chain extends into the subduction, and could result in the tremor activity patterns observed there. Segment 3 that produces lower, longer-recurrence activity can be a smoother, less damaged region, in which less fluid is channeled because of a lower channeling effect, and a lower overall permeability. In this case, the low flux create a longer recurrence time of periods of opening and activity. Although it might explain the observed activity patterns and it seems like the most direct interpretation of the intermittence of a system of hydraulic pressure accumulation and release [START_REF] Mclellan | Margin-wide variations in slab dehydration in Cascadia and their relationship to slow slip[END_REF], this interpretation might be at odds with observations of tidal and dynamic triggering on tremor in patches of highest activity in the end segments of the zone, and not in the buffer zone [START_REF] Miyazawa | Learning from dynamic triggering of low-frequency tremor in subduction zones[END_REF][START_REF] Chao | A Global Search for Triggered Tremor Following the 2011 Mw 9.0 Tohoku Earthquake[END_REF][START_REF] Chao | Triggered tectonic tremor in various types of fault systems of Japan following the 2012 M w 8.6 Sumatra earthquake: TRIGGERED TREMOR IN JAPAN[END_REF][START_REF] Kurihara | Migration of Deep Low-Frequency Tremor Triggered by Teleseismic Earthquakes in the Southwest Japan Subduction Zone[END_REF]. Indeed, if seismicity in the segments is triggered by minute strains on the interface, it suggests that those segments might actually be closer to criticality: valves would be more sensitive to the very small δp changes across them generated by the dynamic strains of tides and teleseismic waves.

It therefore seems that the segmentation could arise from spatial variations of permeability structure, caused by large-scale heterogeneity in the subducting plate topography. In a first order analysis, we can assume that the input flux in fault channels, the background transport properties and the permeability valve characteristics do not significantly change across the strike of the tremor zone. A contrast of valve density and clustering in each segment would therefore explain the segmentation of intermittence.

In this interpretation, we can explain the characteristic activity in the buffer zone, segment 3, as being due to a higher level of interaction between valves, for instance coming from a higher valve density in this region compared to the neighboring ones. In geological terms, this would mean a lower overall permeability in segment 3, but a larger share of the system being dynamic. This could be due to a more homogeneous permeable system, perhaps because of a relatively smooth segment of oceanic plate being subducted in this region. A dense, highly interacting valve network, would produce a remarkably synchronized and periodic activity in both time and space and lower overall activity rate (Figure 9c and 11), and seems that it well describes the seismicity in segment 3.

In this interpretation, the high valve density in the buffer zone would behave as a subductionscale valve for the pressure circulation along-strike.

Scope and limitations of model geometry

In the previous section, we used our results on space-averaged intermittence of activity in along-dip channels as representative of the activity in wide segments along the strike of the subduction. This approach assumes that the behavior of each segment can be explained by along-dip dynamics of fluid pressure, and that it is isolated from its neighbors, as we neglect the connectivity along-strike, within and between segments. As tremor zones often extend much farther along-strike than along-dip (e.g. Cascadia and Nankai), and as tremor exhibit migrations of activity along-strike, and finally, as hydraulic and solid stresses should propagate in all directions, our approach simplifies the problem greatly.

However, we argue that the dynamics in the segments we defined in Shikoku can reasonably be collapsed in isolated, along-dip channels for our analysis, which focuses on time and space averaged intermittence behavior. Indeed, for the most part, episodes of activity mostly span the whole along-strike extension of each segment, and they only occasionnally communicate.

Eurasian Plate

Philipine Sea Plate 1a is shown with the white box. The three segments defined using the activity are indicated by the numbered circles and dotted lines that divide the tremor zone. We single out rough topographic features on the subducting Philippine Sea plate: the Kinan seamount chain, and a subducting seamount [START_REF] Yamazaki | Subducting seamounts and deformation of overriding forearc wedges around Japan[END_REF]. Because of their alignment with the tremor zone, and their consistent spatial scale, the presence such topographic features could be responsible for the segmentation of activity witnessed in Shikoku, and displayed as a 2D

histogram on the map (cf Figure 1 for details).

In our model, generation of tremors is occurring in parts of the faults where the Further modeling work is needed to understand how the source/sink geometry, the two-dimensionnal heterogeneity and anisotropy of permeability, and permeability-opening processes could conspire to fuel the complex patterns of migrations of tremor along both the strike and dip of the interface.

Interactions between fluid circulation and fault slip

The larger question our work aims to ask is: is fluid pressure and permeability driving the dynamics of tremor activity or is fault slip. This model does not intend to faithfully reproduce spatiotemporal patterns of tremor activity, or fluid pressure variations in the subduction interface. It is however a framework to reflect about the role of hydraulic stress dynamics with more spatial and temporal complexity. By only looking at the dynamics of the hydraulic system, we show that it can reproduce characteristic behavior of tremor and fluid pressure in the fault zone, and that spatial variations of hydraulic properties can lead to variations in tremor intermittence, thereby providing a mechanism that could shape the observed along-strike segmentation of tremor.

It is also clear that shear slip is associated with tremor generation and that the spatial variations of its dynamics play a role in the segmentation of tremor [START_REF] Nakajima | Tremor activity inhibited by welldrained conditions above a megathrust[END_REF][START_REF] Kano | Strength of tremor patches along deep transition zone of a megathrust[END_REF][START_REF] Cattania | Precursory Slow Slip and Foreshocks on Rough Faults[END_REF]. Fault slip and tremor activity are linked because their occurrences are to some degree correlated in space and time (e.g. [START_REF] Rogers | Episodic Tremor and Slip on the Cascadia Subduction Zone: The Chatter of Silent Slip[END_REF][START_REF] Hall | Peak Tremor Rates Lead Peak Slip Rates During Propagation of Two Large Slow Earthquakes in Cascadia[END_REF]. At the same time, this correlation is not perfect and detailed observations that the slow slip and tremors are not perfectly co-located on the same fault segments (e.g. [START_REF] Kostoglodov | The 2006 slow slip event and nonvolcanic tremor in the Mexican subduction zone[END_REF]. Nevertheless, the tremor activity can be used to detect the geodetic signature of small slow slip transients otherwise buried in the noise [START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes: UNCOVERING SILENT SLIP[END_REF]. In this interpretation, the tremor is generated as brittle asperities rupture when a slip transient occurs on the subduction fault [START_REF] Ando | A slip pulse model with fault heterogeneity for low-frequency earthquakes and tremor along plate interfaces[END_REF]. The seismic characteristics of LFEs -their radiation pattern mostly -are also consistent with shear slip on the subduction interface [START_REF] Ide | Mechanism of deep low frequency earthquakes: Further evidence that deep non-volcanic tremor is generated by shear slip on the plate interface[END_REF][START_REF] Royer | A comparative study of low frequency earthquake templates in northern Cascadia[END_REF], although observations interpreted that way are sometimes ambiguous, and can sometimes fit a single-force mechanism as well [START_REF] Wech | Cascadia tremor polarization evidence for plate interface slip[END_REF][START_REF] Shapiro | Low-Frequency Earthquakes and Pore Pressure Transients in Subduction Zones[END_REF][START_REF] Ohmi | Deep low-frequency earthquakes beneath the focal region of the Mw 6.7 2000 Western Tottori earthquake: DLF BE-NEATH THE FOCAL REGION OF THE[END_REF]. The spatio-temporal patterns of tremor could therefore be influenced by the dynamics of slip in conditions of near-lithostatic fluid pressure and heterogeneity of frictional properties of the fault interface (e.g. [START_REF] Wech | A continuum of stress, strength and slip in the Cascadia subduction zone[END_REF][START_REF] Sweet | Variations in Cascadia Low-Frequency Earthquake Behavior With Downdip Distance[END_REF][START_REF] Luo | Rate-and-State Model Casts New Insight into Episodic Tremor and Slow-slip Variability in Cascadia[END_REF].

The processes of unclogging -or hydrofracture of low permeability barriers -and transient fluid pressure that we describe should interact closely with slip on the fault.

As slip occurs on the fault, be it seismically or aseismically, it generates an extensional regime in parts of the fault zone, allowing nearby extension veins to open and pump fluid into the fault [START_REF] Kotowski | Length scales and types of heterogeneities along the deep subduction interface: Insights from exhumed rocks on Syros Island, Greece[END_REF], thus strongly modifying the local permeability and fluid pressure in the slipping region. Damage and dilatancy in and around slip also create permeability [START_REF] Tenthorey | Cohesive strengthening of fault zones during the interseismic period: An experimental study[END_REF][START_REF] Mitchell | Experimental measurements of permeability evolution during triaxial compression of initially intact crystalline rocks and implications for fluid flow in fault zones[END_REF][START_REF] Im | Cyclic Permeability Evolution During Repose Then Reactivation of Fractures and Faults[END_REF], and can therefore also affect the fluid pressure field.

On the other hand, it is known that high fluid pressure can trigger slip on a fault by lowering the effective fault strength. The most plausible cause of an increase of fluid pressure is an input of fluid in a partly sealed region. In a homogeneous channel, such fluid flow is necessarily associated with a smooth pressure gradient, implying that the high fluid pressures that are required for sliding motion are limited to the neighborhood of the source. Our study shows how heterogeneous and transient permeability behavior in the fault could generate high fluid pressures locally and transiently, in wide portions of the subduction zone. The collective behavior of valves creates a cumulative effect on fluid pressure that translates into multi-scale fluid pressure increases and drops, that could shape the behavior of fault slip from the source to the subduction scale, thus shaping tremor activity, as hinted by the geological record [START_REF] Angiboust | Probing the transition between seismically coupled and decoupled segments along an ancient subduction interface[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][START_REF] Muñoz-Montecinos | Episodic hydrofracturing and large-scale flushing along deep subduction interfaces: Implications for fluid transfer and carbon recycling (Zagros Orogen, southeastern Iran)[END_REF][START_REF] Tarling | A common type of mineralogical banding in serpentine crack-seal veins[END_REF].

We believe the next step in our approach is to understand what scales and observable behaviors of tremor activity, fluid pressure or slip are specific to a fluid-dominated or a deformation-dominated regime, and therefore can help distinguish between them in the field. To do so, experimental, theoretical and modeling work is needed to better describe the coupling between deformation (induced by pressure and/or slip) and permeability, and what physical parameters control which regime the system is in.

Conclusion

The present work should be understood as a simple framework emphasizing the role of fluid pressure transients in the occurrence of tectonic tremor in subduction zones. Our study is based on the premise that clustering and quasi-periodicity of seismic activity should occur when sources synchronize, through interactions. As they synchronize, the activity of numerous sources becomes coherent over large scales of space and time. In order to measure as explicitly as possible the synchronization of sources, we build novel and simple measures of temporal clustering and periodicity in a point-process description of tremor activity. We propose that the interactions at the origin of the observed intermittency of tremor may be mediated by fluid pressure transients in the permeable fault zone. In our description, elementary tremor events (LFEs) occur when permeability valves open within the fault zone, stressed by the incoming metamorphic fluid flux channeled in the fault zone. We find that how close valves are to their activation threshold -resulting from a competition between the intensity of the fluid pressure source and their mechanical strength -and their spatial distribution in the fault zone both control the intermittency of their collective behavior.

In our framework, the highly-clustered, long-period activity of the buffer segment in the Shikoku tremor zone develops when densely-packed, highly-interacting valves activate collectively, thus building subduction-scale valving behavior accompanied by large bursts of activity. The large timescales of activity, and remarkable spatial coherence of the tremor bursts thus emerge from small-scale fluid pressure and permeability cycling in the dynamic fault zone, building up through interactions between dynamic segments.

The segmentation of activity in the tremor zone seems directly linked with the topography of the subducting Philippine Sea plate. While the rough terrains in the extremity of the tremor zone should favor less critical valves and a rougher, sparser distribution of valves, the smoothness of the slab surface in the buffer segment produces a dynamics adequately represented by a dense valve system. Our model therefore provides a simple, physically-based mechanism to describe the influence of hydrological processes on shaping tremor patterns in both space and time. More work is needed to refine the physical description of hydromechanical processes in this framework, in particular the interplay between the slip, the fluid pressure and the permeability. However, we argue that including dynamic hydrological processes is crucial to understand tremor and microseismicity patterns in general. Beyond the subduction zone setting, we therefore expect the permeability-valve framework has a broader reach to interpret how unsteady fluid circulation processes shape the dynamics of a wide range of geologic plumbing systems [START_REF] Journeau | Seismic tremor reveals active transcrustal magmatic system beneath Kamchatka volcanoes[END_REF][START_REF] Wech | Deep long-period earthquakes generated by second boiling beneath Mauna Kea volcano[END_REF][START_REF] Ross | 3D fault architecture controls the dynamism of earthquake swarms[END_REF][START_REF] Gosselin | Seismic evidence for megathrust fault-valve behavior during episodic tremor and slip[END_REF][START_REF] Materna | Dynamically Triggered Changes of Plate Interface Coupling in Southern Cascadia[END_REF].

Open Research Section

The catalogs of tremor and LFE used in this study are available at https://doi .org/10.31905/UOQ9LVHZ for Shikoku, Japan (Poiata, Vilotte, Shapiro, Supino, & Obara, 2021;Poiata, Vilotte, Shapiro, Obara, & Supino, 2021), in the supporting information of [START_REF] Husker | Characteristic Tectonic Tremor Activity Observed Over Multiple Slow Slip Cycles in the Mexican Subduction Zone[END_REF] for Guerrero and Oaxaca (Mexico), and at https://pnsn.org/ tremor for Cascadia [START_REF] Wech | Cataloging Tectonic Tremor Energy Radiation in the Cascadia Subduction Zone[END_REF]. The code -still in active development -for the model is available at https://github.com/gfarge/PPvalves.

  a year. It is the most temporally clustered activity of the three regions. In other subduction zones, similar variations of temporal clustering and recurrence timescales can be observed, showing a clear along-strike segmentation of tremor activity. In the supplementary materials, we show time series of tremor activity rate in the Central American subduction zone (Figures S1, tremor catalog from Husker et al. (2019)) and along the Cascadia subduction zone (Figure S2, tremor catalog from Wech (2021)). In Shikoku, Guerrero, and Cascadia, the segments are of geological scale, spanning tens to hundreds of kilometers along the strike of the subduction.

Figure 1 :

 1 Figure 1: Caption on next page.

Figure 1 .

 1 Figure 1. On previous page. Patterns of temporal clustering and recurrence of LFE activity,

  (b) Event count binned in 3-km wide bins along-strike, for the full period. (c) Intensity of temporal clustering along-strike, measured by the clustering index c in bins of 20 km, with a 10 km overlap. (d) Recurrence timescale of clusters along strike, with the same spatial bins as in (c). Details on the cluster detection and clustering index can be found in Section 2.3. The blue, shaded area around the measurement represents 1-σ interval around the estimated recurrence timescale. A visually estimated recurrence is shown as a dotted blue line. (e) LFE activity along-strike in time, each black dot represents an event. Three major segments are highlighted in blue (zone 1), yellow (zone 2) and red (zone 3). See text for details, and Figure 12 for wider geographical context.

:

  Figure 2. (a) Conceptual representation of fluid circulation in the fault zone. Downdip and below the tremor source region, metamorphic dehydration reactions release fluid, which is channeled in the permeable subduction interface, circulating at high fluid pressures. (b) Schematic representation of the model set up: a 1D permeable channel in which identical valves are distributed (black bars), opening and closing in response to the fluid pressure field. The channel is fed at its base by input fluid flux qin, which is the main control parameter of the experiment. At its top, the channel is connected to a high permeability region such that the fluid pressure is close to the hydrostatic pressure pout = 0 (see text for details on the scaling of physical properties). (c)Each valve opening is associated with a seismic event, a model LFE. The activity is represented in a space-time graph, with time on the x-axis, and distance along the channel on the y-axis.

Figure

  Figure reproduced and modified from Farge et al. (2021).

  old δp break c , they instantly change to an open state, and their permeability switches to k hi . When the accumulated overpressure across a now-open valve diffuses and δp drops below a lower threshold δp clog c , the valve then closes instantly, and the permeability in the segment goes back to k lo . We do not model in details the dynamics either of the opening or closing of permeability. Although simply parameterized, they are based on a consistent modeling of how porous media clogs and unclogs in response to flux, pressure gra-

c,

  δp across a closed valve will eventually go over δp break c , and the valve will open. As the valve opens, the steep pressure gradients diffuses and δp drops across the valve. For the valve to close back, the strength of the flux has to allow δp to drop below the closing threshold δp clog c

)

  In the previous equations, δp break c and δp clog c are defined as the difference of pressure from down-to updip of the valve, hence a different sign convention than equation 2. It is intuitive that δp break c > δp clog c . However, we show in Farge et al. (2021) (Section 4.1) that because both tortuosity and pore aperture change when pores open, the change of permeability is more important than the ratio of the δp threshold, which only depends on pore aperture, therefore k hi /k lo > δp break c

  , the flux is just above q break c . After a rapid increase of δp when the valve closes, the closed valve is close to an equilibrium and it lets fluid seep through it as the flux loads it, at a barely superior rate. δp will relatively slowly reach the breaking threshold δp break c . When the valve eventually opens, the background flux is so low that δp quickly drops to the closing point δp clog c . Overall, the valve spends most of the time closed, close to opening conditions for q in just above q break c . Conversely, if q in is higher, just below q clog c , the valve takes a much longer time to close when open, and opens very quickly when closed: it spends more time open overall. This behavior is reflected in the time averaged properties of the valve cycle as a function of flux, displayed in Figure3g and h. When close to threshold, the valve rarely switches states, and produces few events -openings. When q in is far from the critical values, the valve rapidly loads and relaxes. Because the cycle is faster, the activity rate is higher. Finally, the proportion of the time the valve is open smoothly transitions from a mostly closed state to a mostly open state from q break c to q clog c , matching the static regime beyond those values. This progressive permeability evolution from a clogged to an unclogged state for increasing values of the hydraulic stressing rate has been observed in injection experiments in the lab[START_REF] Candela | Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing[END_REF]. Two other parameters shape the valve cycle: valve width w v and the ratio of closed to open permeability λ = k hi /k lo . Valve width governs how quickly fluid can seep through the valve, and how strong the fluid flux near a valve has to be to impose a given δp across

Figure 3 .

 3 Figure 3. Temporal behavior of a single valve, subjected to different values of the input flux qin, within the range allowing sustained valve activity, q break c

  4b) eventually synchronize, with a short, consistent delay between the first valve's opening (downdip valve) and the second valve's opening (updip valve). The valves open and close synchronously, and events associated with each opening therefore occur in a two-event burst, recurring every valve cycle.As we zoom in on the last opening sequence of the simulation (Figure4c), we observe that the updip valve is affected with a short delay by the downdip valve's opening and the flux pulse that it produces. δp across the closed updip valve suddenly rises due to the fluid influx behind it, and as it is already close to failure, it is sufficient to bring it to open. It should be noted that in[START_REF] Farge | Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface[END_REF] Section 4.2, we have shown that this triggering interaction does not necessarily occur from the downdip valve towards the updip valve. Indeed, the effect of a valve opening updip creates a transiently low pressure in front of the downdip valve, that can trigger its opening.

Figure 4 .

 4 Figure 4. Interactions in a two-valve system. (a) A system with two valves (v1 and v2) spaced by dv = 0.008 = 0.8wv, subjected to a constant flux qin = 0.13, and a constant pressure output pout = 0. (b) Cycles of fluid pressure accumulation and release are shown with the timeseries of the pressure difference δp across v1 and v2. Closing and opening occur when δp reaches a pre-defined threshold. (c) As v1 opens first and when v2 is critically stressed, the pressure transient that its opening produces (in red) triggers the opening of v2, after a short delay.

Figure 5 .

 5 Figure5. Effect of valve spacing on valve-to-valve interaction. (a) Two timescales describe a two-event burst: δtintra is the delay between events in the burst, and δtinter is the recurrence delay of the burst. In (b) and (c), we observe how δtinter (no markers on the line) and δtintra (circle markers) evolve as a function of valve spacing for low (red, dotted line), intermediate (yellow, dash-dot line) and high (green, dashed line) flux values. Both timescales are normalized by δt0, the period of activation of a valve subjected to the same flux. For sufficiently large valve spacing (dv > 5wv), valves do not interact and act as isolated valves. As valves get closer, events get closer, and the bursts are less frequent.

578

  closer to critical opening and closing thresholds, the deviation of δt inter from the refer-579 ence cycle period (grey line) happens for valve spacings that are wider than for the medium 580 flux range (yellow curve). Indeed, when subjected to a medium flux range, valves are fur-581 ther away from criticality most of the time, and they have to be closer neighbors to in-582 teract with a similar intensity as when they are more critical. The same observation can 583

tion 7 .

 7 Indeed, when u → +∞, the distribution tends toward a Dirac in d 0 -valves are regularly spaced, at a distance of d v = d 0 . When u → 0, the Weibull distribution tends to a power law distribution p(d v ) ∝ d -1 v . In this case, most valves are spaced with very small d v , and rarely, a large spacing is drawn, producing a clustered valve distribution. Finally, for u = 1, the Weibull distribution reduces to an exponential distribution p(d v ) = exp(-d v /d 0 )/d 0 , of mean d v = d 0 . In this case, the number of valves in a segment of given length follows a Poisson distribution, and valves are more or less homogeneously distributed across the channel, although at random. As a rule of thumb, this distribution produces clustered point distributions for u < 1, and more homogeneous and regular valve distributions for u ≥ 1. Finally, the choice of d 0 will allow to specify the average valve distance d v , and therefore number of valves N v in the system (equation 8).

  d

Figure 6 .

 6 Figure 6. Random realizations of valve distributions for different number of valves Nv and spatial clustering. Valve distributions are generated by drawing valve distances dv from a Weibull distribution. The valve density, or average valve distance in the domain dv, and therefore the number of valves in the system, are controlled by the theoretical mean of the distribution. The shape parameter u of the Weibull distribution allows to control how clustered or regularly-spaced valves are in space.

644

  Each system thus defined (3840 total) is subjected to different values of input flux q in = 645 0.11-1.04 (11 values), with denser sampling close to the thresholds allowing valve open-646 ing (q break c = 0.1) and closing (q clog c = 1.05). In all 42,240 simulations, the output pres-647

Figure 7 .

 7 Figure 7. Synthetic tremor activity for three values of the input flux qin (a, b, c) in the same, clustered valve system (u = 0.3, Nv = 29, dv = 2wv, system b in Figure 9). Panels (a.i), (b.i) and (c.i) represent the valve distributions along the channel. Panels (a.ii), (b.ii) and (c.ii) display the activity along the channel in time, in a time-space diagram where each dot represents the location and time of an opening event. Panels (a.iii), (b.iii) and (c.iii) show the event count (per bin of δt = 2.5e -4) time series. The colored panels represent a flux scale, on which qin is represented relative to the opening and closing thresholds. The closer to threshold qin is, the more valves are sensitive to interactions, as they spend more times in a given state. This produces a more synchronized, clustered activity, recurring on long periods. The reader should note the dilated time scale in simulation (b).

-Figure 8

 8 Figure 8 Caption on next page.

Figure 8 .

 8 Figure 8. Activity rate, temporal clustering and recurrence timescales for different values of the input flux qin, for all valve distributions. (a) Activity rate, computed as the number of events per valve per δt = 2.5e -4 (scaled time units). (b) Level of activity clustering measured with the clustering index c, for different values of qin. (c) The average recurrence delay for clusters, as a function of qin. The grey dotted line represents the proportion of periodic simulations (see text for explanation) for each value of qin. In panels (a), (b) and (c), the distribution of measured values for each value of qin is visualized with a violin plot, and the median of all measurements for this qin is represented as a cross. The domains where the flux value prohibits valve activity are shown with red (all valves closed) and green (all valves closed) patches. The colored dots show the clustering and recurrence for the activity in a given system, subjected to three values the input flux. Each corresponding time series of activity is displayed in the bottom panels (d), (e) and (f ). The three simulations are performed on a clustered valve distribution (u = 0.3) with Nv = 29 valves (dv = 2wv), system b displayed in Figures 9b and 11e.

Figure 9 .

 9 Figure 9. Synthetic tremor activity for three valve systems with different density and clustering (a, b, c). Panels (a.i), (b.i) and (c.i) represent the valve distributions along the channel.

  summarizes the results. System (a) produces temporally-clustered activity, occurring in bursts, with however less variability than the other two systems. Indeed, System (b)(Figure 9b.i, b.ii, b.iii) shares the same number of valves and average valve spacing as System (a) however the higher spatial clustering of valves creates dense patches of valves in the system, where valves are locally much closer together than the average inter-valve distance d v .The resulting activity is more synchronized, exhibiting clearly separated bursts -of two sizes, in two different regions of the domain. The burst recurrence is quite variable, but is obviously longer than for bursts produced in System (a). Finally, System (c) (Figure9c.i, c.ii, c.iii) is as spatially clustered as System (b), but with more valves. It is the system where on average valves are the closest, and also with the most places where valves are very close locally. The activity it produces is the most clustered of the three: it proceeds in bursts of two sizes, almost without any activity otherwise. They occur almost periodically, with a more constant, longer period than System (b). Strikingly, the overall activity rate in System (c) is lower than in the other two: as in the two-valve system (Section 4.2), the longer-recurrence intervals of activity episodes and stronger clustering act to decrease the overall activity rate of valves in the system.The level of synchronicity of activity increases both as spatial clustering u increases ((a) to (b)) and as valve density increases ((b) to (c))

Figure 10 .

 10 Figure 10. The number of interacting valves N int v is a measure of how much valves are closefrom each other in a valve system, either owing to the density of the system (Nv, dv), or to the spatial clustering (u). As valves get closer, they interact more, and N int v seems to be the best

Figure 11 .

 11 Figure 11. Intensity, clustering and recurrence timescales of activity for varying levels of interaction in complex valve systems. Results for qin = 0.13. (a) Average event count per valve according to the number of interacting valves N int v . (b) Level of activity clustering measured with the clustering index c, as a function of N int v . (c) Timescales of cluster recurrence as a function of N int v . The detected timescales are represented as dots which size scales with the variability of the measured recurrence, estimated as the ratio of the average to the standard deviation of the inter-cluster delays, in the given simulation. The transparency of each data dot is used to convey the density of points in the graph. The dotted line represents the proportion of periodic simulations (see Section 2.3 for details) for each value of N int v . In panels (a), (b) and (c), the circled, colored dots show the values of clustering and recurrence for three valve systems with increasingly strong valve interactions (same qin = 0.13). The three bottom panels (d), (e) and (f ) show the corresponding valve distributions (top) and activity (bottom).

Figure 12 .

 12 Figure 12. Regional context of the Shikoku subduction and tremor zone. The trench is outlined in thick black line, and the depth contours of the slab are shown in thin dotted black line (model from Iwasaki et al. (2015)). The black arrow and wide dotted line show the direction of relative convergence of the Philippine Sea plate, towards the Eurasian plate (fixed). The extent of the map in Figure1ais shown with the white box. The three segments defined using the activ-

  fluid pressure is very close to the lithostatic pressure. Therefore, large-scale fluid pressure gradient and fluid flow in the tremor generating zone is following the lithostatic gradient -from deep metamorphic sources towards shallow sinks. As a consequence the fluid pressure gradient and associated fluid flux are dominantly oriented updip with horizontal components being on average much weaker. Therefore the tremor -assuming it's associated with permeability opening -is mostly associated with upward fluid flow accelerations along-dip. Along-strike migrations of activity thus likely track a pulse of acceleration of updip fluid flow, propagating along the strike of the fault. Such horizontally propagating pulses can be associated with transient local along-strike fluid pressure gradients emerging after valve breaking. This 2D situation can be schematically represented with a series of 1D along-dip channel filled with valves and governed by updip pressure gradients. Breaking one valve triggers rapid local pressure variations inside the respective channel. These variations are then transmitted to a neighbor channel through the along-strike hydraulic connectivity and eventually trigger closely located valves in this new channel. A cascade of such perturbations can propagate near-horizontally and at the same time overall fluid flow is still dominantly oriented upwards along-dip.

Table 1 .

 1 Valve parameters for all simulations presented in this study

		Parameter	Value (scaled)
	w v	Valve width	0.01
	k hi	Open valve permeability	1
	k lo	Closed valve permeability	0.05
	δp

  describes the density-

	636
	clustering (N v -u) space of distributions available to us using this technique. Figures S7
	and Figures S8 in the Supplementary Information file give a more complete illustration

638

of this distribution space.

639

Panels (a.ii), (b.ii) and (c.ii) display the activity along the channel in time, in a time-space diagram where each dot represents the location and time of an opening event. Panels (a.iii), (b.iii) and (c.iii) show event count (per bin of δt = 2.5e -4) time series. The denser and more clustered the valve system, the more clustered and periodic the activity, and the longer the timescales of recurrence. In other words, the closer the valves get in the system, the more synchronized and long-period the resulting activity seems to become. -24-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). G. Farge gratefully acknowledges financial support for his research at MIT by the Fulbright U.S. Doctoral Program, which is sponsored by the U.S. Department of State and the French-American Fulbright Commission. The contents of this publication are solely the responsibility of the author and do not necessarily represent the official views of the Fulbright Program, the Government of the United States, or the French-American Fulbright Commission. Numerical computations were performed on the S-CAPAD/DANTE platform, IPGP, France. The authors thank S. Rondenay, E. Brodsky, C. Cattania, A.

Schubnel, K. Chanard, two anonymous reviewers and the editor for constructive discussion that led to significant improvement of this manuscript.