### 3.3.6 Numerical implementation

Solutions for the pressure field in time were obtained using a Crank-Nicholson finite difference scheme, with a first order, centered discretization of the second order space derivatives of pressure (Press et al., 2007, section 20.2, pp 1045–1048). The scheme is implicit, centered and second order in time, and is unconditionally stable. Valves are either in an open or a closed state, as described in section 3.3.3 above, and when the pressure differential at their boundary reaches a threshold value, their open/closed state is changed accordingly. Each opening of a valve is counted as a seismic event.

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

## 3.4 Time-dependent valve behaviour

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

### 3.4.1 A single valve

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

$$q_{in} = \frac{\rho}{\eta} k_v \frac{\delta p}{w} \qquad (3.13)$$

where $k_v$, the valve permeability, may be equal to $k_{lo}$ or $k_{hi}$, and where $δp$ is again the pressure difference across the valve. If $δp \le δp_c^{break}$ in a closed valve, stuck particles cannot be removed and the valve remains closed. Conversely, if $δp \ge δp_c^{clog}$ in an open valve, particles cannot be captured and the valve
remains open. These conditions can be expressed as follows, using $q_{in}$, the flux a valve is subjected to:

$$q_{in} \leq \frac{\rho}{\eta} k_{lo} \frac{\delta p_{c}^{break}}{w} (= q_{c}^{break}) \text{ for a closed valve,} \quad (3.14)$$

$$q_{in} \geq \frac{\rho}{\eta} k_{hi} \frac{\delta p_{c}^{clog}}{w} (= q_{c}^{clog}) \text{ for an open valve.} \quad (3.15)$$

It introduces two threshold values for the fluid flux, noted $q_c^{break}$ and $q_c^{clog}$, respectively. The two conditions in (3.14) and (3.15) define domains in $(\delta p, k_v)$ space corresponding to a closed valve or an open one in steady-state (Figure 3.4). Table 3.2 lists values for the pressure thresholds and permeabilities, as well as for the flux thresholds used for the valve set up of this study.

Depending on the four valve parameters, there are two different types of solutions. Let us first consider enhanced permeability variations, such that $k_{hi}/k_{lo} > \delta p_c^{break}/\delta p_c^{clog}$. In this case, $q_c^{break} < q_c^{clog}$ and there are no steady-state solutions in the $q_c^{break} - q_c^{clog}$ interval (Figure 3.4a). In this interval, the input flux being higher than $q_c^{break}$, a closed valve will open, and in turn, the input flux being lower than $q_c^{clog}$, the now open valve will eventually close. Therefore, the valve must constantly switch between open and closed states in order to handle the imposed input flux. This oscillatory regime is illustrated schematically in Figure 3.3.

In the other case, called enhanced pressure variations, $k_{hi}/k_{lo} < \delta p_c^{break}/\delta p_c^{clog}$. This implies that $q_c^{break} > q_c^{clog}$ and the domains for a closed valve and an open one overlap in the $q_c^{clog} - q_c^{break}$ interval. In this interval, steady-state solutions are possible for the two valve states (Figure 3.4b). Which one is selected depends on the initial state of the valve. In this case, therefore, there is no valve activity.

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

$$q_c^{break} < q_c^{clog} \iff \frac{\delta p_c^{break}}{\delta p_c^{clog}} < \frac{k_{hi}}{k_{lo}}. \quad (3.16)$$

In a geological context, the tortuosity ratio is small and lies in a restricted range (see section 3.3.3, last paragraph) and in contrast, permeability is expected to change by one order of magnitude (Beckingham, 2017; McDowell-Boyer et al., 1986; Miller & Nur, 2000). Using the argument developed for equation (3.12), where we derive the threshold $\delta p$ ratio as a function of the permeability ratio, we conclude that the instability criterion should be met in geological porous media where erosion/deposition processes are active. Such permeable systems have been reported to self-organize around these critical erosion/deposition conditions (Kudrolli & Clotet, 2016).
<figure><img src="image_42.png" /><figcaption>**Figure 3.4** — Requirements for steady-state valve regimes (equations (3.14) and (3.15)) in $k-\delta p$ space (in log-log scale). The two values for permeability $k$ and threshold values for the valve pressure difference $\delta p$ are shown with vertical and horizontal dashed lines. In this space, iso-flux contours appear in gray lines of slope -1 in log-log scale (-45°). Blue and red segments indicate values of $\delta p$ that are allowed for a closed valve and an open one in steady-state, respectively. Steady-state conditions are achieved in the green and red shaded areas. See explanations in the text. (a) Valve with enhanced permeability variations ($k_{hi}/k_{lo} > \delta p_c^{break}/\delta p_c^{clog}$). (b) Valve with enhanced pressure variations ($k_{hi}/k_{lo} < \delta p_c^{break}/\delta p_c^{clog}$).</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Figure 3.5** — Influence of the input flux $q_{in}$ on the activity of a single valve. (a) Three intervals of input flux $q_{in}$ are delimited by two threshold values $q_c^{clog}$ and $q_c^{break}$. (b) Examples of activity rate evolution in time (in black) for each regime. Activity rate is taken as the number of events per $\delta t = 0.01$, scaled unit. An event corresponds to the opening of a valve. Valve parameters are listed in Table 3.2.</figcaption></figure>
Figure 3.5 shows results for three values of input flux, corresponding to steady-state with a closed valve ($q_{red} < q_c^{break}$), a permanent intermittent regime with an unstable valve that switches between closed and open states ($q_c^{break} < q_{yellow} < q_c^{clog}$), and finally steady-state with an open valve ($q_{green} > q_c^{clog}$). In the first case, $q_{red}$ is too low to allow the valve to open at all, no activity is recorded. In the last case, $q_{green}$ is so high that after an initial transient during which the valve opens and closes to accomodate the incoming flux, the valve ends up open, submitted to a flux that does not allow it to close. In the intermediate case, $q_{yellow}$ is both high and low enough to allow the valve to open when closed, and to close when open. After an initial transient, the valve activates in a very regular fashion.

### 3.4.2 Valve-valve interaction via pore pressure diffusion

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

<figure><img src="image_44.png" /><figcaption>**Figure 3.6** — Interactions between three neighboring valves. (a) Transient evolution of pore pressure from an initial equilibrium profile with all valves closed towards a profile with valves v1 and v3 closed and valve v2 open. Panels (b) and (c): pore pressure difference evolution across valves v1 and v3.</figcaption></figure>
## 3.5 The collective behavior of interacting valves

We now explore systematically a model involving 29 randomly distributed valves (Figure 3.2b). With such a valve density, there are patches with valves that are closer together than on average, thus strengthening interactions that may generate spatial patterns of activity. All the valves are identical. We take the same permeability and threshold $\delta p$ as specified in Table 3.2, implying that the threshold flux $q_c^{clog}$ and $q_c^{break}$ are the same as in previous sections (equations (3.14) and (3.15)). In fault-zones, we expect that conditions for enhanced permeability variations are met (see section 3.4.1, equation (3.16) in particular) and hence focus on the case in which $q_c^{break} < q_c^{clog}$.

### 3.5.1 Characteristics of valve activity in the permanent regime

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

$$k_{eq} = \frac{L}{(L - N_{cl}w)/k_{hi} + N_{cl}w/k_{lo}} \quad (3.17)$$

where $N_{cl}$ is the number of closed valves, and $L$ is the channel length. For any of the three regimes, activity rate, $k_{eq}$ and the pressure difference across the channel $\Delta p$ reach a dynamic equilibrium around a stable value after an initial transient, thus defining a permanent regime (Figures 3.7 and 3.8).

The analysis of section 3.4.1 for a single valve does not depend on the presence of neighboring valves, and holds for a channel with multiple valves. The three permanent regimes for the 29-valve system are illustrated in Figure 3.9: one with intermittent valve activity and two quiescent ones. Valve activity occurs when the flux is within the $q_c^{break}-q_c^{clog}$ range. Within that range, we find that the activity style changes as a function of the input flux (Figure 3.10). Near the two ends of that range, activity is episodic or time-clustered (see also Figures 3.7 and 3.8), and almost periodic for the highest values of $q_{in}$ (as shown later on, Figures 3.11 and 3.12). For intermediate flux values in the middle of the flux range, activity is random and Poisson-like, with a seemingly constant overall rate. This result holds for the valve distribution displayed in Figure 3.2b, and for 30 other random distributions of 29-valves with the same characteristics (parameters in Table 3.2). This will be discussed in more detail below.

Here, the catalog of synthetic events in space and time is analyzed in the
<figure><img src="image_45.png" /><figcaption>**Figure 3.7** — Time evolution of activity for a channel with 29 randomly distributed valves and a low input flux $q_{in} = 0.16$, just above $q_c^{break}$. (a) Fluid pressure differential across the channel $\Delta p$, (b) equivalent permeability $k_{eq}$ and (c) activity rate (events per $\delta t = 4 \times 10^{-4}$ scaled unit). Raw data in grey and smoothed time series in colors (rolling mean with a window of length $\tau = 0.03$). (d) Time-space diagram of events occurrence within the channel. All valves are initially closed. A permanent regime with intermittent activity is reached after a short initial transient. Valve distribution represented in Figure 3.2, valve parameters are listed in Table 3.2.</figcaption></figure>
<figure><img src="image_46.png" /><figcaption>**Figure 3.8** — Time evolution of activity for a channel with 29 randomly distributed valves and a high input flux $q_{in} = 0.81$, just below $q_c^{clog}$. (a) Fluid pressure differential across the channel $\Delta p$, (b) equivalent permeability $k_{eq}$ and (c) activity rate (events per $\delta t = 4 \times 10^{-4}$ scaled unit). Raw data in grey and smoothed time series in colors (rolling mean with a window of length $\tau = 0.03$). (d) Time-space diagram of events occurrence within the channel. All valves are initially closed. A permanent regime with intermittent activity is reached after a short initial transient. Valve distribution represented in Figure 3.2, valve parameters are listed in Table 3.2.</figcaption></figure>
form of a two-dimensional point process ($t_i, x_i$) illustrated in time-location
graphs (Figures 3.7d, 3.8d, 3.11 and 3.12c). To characterize the level of time
clustering (episodicity) of synthetic events, we follow Frank et al. (2016) and
use statistical tools for the analysis of temporal point processes (Lowen & Te-
ich, 2005). The main idea is to compute the spectrum of the event count signal
autocorrelations, and to estimate the power-law exponent $\gamma$ of its high fre-
quency decay. A flat spectrum with $\gamma \approx 0$ corresponds to a nearly Poissonian
process, whereas $\gamma$ values that are significantly larger than 0 are indicative of
time clustering.

As shown in Figure 3.10a, activity is Poisson-like for intermediate values
of the input flux, such that $\gamma \approx 0$. For values of the input flux at both ends of
the $q_c^{break} - q_c^{clog}$ range, activity is characterized by high values of $\gamma$, indicating
time clustering. Time-series of the activity rate reveal episodic bursts (Fig-
ure 3.10b). When the input flux is low (close to $q_c^{break}$), activity appears to be
close to scale-invariant clustering with a burst recurrence rate that depends
on the time scale of the measurement. For higher values of the input flux
close to the upper bound $q_c^{clog}$, activity proceeds in quasi periodic bursts, with
periods that seem to increase with increasing $q_{in}$.

As the input flux is increased from the lower threshold value $q_c^{break}$ to the
upper one $q_c^{clog}$, the channel is increasingly open on average. This may be
assessed from the average proportion of open valves which increases as a
function of $q_{in}$ (Figure 3.10c). The proportion of open valves dictates the evo-
lution of $k_{eq}$, the bulk permeability of the channel as a whole. For a low input
flux ($q_{in} \approx q_c^{break}$), valves are closed most of the time and activity bursts are
linked to phases of opening. In the other extreme, for $q_{in} \approx q_c^{clog}$, permeabil-
ity is open most of the time. Swarms of events then occur mainly when the
channel shuts down temporarily and valve pressures increase until unclog-
ging conditions are met.

**3.5.2 Cascades and migrations of synthetic activity**

In the random distribution of valves that was adopted here, there are three
patches where valves are closer to one another than elsewhere, around $x =$
0.1, $x =$ 0.5 and $x =$ 0.85 (Figure 3.11). In patches of closely packed valves,
the constructive interaction we describe in Figure 3.6 is responsible for the
occurrence of rapid cascades of events. These migrations may proceed in both
downdip and updip directions, apparently with similar propagation rates.
Obviously, the interaction is weaker and slower for valves that are separated
by larger distances. Thus, the rates of cascading and event migration depend
on the valve distribution.

Figure 3.11 depicts typical synthetic activity patterns. The input flux is
<figure><img src="image_47.png" /><figcaption>**Figure 3.9** — Influence of the input flux $q_{in}$ in a representative system with 29, randomly distributed valves. (a) Three intervals of input flux $q_{in}$ are de-limited by two threshold values $q_c^{clog}$ and $q_c^{break}$. (b) Examples of activity rate (number of events per $\delta t = 5 \times 10^{-4}$, scaled unit, in black) and proportion of open valves (in cyan) as a function of time, starting at $t = 0$ with all valves closed. Valve distribution represented in Figure 3.2, valve parameters are listed in Table 3.2.</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>**Figure 3.10** — Different permanent activity styles as a function of the input fluid flux, averaged for 30 different 29-valve systems. (a) In white, the level of time clustering $\gamma$: the spectral slope of autocorrelation of event counts (see text) ; in gray, mean activity rate: number of events per unit time, both as a function of input flux $q_{in}$. For $\gamma$ close to 0, activity is Poisson-like with no time clustering. The higher $\gamma$ is, the more clustered and episodic activity is. (b) Activity rate (in black) and proportion of open valves (in cyan) as a function of time for particular simulations within the different regimes. (c) Time-averaged proportion of open valves in simulations, as a function of input flux. In (a) and (c), each dot represent the average value of the variable over 30 different 29-valve systems, at the given flux. Error bars are present in both panels, representing the standard deviation of this measurement over the 30 systems.</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>**Figure 3.11** — Synthetic tremor activity patterns. (a) Distribution of valves in the permeable channel. This distribution was generated randomly and contains three large clusters. (b) Activity patterns in space and time. Black dots stand for valve opening events. The lower valve cluster ($x = 0 - 0.18$) is active almost continuously. The middle and upper clusters ($x = 0.2 - 0.65$ and $x = 0.7 - 1$) is activated sequentially in bursts. Within a burst, activity is intermittent too. (c) focus on the activity in the middle valve cluster showing fast, short-scale migrations in both the updip and downdip directions. Valve distribution represented in Figure 3.2, valve parameters are listed in Table 3.2.</figcaption></figure>

close to the upper threshold $q_c^{clog}$, $q_{in} = 0.93$. Figure 3.11c shows nearly deterministic patterns of collective behavior at both short time and spatial scales. Cascades of triggering/opening events are associated with migrations of activity that show up as oblique dot alignments in the time-location chart. For this run as well as others, the migration rate of cascades takes similar dimensionless values $v \approx 4 - 6 \times 10^2$, with faster velocities (up to $v \approx 2 \times 10^3$ in dimensionless units) in the densest valve patches.

The fastest migrations appear to start and stop at the edges of densely populated patches. In such a patch, valves are close to one another and act in tandem. This shows up as an activity cluster which hosts its own internal migration patterns. Larger scale valve patches are activated sequentially in the updip direction. These large scale migrations proceed at slower rates than those within a single valve patch (Figure 3.11b). These large-scale migrations propagate at $V \approx 22 \times 10^2$ in dimensionless values from one end of the channel to the other.
<figure><img src="image_50.png" /><figcaption>**Figure 3.12** — Cyclic variations of activity for $q_{in} = 0.93$. (a) Pore pressure difference across the whole channel $\Delta p$, (b) channel bulk permeability $k_{eq}$, and (c) time-location chart of events. (d) Activity in $(\Delta p, k_{eq})$ space, showing hysteretical cycles like those of a single valve. Valve distribution represented in Figure 3.2, valve parameters are listed in Table 3.2.</figcaption></figure>

### 3.5.3 Large-scale valving behavior

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

### 3.5.4 Random valve distributions

In order to evaluate the impact of the valve distribution on the previous results, we have generated 30 other random distributions of the same 29 valves
and kept all other parameters at the same values. In all these distributions, the inter-valve distance varies randomly but there are always a number of patches where valves are closer to one another than on average. By construction, the two threshold values of the fluid flux are the same. For all these distributions, the style of activity depends on the fluid flux in exactly the same manner as for the example valve distribution represented in Figure 3.2b, which activity is presented in Figures 3.7, 3.8, 3.9, 3.10b, 3.11 and 3.12. For the 30 valve distributions tested, Figure 3.10a and c show the averaged values of the clustering exponent $\gamma$, the mean activity rate and the time-averaged proportion of open valves at each value of the input flux. We observe similar spatio-temporal patterns within individual patches and at the scale of the inter-patch distance as shown in section 3.5.1 for the example valve distribution (Figures 3.11 and 3.12). We conclude that characteristic patterns of synthetic activity do not depend on the specific valve distribution, but should rather depend on valve clustering in space and the average inter-valve distance in the simulated channel.

## 3.6 Discussion

The key feature of the model is that pore-pressure diffusion induced by the opening of a valve triggers neighboring valve openings. In spite of the model simplicity, the generated activity exhibits complex patterns of migration and time clustering, comparable to seismic observations. An important result is that, for a given valve distribution, *i.e.* a given permeability structure of the channel, the style of activity depends on the input fluid flux.

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

### 3.6.1 Comparison with LFE activity patterns in Guerrero

We first observe that activity migrates in both the updip and downdip directions within the sweet spot patch at rates that are faster than that from
the sweet spot to the transient zone (Figure 3.1), as in our simulations. On the longest timescale, LFE activity in Guerrero appears time-clustered (Frank et al., 2016, 2014). Activity bursts recur in episodic fashion and are mostly concentrated in the downdip sweet spot patch. At times, during the most vigorous episodes, activity crosses over into the less active updip transient zone. Our model generates a similar pattern when the input flux $q_{in}$ is near threshold values. At the lower end of the dimensionless operating flux range (i.e. $q_c^{break} - q_c^{clog}$), the channel is mostly closed and activity emerges when it opens up in episodic fashion. In contrast, at the upper end of the flux range, the channel is mostly open and activity bursts are concentrated at the times of partial channel closure. Therefore, two different sets of values for valve opening/closing $\delta p$ thresholds, valve permeability and $q_{in}$ could correspond to the Guerrero activity. The first situation is reminiscent of a fault-valve mechanism, as high fault zone permeability correlates in time with seismicity. Several studies interpret tremor activity in Guerrero using this correlation (Cruz-Atienza et al., 2018; Frank, Shapiro, et al., 2015). In the framework of our model, however, the second situation seems as likely as the first one. We thus argue that seismicity bursts may not be systematically linked to an increase of fault permeability and subsequent fluid pressure surge, but may also be due to episodic phases of clogging in the fault, leading to transient fluid pressure build-up and release.

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

### 3.6.2 Scaling numerical results to LFE activity at Guerrero

The main parameter to constrain is the fault-zone diffusivity $D = k/\phi\beta\eta$. Using values from the literature, we end up with a wide range which is not useful. The channel porosity is $\phi = 0.01–0.05$ (Peacock et al., 2011), the fluid-rock effective compressibility is $\beta = 10^{-10}–10^{-9}$ Pa$^{-1}$ (Shvab & Sadus, 2015; Wibberley, 2002), and the dynamic viscosity of supercritical water is $\eta = 10^{-4}–10^{-3}$ Pa.s$^{-1}$ (Shvab & Sadus, 2015). Values for the background permeability $k$ are very sensitive to the measurement method and, above all, to the spatial scale of the measurement. Here we do not consider laboratory-
scale determinations on drill-core samples and focus on measurements that
capture the large-scale channel permeability. This leads to permeability val-
ues that can be as small as $10^{-18}$ m² and as large as $10^{-11}$ m² (Doan et al.,
2006; Frank, Shapiro, et al., 2015; Hendriyana, 2021; Saffer, 2015). With all
these values, we find that diffusivity $D$ lies in a very large $10^{-4}-10^5$ m²·s⁻¹
range.

In order to scale model results to the Guerrero, Mexico, conditions, we
need values for the characteristic spatial scale $\mathcal{X}$, fault-zone diffusivity $D$, and
the characteristic time scale $\mathcal{T}$. For this diffusion process, $\mathcal{T} = \mathcal{X}^2/D$. For the
Guerrero tremor zone, $\mathcal{X} = 100$ km. Using $D = 10^2$ m²·s⁻¹, consistent with
other estimates based on tremor migrations (Hendriyana, 2021), we obtain a
characteristic time scale $\mathcal{T} = 10^8$ s ≈ 3 years. Using values of the hydraulic
parameters that are quoted above, a diffusivity of $D = 10^2$ m²·s⁻¹ implies a
permeability of $k = 10^{-14}-10^{-12}$ m².

We next use the simulated activity in Figure 3.11 as a generic case. Such
a quasi-periodic behavior allows an unambiguous determination of the char-
acteristic recurrence time between large activity bursts, $T_r \approx 0.05$, and their
characteristic duration, $\tau \approx 0.015$, in dimensionless units. Both values are
representative of highly time-clustered activity, which occur for both low and
high fluid fluxes, and in a large number of 29-valve distributions. Scaling
variables back to physical units, we find that $T_r \approx 3$ months and $\tau \approx 2$ weeks,
which are close to the estimates of Frank and Brodsky (2019) for Guerrero:
$T_r^G = 2-3$ months, $\tau^G = 5-20$ days. In physical units, the migration veloci-
ties described in section 3.5.2 scale up to $v = 15-20$ km.hr⁻¹ for short, rapid
migrations, and to $V = 7$ km.hr⁻¹ for the largest migrations that proceed
from bottom to top of the channel. In Guerrero, short migrations propagate
indeed at velocities of a few tens of kilometers per hour (Cruz-Atienza et al.,
2018; Frank et al., 2014), but the larger scale migrations are slower than our
model estimate, with a velocity of about 1 km.day⁻¹. This may be attributed
to the distance that separates valve patches which has a large impact on the
propagation time in a diffusion process.

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

Last, but not least, we evaluate pore pressure and fluid flux values. For
the fluid flux scale, *Q*, we use the estimate of van Keken et al. (2011) for the Cocos subduction. The fluid flux per unit subduction zone width is $Q_T = 10 \text{ Tg.Myr}^{-1}\text{.m}^{-1}$. Using the thickness of the subduction fault zone inferred from geological observations (Angiboust et al., 2015; Tarling et al., 2019), $h = 200 \text{ m}$, we obtain a mass flux per unit area equal to $q_T = Q_T/h = 3 \times 10^{-7} \text{ kg.s}^{-1}\text{.m}^{-2}$.
This estimate is an average value for the whole subduction zone and the actual fluid flow should be focused into channels at several points along the strike of the fault (Ague, 2014; Piccoli et al., 2021). Such a channeling effect may increase the fluid flux by up to one order of magnitude. With this caveat in mind, we assume that the lower fluid flux threshold $q_c^{\text{break}}$, which is equal to 0.1 in dimensionless units in our simulations, is of the same order of magnitude as the dehydration metamorphic flux. We deduce that $q_c^{\text{break}} \approx 3 \times 10^{-7} - 3 \times 10^{-6} \text{ kg.s}^{-1}\text{.m}^{-2}$. We then find that the characteristic pressure scale $P = \eta/\rho k$ $Q \times \mathcal{X} = 30 \text{ MPa} - 30 \text{ kPa}$, using the parameter ranges previously estimated.

The pressure difference between the two ends of the channel $\Delta p$ is close to 1 dimensionless flux unit most of the time. Therefore, in physical units, $\Delta p \approx 30 \text{ MPa}–30 \text{ kPa}$. Each closed valve gets unclogged for $\delta p = 0.04P = 1 \text{ MPa}–1 \text{ kPa}$. The magnitude of the force that is applied to the plug walls depends on area. For a valve in a channel with a $10^3 \text{ m}^2$ cross-section, N. M. Shapiro et al. (2018) were able to match the amplitude of seismic waves with a pressure drop of 5 MPa. Our channel might have a larger cross-section, and hence allow for large forces with comparatively lower pressure drops. The fact that these values are of the same order of magnitude demonstrates a certain self-consistency of our analysis, backed by independent observations of seismic patterns, geological flux measurements, and seismic characteristics of LFEs.

### 3.6.3 Perspectives for diagnosing hydraulic conditions in fault zones

According to our model, the rate and style of seismic activity depends on the hydraulic regime of the permeable channel. Provided that tremor and LFEs are seismic expressions of valve-controlled pore pressure variations, tremor patterns can be used to investigate the hydraulics of fault zones, and eventually other natural systems such as hydrothermal fields and volcanoes where similar seismogenic mechanisms have been invoked (Honda & Yomogida, 1993)).

For given valve parameters that meet the condition for instability ($\delta p_c^{\text{break}} / \delta p_c^{\text{clog}} < k_{hi} / k_{lo}$, equation 3.16), the activity regime depends on a dimensionless input rate of fluid into the channel. This dimensionless input rate is calculated with respect to a Darcy flux for some reference permeability value (chosen to be
$k_{hi}$ here). For application to natural conditions, the same dimensionless input rate can be obtained for different pairs of fluid flux and permeability values. Conversely, the same fluid flux may be responsible for different activity styles in different parts of the subduction interface depending on the local permeability value. In addition, time changes of the fluid flux, due for example to the subduction of more or less hydrated parts of an oceanic plate, may induce changes of activity.

The rate of fluid release by metamorphic reactions depends on the subduction rate and on the type of material that dehydrates. It cannot be the same everywhere as, for example, subduction rates vary by almost one order of magnitude. Therefore, comparing seismicity patterns in different subduction zones is likely to prove very useful. In particular, specific attention should be paid to subduction zone segments that do not seem to generate LFE activity (e.g. Bocchini et al., 2021). As shown in this paper, this may be due to a fluid flux that is either very high or very small.

### 3.6.4 Perspectives for exploring model configurations

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

In all the cases reported here, all valves are identical, with the same width $w$, the same low permeability value $k_{lo}$ in a state of clogging and the same values of opening/closing pressure thresholds $\delta p_c^{break}$ and $\delta p_c^{clog}$. In natural conditions, the hydraulic and mechanical properties of the permeable fault zone should depend on the ambient pressure and hence on depth. The depth dependence of valve properties may explain — within the framework of our model — why tremor bursts recur more frequently at the deeper end of LFE areas (Frank et al., 2014; Wech & Creager, 2011). As pressure and temperature increase with increasing depth, both the background and valve perme-
abilities should decrease under pressure. Depending on how permeability
changes with depth, a flux that is near-critical for updip valves may well be
in the intermediate range for the downdip ones. This may account for the
behavior of the Guerrero tremor area, where activity is much more frequent
and continuous downdip than updip. With lower permeability, migrations
should also be slower. With that in mind, the fast migration rates observed in
the sweet spot in Guerrero would imply that valve density might increase in
this depth range.

Model parameters that must be explored include the valve characteristics (length, permeability in a state of clogging), the distribution of valves and the total number of valves. However, the valve instability criterion (equation 3.16) constrains how permeability and pressure threshold values can vary. As discussed above, we expect that this condition is met in the pore-pore throat configuration that is relevant to many natural permeable systems (Beckingham, 2017; McDowell-Boyer et al., 1986). Ultimately, a more realistic description of the process should include time-dependent opening and closure processes, but very few constraints on the associated time-scales are available.

We have shown in section 3.5.3 that neighboring valves may act as a single
macro-valve in some conditions (Figure 3.12), and this is worthy of a dedi-
cated investigation. Determining which valve spacing and hydraulic condi-
tions enable this behavior is likely to provide key insight into the interaction
mechanisms that generate large bursts of activity. One could also allow for
some random distribution of valve properties. However, we may draw from
our current understanding and predict that the end result would be that only
a subset of valves get activated for a given flux. Other boundary conditions
than those used here are also worth investigating. For example, the subduc-
tion channel may not be connected to a very large permeability medium and
allowing for an output flux that depends on pressure may lead to large-scale
cyclic activity. In addition, dehydration surges may occur due for example
to kinetic constraints on metamorphic reactions or to highly hydrated por-
tions of the downgoing plate crossing the phase-change boundary. Ongoing
and future work with similar models will investigate these mechanisms and
further explore parametric controls on valve behavior.

**3.6.5 Coupling between hydraulic transients and fault slip**

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

## 3.7 Conclusion

In this work, we have developed a model to explain the tremor patterns in subduction zones as symptoms of rapid, strong, localized pressure transients emerging from intermittent fluid transport in a permeable fault interface. We model fluid transport in a 1D permeable channel, fed by a constant input flux. Based on field observations in faults, and laboratory experiments in porous media, we have argued that clogging and un-clogging events occur in the permeable channel, causing large changes of permeability. Based on simple physical and geological arguments, we have shown that such events can alternate in a simple valve mechanism. The strong pressure transient that occurs when a valve opens is able to generate a seismic event, and we have therefore considered valve activity as synthetic seismic activity.

We have found that the behavior of a permeable channel containing valves depends on two valve parameters (pressure ratio $\delta p_c^{break}/\delta p_c^{clog}$ and permeability ratio $k_{hi}/k_{lo}$). In one case ($\delta p_c^{break}/\delta p_c^{clog} > k_{hi}/k_{lo}$), the channel is in a stable steady-state regime with no valve activity for all values of the input flux. In the other case ($\delta p_c^{break}/\delta p_c^{clog} < k_{hi}/k_{lo}$), the channel is unstable with valves that open and close repeatedly if the input flux is within a well-defined operating range. We have shown that only the latter is relevant to geological systems. Therefore, in this framework, the value of the fluid flux going through a fault zone is a major control on the observed seismicity. For values of the flux outside the operating range, the fault zone is permanently closed or open. For values within the operating range, permeability valves opening and closing in cycles generate permanent activity and intermittent fluid
transport. In this permanently unstable regime, spatiotemporal patterns of
modeled tremor emerge from valve-to-valve interactions, and are very simi-
lar to those observed in many subduction zones.

Within the operating range for permanent activity, the input flux controls
the dynamics of the permeable system. The system is mostly closed for low
flux values and mostly open for high values, with time-clustered bursts of
activity. For intermediate flux values within the range, the system is mostly
open, and closes in quasi-periodic bursts of activity. Furthermore, the activ-
ity rate is highest and almost constant through time. Using this model and
more elaborate versions of it, one may hope to use microseismicity patterns
to probe the hydraulics of permeable fault zones and how they change with
time.
# APPENDICES TO CHAPTER 3

## Contents

3.A Numerical accuracy and convergence tests . . . . . . . . . . . . 101

3.B Derivation of fracture aperture from permeability values . . . . 104

## 3.A Numerical accuracy and convergence tests

### Analytical derivation of the diffusion of a pressure ramp

In order to test the accuracy and convergence of our numerical solver for the diffusion equation, we compare the numerical and analytical solutions of the diffusion of a fluid pressure profile that could have been accumulated across a closed valve: a pressure linear ramp from 1 updip to 0 downdip, between $x_1$ and $x_2$:

$$p_0(x) = \mathcal{H}(x_1 - x) + \frac{x_2 - x}{x_2 - x_1} \times (\mathcal{H}(x - x_1) - \mathcal{H}(x - x_2)) \quad (3.18)$$

where $\mathcal{H}$ is the Heaviside function. We choose to test our numerical solver in the case the diffusion of this particular pressure profile as it is analog to the rapid diffusion event at the opening of a valve, and because of the discontinuities in pressure gradient, should be the major source of numerical errors in our solution.

As the diffusion equation is a linear partial differential equation, we express the analytical solution of this problem as the sum of the solution in time of the diffusion of a Heaviside step located in $x_1$: $S_0(x) = \mathcal{H}(x_1 - x)$, and of a limited ramp $R_0$ between $x_1$ and $x_2$: $R_0(x) = \frac{x_2 - x}{x_2 - x_1} (\mathcal{H}(x - x_1) - \mathcal{H}(x - x_2))$.

For a given pressure field $p(x, t)$ defined for $x$ and $t$ in $[-\infty, +\infty]$, following a 1D diffusion equation with diffusivity $D$:

$$\partial_t p = D \partial_x^2 p . \quad (3.19)$$
Using the Fourier transform over the space domain, $\mathcal{F}(p(x, t)) = \tilde{p}(\xi, t)$, the diffusion equation transforms to a first order differential equation.

$$ \partial_t \tilde{p}(\xi, t) = -D\xi^2 \tilde{p}(\xi, t) \qquad (3.20) $$

with solution $\tilde{p}(\xi, t) = \tilde{p}_0(\xi) e^{-D\xi^2t}$. Transforming back:

$$ p(x, t) = p_0(x) * K(x, t) \qquad (3.21) $$

where $p_0(x) = p(x, t = 0)$. Using the inverse Fourier transform, we get $K(x, t) = \frac{1}{\sqrt{4\pi Dt}}e^{-\frac{x^2}{4Dt}}$, the diffusion kernel.

Any transient solution can therefore be computed through the convolution of the initial profile with the diffusion kernel $K(x, t)$. We apply this method to obtain the analytical solution of the diffusion in space and time of the ramp $p_0(x) = R_0(x) + S_0(x)$.

## Results of accuracy and convergence tests

Figure 3.A.1 summarizes the results of our convergence and accuracy tests. Both the accuracy and convergence of the numerical diffusion solver are tested by comparing the numerical and analytical solutions for the diffusion of the pressure ramp described above. The error is given by the difference between solutions. The discontinuities of pressure gradient in $x_1$ and $x_2$ are responsible for the main error signal, but the physical diffusion smoothes rapidly smoothes it in space. At any given time, the error stays below $4.10^{-3}$ pressure units, for a step of height 1. The RMS error around the valve is at least one order of magnitude lower, demonstrating an overall accurate and stable numerical method.

Convergence is defined as the increase in accuracy of a numerical solution when the discretization it is computed on is refined. Figure 3.A.1 demonstrates the convergence of our numerical scheme in our test case: at three given times after the start of the diffusion, it shows that accuracy increases with finer discretization, or higher number of space steps $N_x$. As the spatial discretization is refined, the temporal discretization follows using $\delta t = h^2/2D$ to define the time step $\delta t$ as a function of the space step $h$.
<figure><img src="image_51.png" /><figcaption>a Numerical v. theoretical diffusion - $N_x = 500$</figcaption></figure>

c

<figure><img src="image_52.png" /><figcaption>b Error $p_{num}$ - $p_{th}$ in space - $N_x = 500$</figcaption></figure>

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

d

<figure><img src="image_54.png" /><figcaption>Convergence test: RMS error at three times v. $N_x$</figcaption></figure>

**Figure 3.A.1** — Accuracy and convergence tests, based on the diffusion of a pressure ramp in time and space. (a) Comparison between analytical solution in black dashed line and numerical solution, colored through time. Last time $t = 4.10^{-4}$. (b) Error between numerical and theoretical solutions along the domain, at three given times. (c) Evolution in time of the maximum of the error and root-mean-squared (RMS) error over the domain. (d) RMS error over the domain, at three times during the diffusion of the pressure ramp, as a function of the number of space steps. Convergence of our numerical scheme is shown demonstrated by the increase in accuracy as the space-time discretization is refined.
## 3.B Derivation of fracture aperture from permeability values

Using the value of the open permeability of the channel $k = 10^{-14}-10^{-12}m^2$ obtained in Section 3.6.2, and assuming that the permeability is due to fractures oriented in roughly along the plate interface, we can deduce the aperture of the fractures a. Norton and Knapp (1977) formulated the laminar fluid flux into parallel fractures, and found that it can reduce to a Darcy-like form, with a *fracture permeability*

$$k_{fr} = \frac{(2a)^3 n_{fr}}{12} \qquad (3.22)$$

where $n_{fr}$ is the frequency of fractures per unit of length across the parallel fracture planes. In a well-preserved outcrop of a subduction channel in the Zagros suture zone (Iran), Muñoz-Montecinos et al. (2021) estimates that around 8% of the surface of the outcrop — both serpentinite matrix and blueschist block — is composed of veins. Assuming we can use this surfacic proportion to directly compute the frequency of fractures in the intact, *in situ* medium, the proportion of fractures can be written $p_{fr} = 0.08 = n_{fr} \times a$, therefore:

$$k_{fr} = \frac{2a^2p_{fr}}{3}. \qquad (3.23)$$

We can then use the values of permeabilities mentioned above to find a range of fracture aperture in the medium:

$$a = \sqrt{\frac{3k_{fr}}{2p_{fr}}} \approx 4 \times 10^{-7} - 4 \times 10^{-6}m^2. \qquad (3.24)$$

The fractures in such a medium therefore around the micrometer scale, in the same range as estimated for crack-seal fractures seen in metamorphic rocks (Etheridge et al., 1984; Ramsay, 1980).
<figure><img src="image_55.png" /></figure>
<figure><img src="image_56.png" /></figure>
# CHAPTER 4

## EMERGENCE AND CONTROL OF TREMOR IN-TERMITTENCE AND PERIODICITY

### Contents

<table><tr><td>4.1</td><td>The intermittence of tremor in subduction zones</td><td>108</td></tr><tr><td>4.2</td><td>Characterizing clustering intensity and recurrence timescales</td><td>113</td></tr><tr><td>4.2.1</td><td>Point process description</td><td>114</td></tr><tr><td>4.2.2</td><td>Characterizing clustering intensity</td><td>114</td></tr><tr><td>4.2.3</td><td>Measuring recurrence timescales</td><td>116</td></tr><tr><td>4.3</td><td>A valve model to describe fault zone permeability</td><td>119</td></tr><tr><td>4.3.1</td><td>The subduction's hydraulics as a valve system</td><td>122</td></tr><tr><td>4.4</td><td>Characterizing source interaction in a two-valve system</td><td>124</td></tr><tr><td>4.4.1</td><td>Valves interact through pressure transients</td><td>124</td></tr><tr><td>4.4.2</td><td>Valve spacing controls activity synchronization</td><td>125</td></tr><tr><td>4.5</td><td>Effects of valve proximity in complex valve systems</td><td>128</td></tr><tr><td>4.5.1</td><td>Statistical description of valve distributions</td><td>128</td></tr><tr><td>4.5.2</td><td>Valve distribution control on activity</td><td>129</td></tr><tr><td>4.5.3</td><td>Parameterizing valve proximity</td><td>131</td></tr><tr><td>4.6</td><td>Emergence and variability of synchronization in complex valve systems</td><td>133</td></tr><tr><td>4.6.1</td><td>Patterns of activity clustering with high valve proximity</td><td>133</td></tr><tr><td>4.6.2</td><td>Valve proximity controls recurrence timescales and periodicity of activity</td><td>135</td></tr><tr><td>4.6.3</td><td>Input flux control on criticality and emergence of synchronization</td><td>136</td></tr><tr><td>4.7</td><td>Discussion</td><td>138</td></tr><tr><td>4.7.1</td><td>From observed to simulated tremor activity patterns</td><td>138</td></tr><tr><td>4.7.2</td><td>The hydraulic segmentation of subduction faults</td><td>141</td></tr></table>
<table><tr><td>4.7.3</td><td>Hydrological processes shaping tremor intermittence in Shikoku</td><td>143</td></tr><tr><td>4.8</td><td>Conclusion</td><td>146</td></tr></table>

## 4.1 The intermittence of tremor in subduction zones

In subduction zones where it is detected, tremor and low-frequency earthquake (LFE) activity clearly appears to be organized in time. It is intermittent, proceeding in periods of quiescence and episodes of intense activity often called *bursts*, or *episodes* (Brudzinski et al., 2010; Frank et al., 2014; Idehara et al., 2014; Obara, 2002; Rogers, 2003; Shelly et al., 2007b). The strongest episodes of tremor or LFE activity are associated with slow slip events (SSE) on the fault interface (Frank, Radiguet, et al., 2015; Obara et al., 2004; Rogers, 2003; Shelly et al., 2006) and episodes of strong fluid pressure change in the vicinity of the source region of the activity (Frank, Shapiro, et al., 2015; Gosselin et al., 2020; Nakajima & Uchida, 2018; Tanaka et al., 2018; Warren-Smith et al., 2019). The cycle of mechanical and hydraulic stress accumulation and release seems to be the main control on the intermittence of the episodic tremor and slip phenomenon. As fluid is released at depth by dehydration reactions in the subducted oceanic crust (R. N. Anderson et al., 1976), regions of low permeability above the slab can seal fluid around the interface, leading to a high fluid pressure anomaly, observed by geophysical measurements (Calvert et al., 2011; Shelly et al., 2006; Wannamaker et al., 2014) and in exhumed paleo-subduction channels (Angiboust et al., 2015; Behr & Bürgmann, 2021). The locally high fluid pressure leads to decreased effective normal stress on the fault, and therefore lower frictional strength of the plate interface, and slip can nucleate with minor changes of the stress field. As the slip front propagates, fueled and stabilized by the high fluid pressures (Segall et al., 2010), the associated fracturation should open pathways for the fluid to diffuse, along the fault interface or to the overriding plate, towards the surface (Frank, Shapiro, et al., 2015). Fluid pressure in the source region therefore drops as fluid diffuses away in the newly created permeability. As the permeability around the fault heal, a new cycle of pressurization happens again (R. H. Sibson, 2013). This cycle of slip, fluid pressure and strength of the fault interface has been coined *fault-valving* by R. Sibson (1992). The recurrence behavior of large bursts of tremor should therefore be controlled by the dynamics of both the tectonic loading and strength on the interface, and the fluid pressure accumulation and release in the subduction shear zone. Mon-
itoring the characteristics of intermittence of tremor activity on the fault thus
provides insight on the hydrological and mechanical state of the fault.

The intermittence of seismic activity can be analyzed using statistical tools
describing the temporal clustering of events (Beaucé et al., 2019, 2022; Frank
et al., 2016; Idehara et al., 2014; Poiata et al., 2021). In order to characterize
it, one can look at the extent with which events occur within clusters com-
pared to the background activity, the distribution size and duration of those
clusters, the recurrence delays between them. There is a fundamental differ-
ence in the way tremor and LFEs are detected and cataloged. Tremor is a
long, non-impulsive signal, and although its source can be located, it cannot
be reduced to one or several events, its activity is therefore often cataloged by
identifying periods of activity instead of individual events (e.g. A. Husker
et al., 2019; Wech & Creager, 2008). Low-frequency earthquakes however
can perfectly be described as a point process: they are individual, repeat-
ing, localized events. They are often considered to be elementary tremor
events (Shelly et al., 2006), a direct proxy for tremor activity. In order to
characterize the activity of tremor using point process statistical tools, we
will work with LFEs catalogs, considering they are just another observation-
nal representation of the same phenomenon. While in some locations large
bursts of tremor activity recur with a visible and well-defined recurrence
timescale, quasi-periodically (Brudzinski & Allen, 2007), in other locations,
individual clusters and their recurrence timescales are harder to define be-
cause the temporal patterns of activity appear to be scale-free (Frank et al.,
2016; Poiata et al., 2021). These different behaviors can be understood statis-
tically as a deviation from a random and independent occurence of events,
as event occurence depends on the occurence of events closeby. Clustering
should therefore emerge from a process that allows seismic sources to inter-
act and synchronize. The stronger the synchronization, the more clustered in
time tremor activity is, thus producing more isolated bursts of activity. With
this framework in mind, we think of perfectly periodic burst recurrence as an
extreme case of time-clustering: as tremor sources become perfectly synchro-
nized, all sources of a region will activate together, producing well-defined
bursts of a characteristic size, that should recur with a periodicity linked to
the specific process of energy accumulation and release that drives tremor
activity.

The degree of activity synchronization clearly varies between subduction
zones, but it is also strikingly segmented along the strike of subduction zones.
In the Shikoku tremor zone in Japan, Poiata et al. (2021) shows that seven
along-strike segments can be identified in the active tremor region, each char-
acterized by a characteristic temporal clustering of LFE activity. In Figure 4.1,
using their LFE catalog, we reproduce their observation by measuring the ex-
<figure><img src="image_57.png" /><figcaption>Figure 4.1: Caption on next page.</figcaption></figure>
**Figure 4.1** — On previous page. Patterns of temporal clustering and recurrence of LFE activity, in Shikoku, Japan (2013–2017). (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 et al. (2021). Depth contours of the slab from Slab2 (Hayes, 2018). The dotted white line indicates a linear fit to the event distribution in space, used as the along-strike direction. (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$ ($c = 0$, no clustering, $c = 1$, all events occur in clusters, see Section 4.2 for details). Bins of 20 km, with a 10 km overlap. (d) Recurrence timescale of clusters along strike, with the same spatial bins as in (c). Clusters of events are detected and their recurrence measured with the algorithm described in Section 4.2. The blue, shaded area around the measurement represents $1-\sigma$ interval around the estimated recurrence timescale. The blue dotted line shows a visually estimated recurrence, based on counting the number of high-activity episodes in the observation period. (e) LFE activity along-strike in time, each black dot represents an event. Activity is clearly segmented along strike, in three major segments (see text for details).

tent to which events cluster in bursts (clustering index $c$) and the timescale of recurrence of those bursts, in bins of 20 km along the strike of the subduction, with 10 km overlap. More details on the method can be found in Section 4.2. In Figure 4.1a and b, the density of events in space reveals a clear segmentation along-strike. The activity in time and space (Fig 4.1e) and the measured clustering intensity and recurrence timescale (Fig 4.1b and c) allow us to outline three major behavior along-strike. In the western zone (distance along-strike < 125–130 km), activity is characterized by a relatively low level of clustering, a high background rate of events. Bursts are frequent, and recur on timescales of 2–3 months, although with a strongly variable delay. In the easternmost zone (distance along strike > 180–190 km), activity is much 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 zones, the activity is least intense, with a very low background rate, but large transients of activity occasionnaly cross it, every 6 months to slightly less than a year. It is the most temporally clustered activity of the three regions. In the few zones where tremor or LFEs have been cataloged on a wide enough extent along strike, their activity is similarly segmented (Brudzinski & Allen, 2007; A. Husker et al., 2019). In Figures 4.A.2 and 4.A.1, we display the tremor activity in the Cascadia (tremor catalog from Wech & Creager, 2008) and Central American subduction zones (tremor log from A. Husker et al., 2019). In both cases, segments of the subduction zone can be clearly delineated based on the characteristics of the tremor activity they produce: the level with which it is concentrated in bursts and how frequently those recur.
The observed segmentation of LFE and tremor activity patterns along the subduction interface is thought to take its roots in structural heterogeneity and spatially organized dynamic processes in the fault zone (e.g. Poiata et al., 2021). On the one hand, the structural properties of the subduction interface, and heterogeneities thereof, crucially determine the distribution of stress and strength on frictional interfaces in the shear zone, and shape how fluid ac- cumulates or circulates in and through the subduction interface. The local geology of the overriding plate (Brudzinski & Allen, 2007), the local dip of the slab, which controls the temperature on the interface (Maury et al., 2018), and the topography of the seafloor on the subducting slab (Ide, 2010) are key parameters controlling the state of stress and therefore the frictional and hy- draulic behavior of the interface. These structural characteristics of the slab shape heterogeneous distributions of tremor activity as they affect the dis- tribution of frictional weakness (Cattania & Segall, 2021; Kano et al., 2018; Nakajima & Hasegawa, 2016) and the conditions of fluid circulation in and around the interface: the dehydration rate, how much fluid escapes to the overriding plate or is trapped and channeled along it (Audet & Bürgmann, 2014; Gosselin et al., 2020; Halpaap et al., 2019). On the other hand, dynamic processes like damaging and healing of permeability or friction, or fluid chan- neling processes can produce heterogeneous patterns of activity on the fault by self-organizing spatially (Halpaap et al., 2019). The memory of previous processes would in this case be the origin of the apparent segmentation of the fault zone, as opposed to intrinsic properties of the components of the subduction system.

Several studies have begun to explore the role of the spatial heterogene-
ity of frictional properties to explain the different intermittence behavior ob-
served along the strike and dip of the subduction fault (Luo & Liu, 2019;
Wech & Creager, 2011). Few studies include permeability and fluid pres-
sure variations to study the role of the hydrological system on shaping pat-
terns of tremor activity (e.g. Cruz-Atienza et al., 2018). However, as evi-
denced in the geological record of subduction channel processes, fluid pres-
sure and permeability are spatially heterogeneous and temporally variable
in the tremor source region, through processes of vein opening and heal-
ing (Platt et al., 2018; Taetz et al., 2018; Tarling et al., 2021), hydrofractura-
tion (Muñoz-Montecinos et al., 2021) and pressure cycling (Angiboust et al.,
2015; Kotowski & Behr, 2019). These strong and sudden variations of perme-
ability and pressure are thought to be closely associated with the generation
of seismic waves that are then recorded as tremor or LFEs (Kotowski & Behr,
2019; Muñoz-Montecinos et al., 2022; N. M. Shapiro et al., 2018). In Chapter 3,
we have elaborated a model of how fluid pressure transients and variable
fluid transport properties of the fault zone can generate complex patterns of
seismicity — temporal clustering, periodic behavior, migrations. We have shown that in this system, the emergence of patterns of activity is based on interactions through fluid pressure transients between elementary permeability valves, which simulate the local changes of permeability and fluid pressure that can be associated with the source processes of LFEs. In this chapter, we will demonstrate how temporal clustering emerges from spatial clustering of our model seismic sources, how hydraulic interactions between tremor sources should be controlled by the local structure of the permeable network in the fault. In turn, we will discuss how the permeability structure and therefore spatial source distribution can be shaped by larger scale structural features of the subduction system, as the topography of the subducted oceanic plate.

## 4.2 Characterizing clustering intensity and recurrence timescales

In this study, we aim at analyzing seismicity through the prism of how events cluster in time, on which timescales those clusters recur, and how regularly they do so. 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. Those measures 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. Beaucé et al., 2019) — regarding which specific patterns the measured timescales are related to. In proceeding this way, we 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 (Frank et al., 2016; Idehara et al., 2014; Poiata et al., 2021), or intermittent microseismic activity in general (Beaucé et al., 2019, 2022), where scale-free behavior is often observed. In our case however, it will allow us to assess the degree of synchronization that can be found in activity, assuming that periodicity emerges from an extreme case of synchronization of sources.

The measures and algorithm we will describe in this section 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 simple understanding of the characteristics of synthetic seismicity for large batches of simulations (several tens of thousands). Figures 4.2 and 4.3 describe the two steps of this analysis.
### 4.2.1 Point process description

Based on Lowen and Teich (2005), we introduce elements of the point process formalism that will be used in the rest of the study. In the point process formalism, an event catalog — earthquakes, LFEs, valve openings — can be represented by the function $N_{ev}(t)$, a monotonically increasing function of time $t$, that increases in steps of 1, when the time is equal to an event time $t = t_i$. The point process can thus be formally written as $dN_{ev}(t)$. However, a more practical and as complete description of the process is to use the series of event times

$$ \{t_i\} \text{ for } 1 \le i \le N_{ev}^{tot}, \qquad (4.1) $$

where $N_{ev}^{tot}$ is the total number of events considered in the analysed subset. One can also work with the *inter-event times* $\Delta t_i$, defined as the time after the previous event

$$ \{\Delta t_i\} = t_i - t_{i-1} \text{ for } i > 1. \qquad (4.2) $$

The *event count time series* $C_k(\delta t)$, or *sequence of counts*, is another representation of the point process. It counts the number of events in successive and contiguous bins of duration $\delta t$, between time $t = k\delta t$ and $t = (k + 1)\delta t$:

$$ C_k(\delta t) = N_{ev}[(k + 1)\delta t] - N_{ev}[k\delta t] \qquad (4.3) $$

### 4.2.2 Characterizing clustering intensity

The objective of this analysis is to quantify the extent to which events occur within clusters, rather than as a more homogeneous process in time. We detect clusters as peaks in the event count time series (see equation 4.3): a count of the number of events in bins of length $\delta t$. In our parametrization, peaks of activity are defined as all consecutive bins in which the event count is higher than a given threshold. In order to select more intense periods of activity, peaks are defined as periods during which activity is more intense than 66% of the time, and should therefore cover 34% of the analysed period (Fig 4.2, 4.B.1a and b, 4.B.2a and b): the 34% most active $\delta t$-bins. We can define the peak cover $f_{peaks}$ as the fraction of the time covered by peaks, compared to the total time of the analyzed period. Because of the discrete nature of the event counting process, the fraction of the time that is covered by peaks $f_{peaks}$ is never exactly 34% of the time:

$$ f_{peaks} = T^{peaks}/T^{tot} \approx 0.34 \qquad (4.4) $$

where $T^{peaks}$ is the cumulative duration of peaks in the time series, and $T^{tot}$ the total duration of the time series.
Characterizing the temporal clustering of activity

<figure><img src="image_58.png" /><figcaption>**Figure 4.2** — Clustering analysis of (a), (b) a simulated, highly clustered activity in a clustered valve system with $N_v = 43$ valves ($\bar{d}_v = 1w_v$), $u = 0.3$, under a input flux of $q_{in} = 0.16$, (b), (c) a simulated activity without temporal clustering, in a Poissonian valve system with $N_v = 29$ valves ($\bar{d}_v = 2w_v$), $u = 1.0$, under a input flux of $q_{in} = 0.16$, and (e), (f) a Poisson process, with event rate $r = 30$ events per $\delta t = 3 \times 10^{-4}$ scaled time units. In (a), (c), and (e), the event count time series are represented in black with the threshold defining peaks as the red horizontal line, corresponding to the 66% quantile of the event count distribution, represented in (b), (d) and (f). The fraction of time that peaks represent $f_{peaks}$, the proportion of events within peaks $p_{ev}^{peaks}$ and clustering index $c$ for each case are listed on the right.</figcaption></figure>
Although peaks always cover about a third of the studied period, in a clustered time series, most events — more than a third — will occur during peaks, that is 34% of the time. To discriminate between a clustered and not clustered activity, we will therefore work with the proportion of events that occur within peaks, $p_{ev}^{peaks}$. When the activity can be described as a Poisson point process, or if the event count is normally distributed around an average value, the events are as much in “peaks” that in the “background”, as defined by the 34%-threshold (Fig 4.2c, d, e and f). The number of events in peaks $N_{ev}^{peaks}$ should therefore sum to about a third of events, as peaks are defined the upper tritle of event counts. On the other hand, if the events cluster in time, events should concentrate in peaks. Peaks still cover about 34% of the time series, but most event (at least more than 34% of them) occur in peaks (e.g Fig 4.2a and b). A simple measure of clustering can therefore be constructed by mapping linearly the proportion of events in peaks $p_{ev}^{peaks}$ to 0 – 1. We thus define the clustering index c: 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.

$$c = \frac{p_{ev}^{peaks} - f_{peaks}}{1 - f_{peaks}} \qquad (4.5)$$

where $f_{peaks}$ is the fraction of the studied period that is covered with peaks, that is always approximately $f_{peaks} \approx 0.34$. In Figure 4.2, three cases show values of the clustering index c for two simulations (highly and not temporally clustered) and a Poisson process. We will use the clustering index c as a measure of clustering intensity in observed LFE catalogs (Figure 4.1) and in simulated event catalogs (Section 4.5 and 4.6, Figures 4.7 and 4.9).

### 4.2.3 Measuring recurrence timescales

In both real cases (Figure 4.1) and simulations (Figure 4.7), clusters of seismic events can exhibit identifiable, characteristic timescales of recurrence. Our aim is to quantify these recurrence timescales, and estimate how variable they are throughout the studied period. To do so, we design an automated analysis that explicitly identifies clusters of events based only on the proximity of their occurrence time $t_i$, classifies clusters based on the number of events they are made of, and finally measures the recurrence of clusters of a same class of size. We will perform this analysis only in cases when the activity is clustered: when more than 50% of events occur in peaks ($c > 0.25$) — peaks are defined as the 34% most active periods of the studied period, see above Subsection 4.2.2. Figure 4.2c to f show two examples of simulations that exhibit no clustering of activity, on which no timescale analysis is therefore per-
formed.

The first step of the recurrence analysis is to identify clusters of events, based on their proximity in time. In order to de-noise the event time time-series, we start by removing all events outside of peaks. We then regroup the remaining events into clusters, by performing a density-based clustering algorithm (DBSCAN, (Ester et al., 1996), scikit-learn implementation (Pedregosa et al., 2011)) on their occurence times $t_i$. The clustering algorithm labels events as isolated or belonging to a cluster based on two parameters: the neighborhood radius (a delay in time, in our case) $\epsilon$ and the minimum number of events $N_{min}$ that have to be at a delay of $\epsilon$ from one another so that they are considered in the same cluster. Those two parameters therefore provide an estimation of the minimum density of events needed around an event for it to be in a cluster. By iterating from event to event in a dense neighborhood, the algorithm identifies which events are within clusters, and which are isolated. We choose $\epsilon = \max_i(\Delta t_i)/4$, a fourth of the maximum inter-event delay measured in the de-noised event times, and $N_{min} = 10$ for simulated seismicity ($N_{min} = 300$ for the Shikoku case, in Figure 4.1), based on the total number of events in each case. At the end of this step, all events are affiliated to a cluster (labelled) $k$, or identified as isolated (noise). In order to avoid windowing artifacts, we discard the events of the first and last identified clusters as noise. In the following paragraphs. Figure 4.3a and b depict this first step of the recurrence analysis.

In some cases, the detected clusters can be grouped into identifiable classes of size, and in each class, the clusters recur with a specific delay. In this second step, our objective is therefore to group event clusters based on their sizes $S_k$, defined as the number of events that belong to it. As cluster sizes sometimes span orders of magnitudes, we perform the DBSCAN clustering algorithm on the logarithm base 10 of cluster sizes. For the DBSCAN parameters, we choose the neighborhood distance $\epsilon = 0.15 * \log_{10}(S_{max}/S_{min})$ where $S_{max}$ and $S_{min}$ are the maximum and minimum cluster sizes in the analyzed case, and $N_{min} = 2$. Those parameters successfully group clusters into classes of different characteristic sizes, while allowing for variability in size within a given class of clusters. At the end of this step, all clusters are labelled as belonging to one or several classes, based on their size. It is more common that one size of cluster only is identified, but in some cases, two (Fig 4.3) or more rarely three classes can be identified with certainty. In the example simulation shown in Figure 4.3b and c, the algorithm detected a class of large clusters (in hot colors) and a class of small clusters (in cold colors).

The final step is to measure the recurrence delay between clusters. Within a class of clusters, we measure the delays $\Delta T_k$ between the centroids $T_k$ of contiguous clusters (defined as the average event occurence time within that
Characterization of event clusters recurrence

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

<figure><img src="image_60.png" /><figcaption>**Figure 4.3** — Recurrence analysis of a simulated, highly clustered activity in a clustered valve system with $N_v = 43$ valves ($d_v = 1w_v$), $u = 0.3$, under a input flux of $q_{in} = 0.16$ (activity count in Figure 4.2a). In a first step, the algorithm detects clusters of events using the occurrence times. (a) Time-dip representation of activity, each point corresponds to an event. Events are colored according to which cluster they are labelled into. The grey dots are labelled as noise, either in the de-noising step, or by the clustering algorithm. The clustering algorithm is performed only on occurrence times, represented as horizontal lines in (b), with the same color-coding as in (a). The centroid of each cluster $T_k$ is computed as the average occurrence time of events that belong to each one of them. In a second step, clusters of events are grouped into classes, based on their sizes. Two sizes of clusters are detected for this simulation, small clusters in blue shades, and large clusters in orange to red shades. (c) Time-dip representation of events, with events colored according to the cluster and class they belong to. The delays $\Delta T_k$ between cluster centroids of each class are indicated with arrows. (d) The size-duration distribution of clusters of events. The classifying step only considers the size of clusters.</figcaption></figure>
cluster). As it appears that both in simulations and real cases the smaller clusters occur more frequently than, and therefore in between, larger clusters, contiguous clusters are considered to be clusters of the same size with no larger clusters occurring in between. In practice, the delays between two clusters is measured only when there are no larger cluster in between (see Figure 4.3c). The average delay mean$_k(ΔT_k)$ between clusters of a class gives an estimation of the recurrence timescale for this class, and the standard deviation std$_k(ΔT_k)$ allows to estimate how variable it is. The ratio of the standard deviation to the average recurrence timescale std$_k(ΔT_k)/mean$_k(ΔT_k)$, called the coefficient of variation, gives an estimation of how regular the recurrence is, or how periodic it is.

Figures 4.B.1 and 4.B.2 provide additional examples of the clustering and recurrence analysis for the highly clustered, loosely periodic activity of segment 6 of the Shikoku tremor zone (Poiata et al., 2021) and finally an intermediate case, where the simulated seismicity is weakly clustered, and the recurrence is measured with a strong variability.

## 4.3 A valve model to describe fault zone permeability

As detailed in Chapter 3, Section 3.3, we describe fluid circulation along the permeable fault interface as being focused in channels of higher permeability than the surrounding medium. The channels are one-dimensional, and their permeability is both heterogeneous and dynamic. In this description, the fluid pressure $p(x, t)$ evolves according to a diffusion equation (Equation 3.5), with a diffusivity

$$D(x, t) = \frac{k(x, t)}{\eta\beta\phi}$$

where $η$ is the fluid viscosity, $β$ a fluid-matrix compressibility, $φ$ the fault zone permeability, and $k(x, t)$ the permeability of the fault zone. Because we dimensionalize the system, the main control parameter on diffusivity is $k(x, t)$.

Permeability $k(x, t)$ is described in a discrete, binary manner: short segments of the channel have a lower permeability $k_{lo}$ than the surrounding $k_{hi}$, and can switch between a closed, low-permeability $k_{lo}$ state, to an open, high-permeability $k_{hi} = 20k_{lo}$ state, in response to the pressure differential $δp$ that builds up across the segment. Those dynamic permeability segments are called *valves* and mimic the dynamics of clogging and unclogging of pore throats by colloids or cataclastic fines. In our model, the rapid changes of fluid pressure that follow unclogging of a segment can act as a seismic source or trigger one, and each opening of a valve can be considered an elementary
seismic event. For a detailed description of the clogging/unclogging mechanisms observed in faults, their link to seismicity, and how we model them, the reader can refer to Chapter 3, Section 3.3. In this section, we will succinctly describe the valve mechanism as it is built into our model, in order to explain which valve parameters control its opening/closing cycle, and how.

In our description, the switch between the closed and open state (and vice versa) occurs instantaneously, when $\delta p$ reaches a given threshold. If the valve is closed, it impedes fluid flow around it and $\delta p$ can increase above the threshold $\delta p_c^{break}$. In this case, the permeability unclogs and the valve breaks open, its permeability increases from $k_{lo}$ to $k_{hi}$. On the other hand, if the valve is open and pressure diffusion across the valve lowers $\delta p$ to $\delta p_c^{clog}$, the permeability clogs, and the valve closes: valve permeability drops from $k_{hi}$ to $k_{lo}$.

The fluid input rate $q_{in}$ into the system is an essential control parameter that governs if the valve can break open when closed, and close when open. In Chapter 3, we show that two threshold flux values can be derived from the threshold $\delta p$ using the Darcy law for the flux through a valve in a permanent regime, and that the value of $q_{in}$ relative to those threshold controls if the valve stays closed, open, or opens and closes in cycles when submitted to a constant flux $q_{in}$. The flux above which a closed valve will eventually close is:

$$q_c^{break} = \frac{\rho}{\eta} k_{lo} \frac{\delta p_c^{break}}{w_v}, \qquad (4.6)$$

where $\rho$ and $\eta$ are the density and viscosity of the fluid. The flux below which an open valve will eventually open is:

$$q_c^{clog} = \frac{\rho}{\eta} k_{hi} \frac{\delta p_c^{clog}}{w_v}. \qquad (4.7)$$

In Chapter 3, we demonstrate that in geological media $q_c^{break} < q_c^{clog}$. When $q_c^{break} < q_{in} < q_c^{clog}$, a closed valve will open because $q_{in}$ is higher than the breaking threshold $q_c^{break}$, and as $q_{in}$ is lower than the clogging threshold $q_c^{clog}$, the now open valve will eventually close. In this case, the valve will be permanently unsteady, opening when it is loaded enough by the input flux, and closing when the fluid pressure differential across it has diffused to a low enough value to allow it to close, and repeating this sequence in cycles.

The closer the input flux $q_{in}$ is to either threshold values within the unsteady regime range, the closer the valve will be to a permanently closed or open state. For instance, if the flux is just above $q_c^{break}$, $\delta p$ will very slowly reach the breaking threshold $\delta p_c^{break}$, as the loading speed is comparable with the diffusion speed, allowing fluid to sip through the valve closed permeability at the same time as it slowly loads. When the valve opens, the background flux is so low that the closing threshold is quickly reached. Overall, the valve
<table><caption>**Table 4.1** — Valve parameters for all simulations presented in this chapter</caption><thead><tr><td>Parameter</td><td>Value (scaled)</td></tr></thead><tbody><tr><td>w</td><td>Valve width</td><td>0.01</td></tr><tr><td>k<sub>hi</sub></td><td>Open valve permeability</td><td>1</td></tr><tr><td>k<sub>lo</sub></td><td>Closed valve permeability</td><td>0.05</td></tr><tr><td>δp<sup>break</sup><sub>c</sub></td><td>Threshold p difference for opening</td><td>0.02</td></tr><tr><td>δp<sup>clog</sup><sub>c</sub></td><td>Threshold p difference for closing</td><td>0.0105</td></tr><tr><td>q<sup>break</sup><sub>c</sub></td><td>Flux above which a closed valve can open</td><td>0.1</td></tr><tr><td>q<sup>clog</sup><sub>c</sub></td><td>Flux below which an open valve can close</td><td>1.05</td></tr></tbody></table>

spends a longer time in the closed state than in the open state when $q_{in}$ is just above $q_{c}^{break}$. Conversely, if $q_{in}$ is high, just below $q_{c}^{clog}$, the valve takes a much longer time to close when open, and opens very quickly when closed: it spends more time open overall.

Values for valve parameters used throughout this Chapter are listed in Table 4.1. In order to test denser valve distributions than in Chapter 3, we take a valve width $w_v$ twice smaller. Valve width is a critical parameter, as it governs how quickly fluid can sip through the valve, and how strong the fluid flux near a valve has to be to impose a given $δp$ across 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.

Another dimensionless parameter governing the dynamics of a valve is the ratio of closed to open permeability $k_{lo}/k_{hi}$. As permeability is the only transport property that is variable in space, the value of $k_{lo}/k_{hi}$ describes how quickly fluid pressure diffuses inside a valve relative to outside of it in the high permeability channel. When this ratio is very low, a closed valve behaves like a barrier to fluid diffusion: the effects of the pressure variations outside the valve will dominate the dynamics, and not the diffusion through the valve. In the case where a valve is isolated and the input fluid rate is low and very close to the opening threshold value $q_{c}^{break}$, the pressure profile through around and through the closed valve will have time to slowly evolve towards a stable, low flow of fluid through the valve, before eventually opening if $q_{in} > q_{c}^{break}$. In this case, the diffusion through the valve can be considered an important factor, and the exact value of the ratio matters. However, when several valves are in close proximity, the variations of pressure due to the opening of a neighboring valve will dominate the dynamics. For this reason, we can consider that with $k_{lo}/k_{hi} = 1/20$, closed valves virtu-
ally behave as barriers to the fluid diffusion, and the exact value of this ratio essentially does not influence valve dynamics.

### 4.3.1 The subduction's hydraulics as a valve system

In the following paragraphs, we describe how to interpret the main parameters controlling seismicity in our model in geological terms, as properties of the hydraulic system in the subduction interface: the criticality of valves, that is the value of the input flux $q_{in}$ relative to the critical opening and closing fluxes $q_c^{break}$ and $q_c^{clog}$, and the number valves $N_v$ and their distribution.

Valves should be thought of as segments of permeable channels within the fault interface. Their relatively low permeability could emerge from a local reduction in the aperture of fractures or constriction of the cross section of the permeable channel, reducing the options for the fluid to circulate. Because of this, those segments are prone to clog with the cataclastic fines or colloids that fluid movements carry. The number and proximity of those segments therefore characterizes the heterogeneity or coarseness of the transport properties along the channel in the subduction interface, on a short scale of around 100 m length. As the system is adimensionalized and valves are kept a fixed width in our study, the number of valves $N_v$ represents the proportion of the system that is dynamic, and is thus directly linked with the proportion and scale of the dynamic segments of the permeable system. In reality, the hydraulic structure of the permeable system in the fault interface is at the origin of the heterogeneity of permeability and dynamic properties. The size distribution of impermeable block in the interface, the quality and grain size of the matrix between larger blocks both govern the permeable section and number of available pathways. Those properties are a product of the amount of subducted sediments on top of the oceanic crust, and the fracturation of the upper oceanic crust and overriding plate. More fracturation could for instance result in a larger availability of pathways for the fluid to percolate in or out of the fault interface (Ide, 2010). It would be modelled with a lower number of low-permeability valve segments in the channels along the interface $N_v$. Valve segments are a short scale heterogeneity, due to fracturing or composition of the interface, which are shaped by larger scale (10–100 km) structural properties of the slab and overriding plate.

Such hydraulic properties like the fracture aperture and size distribution of the finer material in the fault zone is also a fundamental control on how prone to clogging and unclogging the permeability is. Theoretical and experimental studies of clogging in model porous media can help understand how the relative values of the input flux $q_{in}$ and thresholds of opening $q_c^{break}$ and closing $q_c^{clog}$ relate to the structure of the permeable system, and the strength
of the fluid source in the fault zone. For a single constriction of diameter $D$, the ratio $D/d$ of permeable aperture to particle size $d$ to permeable aperture governs how permanently clogged, intermittently clogged or fully open the permeable constriction is (Souzy et al., 2020). When the constrictions are smaller than 2.5 times the size of particles ($D/d < 2.5$), the constriction is fully open and permanently clogged. When the constriction is more than 5 times larger than the particle size ($D/d > 5$), the constriction is fully open, and particles flow through it without clogging it. Finally when $2.5 < D/d < 5$, the constriction clogs and unclogs intermittently, continuously transitioning from a mostly clogged state for relatively large particles compared to the constriction ($D/d \approx 2.5$), to a mostly open state when particles are small compared to the constriction ($D/d \approx 5$). The lower and upper values of the $D/d$ ratio are very similar to the threshold values of flux that allow valve to open or close, as they delimit a range of conditions in which the permeability can be dynamically opened or closed. In reality, the particles carried by the flow in the subduction can be micronic to submillimetric cataclastic fines (Han & Hirose, 2012), or colloid aggregates (Amagai et al., 2019), of micronic scales. In Chapter 3 Section 3.6.2, we found a base permeability for the channel of $k_{open} 10^{-14}-10^{-12}m^2$, which corresponds to a fracture aperture of about 10 microns (see Appendix 3.B for details on the calculation). Only the finer particles can act in the valving mechanism, but the scales are arguably right to allow for the permeability to be dynamic.

The length of the clogged segment is also directly linked with the pressure it is necessary to apply to unclog the segment: the longer the clog, the higher the pressure needs to be to unclog it (Jäger et al., 2018). Finally, it has been experimentally shown that the input flux of a particle-laden fluid into a porous medium controls the intermittence of its clogging, displaying clear regimes for which at low flux values, the medium remains clogged, and progressively and intermittently unclogs as the fluid input rate increases, until it is fully open (Bianchi, Thielmann, et al., 2018). This behavior is once again very similar to the regime diagram drawn in Chapter 3, Figure 3.5. In subduction zones, the fluid flux into the dynamic permeable channels is controlled by the dehydration rate (source) and amount of fluid that escapes to the overriding plate (sink). The value of $q_{in}$ relative to the critical fluxes that allow opening and closing of valves in our model is thus a very direct description of those phenomena described in porous media, in which clogging and unclogging are controlled by a competition between the hydraulic stressing rate of the input fluid source, and the weakness and healing capabilities of the low aperture segments, resulting from the relative size and shape of the pathways and fines that can clog them.
## 4.4 Characterizing source interaction in a two-valve system

### 4.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$, expressed as a function of $w_v$ (Fig 4.4a).

In Figure 4.4, we describe a simulation in which a system with two valves at a distance of $d_v = 0.8w_v$ is submitted to a low flux $q_{in} = 0.13$ — close to $q_c^{break}$, the threshold flux above which closed valves will eventually open. Independently of the initial conditions, we observe that the cycles of pressurization and release of both valves — as seen through $\delta p$ the pressure difference across each of them — eventually synchronize (Fig 4.4b). When the permanent regime is reached, the two valves open and close synchronously, with a short delay between the first valve's opening (downdip valve) and the second valve's opening (updip valve). The 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 (Fig 4.4c), we observe that the updip valve is affected with a short delay by the downdip valve opening and the fluid flux pulse that it triggers. $\delta 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 Chapter 3 Section 3.4.2, we have shown that this triggering interaction does not necessarily occur from the downdip valve towards the updip valve, but that in a two-valve system, the effect of a valve opening updip creates a transiently low pressure in front of the downdip valve, that can trigger its opening.

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 slower and weaker the interaction between the valves. In other words, the distance between valves $d_v$ should be an essential control on the interaction strength and synchronization in this two-valve system.

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 $\delta p$ that will make them change state — open or close. The amount of time
<figure><img src="image_61.png" /><figcaption>**Figure 4.4** — Interactions in a two-valve system. (a) A system with two valves (v1 and v2) spaced by $d_v = 0.008 = 0.8w_v$, submitted to a constant flux $q_{in} = 0.13$, and a constant pressure output $p_{out} = 0$. (b) Cycles of fluid pressure accumulation and release pressure are shown with the time-series of the pressure difference $\delta p$ across v1 and v2. Closing and opening events are indicated, and occur when $\delta p$ reaches a pre-defined threshold. (c) As v1 opens first and when v2 is critically stressed, the pressure transient its opening produces (in red) triggers the opening of v2, with a short delay.</figcaption></figure>

$\delta p$ 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_c^{break}$ or $q_c^{clog}$, the closer to thresholds the valves are. In Figure 4.4b and c, the input flux is $q_{in} = 0.13 \approx q_c^{break} = 0.1$ and it is visible that both v1 and v2 tangentially approach the threshold $\delta p$ for breaking. 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.

### 4.4.2 Valve spacing controls activity synchronization

In order to characterize the influence of the interval valve 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 submitted 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 4.5, we display the results of this experiment for three domains of $q_{in}$ which share similar behavior, and for the full range of interval valve distances.

As the system is noise-free and the valves are identical, during the permanent regime of all simulations, the events occur in identical, repeated se-
<figure><img src="image_62.png" /><figcaption>**Figure 4.5** — Effect of valve spacing on valve-to-valve interaction. (a) Two timescales describe a two-event burst: $\delta t_{\text{intra}}$ is the delay between events in the burst, and $\delta t_{\text{inter}}$ is the recurrence delay of the burst. In (b) and (c), we observe how $\delta t_{\text{inter}}$ (no markers on the line) and $\delta t_{\text{intra}}$ (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 $\delta t_0$, the period of activation of a valve submitted to the same flux. For sufficiently large valve spacing, valves do not interact and act as isolated valves. As valves get closer, events get closer, and the bursts are less frequent.</figcaption></figure>
quences (Fig 4.5). Two timescales therefore characterize the activity: $\delta t_{intra}$
the time interval between the first and second opening of the two-event burst,
$\delta 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 4.5b and
c, we show the evolution of these two timescales. The first observation is
that at large distances ($d_v > 5w_v$, Fig 4.5b), bursts — therefore events for each
valves — occur with a recurrence corresponding to an isolated valve: the cy-
cle of either valve is not affected by the other valve. In Figure 4.4b and c, both
$\delta t_{intra}$ and $\delta t_{inter}$ are normalized by $\delta t_0$ the period of the cycle of an isolated,
identical valve submitted to the same flux.

When the distance between the valves gets shorter, the interaction be-
tween valves becomes stronger, producing two effects on the valves’ cycle.
First, in Figure 4.5c, we observe that the closer the valves are (the lower $d_v$),
the shorter $\delta t_{intra}$ is. As valves get closer, the triggering transient travels faster
from one valve to the other, and therefore, the events are closer in time. Sec-
ond, in Figure 4.5b, we observe that as valves get very close ($d_v < 0.5w_v$ for
$q_{in} = 0.11$), $\delta t_{inter}$ increases for lower and lower $d_v$. This means that bursts
recur with a longer and longer delay, and that overall, valves produce less
events per unit of time. As the inter-event time is so short for such close prox-
imity between valves, the two valves together can be considered as a *macro-
valve*. Interestingly, in Figure4.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 in sys-
tems of many strongly-interacting valves in the previous chapter (Chapter 3,
section 3.5.3).

To summarize, a valve interact with its neighbor through pressure tran-
sients, and the effect gets stronger as the valves are closer together. This in-
teraction has two effects: 1) valves synchronize, their activation (opening) is
closer in time the closer the valves are together, 2) the closer the valves, the
longer the period of their cycle, and therefore the longer the recurrence pe-
riod of the two-event bursts. In the extreme case where the valves are almost
touching, the synchronization is so strong that they can be considered as a
twice wider macro-valve, with a longer period of activation than any of the
single valve that it is made of.

Although this experiment was mainly meant to test the interaction strength
based on inter-valve distance, a comparison of the curves corresponding to
each input flux range in Figure 4.5b and c allows to assess the effect of the
$q_{in}$ on the interaction strength. In Fig 4.5b, we see that for low and high $q_{in}$
ranges (red and green curves), which are closer to critical opening and closing
thresholds, the deviation of $\delta t_{inter}$ from the reference cycle period (grey line)
happens for valve spacings that are wider than for the medium flux range
(yellow curve). Indeed, when submitted to a medium flux range, valves are further away from criticality most of the time, and they have to be closer neighbors to interact with a similar intensity as when they are closer to criticality. The same observation can be made for the delay between events in bursts $\delta t_{intra}$: shorter valve spacings are necessary 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.

In reality, a permeable channel in the fault zone should have many of those valves, distributed according to complex spatial distributions. The elementary interaction we describe here will be key to explain how a given distribution of valve, and therefore structure of the permeable channel, is conducive to the emergence of intermittence, clustered seismicity, or even quasi- or perfect periodicity of bursts.

## 4.5 Effects of valve proximity in complex valve systems

### 4.5.1 Statistical description of valve distributions

In order to represent the heterogeneity and dynamic of permeability of a fault channel, we design valve systems with a number of identical valves of width $w_v$, distributed in more or less regular fashion. 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 \left(\frac{d_v}{d_0}\right)^{u-1} \times \exp((-d_v/d_0)^u), \quad (4.8)$$

and its mean

$$\bar{d}_v = d_0 \Gamma \left( 1 + \frac{1}{u} \right). \quad (4.9)$$

where $d_0$ is the scale parameter, $u$ the shape parameter of the distribution, and $\Gamma$ is the gamma function.

We use the Weibull probability distribution because it allows for a range of different distributions within a single and continuous framework. Depending on the value of the shape parameter $u$ in equation 4.8, it can generate valve distributions that go continuously from perfectly regular to highly clustered spatially. Indeed, when $u \to +\infty$, the distribution tends toward a Dirac in $d_0$, for which all valves are separated by $d_v = d_0$: it is a regular distribution. When $u \to 0$, the Weibull distribution tends to a power law distribution $p(d_v) \propto d_v^{-1}$. In this case, valves form clusters, most of the time, they are spaced with very small $d_v$, and rarely, a large spacing is drawn, separating
Parameterizing valve distributions

<figure><img src="image_63.png" /><figcaption>**Figure 4.6** — Parameterizing valve distributions. Valve distributions are created by drawing valve distances $d_v$ using a Weibull distribution. Using the theoretical mean, we can control the valve density, or average valve distance in the domain $\bar{d}_v$. The shape parameter $u$ allows to control how clustered or regularly-spaced valves are in space.</figcaption></figure>

valve clusters. Finally, for $u = 1$, the Weibull distribution reduces to an exponential distribution $p(d_v) = \exp(-d_v/d_0)/d_0$, of mean $\bar{d}_v = d_0$. In this case, the number of valves in a segment of given length follows a Poisson distribution, valves are more or less homogeneously distributed over the system, even if 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 \ge 1$.

We construct a valve distribution by specifying a buffer zone on each sides of the system, common for all distributions, then by drawing inter-valve distances using the chosen distribution ($u$ and chosen $\bar{d}_v$), and then distributing them on the regular space grid, by simply rounding down the interval valve distance to the closest discrete distance possible. We also ensure that the spacing of the last valve to the one before that is not too far from the target interval valve distance, so that all valves do not end up near the input for the most clustered valve distributions. Figure 4.6 describes the density-clustering ($N_v-u$) space of distributions available to us using this technique.

4.5.2 Valve distribution control on activity

In Section 4.4.2, we showed using a simplified system that the interaction between two valves gets stronger and faster as they are closer together. As
<figure><img src="image_64.png" /></figure>

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

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

**Figure 4.7** — Synthetic tremor activity for three different valve systems: (a) a "Poisson" system (Weibull $u = 1$, $N_v = 29$ valves), (b) a clustered distribution ($u = 0.3$, $N_v = 29$), (c) a denser, clustered distribution ($u = 0.3$, $N_v = 43$). 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 event count (per bin of $\delta 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.
valves get closer, events occur with a shorter and shorter delay after one an-
other, and these two-event bursts recur with a longer and longer period. At
this point, we could make the hypothesis that this effect should persist in
a more complex system: the denser the valves in the system, the stronger
and faster the interactions between them, and therefore, the more clustered
in time and long-period the activity should be.

In order to investigate this, we run a basic test using three different valve
systems, submitted to the same input flux $q_{in} = 0.13$ (near the critical flux
allowing valves to open, $q_c^{break} = 0.1$). The results are presented in Figure 4.7.
System a (Fig 4.7a.i, a.ii, a.iii) has $N_v = 29$ valves, an average inter-valve dis-
tance of $\bar{d}_v = 2w_v$, and a disorder parameter of $u = 1$, corresponding to an
exponential distribution of valve distances, it can be considered rather ho-
mogeneous. The activity it produces is clustered in time, occuring in bursts,
with however less variability than the other two systems. Indeed, System b
(Fig 4.7b.i, b.ii, b.iii) shares the same number of valves and average valve
spacing as System a ($N_v = 29$, $\bar{d}_v = 2w_v$), however the higher spatial clustering of valves ($u = 0.3$) creates dense patches of valves in the system, where
valves are locally much closer together than the average $d_v$. The resulting
activity is more intermittent, it proceeds in 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 (Fig 4.7c.i, c.ii, c.iii) is as spatially clustered as System b
($u = 0.3$), but with more valves ($N_v = 43$, $\bar{d}_v = 1w_v$). It is the system where
on average valves are the closest, and also with the most places where locally
valves are very close. The activity it produces is the most clustered of the
three: it proceeds in bursts of two sizes, almost without any activity other-
wise. They occur almost periodically, with a more constant, longer period
than System b.

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*<sub>v</sub>, lower *d*<sub>v</sub>), the more synchronized their activity is: they all activate at the
same time during bursts, and are inactive outside of those episodes. In ad-
dition to it, the activity seems to also display longer and longer recurrence
timescales as the valves are closer and closer locally in the system.

4.5.3 Parameterizing valve proximity

In order to systematically confirm that the closer the valves get locally, the
more they are interacting, and therefore the more clustered and periodic the
activity gets, we design a series of simulations aimed at understanding how
activity is shaped by the valve distribution, representing the heterogeneity of
<figure><img src="image_67.png" /><figcaption>**Figure 4.8** — The number of interacting valves $N_v^{int}$ is a measure of how much valves are close from each other in a valve system, either owing to the density of the system ($N_v$, $\bar{d}_v$), or to the spatial clustering ($u$). As valves get closer, they interact more, and $N_v^{int}$ should allow to capture the consequences in the most straightforward way.</figcaption></figure>

transport properties in the fault channel. We investigate the effects of two parameters describing the valve distribution: the disorder parameter $u$ of the Weibull probability distribution, describing how clustered valves are in space, and the number of valves $N_v$ (bijectively linked to the average inter-valve spacing $\bar{d}_v$), describing how dense the valve distribution is.

We run more than 42,000 simulations, for distributions described by $N_v = 9–71$ ($\bar{d}_v = 9.5w_v–0.2w_v$, 16 values), $u = 0.2–9$ (8 values). For each ($N_v, u$) (128 theoretical distributions), we draw 30 valve systems to average out the effects of specific valve arrangements on the overall behavior. In all systems, valves have identical width, closed permeability and opening/closing $\delta p$ thresholds (Table 4.1). Each system thus defined (3840 total) is submitted to different values of input flux $q_{in} = 0.11–1.04$ (11 values), with a focus on values of the flux close to the thresholds allowing valve opening ($q_c^{break} = 0.1$) and closing ($q_c^{clog} = 1.05$).

As we showed that valves interact more strongly when they are closer, both a higher level of spatial clustering in the system and a higher valve density put more valves at interacting distance. Systems with higher $N_v$ (lower $\bar{d}_v$) and lower $u$ (more clustered systems) should therefore produce a more clustered activity, allowing the emergence of characteristic recurrence scales. In order to highlight this effect, we will represent the results along a single parameter describing the number of interacting valves in a system: $N_v^{int}$, which is 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 two-valve experiments described in Section 4.4: in Figure 4.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_v^{int}$, we capture the effects of the local proximity of valves in
the system, either due to spatial clustering or overall density (see next Sec-
tion). Figure 4.8 shows systems with different $N_v^{int}$, showing 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.

In the next section, we describe the results of those simulations: how a
higher number of interacting valves $N_v^{int}$ and specific values of the input flux
$q_{in}$ in the system allow for the emergence of temporally clustered and periodic
activity.

## 4.6 Emergence and variability of synchronization in complex valve systems

For all simulations, we measure:

* the activity rate, measured by the average event count per valve in the time period,

* 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),

* and the recurrence timescales of clusters in cases when the activity is clustered in time.

The details of how the latter two measures are computed can be found in
Section 4.2. In Figure 4.9, the results are displayed as a function of the number
of valves at interacting distance, $N_v^{int}$. The results displayed in Figure 4.9a, b
and c are for all valve distributions, but limited to a low value of flux $q_{in} =$
0.13, close to the critical value $q_c^{break} = 0.1$. Three examples of activity are
shown to illustrate the trends visible in the Figure (Fig 4.9d, e and f).

### 4.6.1 Patterns of activity clustering with high valve proximity

We first observe that as $N_v^{int}$ increases, more valves interact in the system,
and they produce less events per valve (Fig 4.9a). This is consistent with the
observation made in Figure 4.5, where the period of valve cycles in a two-
valve system is longer and longer as they are closer together, and interactions
make them act synchronously, as larger, macro-valves.
<figure><img src="image_68.png" /><figcaption>Figure 4.9: Caption on next page</figcaption></figure>
**Figure 4.9** — Intensity, clustering and recurrence timescales of activity for varying levels of interaction in complex valve systems. Results for $q_{in} = 0.13$. (a) Average event count per valve, measuring the activity rate, for varying levels of interactions between valves, parametrized by the number of interacting valves $N_v^{int}$. (b) Level of activity clustering measured with the clustering index $c$, as a function of $N_v^{int}$ (c) The detected timescales of cluster recurrence are represented as dots, as a function of the number of interacting valves. The size of the dot 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 dotted line represents the proportion of periodic simulations (see Section 4.6.2 for details) for each value of $N_v^{int}$. In panels (a), (b) and (c), the circled, colored dots show the values of clustering and recurrence for three valve systems with displaying increasingly strong valve interactions, as shown in (d), (e) and (f). The three bottom panels show the valve distributions for three systems (top) and a short period of activity represented as the event count in time for three systems (bottom), submitted to $q_{in} = 0.13$.

As expected, as the number of interacting valves $N_v^{int}$ increases, the activity is more and more clustered in time (Fig 4.9b). When the system is mostly made of interacting valves (high values of $N_v^{int}$), the activity happens quasi-exclusively in bursts ($c = 1$).

## 4.6.2 Valve proximity controls recurrence timescales and periodicity of activity

The results for the measurements of recurrence timescales are shown in Figure 4.9c. For each simulation that displays a high level of clustering ($c > 0.25$), recurrence scales are measured as the average of delays $\Delta T_k$ between clusters of similar size. The ratio of this average to the standard deviation of the measured delays $mean_k(\Delta T_k)/std_k(\Delta T_k)$ gives an estimate of how stable in time the recurrence is, how periodic the activity is. We consider that when $mean_k(\Delta T_k)/std_k(\Delta T_k) > 2$, activity is nearly periodic. The recurrence timescale and its variability are shown with the blue dots, which size scales with the stability ratio. The proportion of periodic simulations (as previously defined) for systems sharing a given $N_v^{int}$ is shown with the dotted gray lines.

As more valves interact, valve systems tend to produce a more periodic activity, with longer and longer periods. When most of the system is made up of interacting valves ($N_v^{int}$ above 20–30), several size of clusters and recurrence timescales are detected in the simulated activity, both short- ($T \approx 3 \times 10^{-2}$) and long-period ($T \approx 3 \times 10^{-1}$).

The last three panels (Fig 4.9d, e and f) show examples of activity in three systems of increasing $N_v^{int}$ (same simulations as in Fig 4.7).

It is noteworthy that the periods of activity bursts ($5 \times 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⁻⁴ scaled time units for an
isolated valve with qᵢₗ = 0.13). Figure 4.9 shows that the collective behavior
of valves built upon their proximity is key to the emergence of those long
periods.

This behavior is governed by the number of interacting valves in the sys-
tem, but also by the connections between strongly interacting patches of valves
in the system. Indeed, it might seem counterintuitive to see that the *x* = 0.5–1
cluster in Figure 4.7b is larger than any clusters in system c in Figure 4.7c,
and that it nevertheless produces a shorter period. This is because as sys-
tem c is denser overall, the smaller valve clusters behave like a large cluster
of macrovalves, and generates a longer periodicity as the whole system acti-
vates.

The simulations presented in Figure 4.9 span the whole possible space
of valve distributions, but are limited to a single value of the input flux,
$q_{in} = 0.13$. However, we obtain identical trends along $N_v^{int}$ for different values
of the input flux, especially around the critical values $q_c^{break}$ and $q_c^{clog}$, where
valves spend most of the time close to breaking or clogging, and are therefore
very sensitive to small variations of pressure, favoring interactions between
them. When the flux values are further from those thresholds, valve systems
produce an overall less clustered, with much shorter and less stable recur-
rence of event clusters. In the next paragraphs, we describe the influence of
the input flux $q_{in}$ on the characteristic activity that systems produce.

4.6.3 Input flux control on criticality and emergence of syn-
chronization

The experiment described in Section 4.5.3 was designed to test the effects of
the heterogeneity of the system in terms of permeability and dynamic struc-
ture. Activity style, however, also depends on the fluid input rate $q_{in}$. For
all systems, we test 11 values of the fluid input rate that allow for perma-
nent valve activity ($q_{in} = 0.11–1.04$), to have a vision of how flux affects the
synchronization of activity in complex valve systems. Flux values sample
densely the flux regions that should show the strongest variation, closer to
thresholds of opening $q_c^{break} = 0.1$ and closing $q_c^{clog} = 1.05$. In Figure 4.10, we
show measures of activity rate, clustering intensity and recurrence timescales
of bursts for systems in which the valve distribution allows a high level of in-
teraction ($N_v^{int} > 20$). We reproduce and generalize the results that we found
for a Poissonian distribution of $N_v = 29$ valves (30 different sampling with
the same parameters), in Chapter 3, Section 3.5.1, Figure 3.10. The closer $q_{in}$
is to critical values of the flux, the closer valve systems are to staying closed
Characteristics of activity for various input flux $N_v^{int} = 20$

<figure><img src="image_69.png" /><figcaption>Figure 4.10 Caption on next page.</figcaption></figure>
**Figure 4.10** — Clustering and recurrence timescales for different values of the input flux $q_{in}$, restricted to highly interacting valve distributions $N_v^{int} > 20$. (b) Level of activity clustering measured with the clustering index c, for different values of $q_{in}$. (c) The detected timescales of cluster recurrence are represented as dots, as a function of $q_{in}$. The size of the dot 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 dotted line represents the proportion of periodic simulations (see text for explanation) for each value of $q_{in}$. In panels (b) and (c), the domains where the flux value prohibits valve activity are shown with red (all valves closed) and green (all valves closed) patches. The circled, colored dots show the values of clustering and recurrence for three values of increasing values of the input flux, as shown in the three bottom panels. (d) shows the activity for a low, near-critical flux $q_{in} = 0.13$, (e) shows the activity for an intermediate, far from critical input flux $q_{in} = 0.57$, (f) shows the activity for a high flux, near critical $q_{in} = 1.02$. The three simulations are performed with a clustered system ($u = 0.3$) with $N_v = 29$ valves ($\bar{d}_v = 2w_v$), the same system displayed in Figures 4.7b and 4.9e.

$(q_{in} \approx q_c^{break})$ or open ($q_{in} \approx q_c^{clog}$): the average event count per valve is therefore lower (Fig 4.10a), the activity more clustered (Fig 4.10b) and more periodic, with longer periods (Fig 4.10c). When valves in the system witness 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, or with shorter periods.

## 4.7 Discussion

### 4.7.1 From observed to simulated tremor activity patterns

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. Bernaudin & Gueydan, 2018; Cruz-Atienza et al., 2018; Luo & Liu, 2019). In this work, we propose a simple representation of the hydraulic processes in the fault based on elementary fault-valve mechanisms, in order to investigate the role that hydraulic properties 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 (Frank et al., 2016; Idehara et al., 2014; Poiata et al., 2021), the long, quasi-periodic recurrence of tremor bursts (Brudzinski & Allen, 2007; Frank et al., 2014; A. Husker et al., 2019; Poiata et al., 2021), and the spatial disparities of those two phenomena along the strike of tremor source regions in subduction zones, at scales of tens to hundreds of kilometers (Brudzinski & Allen, 2007; Poiata et al., 2021). Those
patterns have been interpreted by some authors as resulting from the dynam-
ics of slip in conditions of near-lithostatic fluid pressure and heterogeneity of
frictional properties of the fault interface (e.g. Luo & Liu, 2019; Wech & Crea-
ger, 2011). Our work aims to understand how strong transients of perme-
ability and fluid pressure could be a key element in the origin of intermittent
seismicity patterns, as hinted by the geological record (Angiboust et al., 2015;
Muñoz-Montecinos et al., 2021; Taetz et al., 2018; Tarling et al., 2021).

The model of we developped reproduces such patterns of seismicity thanks
to the interaction between elementary seismic sources, through fluid pres-
sure transients in the permeable fault interface. As a nearby valve activates,
a given valve will be more likely to activate itself, depending on the strength
of the pressure transient that reaches it, and how close it is to activating: how
high and close to the opening threshold the fluid pressure gradient is across
it at the moment when it feels the activation of its neighbor. In a subduction
zone, we assume that many segments of the channels that transport fluid
along the interface behave as valves, and that they should interact. When
those segments are close enough together, they will activate in unison, by
triggering neighbors when activating.

In experiments with complex valve systems, we showed that temporal
clustering of activity and long-timescale periodicities emerge when enough
valves are in close proximity. The temporal clustering and periodicity of ac-
tivity in the model can therefore be understood as resulting from the inter-
action between many neighboring valves that synchronize their activity, be-
having as a whole macro-valve which activity is composed of the collective
and repeated activation of all the valves that constitute it. As more valves are
in close proximity in the system, they form larger and larger spatial clusters
(Fig 4.7b, for x > 0.5), and clusters of clusters (Fig 4.7c). This can either be
due to a high proportion of the system being made of valves (high number
of valves Nv, low inter-valve distance dv), or the valve system being spatially
clustered (parametrized with u, shape parameter of the Weibull distribution).
We defined a proximity index that describes how close valve gets in the sys-
tem: the number of valves at interacting distance Nint. Under the right flux
conditions, valve clusters behave as macro-valves, and we showed that as
the number of interacting valves Nint increases, the valve activity is more
and more temporally clustered, and the observed recurrence delay of activity
bursts gets longer and longer (Figure 4.9). As the system is adimensional-
ized and valves are kept a fixed width, Nint represents the proportion of the
system that is in interaction, and is thus directly linked with the scale of the
heterogeneity of permeability, with the proportion and scale of the dynamic
segments of the permeable system. The spatial scale of the dynamic segment
built out of interacting valves seems to therefore directly control the period-
icity: the larger the macro-valve, the longer its loading and unloading cycle will be, and therefore the longer the recurrence of its activation. We see that it is the scale of the heterogeneity of permeability and dynamic properties of the system that seems to govern the seismic behavior of the system. In our simulations, we also showed that the fluid input rate is an essential control on how temporally clustered and periodic the activity gets in a system. More than the absolute value of the input fluid rate $q_{in}$, it is its value relative to the threshold values of flux that defines the regime the valve is in: permanently closed when the flux is so low that valves cannot unclog $q_{in} < q_c^{break}$, permanently open when the flux is so strong that valves cannot close $q_{in} > q_c^{clog}$, and sustained cycles loading, opening, unloading and closing when the flux allows valves to close and open $q_c^{break} < q_{in} < q_c^{clog}$. When $q_{in}$ is within the regime that allows valve activity, and it is near the threshold values, we showed that valves spend a long time near the critical $\delta p$ that will trigger opening or closing (Figure 4.4), maximizing the chance that a minute change of pressure from a nearby valve activation will activate the given valve. Therefore, the role of the value of the input flux $q_{in}$ in the fault zone relative to valve properties can be seen as a competition between the loading rate and the strength of valves. The closer to threshold valves are, the weaker they are, and therefore the more susceptible they are to be triggered. In summary, values of $q_{in}$ near the threshold values put valves in a sweet spot for them to synchronize, build temporal clustering and periodicities based on interactions. The closer they are to critical conditions, the more sensitive they are to smaller changes of the pressure field around them. In those conditions, the effective size of the macro-valve increases, as the effective interaction distance increases. A more critical hydraulic system in those terms will therefore produce a more clustered and periodic seismicity, with longer recurrence of activity bursts (Figure 4.10).

Tremor activity observed in subduction zones displays high levels of clustering. In the Shikoku tremor zone for instance (Fig 4.1), we measure a clustering index mostly superior to $c = 0.5$, with wide ($\approx 50$ km) regions where $c > 0.8$. Although there is no question that the dynamic of slip on the interface should play a role in shaping the long-wavelength spatial patterns activity, our work provides a simple, alternative and complementary framework to interpret those patterns as resulting from the dynamic and structure of the hydrological system: how strong the fluid pressure sources are relative to the strength of permeability valves, how heterogeneous and dynamic the permeability is in the fault interface. The high levels of clustering $c > 0.5$ in Shikoku can be interpreted in two ways, through the lens of our model:

* The valve system that could model the reactive permeability of the fault zone should allow a high level of interactions $N_v^{int} > 20$, either owing
to its density or spatial clustering of valves, segments of low, dynamic
permeability that can interact together. For the valve widths we showed
in this chapter $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.

* Only near-threshold valves can produce clustering indices on the level measured in Shikoku ($c > 0.5$). The input flux $q_{in}$ should therefore be quite close to the opening threshold $q_c^{break}$ or to closing the threshold $q_c^{clog}$ of valves in the system. For the former, flux is relatively low, the permeability of valves in tremor zone should therefore be nearly shut on average, and opening with slight variations of fluid pressure. In the latter case, the flux is high, the permeability of the system is nearly fully open, and closes at the slightest local drop of fluid flux. It seems less realistic to have a fully connected, high-flow plumbing system at such depths in the fault-zone, where the permeability seems to close instantly after opening (Tarling et al., 2021). We will therefore assume that $q_{in}$ should be low, close to $q_c^{break}$, and valves must be mostly closed, and open only in rare bursts.

In the next sections, we first discuss how those parameters can be linked
with the hydraulic properties of the subduction interface and how the hy-
draulic properties of the subduction interface described by $N_v^{int}$ and $q_{in}$ can
be inherited from the structure of the oceanic plate that is subducting, and
then, how long-wavelength variations of topography of the subducting plate
in Shikoku can result in the segmented seismicity we observe there, by shap-
ing the hydraulic dynamics of the fault on such a scale.

In the rest of this discussion, we will use the link developed previously
in section 4.3.1 to interpret how the hydraulics of the fault zone are shaped
by large-scale heterogeneities of the subducting plate structure.

4.7.2 The hydraulic segmentation of subduction faults

Long-wavelength geophysical observations in subduction zones, of magnetic or gravimetric anomalies for instance, reveal that subduction zones are struc- turaly segmented on a scale of tens to a hundred of kilometers. This segmen- tation seems to emerge mostly from the topographic and internal structure of the subducting oceanic plate (Bassett & Watts, 2015; Blakely et al., 2005; Shillington et al., 2015; K. Wang & Bilek, 2014; Wannamaker et al., 2014), 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 into the subduction interface, and modifies the fluid pressure conditions.
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 along strike in subduction zones, on a very similar scale as the observed topographic features of the incoming seafloor (Bassett & Watts, 2015; Blakely et al., 2005; Shillington et al., 2015; K. Wang & Bilek, 2014). K. Wang and Bilek (2014) and Bassett and Watts (2015) show that the subduction of a topographically rough feature (seamount or ridge) of an oceanic plate causes zones of aseismic slip, or a larger amount of small earthquakes compared to large earthquakes (larger b-values). As the topographic high plows into the overriding plate, it locally scars the fault zone, which becomes more damaged, and weaker frictionally, it also generates high fluid pressure conditions, by eroding and carrying into subduction fluid-filled sediments in the accretionary prism, or simply because of its thicker section of hydrated oceanic crust.

Similarly, the topographic features of the subducting seafloor have a strong impact on the hydraulics of the subduction interface. In a study of the spatial and temporal behavior of tremor in the Shikoku subduction, Ide (2010) interprets the striations formed by areas of high, and frequently recurring tremor activity are due to the scars of subducting seamounts, enhancing fault zone permeability and therefore fluid availability in this zone by channeling more fluid along the interface. Different parameters can impact the flux of fluid channeled through the interface, $q_{in}$. The amount of fluid that is freed and channeled in the interface depends on the importance of the dehydration source, which depends on the volume of crust and sediments that are dehydrating through metamorphic reactions. The thicker the crust, because of a seamount or ridge, the more volume of hydrated crust will dehydrate at depth, increasing $q_{in}$. The local geometry of the plate also changes how important the dehydration rate is deeper along the slab, and can modify $q_{in}$ (Maury et al., 2018). The fluid input rate into channels also depends on the importance of sinks, where the fluid escapes to the overriding plate. It depends on the permeability contrast between the along-dip/along-strike directions in the fault plane and the direction across the fault plane towards the overriding plate. This contrast would be lowered if rough terrain has strongly damaged the overriding plate: more fluid can escape the interface through fractures in the sole of the overriding lower crust or mantle wedge, decreasing in turn the amount of fluid circulating in the channel, and thus the value of $q_{in}$.

Additionally, subducting topographic highs also have a strong impact on the structure of the permeable system, and therefore on the number of valving segments $N_v$, their proximity $N_v^{int}$, and how close to clogging/unclogging
threshold they are, that is how close $q_{in}$ is to the critical fluxes for opening and closing of valves. A subducting seamount might for instance strongly damage the neighboring of the fault, increasing the number of fluid pathways around the faults, therefore lowering the overall number of valves. The damaging and milling process in the fault could also generate more fines, and favor clogging, thus increasing the value of $q_c^{break}$ higher and closer to the local fluid flux $q_{in}$. Finally, the criticality of valves can depend on how much silica colloids can precipitate in the fault zone and clog the permeability, which critically depends on the temperature, and therefore on the local geometry of the subduction (Audet & Bürgmann, 2014). Contrary to the input fluid rate, it is harder to map the actual effects of fracturing, block size distribution and sediment content in the fault zone onto the criticality of valves or their spatial distribution. More in situ and geologic studies of the architecture and dynamics of the fault zone are needed to understand how dynamic the permeability in the interface is, and how it relates to the local structure of the subducting plate. We can however use the patterns of seismicity that we observe in Shikoku to link them with structural features of the subducting plate, in order to bring some insight on how the hydraulic system is structured, and how it shapes seismicity in turns.

### 4.7.3 Hydrological processes shaping tremor intermittence in Shikoku

From Figure 4.1, the Shikoku tremor zone can schematically be divided into three segments based on its activity. Two segments (segment 1 and 3) at each along-strike extremity ($d < 130$ km, and $d > 200$ km, Fig 4.1, Fig 4.11), characterized by a high level of clustering and relatively short timescales of recurrence, and a buffer segment in between, segment 2, characterized by low, almost exclusively clustered activity, with a longer period. In Figure 4.11, we see that the along-strike extent of those segments and spots of high tremor activity within them is quite similar to the size of seamounts on the seafloor of the Philippine plate, or the ones observed in the accretionary prism (Yamazaki & Okamura, 1989), of 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.

In the next paragraphs, we will interpret the characteristic intermittence observed in those segments of tremor zone first as being mainly governed by criticality of the permeability valves in each segment, and then as a symptom of the valve structure of the interface primordially.
### 4.7.3.1 Criticality conditions

We can interpret the segmentation of intermittence in Shikoku as a symptom of how close valves are to threshold. For simplicity, we consider that flow conditions are such that the metamorphic flux is relatively low, just above to the threshold for which valves can simply bear the pressure gradient it imposes on them, and not break: $q_{in} \approx q_c^{break}$, and that transport properties and valve distribution do not vary significantly along-strike.

If the valves are relatively similar across the whole zone, the lower activity, higher temporal clustering and clearer-defined, longer periodicities in the buffer segment (segment 2, $130 \text{ km} < d < 200 \text{ km}$ along-strike) can be due to a relatively lower metamorphic flux in the region compared to the neighboring regions 1 and 3. 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 2 would therefore produce a lower seismicity, and a more clustered, long-period activity than the end segments 1 and 3, as valves would be closer to threshold. In the latter segments ($d < 130 \text{ km}$ and $d > 200 \text{ km}$ along-strike), the permeability valves would be further from threshold because of a higher fluid flow, open and close more frequently, with less possibility of interaction, thus reproducing the higher, less clustered seismicity we observe.

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 channels — does not vary significantly along-strike. The valves in the buffer, segment 2, could indeed be stronger, with a higher fluid pressure gradients $\delta p_c^{break}/w_v$, 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 3 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 ore several subducting seamounts (Ide, 2010), 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 4.11, the alignment of the Kinan seamount chain with the subducting seamount detected by Yamazaki and Okamura (1989) and the tremor zone could indicate that the chain extends into the subduction, and could result in the tremor activity patterns observed there. Segment 2 that produces lower, longer-reccurence 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, this interpreta-
tion 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 (Chao & Obara, 2016; Chao et al., 2013; Kurihara et al., 2018;
Miyazawa et al., 2008). Indeed, if seisimicity in the segments is triggered by
minute strains on the interface, it suggests that those segments might actu-
ally be closer to criticality: valves would be more sensitive to the very small
$\delta p$ changes across them generated by the dynamic strains of tides and tele-
seismic waves.

4.7.3.2 Permeability structure

If we make the simplification that the input flux in fault channels, the back-
ground transport properties and the permeability valve characteristics all do
not significantly change across the strike of the tremor zone, a variation of
valve density and clustering in each segment can still explain the segmenta-
tion of intermittence. As we showed in Section 4.5.3, the closer the valves
locally get in the system — either due to high spatial clustering, or high over-
all density —, the stronger the interactions between them. With strong in-
teractions come a highly clustered activity in time, with well-defined, long
recurrence timescales of clusters of activity. Also, when the interaction be-
tween valves is stronger, the activity rate of each valve is lower than when it
is further away from its neighbors. We can therefore interpret the character-
istic activity in the buffer zone, segment 2, 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 2, but a larger share of
the system being dynamic. This could be due to a more homogeneous perme-
able system, perhaps because of a relatively smooth segment of oceanic plate
being subducted in this region.

In this interpretation, the buffer zone high valve density would behave
as an hydrological barrier for the pressure circulation along strike. Although
the boundary conditions and geometry we explored for the model are tai-
lored to describe along-dip propagations of seismicity, we can draw from the
behavior of the system to imagine how fluid would circulate along-strike, in a
network of channels that would extend in both the dip and strike directions,
in the whole the tremor zone. Circulations of seismicity along-strike would
therefore occur as pressure is transmitted to neighboring segments through
along-strike channels. As it meets segments with a higher valve density, the
pressure would preferentially escape along-dip, arresting the seismicity mi-
gration. In Figure 4.7a, a dense patch of valves around $x = 0.3$ seems to act
such — a barrier between the neighboring parts of the system — and activates
in a more clustered, less frequent manner, letting the fluid to circulate more
intermittently than in the neighboring regions. This could explain why the
buffer zone seems to act as a barrier to large migrations of seismicity in the
neighboring segments: it opens and lets pressure transients circulate through
more rarely, only when a large enough pressure difference along-strike has
built up, depending on the past activity in the neighboring segments.

4.8 Conclusion

In this study, we have shown how the fluid flow and the structure and dy-
namics of permeability in the subduction fault interface can shape the inter-
mittence of tremor activity, and explain the segmentation of temporal clus-
tering and periodicity of activity observed in Shikoku, Cascadia, or Mex-
ico. We propose that the interaction between elementary tremor sources is
mediated by pressure transients, and that their interactions allow temporal
clustering and well-defined timescales to emerge. We built a model of tran-
sient fluid transport processes in the fault that encapsulates this idea. In our
model, events are generated as fluid flow breaks barriers of permeability open
in a porous channel within the fault zone. The transient of fluid pressure
that result from the sudden opening of permeability can trigger neighboring
sources, if they are close enough in space, and close enough to their rupture
conditions. Those two factors are the main controls on how valves interact in
our system, and we show that they crucially shape the intermittence of the
simulated activity: the closer valves are together locally in the system, the
closer they are to their opening or closing thresholds, the stronger their in-
teraction, and the more temporally clustered, periodic and longer period the
simulated activity is. This framework can bring essential insight into how
hydrological processes in the fault zone shape the observed spatial variation
of tremor activity. Our model provides a simple but physically-based mech-
anism to describe the influence of hydrological processes on shaping tremor
patterns. While more work is needed to refine the physical description of
fluid processes in this framework, and understand how they are connected to
deformation and slip in the subduction shear zone, we argue that including
dynamic hydrological processes in our interpretations of tremor, and micro-
seismicity patterns in general, is essential to understand fault zone phenom-
ena.
<figure><img src="image_70.png" /></figure>

**Figure 4.11** — 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 plate, towards the Eurasian plate (fixed). The extent of the map in Figure 4.1a is shown with the white box. We single out rough topographic features on the subducting Philippine plate: the Kinan seamount chain, and a subducting seamount (Yamazaki & Okamura, 1989). A 2D histogram overlaid on Shikoku island shows the distribution of activity in the tremor zone, which can be separated in three zones characterized by different patterns of activity (cf Figure 4.1 for details).
# APPENDICES TO CHAPTER 4

## Contents

<table><tr><td>4.A Tremor patterns in two other subduction zones</td><td>149</td></tr><tr><td>4.B Additional examples of the clustering analysis</td><td>152</td></tr><tr><td>4.C Full distribution space</td><td>155</td></tr></table>

## 4.A Tremor patterns in two other subduction zones

In subductions where the extent of observations along strike is wide enough, the observed tremor activity is segmented along strike, varying both in temporal clustering and recurrence timescale of bursts.
<figure><img src="image_71.png" /></figure>

**Figure 4.A.1** — Segmentation of tremor intermittence in Mexico (Central American subduction zone, 2009–2018). (a) Map of tremor activity in the subduction zone. Colored patches outline the extent of the source regions (Brudzinski et al., 2010; Frank et al., 2014). The stations used to detect daily tremor activity are shown with black triangles. Depth contours of the slab from Slab2 (Hayes, 2018). The thick black line outlines the trench. Hours of tremor per day in (b) Guerrero, measured on station ARIG, and (c) Oaxaca, measured on station YOIG. Tremor data from A. Husker et al. (2019). On both (b) and (c) the periods shaded in dark grey represent geologically detected slow slip events (A. Husker et al., 2019). The two zones show different intermittence patterns. In Guerrero, large bursts of tremor occur quasi-periodically every 2–3 months, in between the intense tremor episodes occurring with SSEs every 3.5 years. In Oaxaca, inter-SSE bursts are less periodic, recurring every 1–3 months, in between longer tremor episodes associated with large SSEs, roughly every 2.5 years.
<figure><img src="image_72.png" /></figure>

**Figure 4.A.2** — Segmentation of tremor intermittence in Cascadia (2010–2022). (a) Map of tremor activity in the subduction zone, colored according to source density and segmentation (as defined by Brudzinski and Allen (2007)), catalog from Wech (2010); Wech and Creager (2008). Depth contours of the slab from Slab2 (Hayes, 2018). The thick black line outlines the trench. Tremor counts per day under (b) Vancouver Island and Northern Washington, (c) South Washington and Oregon, (d) South Oregon and Northern California. The periodicity estimates are from Brudzinski and Allen (2007). The activity is segmented along strike: the northernmost segment is most periodic, the center segment is less periodic, with a longer timescale of recurrence and higher level of clustering, and the southernmost segment is least clustered, with a shorter recurrence scale of tremor bursts.
4.B Additional examples of the clustering analysis

In this section, we depict how the clustering analysis characterizes the activity
in two additional cases: on observed LFE activity in Shikoku (Fig 4.B.1), and
on a synthetic case, exhibiting weakly clustered activity (Fig 4.B.2), which
demonstrates the limit of our analysis for weakly clustered activity.
<figure><img src="image_73.png" /><figcaption>Clustering analysis: Shikoku tremor zone, segment 6 [Poiata et al, 2021]</figcaption></figure>

**Figure 4.B.1** — Clustering analysis of LFEs in the segment 6 (Poiata et al., 2021) of the Shikoku tremor zone. (a) Event count time series. The threshold defining peaks is represented as the red line, and corresponds to the 66% quantile of the event count distribution, represented in (b). (c) shows the activity in time and space, projected along the dip of the fault. Colors show how events are labelled into clusters. (d) shows the duration and size distribution of clusters identified for the full simulation. Only one type of cluster is identified by the algorithm in this case.
Clustering analysis: $N_v = 25$, $\bar{d}_v = 2.5w_v$, $u = 0.7$, $q_{in} = 0.13$

<figure><img src="image_74.png" /><figcaption>**Figure 4.B.2** — Clustering analysis of a simulated, weakly clustered activity in a clustered valve system with $N_v = 25$ valves ($\bar{d}_v = 2.5w_v$), $u = 0.7$, under a input flux of $q_{in} = 0.13$. (a) Event count time series. The threshold defining peaks is represented as the red line, and corresponds to the 66% quantile of the event count distribution, represented in (b). (c) shows the activity in time and along the channel, colors show how events are labelled into clusters. (d) shows the duration and size distribution of clusters identified for the full simulation. In this case, only one class of clusters is identified by the algorithm, with a high variability of sizes and durations.</figcaption></figure>
4.C Full distribution space

The next two Figures (Fig 4.C.1 and 4.C.2) show examples of a valve distribu-
tion drawn from a Weibull distribution for all $(u, \bar{d}_v)$ explored.
<figure><img src="image_75.png" /><figcaption>Valve distributions, homogeneous to regular ($u \ge 1$)</figcaption></figure>

**Figure 4.C.1** — An example of valve distribution drawn from every ($u, \bar{d}_v$) explored in the simulations, here for Poissonian to regular valve distributions.
<figure><img src="image_76.png" /><figcaption>Valve distributions, clustered to homogeneous ($u < 1$)</figcaption></figure>

**Figure 4.C.2** — An example of valve distribution drawn from every ($u, \bar{d}_v$) explored in the simulations, here for clustered to Poissonian valve distributions.
<figure><img src="image_77.png" /></figure>
<figure><img src="image_78.png" /></figure>
<figure><img src="image_79.png" /></figure>
Part IV

Conclusion
<figure><img src="image_80.png" /></figure>
# CHAPTER 5

## SCOPE, LIMITATIONS AND PERSPECTIVES OF THIS WORK

### Contents

<table><tr><td>5.1 Summary and scope of this work</td><td>163</td></tr><tr><td>5.2 Confronting modeled and observed processes</td><td>167</td></tr><tr><td>5.2.1 Changes of permeability in the geological record</td><td>167</td></tr><tr><td>5.2.2 Interactions between fluid circulation and fault slip</td><td>169</td></tr><tr><td>5.3 Application to volcanic long-period seismicity</td><td>171</td></tr><tr><td>5.4 Beyond the subduction zone</td><td>174</td></tr></table>

To conclude this manuscript, we will first summarize key takeaways from our work. We will then discuss how realistically the physics we base our model of unsteady fluid circulation represents evidence of permeability changes in the fault zone, and how to think about the relationship of the described processes here with slip and deformation in the plate interface. Finally, we will describe how our model could be applied to understand spatio-temporal patterns of volcanic tremor as symptom of unsteady magmatic circulation in volcanoes, before concluding on the more general scope of this model, and what it could contribute to our understanding of geological plumbing systems.

### 5.1 Summary and scope of this work

[1] In Chapter 2, we show that LFEs in Guerrero, Mexico, seem to have a corner frequency, and therefore duration, independent from their seismic mo-
ment, with a characteristic value just between 2–3 Hz. This behavior is con-
sistent with a source process controlled by the timescale of fluid diffusion
in the fault interface. The sudden opening of a permeability barrier in the
fault could trigger a fast fluid pressure transient responsible for generating a
seismic signal consistent with of LFEs’ amplitudes and radiation patterns in
subduction zones.

[2] In Chapter 3, we focus on understanding how intermittent, migrating pat-
terns of LFE activity could emerge from fluid pressure transients in the sub-
duction interface. We build a model of circulation of fluid, channelized in a
permeable conduit along the dip of the fault zone, fed by a constant metamor-
phic source. Based on observations that fluid pressure and permeability are
heterogeneous, vary strongly and cyclically in the fault, we model the perme-
ability of the fault channel in a binary fashion: high-permeability segments
and elementary low-permeability segments. The low permeability segments
behave as valves, opening and closing in response to the local fluid pressure
gradient. This behavior is based on the physics of clogging and unclogging
of a permeable medium conducting a pressurized, particle-laden fluid.

[3] Each sudden opening of a valve segment of the permeable channel gener-
ates a strong fluid pressure transient as fluid rushes through the open valve,
suddenly decreasing the force applied to the matrix by the impeded flow.
This strongly, rapidly varying force can act as a source of seismic waves, and
its radiation pattern is both consistent with the radiation pattern of LFEs, and
nearly indistinguishable from a shear mechanism on the plate interface, with
the geometry of observation of LFEs in subduction zones. Through this ex-
plicit association in our description, the fluid circulation model is fully cou-
pled with processes generating seismicity. The actual process generating the
seismic event is unknown and still debated, and a small shear rupture trigger-
gered by the changing effective normal stress on the fault could also be con-
sistent with our description. However, our focus is to model the interactions
between sources, and we leave this subject open for discussion.

[4] We find that in geological conditions, when submitted to a constant flux,
permeability valves can open and close in cycles if the flux is high enough
that valves can be broken, but low enough that they can heal closed after be-
ing broken. Two critical values of flux delimit this operating range for valving
and seismic activity: the flux value above which a closed valve will eventu-
ally open, and the flux value below which an open valve will necessarily
close. When any valve system is submitted to a fluid flux in this range, it
generates intermittent pressure, flow and permeability conditions in the fault
zone. As seismicity is coupled with the opening of permeability valves, the
intermittent loading and unloading of fluid pressure in the hydraulic sys-
tem produces intermittent seismic behavior. Under a constant source of fluid,
a complex, intermittent behavior of fluid circulation emerges, shaping sus-
tained, permanently intermittent patterns of seismic activity.

[5] As closed valves open when the pressure gradient exceeds a certain thresh-
old across their spatial extent, the flow transient associated with the opening
of a neighboring valve can trigger their opening. This triggering interaction
can act in both up and downdip directions, by sucking the fluid in front of
them if the opening neighbor is in front of the valve, or by pushing more
fluid against it, if the opening neighbor is behind it. Valves, and thus seismic
sources, therefore interact through the fluid pressure field, and can trigger
each other in cascades. This produces fast migrations of activity in the chan-
nel, both in the up and downdip directions.

[6] In a system with randomly distributed valves, complex patterns of activ-
ity can emerge, reproducing characteristic behavior observed in tremor and
LFE activity in subduction zones. As valves interact, their activity synchro-
nizes, and intermittent and quasi-periodic activity can occur. Fast migrations
of activity occur as valves cascade open, reproducing fast tremor bursts. From
the synchronization of their cycles of fluid pressure loading and unloading, a
large-scale valving behavior emerges at the scale of the system. The perme-
ability in the system episodically opens and closes, associated with episodes
of high activity, thereby reproducing the fault-valve behavior observed on
long timescales at the scale of subduction zones. Building on elementary cas-
cades of activity and synchronous cycles of fluid accumulation at the scale of
the system, large scale migrations of activity are produced, and slowly prop-
agate across the whole channel.

[7] In Chapter 4, we explore how the input fluid flux and distribution of the valves in the system control the temporal clustering and timescales of recurrence of activity bursts in our model subduction channel. We start by demonstrating that the strength of interactions in a two valve system depends on the distance between the valves and the flux they are submitted to. The former controls the delay and strength of the triggering pressure transient that travels from the opening valve to its neighbor. The latter controls how close to breaking or critical the neighbor is in average, when it receives the triggering transient from the opening valve. The closer and more critical the valves are, the stronger the interactions are, and therefore the stronger their synchronization.

[8] In a systematic study of how valve distributions affect the simulated seismic activity, we show that the number of valves at interacting distance in the system is a key factor shaping the intermittence of activity. The more valves are interacting in the system, the more clustered and periodic the activity gets, and the longer the observed timescales of recurrence of the bursts. The quantity of interacting valves describes how dynamic and heterogeneous the fault
is, in terms of its permeability. It is a fundamental characteristic of the sub-
duction channel, and should be considered a major factor shaping the activity
of LFE and tremor seismicity.

[9] The same result is obtained when flux is set to values near the critical values bounding the operating range for seismic activity. Under such flux conditions, valves spend most of their time closed, and very close to breaking open when the flux is low, or open, near closing conditions when the flux is high. In those cases, valve can change state when they feel minute changes of pressure, and are therefore strongly interacting with their neighbors. Their high criticality produces highly clustered, periodic activity, recurring on long timescales. The input fluid rate into the fault zone is controlled by the kinetic of the dehydration in the subducting oceanic plate. Other controls include channeling effects that concentrate fluid flow as a result of impermeable sealing of the plate interface, or on the contrary, pervasive damaging of the overriding crust allowing fluid to escape the fault zone. We argue that the relationship between sources and sinks of fluid in the fault zone is a second key factor to take into account when understanding how fluid flow processes could shape tremor activity.

[10] A central result of this study is the emergence of long periodicities of activity in a system where isolated sources recur orders of magnitude more frequently. The emerging collective behavior of valves is due to their collective effect on the pressure field. As they interact, they synchronize, and the cumulated effect of their closing shapes cyclic accumulation and releases of pressure on a large scale within the channel. The buildup and release of pressure on such spatial scales is long, as large portions of the fault alternatively fill up and empty, generating the long-term intermittence that we observe when valves interact strongly.

[11] The spatial variations of source and sinks of fluids, the channeling ef-
ficiency, dynamics and heterogeneity of the permeable channel can help in-
terpret the segmentation of seismicity in subduction zones as a result of the
characteristic of fluid circulation in each segment. Segments where the activ-
ity is most clustered could for instance be the theatre of clearly intermittent
fluid circulation, due to a low fluid input rate into the fault zone, or high
number of valves, that collectively block the fluid flux in the fault, sealing
the interface and generating high fluid pressure downdip of the tremor zone,
and then collectively breaking, allowing the fluid to episodically escape and
producing episodes of high low-frequency seismicity.

[12] The characteristic lengthscale of the segmentation seems to match the size of geological features of the subducting plate, that could influence the hydraulics of the system when such structures reach the depth of the tremor-genic zone. A subducting seamount could cause strong damaging of the fault
and overriding plate, resulting in a lower channeling of fluid in the interface, or in weaker, less critical valves. A locally thickened subducted crust could for instance increase the amount of dehydration, and therefore increase the fluid flow in the fault, generating less critical conditions, thus resulting in a less intermittent seismicity.

The model we have built provides a simple framework to understand how from small scale processes of fluid pressure cycling generated by dynamic permeability, fluid pressure transients can emerge in a subduction channel on a much longer and larger scale, and shape the activity of tremor. In our description, the smaller scale hydraulic processes are directly linked with the activity of elementary sources of tremor, and the interactions between these sources through the pressure field allow them to synchronize. It links in one framework the source of LFEs and their intermittent, migrating activity patterns. We believe this work can help understand better how the complexity of fluid flow processes in the subduction — especially small-scale valving behavior — can shape the activity of tremor, from source to swarm. Although our model only depicts simplified fluid flow mechanisms, it can be used to complexify our understanding of fluid flow in subduction zones by including more dynamic and heterogeneous phenomena, and diagnose their activity using the characteristics of intermittence of tremor or LFE activity in a subduction zone.

## 5.2 Confronting modeled and observed processes

### 5.2.1 Changes of permeability in the geological record

The valve mechanism we design for opening and clogging permeability in the fault is based on a simplified model of how pores clog and unclog in a permeable medium channeling a particle-laden fluid. The conditionally intermittent flow through the valve we obtain in our model (Chapter 3, Section 3.4.1, Figures 3.5 and 3.9) is very similar to the one observed in simulations and experiments of filtration of a fluid carrying particles (Bianchi, Thielmann, et al., 2018; Jäger et al., 2017b; Souzy et al., 2020). However, this mechanism has only been observed in the lab (Candela et al., 2014), in shallow crustal faults (Manga et al., 2012) and oil reservoirs (Sahimi et al., 1999), and there exists no direct evidence of its existence deep in the subduction fault zone where tremor occurs. One of the only ways to find evidence for this process would be to find its signature in the geological record, but such short-lived, transient processes would most surely leave no trace in exhumed subduction channel, either because they are so short and transient that their signature would be quickly overprinted by more important processes, or sim-
ply because they do not occur. In this section, we discuss how relevant this
mechanism and its physical description is to describe the dynamics of perme-
ability in the fault zone, by comparing it to the geological record of processes
associated with fluid flow and tremor in exhumed subduction channels.

Our model of fluid circulation in the fault assumes that fluid is channeled along a permeable conduit in the fault, which is segmented between pockets of highly connected fluids and low permeability barriers in between these pockets. The high permeability segments that form the pockets of fluid can connect when barriers to the flow break, letting the fluid circulate, channeled through the fault zone. The geological evidence seems quite unambiguous that fluid should be channeled from deep and distributed metamorphic sources along the plate interface (Angiboust et al., 2014; Piccoli et al., 2021; Saffer, 2015). In fossilized fault zones, there even is evidence of episodic fluid flow when veins suddenly interconnect domains of overpressurized fluid, forming highly connected networks that channel fluid from external sources (Taetz et al., 2018). In this system, veins would eventually be cemented and permeability would shut down, remodeling the hydraulic network. It is not obvious that as in our model, permeability would close down and reopen in cycles in the same spot, and that valves should remain stable in time.

Evidence of metastomatism around the larger veins indicates that they
should remain open for periods of months. In a network of such veins, cat-
aclastic fines or silica colloids could be transported by the flow and clog the
permeability, reproducing a valving behavior similar to the one we model.
Even if this mechanism really occured, it would occur in micrometer to mil-
limeter wide fractures, and the length of the clogs should also be on the same
order of magnitude. In order to generate pressure transients strong and rapid
enough to be responsible for hydrofracturation (Muñoz-Montecinos et al.,
2021) or the generation of an LFE-like seismic signal (N. M. Shapiro et al.,
2018), many strong unclogging events should occur rapidly in cascade in a
wide area. Even with large clogs in wide fractures, it seems hard for a syn-
chronous opening at many valves to bridge the gap between the millimet-
eric scale (at most) of the process, and the width (~ 300 m), surface area (~
1000 m²), and pressure drop (5 MPa) across the valve assumed by N. M. Shapiro
et al. (2018) to account for the seismic characteristics of LFEs. In lab experi-
ments, Bianchi, Thielmann, et al. (2018) have shown that the pressure jumps
detected as the medium clogs and unclogs follow a power-law size-frequency
distribution. Assuming it would be the case in the fault zone, and that this
law is valid for pressure jumps that could generate LFEs, small uncloggings
should occur orders of magnitudes more frequently than large collective un-
cloggings that could generate LFEs. Following this line of reasoning, this
phenomena would have to be a pervasive feature of the hydraulic system of subduction zone to explain the observed cyclic behavior of fluid flow, pressure and seismicity in the subduction zone. Although such processes are perhaps unlikely to be responsible for generating LFEs only through fluid pressure transients, their presence in the fault zone would still generate transient, cyclic and migrating fluid pressure as our results show. This might be sufficient to shape complex tremor activity by triggering seismic sources, by transiently lowering the effective stress on near-rupture asperities in the fault.

Exhumed subduction channels tend to show large changes of permeability in the fault zone, which seem to mainly come from mechanisms of hydrofracturation (Muñoz-Montecinos et al., 2021) or vein opening (Platt et al., 2018; Taetz et al., 2018; Tarling et al., 2021). After opening rapidly, the permeability healing due to processes of crystallization in the created fractures (Tarling et al., 2021), or ductile deformation sealing and isolating the fluid-filled, overpressurized cracks (Gratier & Favreau, 2003). The physical description of our model arguably approximates such behavior. Hydrofracturation occurs because of a strong, localized overpressure that overcomes the cohesion of the material in the fault, creating fractures for the fluid to circulate. The opening of valves in our model is quite close to this behavior, as it is controlled by the pressure gradient across low-permeability barriers. However, the mechanical and chemical sealing processes should be mainly pressure and time dependent, and depend only weakly on the pressure gradient. A model that described opening and sealing of elementary fault valves in such a way would have much in common with the physical description we make of the interaction between fluid pressure and permeability. It should reproduce strong and cyclic fluid pressure transients as hydrofracturation and sealing alternate, cascades of permeability openings as pressure diffuses in the overpressurized fault zone allowing near-breaking areas to interact, and be potentially more suited to generate single-force LFEs, because of the larger space and pressure scales involved. A central concept worth exploring would be the role of cohesion and damage in the fault, as it might provide a way to link the mechanics of failure and healing of the rock and the evolution of its permeability. Such a description might even be suitable for a more continuous, less binary description of the system's dynamics, with continuously evolving, continuously defined permeability in the fault, allowing to explore how the hydraulic structure of the system could evolve in time and space.

### 5.2.2 Interactions between fluid circulation and fault slip

The simple permeability valve model we have developed in this work was aimed at exploring how transient fluid pressure and permeability in the fault
zone could interact, and how they could affect and shape tremor activity. Because our goal is to understand the role of the hydraulic system in a simple system, fault slip is not at all present in the model, not even implicitly, although it seems clear that it should at least be partly implicated in generating tremor. Fault slip and LFE activity are intimately linked, since they occur in the same places at the same time (e.g. Rogers, 2003), and LFE activity can even be used to detect the geodetic signature of small slow slip transients otherwise buried in the noise (Frank, Radiguet, et al., 2015). In this interpretation, tremor is generated as brittle asperities rupture when a slip transient occurs on the subduction fault (Ando et al., 2010). The seismic characteristics of LFEs — their radiation pattern mostly — are also consistent with shear slip on the subduction interface (Ide, Shelly, & Beroza, 2007), although observations interpreted that way are sometimes ambiguous, and could better fit a single-force mechanism (Wech & Creager, 2007), as noted by Ohmi and Obara (2002).

The processes of unclogging — or hydrofracturation 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 the fault zone, allowing extension veins to open and pump fluid into the fault (Kotowski & Behr, 2019), thus strongly modifying the local permeability and fluid pressure in the slipping region. Damaging and dilantancy in and around a slipping portion of the fault also create permeability (Im et al., 2019b; Mitchell & Faulkner, 2008; Tenthorey & Cox, 2006), 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 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 proximity 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 parts of the subduction zone. The collective behavior of valves creates a cumulated effect on fluid pressure that translate into multi-scale fluid pressure increases and drops, that could shape the behavior of fault slip from the source to the subduction scale. Coupling a description of the dynamic hydraulic processes we describe with slip would however require to include a rock mechanics framework, to depict how slip and unclogging/hydrofracturation affect the cohesion, friction and permeability in the fault, and how fluid pressure affects slip in the fault in return. Although it is still an open research question to understand exactly how fluid transport processes and slip on a fault interact, a coupled
model could for instance help assess the time and space scale at which each
process could play the predominant role in shaping the activity of tremor.

**5.3 Application to volcanic long-period seismicity**

The first detections of low-frequency earthquakes and tremor in tectonic set-
tings were quite surprising for researchers who first found them, as such seis-
mic signals had only been detected in volcanoes before (Obara, 2002). Active
magmatic systems produce a very diversified seismic activity, which can be
divided into many different classes of events based on characteristic of their
signal and the processes they are thought to be linked with (B. A. Chouet &
Matoza, 2013). One of the main classes — quite similar to earthquakes de-
tected in faults — is composed of events called volcano-tectonic (VT) earth-
quakes. VT events occur from small ruptures as magma propagation changes
stress in the rock around the magmatic system (Roman & Cashman, 2006).
Another main class of signals is the long-period (LP) seismicity, composed of
LP earthquakes and tremor, which are similar to non-volcanic low-frequency
earthquakes and tremor. LP earthquakes and tremor in magmatic systems
have a relatively low characteristic frequency compared to VT events, be-
tween 1–5 Hz usually (B. A. Chouet, 1996). They seem to relate closely with
magma or gass movements in the plumbing system of volcanoes, indicating
transient flow events and rapid pressure systems as fluid is transported in
the system (B. A. Chouet, 1996). Their activity seems to be a consistent in-
dicator of an impeding eruption, or at least of transient fluid movements in
the magmatic systems. Indeed, their activity often migrates back and forth
from potential reservoirs in the crust towards the surface, before an erup-
tion (N. M. Shapiro, Droznin, Droznina, Senyukov, et al., 2017). Models of
their source all necessitate strong changes of pressure in the magmatic fluid.
To list a few, LP signals can be generated by fast flow of magmatic fluid
through constricted pathways of the magmatic systems (Benson et al., 2008;
Honda & Yomogida, 1993; Julian, 1994; Ukawa & Ohtake, 1987), by the res-
onance of fluid filled cracks after a strong in- or outflux of fluid (B. Chouet,
1986), or by the sudden nucleation, inflation and resonance of bubbles of gas
as pressure drops because of accelerating flow in the magmatic system (Mel-
nik et al., 2020b).
Spatio-temporal patterns of long-period tremor in magmatic systems

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

**Figure 5.1** — Migrations of volcanic long-period tremor in the Klyuchevskoy volcanic group (KVG), in Kachatka, Russia. (a) Map view of the KVG, black triangles showing the position of seismic stations. The yellow dotted line represents the location of the profile in panel (b). Map taken from N. M. Shapiro, Droznin, Droznina, Senyukov, et al. (2017). Red and green stars show locations of past eruptive sites. In panel (b), we show the location of detected tremor under the KVG. The grey scale shades indicate spatial density of tremor sources projected on the fault plane indicate as the yellow dashed line in (a). The panel is still from an animated catalog in Journeau et al. (2022), where each dot corresponds to location emitting tremor. Dots fade in time: yellow are most recent, while orange, red and then black indicate events further and further in the past. A downdip migration is therefore visible: starting about half a day before the 24th Dec. 2015 at 20 km depth, and migrating downdip, where it locates around 45–50 km depth on the 24th Dec 2015. Panels (c) and (d) show the spatio-temporal patterns of activity in a time-depth representation. Each black dot represents a tremor detected at depth. The red hues represent the spatial likelihood function plotted in time, which account for the diffuse nature of tremor sources and the uncertainty in location of the source at a given time. The arrows and velocity labels highlight migrating behavior of tremor, in both the upward and downward direction. Both panels are taken from Journeau et al. (2022).
The migrations of LP activity sometimes span the whole extent of the mag-
matic system from deep sources to the surface, traveling in both upward and
downward pulses of relatively constant velocities around a few to a few tens
of kilometers per hour (Journeau et al., 2022; N. M. Shapiro, Droznin, Drozn-
ina, Senyukov, et al., 2017). In Figure 5.1, we display examples of such mi-
grations under the Klyuchevskoy volcanic group. The speeds, directions and
extent of those migrations seem to be inconsistent with the flow of fluid in
the system, which would for one mainly travel upwards. They are there-
fore linked with acceleration or deceleration of the flow in the magmatic sys-
tem, which cause propagating transients of pressure triggering seismic vibra-
tions along the way (B. A. Chouet, 1996; Journeau et al., 2022; N. M. Shapiro,
Droznin, Droznina, Senyukov, et al., 2017). The mechanism of transmission
of pressure is not well established, as many questions remain regarding the
composition and structure of an active magmatic system, and how magmatic
fluids flows within it.

A common representation of the magmatic system beneath volcanoes is
that magma storage occurs in open, long-lived chambers shallow in the crust,
fed by deeper, similar reservoirs through dike-like fractures. Geological ev-
idence, geophysical measurements and modeling thereof seem to point to-
wards a much more distributed system of storage and circulation of magma
and gasses, from deep mantellic sources to volcanic vents at the surface. Mag-
matic systems seem to be organized in wide trans-crustal structures, built
from the remnants of repeated injections of magma forming sheet-like struc-
tures (sills), that slowly cool down and are heated back up at each reinjec-
tion (Annen et al., 2006). The combined effects of magmatic fluid advection
and heat transfer forms a complex magmatic "mush" system, in which the
melt, crystalline and gas fraction are heterogeneous across the depth of the
reservoir, forming an intermittently connected permeable system for mag-
matic fluids to flow in (Cashman et al., 2017; Sparks et al., 2019). This sys-
tem is also dynamically remodeled, through processes of injection pressuriz-
ing melt-rich portion of the systems, which can then break low-melt-fraction
mush seals, triggering sudden flow transient that ascends through the sys-
tem (Sparks et al., 2019).

The heterogeneity and dynamics of this system could be modeled very
directly in the framework we developed in this thesis, to describe how in-
termittent magmatic transport emerges in this system, generating punctual
eruptions, whereas it is fed by an *a priori* constant source of magma deep in
the mantle. LP seismicity in volcanoes is clearly linked with rapid flow ac-
celeration in the magmatic system, which seems to be linked with sudden
opening of permeability barriers in the description of Sparks et al. (2019).
Our malve model of the permeability changes, how they relate to changes of
pressure, and their association with the activation of a seismic source, might therefore be a simple and powerful way to represent elementary flow processes and their relationship to seismicity in the volcanic plumbing system. It could provide a direct interpretation of seismicity migrations as acceleration and decelerations in the flow of magmatic fluid throughout the system. The seismicity can thus be used to try and delineate the active parts of the system, to understand better the structure and dynamics of the magmatic plumbing system and how it conducts fluid towards the surface.

More precisely, even, our model provides an interpretation for how clustered in time the seismicity gets depending on the strength of the fluid input into the permeable system. Because it changes the average criticality of valves in time, a change of input flux in the system changes the strength of interaction between valves, and therefore their clustering and recurrence behavior. For a progressive reactivation of an initially inactive system, modeled by a vertical channel with valves distributed within it, the progressive increase of flux from low to high values would trigger first a low, extremely clustered seismicity as the system opens intermittently, remaining mainly closed. If the system was submitted to an increasing flux, it would further open, due to the increase of pressure gradient in the fluid, and an increasingly constant, less clustered activity would occur. Such patterns are observed before eruptions in one of the largest volcanic group on the planet, the Klyuchevskoy Volcanic Group (KVG) (Frank et al., 2018). In such cases, our model could provide a straightforward interpretation in terms of fluid flux variations in each part of the system, and — provided it is correctly scaled — provide estimates of the mass of fluid transported throughout the system, or transport properties of the latter, opening promising alleys for volcano monitoring.

Even though these prospects are encouraging, many questions have to be addressed before directly applying our model to volcanoes, in particular regarding the rheology of the fluid and permeability dynamics. The nature and rheology of the fluid (gas or melt?), the strong thermal effects on fluid viscosity, melt-fraction and therefore permeability, the dependence of fluid flow processes on depth or the stage of maturation of the plumbing system are examples of specificities of the magmatic case that need thorough consideration, in order to estimate how applicable and precise the valve system description can be to interpret patterns of seismicity as symptoms of fluid transport processes in volcanoes.

## 5.4 Beyond the subduction zone

Even though the model we presented is a simplified and non-exhaustive description of the dynamics of permeability and its connection with the source
processes of LFEs and tremor in the subduction zone, it encapsulates key and
realistic aspects of how fluid flow seems to occur in the subduction fault and
shape seismicity, as shown by geological and geophysical observations:

* Fluid pressure and its gradient should have a strong effect on permeability.

* The opening of permeability should be associated with — caused by, or causing — the generation of seismic waves.

* Neighboring sources should therefore interact through rapid pressure transients in the permeable subduction channel.

It provides a useful conceptual framework to understand how heterogeneous,
intermittent fluid pressure conditions can emerge in a fault. Our model of
how the seismic source is triggered by the fluid pressure field and affects it
in return allows to conceptualize how elementary sources can interact in an
hydraulic system, and build up patterns of seismicity.

Beyond the subduction zone, this framework has a broader reach to inter-
pret how fluid unsteadily flows in geologic plumbing systems, magmatic (Journeau
et al., 2022) and hydrothermal systems (Wech et al., 2020), permeable reser-
voirs (Sahimi et al., 1999), faults (Ross et al., 2020), using the seismicity they
generate and the hydraulic parameters that can be observed. The permeabil-
ity valve mechanism is a powerful yet simple description that lets one think
of fluid transport as a series of fluid pressure accumulation and release, that
shape intermittent transport from a generally constant fluid source towards
sinks at the surface, or less active, more diffuse reservoirs of fluid. Transient
patterns of fluid flow in geologic systems can help understand seismicity pat-
terns as symptoms of intermittent transport properties and fluid flow in the
fault that connect seismic sources, controlling their activation and interaction
in time and space. Regardless of the exact processes that trigger the activa-
tion of sources or exactly how permeability barriers break, it is a promising
framework to explain the role of unsteady fluid circulation processes in shap-
ing the dynamics of a wide range of systems (Gosselin et al., 2020; Materna et
al., 2019; Ross et al., 2020; Wech et al., 2020).
<figure><img src="image_82.png" /></figure>
REFERENCES

Abercrombie, R. E. (1995, December). Earthquake source scaling relationships from -1 to 5 $M_L$ using seismograms recorded at 2.5-km depth. *Journal of Geophysical Research: Solid Earth, 100*(B12), 24015-24036. https://doi.org/10.1029/95JB02397
53, 56, 58

Agbangla, G. C., Bacchin, P., & Climent, E. (2014). Collective dynamics of flowing colloids during pore clogging. *Soft Matter, 10*(33), 6303-6315.
https://doi.org/10.1039/C4SM00869C
74

Ague, J. J. (2014, January). 4.6 - Fluid Flow in the Deep Crust. In H. D. Holland & K. K. Turekian (Eds.), *Treatise on Geochemistry (Second Edition)* (pp. 203–247). Oxford: Elsevier. https://doi.org/10.1016/B978-0-08-095975-7.00306-5
95

Aki, K., & Richards. (2002). *Quantitative Seismology* (Second Edition ed., Vol. 1). University Science Books.
32

Allmann, B. P., & Shearer, P. M. (2009, January). Global variations of stress drop for moderate to large earthquakes: GLOBAL STRESS DROP VARIATIONS. *Journal of Geophysical Research: Solid Earth, 114*(B1). https://doi.org/10.1029/2008JB005821
26, 35, 36, 42, 43, 57

Amagai, T., Okamoto, A., Niibe, T., Hirano, N., Motomiya, K., & Tsuchiya, N. (2019, December). Silica nanoparticles produced by explosive flash vaporization during earthquakes. *Scientific Reports*, *9*(1), 9738.
https://doi.org/10.1038/s41598-019-46320-7
18, 74, 123

Anderson, J. G., Das, S., Boatwright, J., & Scholz, C. H. (1986). Implication of attenuation for studies of the earthquake source. *Earthquake Source Mechanics*, *37*, 311–318.
56

Anderson, J. G., & Hough, S. E. (1984). A model for the shape of the Fourier amplitude spectrum of acceleration at high frequencies. *Bulletin of the Seismological Society of America*, 74(5), 1969–1993.
53, 57

Anderson, R. N., Uyeda, S., & Miyashiro, A. (1976, February). Geophysical and Geochemical Constraints at Converging Plate Boundaries–Part I: Dehydration in the Downgoing Slab. *Geophysical Journal International*, 44(2), 333–357.
https://doi.org/10.1111/j.1365-246X.1976.tb03660.x
108

Ando, R., Nakata, R., & Hori, T. (2010). A slip pulse model with fault heterogeneity for low-frequency earthquakes and tremor along plate interfaces. *Geophysical Research Letters*, 37(10). https://doi.org/10.1029/2010GL043056
4, 8, 14, 170
Angiboust, S., Kirsch, J., Oncken, O., Glodny, J., Monié, P., & Rybacki, E. (2015, June). Probing the transition between seismically coupled and decoupled segments along an ancient subduction interface. *Geochemistry, Geophysics, Geosystems*, 16(6), 1905–1922. https://doi.org/10.1002/2015GC005776
16, 65, 66, 95, 108, 112, 139

Angiboust, S., Pettke, T., De Hoog, J. C. M., Caron, B., & Oncken, O. (2014, May). Channelized Fluid Flow and Eclogite-facies Metasomatism along the Subduction Shear Zone. *Journal of Petrology*, 55(5), 883–916.
https://doi.org/10.1093/petrology/egu010
66, 168

Angiboust, S., Wolf, S., Burov, E., Agard, P., & Yamato, P. (2012). Effect of fluid circulation on subduction interface tectonic processes: Insights from thermo-mechanical numerical modelling. *Earth and Planetary Science Letters*, 357–358, 238–248.
https://doi.org/10.1016/j.epsl.2012.09.012
45

Annen, C., Blundy, J. D., & Sparks, R. S. J. (2006, March). The Genesis of Intermediate and Silicic Magmas in Deep Crustal Hot Zones. *Journal of Petrology*, 47(3), 505–539.
https://doi.org/10.1093/petrology/egi084
173

Aochi, H., Poisson, B., Toussaint, R., Rachez, X., & Schmittbuhl, J. (2013, December). Self-induced seismicity due to fluid circulation along faults. *Geophysical Journal International*, 196(3), 1544–1563. https://doi.org/10.1093/gji/ggt356
67

Archuleta, R. J., Cranswick, E., Mueller, C., & Spudich, P. (1982, June). Source parameters of the 1980 Mammoth Lakes, California, earthquake sequence. *Journal of Geophysical Research: Solid Earth*, 87(B6), 4595–4607. https://doi.org/10.1029/JB087iB06p04595
56

Aristide, A. (2020, June). *Modeling Fluid Migration in Subduction Zones*. http://eos.org/science-updates/modeling-fluid-migration-in-subduction-zones.
16

Atkinson, G. M. (1995, October). Attenuation and source parameters of earthquakes in the Cascadia region. *Bulletin of the Seismological Society of America*, 85(5), 1327–1342.
https://doi.org/10.1785/BSSA0850051327
54

Audet, P., & Bürgmann, R. (2014, June). Possible control of subduction zone slow-earthquake periodicity by silica enrichment. *Nature*, 510(7505), 389–392.
https://doi.org/10.1038/nature13391
112, 143

Audet, P., & Kim, Y. (2016, February). Teleseismic constraints on the geological environment of deep episodic slow earthquakes in subduction zone forearcs: A review. *Tectonophysics*, 670, 1–15. https://doi.org/10.1016/j.tecto.2016.01.005
7, 64, 65

Auradou, H., Drazer, G., Hulin, J. P., & Koplik, J. (2005, September). Permeability anisotropy induced by the shear displacement of rough fracture walls: PERMEABILITY ANISOTROPY. *Water Resources Research*, 41(9).
https://doi.org/10.1029/2005WR003938
67

Baltay, A. S., & Beroza, G. C. (2013, December). Ground-motion prediction from tremor. *Geophysical Research Letters*, 40(24), 6340–6345.
https://doi.org/10.1002/2013GL058506
56

Barbot, S. (2019). Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault. *Tectonophysics*, 768, 228171.
https://doi.org/10.1016/j.tecto.2019.228171
65

Bassett, D., & Watts, A. B. (2015, May). Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief: CRUSTAL STRUCTURE AND SEISMICITY: 1. SUBDUCTING RELIEF. *Geochemistry, Geophysics, Geosystems*, 16(5), 1508–1540. https://doi.org/10.1002/2014GC005684
141, 142

Beaucé, E., Frank, W. B., Paul, A., Campillo, M., & Hilst, R. D. (2019, November). Systematic Detection of Clustered Seismicity Beneath the Southwestern Alps. *Journal of Geophysical Research: Solid Earth*, 124(11), 11531–11548.
https://doi.org/10.1029/2019JB018110
71, 109, 113

Beaucé, E., van der Hilst, R. D., & Campillo, M. (2022, March). *Microseismic Constraints on the State of the North Anatolian Fault Thirteen Years after the 1999 M7.4 Izmit Earthquake* [Preprint]. http://www.essoar.org/doi/10.1002/essoar.10510912.1. Earth and Space Science Open Archive. https://doi.org/10.1002/essoar.10510912.1
109, 113

Beckingham, L. E. (2017). Evaluation of macroscopic porosity-permeability relationships in heterogeneous mineral dissolution and precipitation scenarios. *Water Resources Research*, 53(12), 10217–10230. https://doi.org/10.1002/2017WR021306
74, 77, 81, 97

Bedford, J. R., Moreno, M., Deng, Z., Oncken, O., Schurr, B., John, T., ... Bevis, M. (2020, April). Months-long thousand-kilometre-scale wobbling before great subduction earthquakes. *Nature*, 580(7805), 628–635. https://doi.org/10.1038/s41586-020-2212-1
4

Beeler, N. M., Thomas, A., Bürgmann, R., & Shelly, D. (2013, November). Inferring fault rheology from low-frequency earthquakes on the San Andreas. *Journal of Geophysical Research: Solid Earth*, 118(11), 5976–5990. https://doi.org/10.1002/2013JB010118
26, 45, 65

Behr, W. M., & Bürgmann, R. (2021, March). What's down there? The structures, materials and environment of deep-seated slow slip and tremor. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 379(2193), 20200218.
https://doi.org/10.1098/rsta.2020.0218
16, 66, 108

Ben-Zion, Y. (2012, May). Episodic tremor and slip on a frictional interface with critical zero weakening in elastic solid. *Geophysical Journal International*, 189(2), 1159–1168.
https://doi.org/10.1111/j.1365-246X.2012.05422.x
65

Benson, P. M., Austria, D. C., Gehne, S., Butcher, E., Harnett, C. E., Fazio, M., ... Tomas, R. (2020, December). Laboratory simulations of fluid-induced seismicity, hydraulic fracture, and fluid flow. *Geomechanics for Energy and the Environment*, 24, 100169.
https://doi.org/10.1016/j.gete.2019.100169
9

Benson, P. M., Vinciguerra, S., Meredith, P. G., & Young, R. P. (2008, October). Laboratory Simulation of Volcano Seismicity. *Science*, 322(5899), 249–252.
https://doi.org/10.1126/science.1161927
2, 9, 66, 77, 171

Bernaudin, M., & Gueydan, F. (2018, April). Episodic Tremor and Slip Explained by Fluid-Enhanced Microfracturing and Sealing. *Geophysical Research Letters*, *45*(8), 3471–3480. https://doi.org/10.1029/2018GL077586
66, 138

Beroza, G. C., & Ide, S. (2011, May). Slow Earthquakes and Nonvolcanic Tremor. *Annual Review of Earth and Planetary Sciences*, *39*(1), 271–296.
https://doi.org/10.1146/annurev-earth-040809-152531
25, 64

Bhattacharya, P., & Viesca, R. C. (2019, May). Fluid-induced aseismic fault slip outpaces pore-fluid migration. *Science*, *364*(6439), 464–468.
https://doi.org/10.1126/science.aaw7354
3

Bianchi, F., Thielmann, M., de Arcangelis, L., & Herrmann, H. J. (2018, January). Critical Bursts in Filtration. *Physical Review Letters*, *120*(3), 034503.
https://doi.org/10.1103/PhysRevLett.120.034503
17, 18, 67, 74, 123, 167, 168

Bianchi, F., Wittel, F. K., Thielmann, M., Trtik, P., & Herrmann, H. J. (2018, March). Tomographic Study of Internal Erosion of Particle Flows in Porous Media. *Transport in Porous Media*, *122*(1), 169–184. https://doi.org/10.1007/s11242-017-0996-8
74

Blakely, R. J., Brocher, T. M., & Wells, R. E. (2005). Subduction-zone magnetic anomalies and implications for hydrated forearc mantle. *Geology*, *33*(6), 445.
https://doi.org/10.1130/G21447.1
141, 142

Bocchini, G. M., Martínez-Garzón, P., Harrington, R. M., & Bohnhoff, M. (2021). Does Deep Tectonic Tremor Occur in the Central-Eastern Mediterranean Basin? *Journal of Geophysical Research: Solid Earth*, *126*(1), 2020JB020448.
https://doi.org/10.1029/2020JB020448
96

Bostock, M. G., Royer, A. A., Hearn, E. H., & Peacock, S. M. (2012, November). Low frequency earthquakes below southern Vancouver Island. *Geochemistry, Geophysics, Geosystems*, *13*(11), n/a-n/a. https://doi.org/10.1029/2012GC004391
5, 45, 64

Bostock, M. G., Thomas, A. M., Rubin, A. M., & Christensen, N. I. (2017, January). On corner frequencies, attenuation, and low-frequency earthquakes. *Journal of Geophysical Research: Solid Earth*, *122*(1), 543–557. https://doi.org/10.1002/2016JB013405
7, 45, 58

Bostock, M. G., Thomas, A. M., Savard, G., Chuang, L., & Rubin, A. M. (2015). Magnitudes and moment-duration scaling of low-frequency earthquakes beneath southern Vancouver Island. *Journal of Geophysical Research: Solid Earth*, *120*(9), 6329–6350.
https://doi.org/10.1002/2015JB012195
5, 8, 10, 12, 26, 27, 35, 36, 37, 38, 42, 43, 66, 68

Bouchon, M. (2003). A review of the discrete wavenumber method. In *Seismic motion, lithospheric structures, earthquake and volcanic sources: The keiiti aki volume* (pp. 445–465). Springer.
33

Brodsky, E. E. (2019, May). The importance of studying small earthquakes. *Science*,