# Sounds of geological plumbing systems - How transient fluid circulation processes in faults and volcanoes shape sources and patterns of microseismicity

Gaspard Farge

► To cite this version:

Gaspard Farge. Sounds of geological plumbing systems - How transient fluid circulation processes in faults and volcanoes shape sources and patterns of microseismicity. Geophysics [physics.geo-ph]. Université Paris Cité; Institut de Physique du Globe de Paris (IPGP), 2022. English. NNT: tel-03816151

HAL Id: tel-03816151

https://cnrs.hal.science/tel-03816151

Submitted on 15 Oct 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
UNIVERSITÉ PARIS CITÉ
INSTITUT DE PHYSIQUE DU GLOBE DE PARIS
École doctorale STEP'UP, ED 560
Équipe de Dynamique des Fluides Géologiques

SOUNDS OF GEOLOGICAL PLUMBING
SYSTEMS

How transient fluid circulation processes in faults and volcanoes
shape sources and patterns of microseismicity

GASPARD FARGE

Juillet 2022

Thèse de doctorat de Sciences de la Terre et de l'Environnement, spécialité Géophysique

dirigée par CLAUDE JAUPART et NIKOLAÏ SHAPIRO

Présentée et soutenue publiquement le 7 Septembre 2022, devant un jury composé de:

<table>
  <tr>
    <td>EMILY BRODSKY</td>
    <td>Professeure (UC Santa Cruz)</td>
    <td>Rapportrice</td>
  </tr>
  <tr>
    <td>CAMILLA CATTANIA</td>
    <td>Professeure associée (MIT)</td>
    <td>Examinatrice</td>
  </tr>
  <tr>
    <td>KRISTEL CHANARD</td>
    <td>Chargée de recherche (IGN — UPC, IPGP)</td>
    <td>Examinatrice</td>
  </tr>
  <tr>
    <td>STÉPHANE RONDENAY</td>
    <td>Professeur (Uni. Bergen)</td>
    <td>Rapporteur</td>
  </tr>
  <tr>
    <td>ALEXANDRE SCHUBNEL</td>
    <td>Directeur de recherche (CRNS — ENS)</td>
    <td>Examinateur (Pdt.)</td>
  </tr>
</table>

<table>
  <tr>
    <td>CLAUDE JAUPART</td>
    <td>Professeur (UPC, IPGP)</td>
    <td>Directeur de thèse</td>
  </tr>
  <tr>
    <td>NIKOLAI SHAPIRO</td>
    <td>Directeur de recherche (CNRS — ISTerre)</td>
    <td>Co-directeur de thèse</td>
  </tr>
</table>
<figure><img src="image_2.png" /></figure>
# RÉSUMÉ DE LA THÈSE

L'activité des plus petits séismes détectés dans la Terre est symptomatique des changements infimes qui se produisent sous nos pieds. L'activité microsismique peut notamment être utilisée pour surveiller la façon dont les fluides circulent dans les failles ou les volcans, depuis des sources de fluides en profondeur jusqu'à la surface. Dans les zones de subduction par exemple, des signaux sismiques faibles, intermittents et de basse fréquence appelés tremor, composés de petits événements de basse fréquence (*low-frequency earthquakes*, ou LFEs) sont émis depuis l'interface entre les plaques. Dans ce manuscrit, nous essayons de comprendre comment des processus instables de circulation de fluide dans la zone de faille de subduction peuvent façonner le signal et les motifs spatio-temporels que dessinent le tremor, et comment celui-ci peut donc être utilisé pour diagnostiquer la dynamique du système hydraulique de la faille de subduction.

Dans une première partie, nous caractérisons la source des LFEs à Guerrero, au Mexique, et constatons que leur durée semble indépendante de leur magnitude. Cette observation est incompatible avec la vision classique selon laquelle les séismes devraient être générés par une rupture sur une faille. Elle peut cependant être reproduite par un modèle qui explique la source sismique des LFEs par de rapides changements de pression de fluide dans l'interface de la faille.

La source des LFEs n'est pas le seul indice que de fortes variations de pression de fluide se produisent dans l'interface de subduction. Quand cette zone est exhumée, les affleurements sont marqués par la signature d'ouvertures et de fermetures répétées de la perméabilité de la zone de faille, générant des conditions instables de pression dans le fluide. La dynamique de ce système hydraulique doit contribuer à la complexité des motifs d'activité du tremor, bien qu'elle soit rarement considérée dans les études tentant de les décrire dans leur entièreté.

Dans une deuxième partie, nous élaborons donc un modèle qui permettrait de décrire les processus dynamiques de circulation des fluides dans la faille de subduction, et leur lien avec le tremor. Nous modélisons l'ascension des fluides depuis des sources profondes de déshydratation métamorphique vers la surface, au travers d'un canal perméable dans l'interface de faille. Les segments faiblement perméable du canal entravent le flux, et génère localement de forts gradients de pression. La perméabilité des barrière s'ouvre puis se referme dans un mécanisme de valve, en réponse au gradient local de pression de fluide. Lorsque les valves se rompent, le changement brutal de pression génère un LFE. Les valves voisines communiquent par des transitoires de pression de fluide, qui déclenchent leur ouverture quand une valve à proximité s'ouvre. Cette interaction produit des motifs réalistes de sismicité synthétique, semblables à ceux observés pour le tremor : cascades d'activité, synchronisation et migration de l'activité font apparaissent et se développent sur
plusieurs échelles spatiales et temporelles. Le flux fluide à l'entrée du système con-
trôle la présence et le style d'activité produite par les valves. Nous montrons que
plus les interactions entre valves sont fortes, plus une activité clusterisée dans le
temps apparait, et plus les épisodes d'activité se reproduisent périodiquement. La
force des interactions transportées par le champ de pression fluide est d'autant plus
forte que les valves sont nombreuses et groupées dans le système, et que la valeur
du flux d'entrée les rend critiques — proche des conditions d'ouverture ou de fer-
meture —, et donc susceptible d'interagir. Sur la base de nos résultats, nous émet-
tons l'hypothèse la complexité et des longues périodes de récurrence de l'activité
du tremor pourraient provenir du comportement collectif de valves dans la zone de
faille, qui façonnent des régimes hétérogènes et transitoires de circulation du fluide,
à des échelles spatiales allant de la source à la zone de subduction, et des échelles
temporelles allant de la journée à plusieurs années.

Ce étude pose un cadre conceptuel permettant de comprendre comment des con-
ditions de pression de fluide hétérogènes et intermittentes peuvent émerger dans la
faille de subduction. Notre modèle permet de concevoir l'émergence à grande échelle
de motifs spatio-temporels d'activité grâce à l'interaction entre des sources sismiques
individuelles, en décrivant la manière dont un événement sismique déclenche et est
déclenché par des changements de pression de fluide. Au-delà des zones de sub-
duction, ce cadre de pensée est applicable à d'autres systèmes géologiques de cir-
culation de fluide, comme les réservoirs magmatiques. Le mécanisme de valve de
perméabilité que nous décrivons est une description simplifiée mais efficace, perme-
tant de penser le transport des fluides dans des milieux dynamiques comme étant
composé d'événements d'accumulation et débâcle de fluide. De tels processus façon-
nent l'intermittence du transport de fluide et de l'activité sismique qui en découle,
et leur compréhension nous permet d'utiliser les observations sismiques comme un
diagnostic des conditions hydrauliques en profondeur.

MOTS CLÉS: Subduction, système magmatique, hydrologie, microsis-
micité, tremor, seismes basses fréquence
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>
# SUMMARY OF THE MANUSCRIPT

The activity of the smallest earthquakes that can be detected in the Earth is a symptom of minute changes that occur beneath our feet. Microseismic activity can in particular be used to monitor how fluid flows in geological plumbing systems like faults or volcanoes, from deep sources of fluid to the surface. In subduction zones for instance, the plate interface emits faint, intermittent, low-frequency seismic signals called tremor, composed of small individual low-frequency earthquakes (LFEs). In this manuscript, we try to understand how unsteady fluid circulation processes in the subduction fault zone can shape the source and spatio-temporal patterns of tremor, how tremor can therefore be used to probe dynamic fluid processes at depth.

In a first part, we characterize the source of LFEs in Guerrero, Mexico, and find that their duration seems independent of their magnitude. This observation is inconsistent with the classical view that earthquakes should be generated by a rupture of the rock. It can however be explained by a model of LFE source that involves fluid pressure transients in the fault interface, fast and strong enough to generate seismic waves.

The source of LFEs is not the sole indicator of strong, localized variations of fluid pressure in the fault zone. Geological observations of fossilized subductions indicate that the fault permeability opens and closes episodically, and that the fluid pressure drops and increases as a result. Such widespread unsteady hydraulic conditions in the fault could be a key factor shaping the complexity of spatio-temporal patterns of tremor activity. They are however seldom considered in attempts to describe the processes at the origin of the full spectrum of tremor patterns observed in subduction zones.

In a second part, we therefore design a framework to describe the unsteady processes of fluid circulation in the subduction fault, and how they trigger elementary tremor events. We model fluid ascent from metamorphic dehydration sources through a permeable channel in the fault interface. Low-permeability segments of the channel impede the flux, locally producing strong pressure gradients. In an elementary valving mechanism, the permeability barriers break open and then heal shut in response to the strength of the local fluid pressure gradient. When valves break they produce a simulated LFE, the rapid fluid pressure transient potentially being the source of seismic waves. Valves can trigger each other through pressure transients. Their interaction lead to realistic, tremor-like patterns: cascades, synchronized bursts and migrations of activity, on both short and long time and space scales. We describe how the fluid input rate in the channel allows for activity, and more specifically which conditions best favor the emergence of realistic patterns of activity. We show that stronger interactions, and therefore more clustered, periodic activity can occur when valves are clustered in space, and when flux conditions make
them more critical, and thus more susceptible to interactions. Based on our results, the emergence of large, complex, long-period activity transients could come from the collective behavior of valve like segments in the interface, shaping large scale valving behavior of the fault interface.

This work provides a conceptual framework to understand how heterogeneous, intermittent fluid pressure conditions can emerge in the subduction 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 more general reach to interpret how fluid unsteadily flows in all geologic plumbing systems. The permeability valve mechanism is a powerful yet simple description that allows to think of fluid transport as a series of fluid pressure accumulation and release, that shape intermittent fluid transport, and intermittent seismicity as a symptom.

KEYWORDS: Subduction, magmatic system, hydrology, microseismicity, tremor, low frequency earthquake
<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /></figure>
REMERCIEMENTS

As I am writing those last words in this manuscript — which might be the first,
and even only, words you read —, I realize that the list of all the people I have to
thank for is too long for me to go through here. So I will keep it short and profes-
sional, hopefully.

First I want to very warmly thank my supervisors, Nikolai and Claude, for mak-
ing so much time, being so present and helpful throughout those three years, and
even five years for Nikolai. I want to thank Claude for his kind encouragements
when he saw I was down, and Nikolai for his directivity and friendly realism, which
has helped me regain focus many times I was getting lost trying to follow all the
leads. I want to thank William for the discussions and external point of view he
brought on my work, at a time when I thought I would not regain enthousiasm about
my work, and also for making me feel so welcome in the team on the fourth floor of
the Green's building. I want to thank Camilla too, for accepting to be in my defense
jury of course, but also for the scientific and carreer discussions that made me think
about why to continue further in research, and honestly convinced me. Of course,
I want to thank all members of my jury, for accepting to read my manuscript, and
participating in my defense. I am eager to discuss this work with you.

I think being colleagues is an underrated form of friendship, less engaged than
the classical friendship, but as caring and often much lighter, focusing on simply
making a nice everyday environment for everyone to thrive in. For that I want to
thank everyone in the fantastic Geological Dynamics team, for being the best col-
leagues I could wish for, at lunch, around coffee, growing tomatoes and strawberries
on the balcony, and tasting them for *le goûter*. Sometimes colleagues become the
dearest friends. I want to thank Angèle, Cyril, Émilie, Marie, for being *le sang de la
veine*, vraiment. Farther from Paris, I want to thank Caroline for skating and dancing,
Léonard for walks in the Cambridge night by -10°C, Éric for the Lower Level, and En-
rico and Qingyu for helping me practice my Chinese — I hope you will get that if
you read it, it's been a long time already. Luc has never really been a colleague, but
since we spent so much time on the phone during work hours, he might as well fit in
here. I still hope we will one day collaborate, for the most extraordinary project, on
the creative edge that separates science and art, on our way to find *the truth*. Wow,
I know.

I want to thank my family, my friends for supporting me and everything, but
I suspect most of them will not stumble upon this manuscript. It's ok. If they do
however, I want them to know that I feel very proud of getting there, of graduating
soon and making a job of thinking about how things in the Earth really work, and
that they should know that I love talking to them about it, and that I love them for
listening to me.
<figure><img src="image_7.png" /></figure>
# RÉSUMÉ LONG DE LA THÈSE

L'activité des plus petits séismes détectés dans la Terre est symptomatique des changements infimes qui se produisent sous nos pieds. Dans les dernières années, l'augmentation de la puissance de calcul disponible pour les analyses de données géophysiques a permis une augmentation drastique de la résolution des catalogues de sismicité. Cela a permis de détecter des événements de plus en plus faibles, et ainsi de pouvoir observer les processus lents de libération d'énergie qui se produisent dans les failles et les volcans, qui précédent, mènent ou suivent les événements de libération d'énergie plus catastrophiques qui rythment le cycle de ces systèmes — éruptions et grands séismes. L'activité microsismique peut notamment être utilisée pour surveiller la façon dont les fluides circulent dans les failles ou les volcans, depuis des sources de fluides en profondeur jusqu'à la surface. Dans les zones de faille, la présence de fluides à haute pression a pour effet de favoriser le glissement en lubrifiant la faille, et donc de ruptures sismiques. Les essais de séismes générés par des injections de fluides dans les failles — naturelles ou industrielles — migrent dans l'espace selon des motifs caractéristiques, qui semblent dirigés par la diffusion de la pression au sein du système perméable autour de la zone de faille. Dans les volcans, la migration de phases fluides dans le système magmatique est à la source même d'une partie de la sismicité. Des changements rapides et hétérogènes de pression dans le système perméable peuvent mener à de soudaines accélérations du fluide, qui, en interagissant avec la matrice solide du réservoir, génèrent des ondes sismiques basses fréquence. Ces signaux détectés dans les volcans sont appelés tremor volcanique ou séismes longue période (LP), et migrent rapidement au travers du système magmatique soulignant la structure des conduits actifs et les variations de flux dans le temps et l'espace.

Dans les profondeurs des zones de subductions, on détecte des signaux sismiques similaires: le tremor non-volcanique, un signal long, émergent, de basse fréquence, qui semble composé d'événements impulsifs élémentaires, les séismes basse fréquence (*low-frequency earthquakes* ou LFEs). Ces signaux sont émis depuis l'interface des plaques, et leur activité est caractérisée par une intermittence prononcée. Le tremor se produit en épisodes à structure quasiment fractale, étant eux mêmes composés d'épisodes plus courts séparés de courtes périodes de quiescence. Pendant ces épisodes, l'activité migre sur la faille dans des motifs rappelant la sismicité déclenchée par les accélérations du fluide dans les volcans, en migrations à vitesse constante, rapide pour les épisodes d'activité les plus courts, et plus lent pour les épisodes les plus longs. L'intermittence est caractérisée par une récurrence quasi périodique de l'activité, sur des périodes allant du jours pour les sources individuelles de tremor, aux années pour l'entièreté de la zone sismique. Les caractéristiques de source et les motifs spatio-temporels d'activité du tremor semblent nécessiter l'implication de processus instables de circulation de fluide dans l'interface de faille.

Il existe un corpus grandissant d'observations géologiques et géophysiques qui
semble montrer que les fluides provenant de la déshydratation de la plaque circu-
lent chaotiquement dans la région de production du tremor. Les pressions de fluides
y seraient extrêmes, contrebalançant la pression lithostatique locale, et varieraient
cycliquement, entre épisodes de pressurisation et dépressurisations liés à des ouver-
tures et fermetures soudaines de la perméabilité de la faille, qui alternativement lais-
seraient circuler et bloqueraient le flot. Cette dynamique de la faille semble être ac-
tive à des échelles correspondant à la source sismique (mètres à centaines de mètres),
mais aussi à l'échelle de la zone sismique entière (dizaine à centaines de kilomètres),
et les grands changements de perméabilité, de pression et de flux semblent être cor-
rélés à des épisodes d'intense activité du tremor. Malgré ces observations, la plupart
des efforts de modélisation des processus générant le tremor dans toute sa complex-
ité ne prennent en compte l'hydraulique de la faille que comme un paramètre sta-
tique du modèle, affectant la dynamique du glissement, mais n'étant pas affectée
par elle. Dans ce manuscrit, nous nous attachons donc à comprendre comment des
processus instables de circulation de fluide dans la zone de faille de subduction peu-
vent par eux seuls façonner le signal et les motifs spatio-temporels que dessinent le
tremor, et comment celui-ci peut donc être utilisé pour diagnostiquer la dynamique
du système hydraulique de la faille de subduction.

Dans une première partie, nous caractérisons la source des LFEs à Guerrero, au Mexique. Nous sélectionnons 3498 événements des plus impulsifs et de grande amplitude du catalogue, et déterminons leur moment sismique et fréquence coin par une analyse de leur spectre sur plusieurs stations. Les moments sismiques mesurés correspondent à des magnitudes de moment de M2.2 en moyenne, comprises entre M1.7 et M3.9. Les fréquences coins mesurées sont en moyenne autour de 3 Hz et comprises entre 1 et 6 Hz, plus basses d'un facteur dix environ que les séismes classiques de même magnitude. Nous observons que les fréquences coins ne dépendent pas de la magnitude des événements. Cette observation est incompatible avec la vision classique selon laquelle les séismes devraient être générés par une rupture sur une faille. En effet, pour un séisme qui procède d'une rupture dite auto-similaire, la durée de la rupture augmente plus la magnitude du séisme est grand, et donc la fréquence coin du signal — qui est en théorie proportionnelle à l'inverse de la durée — devrait être décroissante en fonction de la magnitude. Un modèle qui expliquerait la source sismique des LFEs par de rapides changements de pression de fluide dans l'interface de faille pourrait expliquer cette observation. Alors qu'une barrière de basse perméabilité cède dans la faille, le fluide sous pression qui s'est accumulé derrière peut soudainement circuler, produisant un transitoire de pression assez rapide et puissant pour générer un signal similaire à celui des LFEs. L'amplitude — proportionnelle au moment sismique mesuré — de ces événements dépend de l'amplitude du transitoire de pression, et peut donc varier, a priori. Cependant, la durée de la source sismique est contrôlée par la vitesse de diffusion du transitoire de pression. Dans ce modèle, et en accord avec notre observation, les durées et fréquences coins devraient donc avoir une valeur caractéristique dépendant des conditions hydrauliques locales dans la faille, et indépendante de l'amplitude du saut de pression.

La source des LFEs n'est pas le seul indice que de fortes variations de pression
de fluide se produisent dans l'interface de subduction. Quand cette zone est exhumée, les affleurements sont marqués par la signature d'ouvertures et de fermetures répétées de la perméabilité de la zone de faille, générant des conditions instables de pression dans le fluide. La dynamique de ce système hydraulique doit contribuer à la complexité des motifs d'activité du tremor, bien qu'elle soit rarement considérée dans les études tentant de les décrire dans leur entièreté. Dans une deuxième partie, nous élaborons donc un modèle qui permettrait de décrire les processus dynamiques de circulation des fluides dans la faille de subduction, et leur lien avec le tremor. Nous modélisons l'ascension des fluides depuis des sources profondes de déshydratation métamorphique vers la surface, au travers d'un canal perméable dans l'interface de faille. Pour refléter l'hétérogénéité des conditions de transport du fluide dans la faille, la perméabilité du canal est décrite comme étant forte, mais localement, des segments de basse perméabilité agissent comme des barrières au flux de fluide, générant une hétérogénéité du gradient de pression quand le système est soumis à un flux constant. Autour des barrières de basse perméabilité, le gradient de pression est fort. Pour modéliser des processus de débouchage de la perméabilité de l'interface de faille, nous décrétons que si le gradient de pression au travers des barrières de basse perméabilité excède une valeur seuil, la barrière cède, et la perméabilité locale augmente brutalement, simulant un processus d'ouverture. Cette ouverture laisse le flux circuler brutalement, générant un fort transitoire de pression qui peut être à la source d'un LFE. Dans notre description, chaque ouverture de barrière est donc synchrone d'un événement sismique élémentaire. Lorsque le flux au travers de la valve est assez bas, le gradient de pression atteint une deuxième valeur seuil, plus basse que celle permettant l'ouverture, et la perméabilité de la barrière se rebouche. Ce mécanisme d'ouverture et fermeture de la perméabilité de la zone de faille nous permet de décrire l'évolution des conditions de transport du fluide dans la faille de manière intermittente. Chaque barrière se comporte comme une valve, et l'action combinée de nombreuse valve fait émerger des pressions de fluides variables et hétérogènes dans l'interface de faille, comme mis en évidence par de nombreuses observations.

Nous mettons en évidence que l'activité des valves dans ce système est soumis à l'intensité du flux d'entrée. Il existe un domaine de flux délimité par deux valeurs critiques, qui produit nécessairement des gradients de pression déclenchant l'ouverture des valves quand elles sont fermées, et leur fermeture quand elle sont ouverte. Quand le flux d'entrée du système est compris entre ces valeurs critiques, les valves s'ouvrent puis se ferment cycliquement, et de manière permanente. Dans ces conditions, le flux, la pression de fluide et les conditions de transport sont intermittentes dans le système, et la sismicité se produit de concert, reflétant cette intermittence. Par ailleurs, lorsqu'une valve s'ouvre, le transitoire de pression qui est généré peut affecter les valves voisines, et déclencher leur ouverture. Cette interaction permet l'apparition de cascades d'ouvertures de valves — et donc de sismicité — qui se propagent vers l'amont et l'aval du canal. Elles sont le composant élémentaire de motifs plus complexes qui émergent dans le système alors que les interactions entre valves leur permet de synchroniser leur activité. Comme pour le tremor observé,
celui que nous simulons est épisodique, les périodes de haute activité sont elles même intermittentes, la sismicité migre le long du canal à petite et grande échelle, et dans certains cas, l'activité est très régulière, et de longues périodes émergent. Dans ces cas là, la synchronisation de l'activité de nombreuses valves dans le système fait apparaitre un comportement de valve de l'ensemble du système, très similaire à ce qui a été décrit et modélisé pendant les épisodes de tremor et glissement lent dans les zones de subduction.

Dans un dernier chapitre, nous explorons comment la force des interactions entre valve contrôle le style d'intermittence produite par le système, dans quels cas l'activité est relativement constante, et au contraire, dans quels cas elle se synchronise en épisodes intermittents, et se reproduisant à des intervalles plus réguliers. Deux paramètres contrôlent la force des interactions, la distance entre valves, et la valeur du flux d'entrée. À partir d'un système à deux valves, nous montrons que plus les valves sont proches, plus les événements qu'elles produisent sont proches dans le temps, le transitoire de pression qui porte l'interaction ayant moins de distance à parcourir. Dans un système avec de nombreuses valves distribuées dans l'espace, la proximité entre valve — qui résulte soit du nombre de valves dans le système, soit de leur regroupement en clusters dans l'espace — contrôle la force de l'interaction. Plus les valves sont proches en moyenne, plus elles interagissent. De son côté, la valeur du flux d'entrée contrôle la proximité des valves au seuil de rupture, et donc leur susceptibilité à être déclenchées par un transitoire de pression généré par l'ouverture d'une valve voisine. Quand ces deux paramètres permettent aux valves d'interagir forte- ment, l'activité est plus intermittente, et de longues périodes d'activité émergent. En conclusion, nous montrons que l'émergence de l'intermittence caractéristique au tremor doit être due à une interaction forte entre segments de valves dans le système hydraulique, produite soit par leur nombre ou proximité dans le système, soit par leur proximité à la rupture. Sur la base de nos résultats, nous émettons l'hypothèse la complexité et des longues périodes de récurrence de l'activité du tremor pourraient provenir du comportement collectif de valves dans la zone de faille, qui façonnent des régimes hétérogènes et transitoires de circulation du fluide, à des échelles spatiales allant de la source à la zone de subduction, et des échelles temporelles allant de la journée à plusieurs années. Ces caractéristiques du système hydraulique de la zone de subduction doivent être héritées de la structure de l'interface des plaques, qui elle même est modelée par la structure de la plaque plongeante. Notre modèle permet donc une interprétation des conditions hydrauliques dans l'interface de faille sur la base de la caractérisation des motifs de sismicité basse fréquence que l'on peut y détecter.

Ce étude pose un cadre conceptuel permettant de comprendre comment des conditions de pression de fluide hétérogènes et intermittentes peuvent émerger dans la faille de subduction. Notre modèle permet de concevoir l'émergence à grande échelle de motifs spatio-temporels d'activité grâce à l'interaction entre des sources sismiques individuelles, en décrivant la manière dont un événement sismique déclenche et est déclenché par des changements de pression de fluide. Au-delà des zones de subduction, ce cadre de pensée est applicable à d'autres systèmes géologiques de cir-
culation de fluide, comme les réservoirs magmatiques. Le mécanisme de valve de
perméabilité que nous décrivons est une description simplifiée mais efficace, perme-
ttant de penser le transport des fluides dans des milieux dynamiques comme étant
composé d'événements d'accumulation et débâcle de fluide. De tels processus façon-
nent l'intermittence du transport de fluide et de l'activité sismique qui en découle,
et leur compréhension nous permet d'utiliser les observations sismiques comme un
diagnostic des conditions hydrauliques en profondeur.
<figure><img src="image_8.png" /></figure>
# TABLE OF CONTENTS

Résumé de la thèse iii

Summary of the Manuscript vii

Remerciements xi

Résumé long de la thèse xiii

I Introduction xxv

1 Introduction 1

1.1 Microseismic activity probes hydraulic processes in geological systems 1

1.2 The seismic source of low-frequency seismicity 5

1.2.1 Source characteristics of low-frequency events in faults 5

1.2.2 Models of LFE source processes 7

1.3 Spatio-temporal patterns of low-frequency seismicity 11

1.3.1 Migrations of tremor and LFE activity 11

1.3.2 The intermittence of tremor and LFE activity 14

1.4 Fluid circulation in dynamically permeable systems 15

1.5 Motivations for this work 19

II The source of low-frequency earthquakes 23

2 Characterization of the seismic source of LFEs in Guerrero, Mexico 25

2.1 Introduction 25

2.2 Low-frequency earthquakes in Guerrero, Mexico 27

2.3 Selection of LFEs with high-quality waveforms 27

2.4 Measuring source parameters 32

2.4.1 Seismic moment 32

2.4.2 Corner frequency 34

2.5 Results 35

2.5.1 Seismic moment and corner frequency of LFEs 35
<table><tr><td>2.5.2</td><td>Moment-duration scaling</td><td>37</td></tr><tr><td>2.5.3</td><td>Stress-drop variations</td><td>40</td></tr><tr><td>2.5.4</td><td>Moment-frequency distribution of LFEs</td><td>40</td></tr><tr><td>2.6</td><td>Discussion and interpretation of the observed scaling law</td><td>42</td></tr><tr><td>2.6.1</td><td>Scaling laws of earthquakes, slow earthquakes and LFEs</td><td>42</td></tr><tr><td>2.6.2</td><td>Models for magnitude-independent duration of LFEs</td><td>43</td></tr><tr><td>2.7</td><td>Conclusions</td><td>45</td></tr></table>

## Appendices to Chapter 2 47

<table><tr><td>2.A</td><td>Effects of bandpass filtering on events corner frequencies</td><td>47</td></tr><tr><td>2.1.1</td><td>Analysis bandpass filter</td><td>47</td></tr><tr><td>2.1.2</td><td>Detection bandpass filter</td><td>51</td></tr><tr><td>2.B</td><td>Effects of attenuation on corner frequency f<sub>c</sub> measurements</td><td>52</td></tr><tr><td>2.2.1</td><td>Seismic attenuation in the crust of Guerrero</td><td>53</td></tr><tr><td>2.2.2</td><td>Effects of a highly-attenuating fault zone on LFE spectra</td><td>56</td></tr></table>

# III Spatio-temporal patterns of tremor activity 61

## 3 Tremor generation through fast fluid pressure transients in the per-meable subduction interface 63

<table><tr><td>3.1</td><td>Introduction</td><td>64</td></tr><tr><td>3.2</td><td>Tremor patterns observed in subduction zones</td><td>68</td></tr><tr><td>3.2.1</td><td>Spatial Segmentation</td><td>69</td></tr><tr><td>3.2.2</td><td>Tremor Migration</td><td>70</td></tr><tr><td>3.2.3</td><td>Temporal Patterns</td><td>70</td></tr><tr><td>3.2.4</td><td>Tremor activity Patterns in Guerrero, Mexico</td><td>71</td></tr><tr><td>3.3</td><td>Model design</td><td>71</td></tr><tr><td>3.3.1</td><td>Diffusion equation for pore pressure variations</td><td>72</td></tr><tr><td>3.3.2</td><td>A mechanism for large changes of permeability</td><td>73</td></tr><tr><td>3.3.3</td><td>A valve mechanism</td><td>75</td></tr><tr><td>3.3.4</td><td>Valves as elementary seismic sources</td><td>77</td></tr><tr><td>3.3.5</td><td>Permeability values, boundary and initial conditions</td><td>77</td></tr><tr><td>3.3.6</td><td>Numerical implementation</td><td>80</td></tr><tr><td>3.4</td><td>Time-dependent valve behaviour</td><td>80</td></tr><tr><td>3.4.1</td><td>A single valve</td><td>80</td></tr><tr><td>3.4.2</td><td>Valve-valve interaction via pore pressure diffusion</td><td>83</td></tr><tr><td>3.5</td><td>The collective behavior of interacting valves</td><td>84</td></tr><tr><td>3.5.1</td><td>Characteristics of valve activity in the permanent regime</td><td>84</td></tr><tr><td>3.5.2</td><td>Cascades and migrations of synthetic activity</td><td>87</td></tr><tr><td>3.5.3</td><td>Large-scale valving behavior</td><td>91</td></tr><tr><td>3.5.4</td><td>Random valve distributions</td><td>91</td></tr></table>
<table><tr><td>3.6</td><td>Discussion</td><td>92</td></tr><tr><td>3.6.1</td><td>Comparison with LFE activity patterns in Guerrero</td><td>92</td></tr><tr><td>3.6.2</td><td>Scaling numerical results to LFE activity at Guerrero</td><td>93</td></tr><tr><td>3.6.3</td><td>Perspectives for diagnosing hydraulic conditions in fault zones</td><td>95</td></tr><tr><td>3.6.4</td><td>Perspectives for exploring model configurations</td><td>96</td></tr><tr><td>3.6.5</td><td>Coupling between hydraulic transients and fault slip</td><td>97</td></tr><tr><td>3.7</td><td>Conclusion</td><td>98</td></tr></table>

## Appendices to Chapter 3 101

3.A Numerical accuracy and convergence tests 101

3.B Derivation of fracture aperture from permeability values 104

# 4 Emergence and control of tremor intermittence and periodicity 107

4.1 The intermittence of tremor in subduction zones 108

4.2 Characterizing clustering intensity and recurrence timescales 113

4.2.1 Point process description 114

4.2.2 Characterizing clustering intensity 114

4.2.3 Measuring recurrence timescales 116

4.3 A valve model to describe fault zone permeability 119

4.3.1 The subduction's hydraulics as a valve system 122

4.4 Characterizing source interaction in a two-valve system 124

4.4.1 Valves interact through pressure transients 124

4.4.2 Valve spacing controls activity synchronization 125

4.5 Effects of valve proximity in complex valve systems 128

4.5.1 Statistical description of valve distributions 128

4.5.2 Valve distribution control on activity 129

4.5.3 Parameterizing valve proximity 131

4.6 Emergence and variability of synchronization in complex valve systems 133

4.6.1 Patterns of activity clustering with high valve proximity 133

4.6.2 Valve proximity controls recurrence timescales and periodicity of activity 135

4.6.3 Input flux control on criticality and emergence of synchronization 136

4.7 Discussion 138

4.7.1 From observed to simulated tremor activity patterns 138

4.7.2 The hydraulic segmentation of subduction faults 141

4.7.3 Hydrological processes shaping tremor intermittence in Shikoku 143

4.8 Conclusion 146
<table><tr><td><strong>Appendices to Chapter 4</strong></td><td><strong>149</strong></td></tr><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>

# IV Conclusion 161

<table><tr><td><strong>5 Scope, limitations and perspectives of this work</strong></td><td><strong>163</strong></td></tr><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>

# Bibliography 200

# List of Figures 203

# List of Tables 204
<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /></figure>
# Part I

## Introduction
<figure><img src="image_11.png" /></figure>
# CHAPTER 1

## INTRODUCTION

### Contents

<table><tr><td>1.1</td><td>Microseismic activity probes hydraulic processes in geological systems</td><td>1</td></tr><tr><td>1.2</td><td>The seismic source of low-frequency seismicity</td><td>5</td></tr><tr><td>1.2.1</td><td>Source characteristics of low-frequency events in faults</td><td>5</td></tr><tr><td>1.2.2</td><td>Models of LFE source processes</td><td>7</td></tr><tr><td>1.3</td><td>Spatio-temporal patterns of low-frequency seismicity</td><td>11</td></tr><tr><td>1.3.1</td><td>Migrations of tremor and LFE activity</td><td>11</td></tr><tr><td>1.3.2</td><td>The intermittence of tremor and LFE activity</td><td>14</td></tr><tr><td>1.4</td><td>Fluid circulation in dynamically permeable systems</td><td>15</td></tr><tr><td>1.5</td><td>Motivations for this work</td><td>19</td></tr></table>

### 1.1 Microseismic activity probes hydraulic processes in geological systems

In the past few years, the rapid enhancement of computing power, allowing the implementation of computation-heavy algorithms, has led to a leap forward in earthquake detection capabilities. It has drastically lowered the magnitude of the smallest detected events, and thus increased their number. Small earthquakes are much more frequent than larger earthquakes: in most places of the world, earthquakes of a given magnitude ($M_w$1 for instance) are ten times more frequent than earthquakes one magnitude point bigger
(M<sub>w</sub>2) (Gutemberg & Richter, 1944). With lower and lower completeness magnitude — the magnitude above which one can consider that virtually all earthquakes are detected in a catalog —, the catalogs now image fainter, slower, smaller-scale processes in the Earth. The study of small earthquakes (M<sub>w</sub> < 3), or *microseismicity*, now seems key to understand the finer workings of geological systems as they provide a finer and finer observation coverage in both time and space, spanning periods thought quiescent between more intense events (eruptions in volcanoes, large earthquakes on faults) (Brodsky, 2019). The improvement of microseismic catalogs allows to precisely map structures at depth, illuminating fault geometry (Ross et al., 2020), the shape and size of seismic areas deep on the subduction interface (Rubin & Armbruster, 2013), or the geometry of magmatic plumbing systems (Hotovec-Ellis et al., 2018; Journeau et al., 2022). By more finely recording the changes in earthquake activity in time, smaller microseismic data allows to detect more elusive fore-shock activity, and in general map the minute seismic signature of the nucleation processes that precede larger fault ruptures (Ross et al., 2019). As the number of detected small earthquake go up, the time between detected events is also lower and lower, allowing to get finer and finer insight into dynamic processes. The activity of small, low-frequency events on subduction faults has for instance been used as a proxy signal for short, slow, low-amplitude slip events deep on faults, which geodetic signature is otherwise hidden in the noise (Frank, Radiguet, et al., 2015; Rousset et al., 2019).

In general, small seismic events can be used to monitor and characterize
minute stress changes on pre-stressed structures, like faults near the rupture
criterion (Brodsky & van der Elst, 2014), or high-pressure fluid flow in mag-
matic systems (Benson et al., 2008; Journeau et al., 2022). Such stress transfer
is often mediated by elastic, viscous or plastic changes in the rock, but in both
volcanic and tectonic settings, processes of fluid circulation can also trigger
and shape microseismic activity. It has been known for several decades now
that the injection of fluids in the crust — e.g. for wastewater disposal from
industrial oil and gas extraction processes, or heat collection from geother-
mal reservoirs — triggers microseismic swarms of earthquakes (Cornet et al.,
1997; Ellsworth, 2013; Majer et al., 2007; Nicholson & Wesson, 1992; Raleigh
et al., 1972). The evolution of the seismicity rate in those swarms differ from
classical, Omori behavior of seismicity observed in faults (Omori, 1895). As
the seismicity is controlled by the temporary injection of fluids — or processes
that it triggers — and not by a permanent tectonic forcing, the seismicity rate
transiently rises, is maintained high for a time, and decays. Such swarms of
seismicity also occur in places where no fluid is injected in the crust by an-
thropic activity. Because the seismicity in those settings has similar patterns,
studies have linked it to natural processes of circulation of fluids in the crust,
from deep sources towards the surface along permeable pathways, around faults (De Barros et al., 2020; Horálek & Fischer, 2008; Ross et al., 2020), in geothermal systems (Shelly et al., 2013) or in magmatic systems (Hotovec-Ellis et al., 2018).

Two factors can explain the triggering of fluid-induced seismicity: simple triggering by fluid pressure increase, or triggering mediated by aseismic slip on the fault. On the one hand, as fluids diffuse in the permeable crust, the local fluid pressure $p$ increases, countering the overburden normal stress $\sigma$ on fracture planes, and thus lowering the effective normal stress $\sigma_{eff} = \sigma - p$. When the effective normal stress that clamps them shut decreases, faults are brought closer to the failure point, as the tangential stress $\tau$ can overcome the friction on the fault. This explanation is coined the *effective normal stress principle*, and is based on the Coulomb failure criterion, that states that a fault ruptures when

$$ \tau > C + f(\sigma - p), \qquad (1.1) $$

where $C$ is the cohesion of the material in the fault, and $f$ the friction coefficient of the interface. Such a simple interpretation seems to be enough to explain migrations of seismicity that follow a diffusion pattern: as fluid diffuses, it brings small faults to failure and thus triggers earthquakes (e.g. S. A. Shapiro et al., 2002). However, in many cases, the spatio-temporal patterns of activity are more complex, and seem to indicate that the fluid pressure increase on the fault could interact with aseismic slip transients, that in turn trigger sources, by increasing the local tangential stress on seismogenic parts of local faults (De Barros et al., 2020; Ross et al., 2017; S. Wei et al., 2015). Field scale experiments of fluid injection in shallow faults seem to confirm this interpretation (Guglielmi et al., 2015): for natural and anthropic fluid injection in the crust, the fluid pressure source could act only as the trigger, initiating an aseismic slip transient that governs the propagation of seismicity (Bhattacharya & Viesca, 2019). Recent studies of the same experiment have however complexified this vision, showing that permeability can change by nearly two orders of magnitude as the slip and high fluid pressure transient travels across the fault, triggering seismicity as it progresses (Cappa et al., 2022). This evolution in the interpretation of the same dataset highlights that dynamic fluid circulation processes can be a fundamental component in explaining how fluid-related seismicity evolves in time and space, but that they are still poorly understood.

Such an implication of fluid in triggering seismicity does not only have an impact on the spatiotemporal patterns of activity, it also influences the source processes of microseismic events, that can be observed through seismic measurements (Lengliné et al., 2014; Lin et al., 2016). One of the clearest
examples of this is volcanic long-period seismicity: tremor or long-period earthquakes. Since they have been first observed, those phenomena have been linked with magmatic fluid circulation processes (B. A. Chouet, 1996; Koyanagi et al., 1987), because of the specificities of their signal: a sometimes monochromatic, especially low frequency content. When signals similar to volcanic tremor and long-period earthquakes have been detected deep on the interface of subduction zones (Katsumata & Kamaya, 2003; Obara, 2002), it was first suggested that the unusually low frequency of their signal was due either to an unusual rupture process (Ide, Beroza, et al., 2007), or to metamorphic fluid circulation around the fault interface, by analogy to volcanic processes (Katsumata & Kamaya, 2003). Although tectonic (or non-volcanic) tremor and low-frequency earthquakes have then been associated with slip on the interface (Ando et al., 2010; Obara et al., 2004; Rogers, 2003; Shelly et al., 2007a), almost excluding the fluid circulation hypothesis, there is mounting evidence that, as long period seismicity in volcanoes, low frequency seismicity in faults is closely associated with fluid circulation processes that shape their source processes and spatio-temporal patterns of activity.

On subduction faults, the circulation of fluid can be responsible for lubricating the plates' interface, producing deep slow slip transients that can load the plate interface enough to trigger the largest megathrust earthquakes (Bedford et al., 2020). Such hydraulic processes on the subduction fault remain enigmatic, but thanks to its high seismicity rate, the quality of monitoring in subduction zones producing it, and recent improvements in the resolution of catalogs of low-frequency events (Frank et al., 2014; Poiata et al., 2021), low-frequency seismicity makes an exceptional probe of the slow and faint processes leading to such catastrophic events. Therefore, our objective in this thesis will be to further the understanding of the link between fluid flow in the subduction zone and the production of low-frequency seismicity there. More generally, we hope to build a conceptual framework that can help understand better the dynamic of fluid flow in geological systems, and how it relates to seismicity. In the next pages, we describe the characteristics of low-frequency seismicity, and what is known of the hydraulic processes underlying it.
## 1.2 The seismic source of low-frequency seismicity

### 1.2.1 Source characteristics of low-frequency events in faults

In subduction zones, low-frequency seismicity locates mostly along the plate interface, between 30–60 km depth (Bostock et al., 2012; Frank et al., 2014; Shelly et al., 2006), slightly deeper than the “seismogenic zone” of the plate interface that produces megathrust earthquakes. Two major phenomena can be singled out by observing seismic waveforms: *tremor*, a long, non-impulsive signal with a characteristic frequency of about 1–8 Hz, and *low-frequency earthquakes* (LFEs), impulsive signals that can be found within tremor, with clear S-wave arrivals, and a similar frequency content as tremor. The sources of tremor and LFEs can be located to the plate interface as localized, finite patches that repeatedly activate, producing similar waveforms each time (Chestler & Creager, 2017b; Rubin & Armbruster, 2013). Because both phenomena share similar frequential content, colocate and occur at the same time, tremor is thought to be made of swarms of small, elementary LFEs (Shelly et al., 2006). Thanks to their characteristics as impulsive, locatable, repeated sources, LFEs are well-fitted to investigate the physical processes that generate tremor activity. Figure 1.1 displays the characteristic location, waveforms and spectral content of LFEs and tremor in subduction zones.

Those phenomena occur on the plate interface at pressures (about 1 GPa) and temperatures (350–550°C) that mostly prohibit brittle fracture, favoring a more ductile accomodation of the shear stresses due to the slab subduction (Platt et al., 2018). This part of the fault interface is described as a frictional transition, between shallow regions where fault slip is unstable, and can accelerate to seismogenic speeds, and deeper regions where fault slip is stable, continuously and aseismically creeping without accelerating (Lay et al., 2012; Scholz, 1998). Low-frequency earthquakes have unusual source properties, that differentiates them from classical earthquakes. The main characteristic of tremor and LFEs is of course their particular frequency content. For events of a similar magnitude, LFEs have a characteristic frequency on the order of 10 times lower than regular earthquakes (Farge et al., 2020; Shelly et al., 2007b, cf Chapter 2). The spectrum of their body waves is quite similar to an earthquake’s, but shows a much lower corner frequency $f_c$ — the frequency delimiting the high-frequency decay of the spectrum and the low-frequency content — that seems to indicate a longer duration $T \sim 1/f_c$ of the source process. In Cascadia (Bostock et al., 2015) and Guerrero (Farge et al., 2020, cf Chapter 2), their corner frequency seems to be invariant with
<figure><img src="image_12.png" /><figcaption>Tremor and low-frequency earthquakes in Guerrero, Mexico</figcaption></figure>

Figure 1.1: Caption on next page.
**Figure 1.1** — (a) Map of the Guerrero (Mexico) tremor zone. The black dots represent tremor events, and purple shaded areas the slow-slipping parts of the fault. Triangles represent volcanoes, in the transmexican volcanic belt. (b) Tomography image of the A-A' transect, overlaid with LFE detections from Frank et al. (2014), as points colored with the density of source. LFEs cluster around the plate interface, in two zones. The oceanic crust appears with a strong low-velocity anomaly, possibly indicating the presence of high-fluid pressure fluid. (c) Tremor burst signal over a few hours, recorded on one horizontal component of three stations in the Guerrero tremor zone, filtered between 1-2 Hz. Each light red dashed line represent an LFE detection. (d) LFE waveforms, on recorded by the same stations and components as the tremor waveforms, filtered between 1-2 Hz. In red, the part of the signal that is inferred to be the S-wave of the LFE. (e) Stacked amplitude spectra of the displacement for different events detected in the Shikoku tremor region: small and large earthquakes, LFEs, tremor and background noise. See Shelly et al. (2007b), Figure 2 for details. The LFE's signal has a similar low-frequency amplitude as small earthquakes, but are strongly depleted in high frequency. (a) and (b) are adapted from Audet and Kim (2016), (c) and (d) are adapted from Frank and Shapiro (2014), (e) is adapted from Shelly et al. (2007b).

respect to the seismic moment of the events: whatever their seismic moment, LFEs' corner frequency range between 1–5 Hz, with very little dependence on their size. This behavior is at odds with the scaling predicted by a simple model of seismic rupture, which predicts a strong dependence of duration and thus corner frequency on seismic moment: as earthquakes grow bigger, the duration of their source increases, and therefore their corner frequency decreases (Kanamori & Anderson, 1975). This deviation from the self-similarity scaling consistent with classical ruptures is not observed everywhere LFEs are detected. For instance, (Supino et al., 2020) find that in this respect, LFEs do not differ with classical earthquakes. In addition to the observation of a potentially constant duration for LFEs in a given region, some studies find that LFEs could result from a scale limited process, that limits their size and therefore their magnitude around a characteristic value (Chestler & Creager, 2017a), which clearly differs from the scale-independence of earthquake magnitudes observed in most seismic regions (Gutemberg & Richter, 1944).

1.2.2 Models of LFE source processes

Although it could be suggested that the especially low-frequency content
and characteristic duration of LFEs could be due to a strong attenuation of
the medium, that would effectively act as a low-pass filter on the seismic
waves (Bostock et al., 2017; Ide et al., 2003), most studies attribute them to the
specific source processes that generate LFEs.
1.2.2.1 Slip models of LFEs

Quickly after they were discovered, major bursts of tremor and LFEs have
been found to occur at the same time, and colocate with slow slip transient
on the deep parts of the subduction interface (Obara et al., 2004; Rogers, 2003).
As detection methods improved, smaller and smaller tremor episodes were
discovered, linking further and further tremor activity with slow slip on the
interface, to the point that it is now used as a proxy to detect geodetically-
invisible slip transients (Frank, Radiguet, et al., 2015; Rousset et al., 2019).
Because of seismic clues like a radiation pattern consistent with shear slip on
the plate interface (Frank et al., 2014; Ide, Shelly, & Beroza, 2007), LFEs have
since been thought to be due to the rupture of small, localized patches on a
fault which mostly slips slowly and in a stable fashion, but can on the said
patches host accelerating slip, after being loaded by the neighboring fault
slip (Ando et al., 2010; Lay et al., 2012; Luo & Liu, 2019). In this interpreta-
tion, the measured seismic moments and corner frequencies LFEs can be used
to deduce rupture parameters. Most studies find that LFEs could be due to
slow ruptures (10% to 20% the shear wave speed), on patches of radius on
the scale of 100–800 m, with very low stress drops, on the order of 10 kPa (Bo-
stock et al., 2015; Farge et al., 2020; Supino et al., 2020; Thomas et al., 2016,
cf Chapter 3). The apparent independence of the duration of LFEs on their
magnitude could owe to the fact that such ruptures are almost limited to the
slip nucleation size, due to the local frictional properties the tremor source
region (Cattania & Segall, 2019; X. Wei et al., 2021).

1.2.2.2 Fluid flow models of low-frequency seismicity

A few studies have taken inspiration in the proposed mechanisms for long
period tremor and long period events in volcanic systems to try and under-
stand how observed properties of LFEs and tremor could be the result of the
implication of strong fluid flow transients in the LFE seismic source mech-
anism (N. M. Shapiro et al., 2018). One of the principal — yet still debated
— mechanism invoked to explain volcanic long-period events and tremor is
a fast flow transient of magmatic fluid — gas or melt — in constricted path-
ways of the magmatic system. As a pocket of high pressure fluid suddenly
breaches into a lower pressure area, or a differential of pressure suddenly oc-
curs across a constriction in the conduits, due to movements of fluid on either
side for instance, fluid gushes through the constriction, and elastic waves are
generated either by the sudden traction force on the walls of the conduits,
or the development of a vibrating instability in the conduit walls as fluid
passes through (Honda & Yomogida, 1993; Julian, 1994; Ukawa & Ohtake,
1987). Such seismic behavior can be reproduced in the lab. In experiments
of sudden decompression of a fluid that permeates through a fractured rock samples, low frequency acoustic emissions that closely resemble long period seismicity are generated as fluid escapes the tortuous fractures of the rock (Benson et al., 2020, 2008). These fast fluid flow transient should occur as fractures open in a fluid-saturated, heterogeneously pressurized subduction channel, and could play an important role in shaping the source processes and therefore the seismic signature of tremor and LFEs.

Within the fluid-saturated fractures in the subduction interface, fluid should concentrate in discrete, high-porosity, high-pressure pockets (Gold & Soter, 1985). Within pockets the fluid is static, and the pressure gradient with depth should therefore be hydrostatic. As pockets are isolated by collapsed-porosity, low-permeability barriers, there exists strong steps of pressure in between fluid pockets, across low-permeability barriers. When such barriers break in a kind of hydrofracture — a change in pressure in either side of the barrier could for instance enable an accelerating flow event that opens the permeability barrier — the pressure change in the fluid should be strong and fast enough to be transmitted to the rock, and propagate as an elastic wave. N. M. Shapiro et al. (2018) models this event as a single force seismic source: the strong difference of pressure across the permeability barrier results in a force along the gradient of pressure, and as the pockets connect, the fluid pressure gradient suddenly smoothes, the force drops and seismic waves are reverberated in the medium. Figure 1.2 depicts this source mechanism. Such a seismic source radiates in a different way than a shear rupture. However, because LFEs' waveforms are only well observed from directly above, the S-wave radiation pattern from double couple source mechanism corresponding to a shear rupture on the subduction interface is almost undistinguishable from a single force radiation pattern (N. M. Shapiro et al., 2018). Several studies have attempted to characterize the source mechanism of LFEs using the patterns of first impulse, polarity or full waveforms of recorded S-waves at the surface (Ide, Shelly, & Beroza, 2007; Ohmi & Obara, 2002; Wech & Creager, 2007). Such mechanisms seem to always show a strong polarization along an assumed direction of slip, consistent with the tectonic setting, however, the data is quite ambiguous, being as consistent with the radiation pattern of a single force oriented along the updip direction of the subduction interface. In addition to that, the duration of such a source would be controlled by the characteristic duration of the pressure variation across the permeability barrier. When considering that this flow event occurs in a poro-elastic material, one can derive that the duration of the source would be proportionnal to the diffusivity of the fluid in the medium, and independent of the magnitude of the force that generates the seismic waves (N. M. Shapiro et al., 2018). Such a strong, fast fluid pressure transient source would therefore display a simi-
Simple force seismic source mechanism generated by
fast fluid pressure transient

<figure><img src="image_13.png" /><figcaption>**Figure 1.2** — Model of simple-force seismic source generated by a fast fluid pressure transient in the fault interface. (a) Fluid pressure profile along the channel, and its evolution in time. The force on the porous matrix resulting from the fluid pressure gradient is represented with arrows of size proportional to the force amplitude. Time for each curve and arrow is color coded, going from light at the start of the transient, to dark at the end. The dotted line is the fluid pressure profile before opening. The shaded area represent the extent of the barrier. (b) The fluid pressure transient is generated as a permeability barrier breaks open in the permeable fault zone. Because the barrier impedes fluid flow in the fault zone, a strong fluid pressure gradient accumulates across it, imposing a force on the matrix in the direction of the flow. As the barrier breaks, the force drops, and the generated displacement in the elastic rock radiates seismic waves away from the barrier (N. M. Shapiro et al., 2018).</figcaption></figure>

lar scaling of seismic moment with duration (Bostock et al., 2015; Farge et al., 2020, cf Chapter 3): a range of seismic moments, and a characteristic duration of the event that does not depend on its moment, controlled by the local hydraulic properties of the fault zone.

The amount of observations that link LFEs and tremor to slip on the fault is compelling, and one should not think that LFEs occur wholly without the implication of slip processes on the fault. However, the consistency of such a mechanism with seismic observations and the geologic record of chaotic fluid flow in the fault is a hint that such a mechanism could plausibly be involved in the source processes generating LFEs and tremor, and shape their activity.
## 1.3 Spatio-temporal patterns of low-frequency seismicity

### 1.3.1 Migrations of tremor and LFE activity

In contexts of anthropic fluid-induced seismicity, seismic activity is observed to migrate away from the point of fluid injection reaching distances following a clear diffusion pattern: the outtermost extent of the active region extends at a distance *R* from the injection point for a duration of injection *T*, and the two parameters follow *R*² ∝ *T* (S. A. Shapiro, 2015; S. A. Shapiro et al., 2002). The diffusion of fluid pressure in the porous medium alone can explain this pattern. As a transient increase in fluid pressure diffuses away from the source — the injection point — the activity increases where pressure allows small ruptures, by lowering the effective normal stress on the fault enough (cf equation 1.1). In other words, the seismicity first diffuses quickly away from the pressure source, and its progression slows down as the seismicity front is further and further away from the source.

Although the simple diffusive pattern can sometimes be observed for LFEs and tremor (Kato & Nakagawa, 2020; Poiata et al., 2021), they often proceed in steps: a first diffusive migration accelerating and slowing down along the fault, and a second and sometimes more steps, as if the source of fluid was reactivated further along the fault, generating a second diffusion front. Most frequently, the LFE or tremor activity migrates with a constant speed within each burst of activity. The larger migrations mostly propagate along the strike of the subduction zone on tens to hundreds of kilometers, at velocities in the range of 1 to 10 km/day (less than 1 km/hr) (Kao et al., 2009; Poiata et al., 2021) — in Guerrero, the largest LFE episodes propagate along the almost 100 km-long flat slab section, along dip (Frank et al., 2014). Smaller migrations, spanning a few to tens of kilometers, are faster, in the range of 1–100 km/h (Cruz-Atienza et al., 2018; Frank et al., 2014; Ide, 2012). High-resolution seismological networks can even image short, rapid migrations of tremor, reaching over a hundred of km/hr (Ghosh et al., 2010; Shelly et al., 2007a). The length *L* of tremor or LFE bursts seems to be inversely related to their propagation velocity *V*. This behavior in which small scale migrations are much faster than larger scale migrations is reminiscent of a diffusive behavior, for which the propagation rate of a diffusive front follows *V* ∝ 1/*R*, where *R* is the advancement of the diffusive front. The directivity of the migrations seems to depend on their scale: larger, slower migrations fill the tremorgenic zone across dip, and propagate along the strike of the interface (Ide, 2012; Kao et al., 2009; Poiata et al., 2021), while the shorter,
Low-frequency bursts and migrations in Guerrero, Mexico

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

**Figure 1.3** — Bursts and migrations of tremor in Guerrero, Mexico, represented in a time-dip diagram. In (a), (b) and (c), each dot represents a detected LFE, located in the diagram by its occurence time and along-dip distance, from the coast of Mexico. The LFE catalog was compiled by Frank et al. (2014). The three panels show progressive zooms in time and space, highlighting that the activity is intermittent at all time and space scales, and proceeds in slow and fast migrations.

faster migrations mostly proceed along-dip (Cruz-Atienza et al., 2018; Ghosh et al., 2010; Ide, 2012; Shelly et al., 2007a), with no obvious downdip or updip preference. It should be noted that rapid tremor migrations within a large burst seem to propagate in the reverse direction as the main, slower front, a phenomenon coined “rapid tremor reversals” (Bostock et al., 2015; Houston et al., 2011). In Figure 1.3, we display LFE activity patterns characteristic of the Guerrero tremor zone. Activity is intermittent, at all scales: each zoom demonstrates that an apparent burst is in fact constituted of smaller bursts of activity, separated by periods of relative quiescence. The activity migrates along the dip of the fault, in both up and downdip directions, at velocities ranging from tens of kilometers per day to tens of kilometers per hour.

The most salient feature of LFEs migrations is their constant migration speed, and that different velocities seem to coexist in the same region of activity. In a modeling study of fluid-induced swarm seismicity, Dublanchet
and De Barros (2021) reproduces constant speed, slow and fast migrations. In order to describe the seismic behavior of a prestressed fault in which fluid is locally injected, they describe fault slip in a rate-and-state framework. The fault is mostly stable frictionnaly, with small unstable patches that generate the seismic events. Fluid pressure propagates diffusively, and diffusivity is coupled with fault slip, to reproduce the increase of permeability that slip should trigger on the fault. They observe a slowly propagating front of seismicity away from the injection point, that tracks the propagation of the a slow slip front triggered by the lowered effective normal stress. Fast migrations of seismicity are observed within the slipping region of the fault, and act as cascades: each event triggering neighboring ones. Although their results seem to match closely the migrating behavior of LFEs and tremor, they also conclude that the fluid pressure is only involved in triggering the slip around the injection, and pre-stressing the fault to criticality. Because there is evidence that fluid circulation should be a strong factor shaping tremor activity, a models that only rely on slip to explain patterns should be applied with caution.

Pulse-like propagations of fluid pressure in a fault can occur when considering that a small increase in fluid pressure will also dilate pores, and thus increase the permeability of the medium (Rice, 1992). The fluid pressure then propagates as a pulse within a pocket of high porosity and permeability, in a mechanism akin to porosity waves that can occur in a pressurized, deforming porous medium (Skarbek & Rempel, 2016). In a study of rapid tremor migrations in Guerrero, (Cruz-Atienza et al., 2018) have used this framework to interpret constant-speed propagations of tremor along the dip of the subduction as propagating pulses of activity triggered by high-fluid pressure and slip pulses. Variations of permeability in response to fluid pressure can in general shape complex migrations of fluid along permeable structures, like faults. Ross et al. (2020) for instance interpret the different phases of propagation of seismicity during the Cahuilla swarm as a step by step diffusion of pressure within a fault. As a fluid source fills a confined region along the fault, the seismicity progresses within the region with a diffusive profile, until it meets a barrier of permeability. Then, pressure rises in the confined region, and a seismic event breaks the barrier, setting off a new phase of diffusive propagation of seismicity along the fault. When taken together, all those phases appear as an intermittent but almost linear propagation of seismicity along the fault. Taken to a limit, one could imagine that if the confined pockets of fluid are small enough, this barrier behavior could explain fast, linear migrations on a short scale, as activity cascades from breaking barrier to breaking barrier, triggered by the influx of fluid (N. M. Shapiro et al., 2018).
1.3.2 The intermittence of tremor and LFE activity

The final aspect of tremor activity we want to tackle here is perhaps the most
obvious, that tremor occurs in “bursts”. Tremor activity is characteristically
intermittent, occuring in finite episodes of high activity, separated by periods
of low or inexistant activity. When considering that tremor is composed of
the cumulated contribution of many elementary sources producing LFEs, the
intermittent behavior can be thought of as a phenomenon of synchronization
of sources in time: all sources being active or quiescent at the same time. It is a
striking component of this system, and characterizing and modeling it could
bring strong insight into the deformation and hydraulic processes generating
tremor.

A first characteristic of the intermittence of tremor is the quasi-periodicity, or characteristic recurrence times of the largest bursts of tremor. In Casca-dia (Brudzinski & Allen, 2007), Guerrero and Oaxaca (Frank & Brodsky, 2019; A. Husker et al., 2019), or the Nankai tremor zones (Poiata et al., 2021), the largest bursts recur relatively regularly, every several months to a few years, almost periodically. When taking a closer look at the activity, it is often hard to define a characteristic size for such bursts, however. Indeed, in some regions, tremor intermittence seem to be scale-free, proceeding in imbricated bursts of all sizes, with different recurrence times depending on their sizes (Frank et al., 2014; Poiata et al., 2021) (cf Figure 1.3). This scale free process can how-ever be characterized, by evaluating the frequency of occurence of bursts of all sizes (Frank et al., 2016). This reveals that the clustering of events changes in time, outside of bursts, the activity is smoother, occuring at a more constant rate than during large transients, where events cluster strongly, small bursts building up larger, longer bursts (Frank et al., 2016; Poiata et al., 2021). The level with which sources synchronize to build bursts is also variable in space, and delineates a clear segmentation of the tremorgenic zones, along the strike and dip of the interface (Frank et al., 2016; Poiata et al., 2021).

Many of those characteristics can be explained in a framework where tremor is generated by small shear sources triggered by a slow-slip transient on the subduction fault (Ando et al., 2010; Obara et al., 2004; Rogers, 2003), however, few studies model the full complexity of tremor patterns in a fault slip framework (Luo & Liu, 2019), from the intermittent activity of a source to the emergence of large scale patterns of recurrence and migrations. In such studies, time clustering of activity, migrations and long-period of recurrence of activity burst emerge, consistent with the observed recurrence of episodic tremor and slip events in subduction zones (Luo & Liu, 2019). The heterogeneity of the fault zone, in terms of asperity size and effective normal stress conditions can lead to domains of the fault that trigger early or late, interfer-
ring with larger scale behavior, and provoking the emergence of smaller scale bursts, recurring on a shorter timescale (Luo & Liu, 2021). During large slip transients, strong asperities start slipping, and amplify the triggering effect of the slip around neighboring asperities. In such a framework, the sources synchronize more during slip transient, and can reproduce the observed strong temporal clustering behavior during high tremor activity transients (Frank et al., 2016).

Fluid pressure and permeability are often key to explain the recurrence and intermittence patterns of tremor in such framework, by shaping the heterogeneity of the frictional parameters on the fault for instance (Luo & Liu, 2019). However they are often described in a static manner, without interaction with slip. The dynamic implication of fluid pressure should play a strong role, by shaping slip patterns or directly triggering seismic waves, and its contribution to the emergence of tremor patterns cannot be studied if it remains a simple static parameter in the description of processes that shape tremor. To our knowledge, few studies have tackled this problem in subduction zones. A study by Zhu et al. (2020) however present a model of a fault-valving behavior in a crustal strike-slip fault, in which fluid is channeled towards the surface from a deep source. In their model, slip on the fault generates a permeability increase, and healing processes shut the permeability close and increase the friction on the fault. Their description allows for both changes of pressure and permeability on the fault, and the patterns of slip generated demonstrate that the involvement of a fault-valve mechanism can generate complex, intermittent seismic behavior, as small rupture cascade in the fault in swarm-like patterns, coupled with the ascent of high-pressure fluid and opening permeability. Even though it is not reported in the aforementioned study, the accumulation and release of fluid in a sealed subduction interface could influence the period of activation of fault slip and therefore seismicity, in a larger scale fault-valve mechanism (Frank, Shapiro, et al., 2015; R. Sibson, 1992). On a smaller scale, the propagation of fast diffusive fluid pressure transients could outpace slip, and be responsible for cascading activity, and therefore temporal clustering of sources (N. M. Shapiro et al., 2018).

## 1.4 Fluid circulation in dynamically permeable systems

An array of geophysical and geological evidence indicates that fluid flow conditions in the permeable subduction channel are near-lithostatic, heterogeneous, and variable. In active subduction zones, seismic imaging of the subduction interface (Gosselin et al., 2020), stress-orientation derived from the
<figure><img src="image_15.png" /><figcaption>Geological evidence for dynamic fluid circulation in the subduction fault zone</figcaption></figure>

**Figure 1.4** — Geological evidence for transient flow conditions in the subduction fault zone. (a) Evidence for large scale hydrofracturation in the Franciscan Complex, in Catalina Islands, (CA, USA). Picture credit: John Platt, from Aristide (2020). The large, beige magmatic dike cross-cuts parallel quartzic veins, which could be symptoms of repeated episodes of hydrofracturation and catastrophic, overpressured fluid flow in the fault. (b) Hydraulic breccia in a paleo subduction channel, in the Zagros orogen (Iran). The breccia likely comes from hydrofracturation of the blueschist matrix, evidence for overpressurized fluid, and brutal variations of fluid pressure, flow, and permeability of the rock. The picture is from Muñoz-Montecinos et al. (2021), Figure 3E. (c) Veins in the Pouébo eclogite melange, in New Caledonia (France). Type I veins are evidence of progressive concentration of dehydration fluid in the medium. As they interconnect and start forming a network, large type II veins form, and channel fluid from outside the region, as evidenced by the chemical anomalies in the selvage pattern around the veins. The photo is taken from Taetz et al. (2018), where the authors find that the chemical anomaly induced by metastomatic reaction between the external fluid and the host rock can help constrain short residence time of the fuids in the veins, on the order of months or less. This is further evidence of dynamic flow processes in the fault zone.

focal mechanisms of plate-interface earthquakes (Warren-Smith et al., 2019), measures of the seismic attenuation above the slab and gravity anomaly measurements (Tanaka et al., 2018) seem to indicate that the high fluid pressure in the tremor source region drop suddenly when episodic tremor and slip occurs in subduction zones, indicating a sudden breach of the seal that confines the fluid in the interface. In outcrops of paleo-subductions, evidence of hydrofracturation (Muñoz-Montecinos et al., 2021) and cyclic vein-opening and healing (Behr & Bürgmann, 2021; Taetz et al., 2018; Tarling et al., 2021) point towards strongly varying and heterogeneous fluid pressure and fluid transport along the subduction channel (Angiboust et al., 2015).

Different processes can generate such changes of permeability. The most explored mechanism is a coupling with slip, in which slip opens up permeability, through damaging that opens micro-cracks for the fluid to circulate
in (Mitchell & Faulkner, 2008), breaking the cement that formed in the gouge of the fault plane an therefore connecting pores through newly formed pathways (Tenthorey & Cox, 2006), or through a dilatant mechanism in which the fault opens during slip, allowing the fluid to circulate (Cappa et al., 2022; Im et al., 2019b; Segall et al., 2010). In either case, slip pumps fluid into the active portion of the fault, and propagates a high-pressure fluid mass as it progresses along the fault. This type of behavior is what is referred to as a “fault-valve” mechanism. However, slip is not necessary to generate strongly heterogeneous and variable conditions of fluid pressure in the fault zone. Mechanisms that only involve the dynamics of the hydraulic system can be sufficient to generate such patterns. We want to highlight two classes of such processes, porosity waves and clogging-unclogging of the porous medium.

Porosity waves occur in mediums where porosity can change in response to a fluid influx and subsequent increase in fluid pressure. As fluid is pumped into the porous fault zone by dehydration reactions, perhaps episodically, a local increase of fluid pressure will propagate as a non-linear, sometimes solitary, wave in the decreasing direction of the fluid pressure gradient, as it dilates the pore space, thus generating new permeability for the fluid to circulate in. The dilation of pore space can be due to a simple poro-elastic effect (Cruz-Atienza et al., 2018; Rice, 1992), or to viscous creep (Skarbek & Rempel, 2016). Such propagations of pulses of high fluid pressure have been linked with the propagation of tremor and slip episodes in subduction zones (Cruz-Atienza et al., 2018; Skarbek & Rempel, 2016). However, as the resulting waves have to follow the gradient of fluid pressure, they cannot be a straightforward interpretation of the many observations of downdip-propagating tremor activity (e.g. Ghosh et al., 2010).

In geological porous medium, fine particles can partially clog the permeable pathways, and be mobilized by the flow if the traction force of the flow becomes high enough. Such behavior has been largely studied in the lab, in experiments of filtration of a fluid carrying particles through a model fault in a rock (Candela et al., 2014) or model porous medium (Bianchi, Thielmann, et al., 2018). In such systems, particles aggregate and form clogs in thin pore throats between pores, thus reducing the medium permeability, as represented in Figure 1.5. They are then suddenly removed from the constrictions as the pressure gradient increases across the clogged pathway, and the traction forces generated by fluid flow in the low-permeability clog mobilizes them. In experiments of continuous filtration, jumps of fluid pressure and variations of flow spontaneously occur as the medium clogs and unclogs. The power law size-frequency scaling of such pressure jumps is characteristic of avalanche behavior, implying that several pores could clog and unclog in cascades (Bianchi, Thielmann, et al., 2018). Both experiments and numeri-
Porosity-scale clogging and unclogging mechanism

<figure><img src="image_16.png" /><figcaption>**Figure 1.5** — Snapshots of a modelisation of pore-scale clogging (Jäger et al., 2018). The panel in the lower right represents the evolution in time of the pressure *P* behind the constriction. In (1) the fluid freely flows through the constriction, in (2) particles clog the constriction and *P* rises because of the reduced permeability, in (3) the clog is eroded. Flow velocities *u* are represented with the rainbow colors. Figure taken from Jäger et al. (2018).</figcaption></figure>

cal simulations of such processes indicate that the imposed pressure gradient
across the porous sample, fluid input rate into the sample and the particle
size compared to the constrictions control the dynamics of this phenomenon:
if the material stays clogged, if it is permanently unclogged as particles can-
not deposit, or if intermittently clogs and unclogs (Bianchi, Thielmann, et
al., 2018; Jäger et al., 2017b; Souzy et al., 2020). Pressure variations imposed
to the fluid can therefore dynamically change the permeability by triggering
clogging or unclogging of the pores, as evidenced by experiments of cyclical
pressurization of a fluid saturated fracture (Candela et al., 2014). In reality,
such behavior has been observed as passing seismic waves mobilize the fluid
in shallow crustal faults, producing strong changes of permeability (Manga et
al., 2012), and in oil reservoirs, where permeability of the reservoirs oscillates
whereas the injection stays constant, due to deposition and erosion of asphalt
modifying the permeability of the medium (Sahimi et al., 1999).

Such behavior would be extremely complex to observe in outcrops of
paleo-subductions, because high temperatures and pressures would have erased
and cristallized evidence of such pore-scale processes. However, it is plau-
sible that either micrometric cataclastic fines (Han & Hirose, 2012) or sil-
ica colloids (Amagai et al., 2019) can clog the fracture permeability in the
subduction channel, and that such clogging/unclogging processes be active
there (Prasianakis et al., 2017). If the pressure drops resulting from unclogging are strong enough, neighboring clogs could open in cascades, triggered by the flow through the initially breaking clog. Strong, rapid and localized pressure transients would therefore occur along the fault zone, and could propagate in both the up and downdip direction, as unclogging can occur due to an increase of pressure behind a plug, or through a decrease of pressure after a plug. Such rapid changes of pressure can act as a seismic source, or be responsible for triggering one (N. M. Shapiro et al., 2018), and therefore reproduce patterns of intermittent, cyclic and migrating seismicity along the fault, as it clogs and unclogs intermittently. It is not obvious how permeability-scale clogging and unclogging could relate with the evidence of hydrofracturation of the fault zone, but if many plugs can unplug at the same time, it would locally generate strong gradient of pressure that could be responsible for hydrofracturation in the fault. Alternatively, the clogging/unclogging description could be a convenient physical approximation for hydrofracturation, during which the cohesion of the fault zone material is overcome by the stress exerted by locally overpressurized fluid.

These processes are seldom taken into consideration when trying to understand tremor and LFE activity, and therefore their implication is not fully understood. However, because they can be responsible for strong, rapid, localized transients of pressure and permeability, and cyclic and migrating behavior of those transients, they could be a major force in shaping the source and spatio-temporal patterns of tremor activity.

**1.5 Motivations for this work**

As developed in the previous paragraphs, many aspects of low-frequency seismicity in subduction zones seem to involve complex fluid flow dynamics in the fault zone, from the seismogenic source processes to the processes responsible for the spatio-temporal patterns of activity. Most models that attempt to explain either the source or patterns however mostly focus on the implication of fault slip, often with a narrow range of interaction with fluid pressure or permeability, which does not allow to explore how their dynamics could shape the activity. In the subduction fault zone however, there is pervasive evidence that permeability is strongly dynamic, and that in itself, the dynamics of the hydraulic system are a dominant factor shaping the activity of tremor. In this thesis, our objective is therefore to build a framework that would help us evaluate this impact, in order to answer two complementary questions:

* From observations of the source and spatio-temporal patterns of low-frequency seismicity in subduction zones, can we get more insight into
how fluid dynamically circulates in the fault?

* How does the implication of transient opening and closing of the fault permeability could shape the source, intermittence and migrations of seismicity?

We will focus on the characteristics of low-frequency seismicity in subduction zones highlighted in this introduction: the unusual moment-duration scaling of LFEs, the intermittence of its activity, characterized by highly synchronized, migrating bursts recurring on long timescales, with quasi-periodic regularity. In order to reproduce such patterns and to understand how the hydraulic structure and dynamics in the fault controls the seismicity, we will build a model that links fluid circulation in the fault with source processes of LFEs, to describe how sources can interact through the hydraulic system, thus building emergent patterns of activity. This model will focus on the strong changes of fluid pressure and permeability observed in active and fossilized subduction faults, in order to build a framework that allows to explore their role in shaping the seismicity we observe. It will connect the hydraulics of the system with processes shaping both the source and patterns of tremor and LFE activity, and produce catalogs that can be directly compared with the seismicity observed in subduction zones.

Although we will argue that this framework is more general, we will focus on understanding low-frequency seismicity in subduction zones. Because subduction zones produce megathrust earthquakes that can grow to be larger than $M_{w9}$, thus a major threat to the region, they are highly monitored. The long and wide coverage of seismic observations thus provided can help image long-term processes with a high resolution. The quality of this data is precious to understand further the mechanics of tectonic and hydraulic stress accumulation in the fault, and our work will therefore participate in building a more detailed and thorough understanding of the earthquake cycle in these faults. More generally, by directly connecting unsteady fluid circulation processes and seismic activity in a simple model, we provide a framework that helps conceptualize how fluid flows in geological systems, how it shapes seismicity, and in return, what the seismicity can tell us about fluid flow processes. Further understanding microseismicity as a symptom of dynamic flow processes is precious to monitor hazardous systems, such as industrial injection of fluids in the crust that can cause intermediate-size earthquakes, or the build-up towards eruptions in volcanic systems.

This thesis is organized as follows. In a first part, we characterize the source parameters of LFEs in Guerrero, Mexico, in order to further understand their source mechanism, and the possible involvement of transient fluid processes in generating them. In a second part, we elaborate a model of how fluid transport properties should change in response to fluid pressure in the
fault, and how such changes can be responsible for seismicity. We then explore the activity regime of this system, in order to find the hydraulic conditions necessary for a realistic, sustained but intermittent seismicity to occur, under constant fluid flow conditions akin to the constant dehydration in subduction zones. After highlighting intermittence and migration patterns that seem to reproduce the behavior of tremor and LFE activity in subduction zones, we proceed to evaluate how hydraulic interactions in the fault zone shape the temporal clustering and periodicity of the simulated activity. This last part of our work allows us to discuss how the structure of the plate interface should play a role in shaping the hydraulic system of subduction, which then affects how seismicity manifests. Before concluding, we take a final look at the consistency of the modeled processes with the geological observations of fossilized fluid flow processes in subduction zones, and discuss the relevance of this framework to explain volcanic low-frequency seismicity as a symptom of fluid circulation in the plumbing system.
<figure><img src="image_17.png" /></figure>
# Part II

## The source of low-frequency earthquakes
<figure><img src="image_18.png" /></figure>
# CHAPTER 2

## CHARACTERIZATION OF THE SEISMIC SOURCE OF LFEs IN GUERRERO, MEXICO

### Contents

<table><tr><td>2.1 Introduction</td><td>25</td></tr><tr><td>2.2 Low-frequency earthquakes in Guerrero, Mexico</td><td>27</td></tr><tr><td>2.3 Selection of LFEs with high-quality waveforms</td><td>27</td></tr><tr><td>2.4 Measuring source parameters</td><td>32</td></tr><tr><td>2.4.1 Seismic moment</td><td>32</td></tr><tr><td>2.4.2 Corner frequency</td><td>34</td></tr><tr><td>2.5 Results</td><td>35</td></tr><tr><td>2.5.1 Seismic moment and corner frequency of LFEs</td><td>35</td></tr><tr><td>2.5.2 Moment-duration scaling</td><td>37</td></tr><tr><td>2.5.3 Stress-drop variations</td><td>40</td></tr><tr><td>2.5.4 Moment-frequency distribution of LFEs</td><td>40</td></tr><tr><td>2.6 Discussion and interpretation of the observed scaling law</td><td>42</td></tr><tr><td>2.6.1 Scaling laws of earthquakes, slow earthquakes and LFEs</td><td>42</td></tr><tr><td>2.6.2 Models for magnitude-independent duration of LFEs</td><td>43</td></tr><tr><td>2.7 Conclusions</td><td>45</td></tr></table>

### 2.1 Introduction

LFEs are a member of the "slow earthquake" class of phenomena, regrouping events which seem to be associated with slow slip processes on faults (Beroza & Ide, 2011). Those events exhibit a wide range of durations (from slightly less than 1 s for LFEs, up to nearly a year for slow slip events) and
magnitudes (from around $M_w1$ up to $M_w7.5$). When accounted as a one fam-
ily of fault slip processes, their seismic moment $M_0$ appears to be directly
proportional to their duration $T$ (Ide, Beroza, et al., 2007), whereas regular
earthquakes align along a $M_0 \propto T^3$ scaling law (Allmann & Shearer, 2009;
Houston, 2001; Prieto et al., 2004) due to the self-similarity of the rupture pro-
cess (Kanamori & Anderson, 1975). However, considering slow-earthquakes
as one entity relies on bridging two different scales of observations: geodetic
observations of large-scale, long-term (sub-daily at minima) slow slip events
and seismic observations of short term events as LFEs, and longer events (30-
60 s) called very-low-frequency earthquakes (VLFEs). There seems to be evi-
dence that the shortest slow slip transients measured to this day align along a
self-similar moment-duration scaling $M_0 \propto T^3$ (Frank & Brodsky, 2019), and
that only at the largest scales it would transition to a $M_0 \propto T$ scaling law
(Gomberg et al., 2016). This observation isolates LFEs and VLFEs from the
bulk of geodetic observations of slow slip. It brings into question whether
they are simply smaller scale manifestations of slow slip, or they have a more
complex coupling relationship with it. In particular, LFEs co-occur with slow
slip events, and investigating their properties provides key insight to under-
standing fault state, geometry of activation and overall rupture process dur-
ing slow slip activity (Beeler et al., 2013; Chestler & Creager, 2017a; Frank,
Shapiro, et al., 2015; Houston et al., 2011, e.g.).

In this study, we investigate the moment-duration relationship for LFEs
in Guerrero, Mexico. As for regular earthquake, evaluating how moment re-
lease scales with source duration is key to understand the dynamics of the
process generating LFEs. Similar studies have already been carried out in
regions witnessing tremor and slow slip. In Cascadia, LFEs source duration
has been found to be very weakly dependent on their seismic moment, scal-
ing along $M_0 \propto T^{10}$ (Bostock et al., 2015). On the other hand, in Nankai, LFEs
exhibit a self-similar behavior, scaling along $M_0 \propto T^3$ (Supino et al., 2020). To
our knowledge this is the first published work to measure source parameters
for LFEs in Guerrero, Mexico. We start with a large LFE catalog compiled
by Frank et al. (2014) and select those whose waveforms have a sufficiently
high quality to measure their seismic moments and corner frequencies. The
latter are estimated with two different methods to check the robustness of
the obtained results. Similar to the study of Bostock et al. (2015) for LFEs in
Cascadia, we find that for the LFEs in Guerrero the corner frequency is very
weakly dependent on the seismic moment.
## 2.2 Low-frequency earthquakes in Guerrero, Mexico

Slow-slip events (SSE), tectonic tremor, and low-frequency earthquakes have been observed in the state of Guerrero, Mexico, in the subduction zone of the Cocos plate under the North American plate (Frank et al., 2013; Kostoglakov et al., 2010; Payero et al., 2008). They have been detected in the Guerrero seismic gap — a corridor extending from Acapulco to Mexico City — that has not witnessed major earthquakes ($M_w > 7$) for more than 100 years (Kostoglakov & Pacheco, 1999). This longer recurrence time of major earthquakes compared to adjacent regions can be related to the large SSEs that accommodate a significant fraction of the plates convergence in the Guerrero seismic gap (Radiguet et al., 2012a).

We use the data collected during the Meso-American Seismic Experiment (MASE) (A. Husker et al., 2008; Perez-Campos et al., 2008), between 1 January 2005 and 15 April 2007. Following the work of Frank and Shapiro (2014), we selected a subset of 10 stations out of the 100 total based on their high signal-to-noise ratio (SNR), the continuity in time of their records, and their proximity to the LFE activity (see Figure 2.1). Frank et al. (2014) used this dataset and the method described in Frank and Shapiro (2014) to compile a large catalog of 1.8 million LFEs grouped in 1120 families of repeating events. We use this catalog as a starting point in our study. It should be noted that the events in the catalog are detected by a matched-filter search using template events, filtered between 1 and 2 Hz. This band-pass is narrower than the range most studies of LFEs using matched-filter search have used up to now, usually 1–8 Hz (e.g Bostock et al., 2015; Chamberlain et al., 2014).

Nonetheless, we show in appendix 2.A that performing the matched-filter search in a narrow band does not bias the event detection and measurements of events’ seismic moments and corner frequencies. Within the range of measured seismic moments, events with a corner frequency higher than 1 Hz do not have a lower probability of being detected because of the band-limited nature of detection.

## 2.3 Selection of LFEs with high-quality waveforms

The catalog of Frank et al. (2014) has been created on the network-based template matching detection method (Gibbons & Ringdal, 2006, for instance). The advantage of this method is that it detects many events with very weak signals hidden in the noise. At the same time, the signal-to-noise ratio for most LFEs in the final catalog is too low to be suited for a determination
<figure><img src="image_19.png" /></figure>

**Figure 2.1** — Spatial and temporal distribution of the selected LFEs. Selected LFEs are plotted as white and blue dots, the bulk catalog of Frank et al. (2014) is in gray. Subduction interface depth contours are displayed in black (Kim et al., 2010). The focal mechanism corresponds to an average mechanism computed from cataloged LFE detections (Frank et al., 2013). (c) shows the time-dip distribution of the events.
of their source parameters. Therefore, we selected LFEs with high-quality waveforms that could be used to robustly measure their seismic moments and corner frequencies.

Some single-component waveforms of catalog detections are contaminated with very high amplitude spikes, instrumental noise or waves from strong earthquakes. In these cases, useful signals are masked and corresponding waveforms cannot be used for determination of the LFE source parameters. We remove those seismograms from the analysis. Determination of seismic moment requires knowing the earthquake source location. However, LFE family locations determined by Frank et al. (2014) contain a strong uncertainty on position in the direction perpendicular to the station array. This uncertainty is inherent to the linear observation geometry of the MASE network. Consequently, a few families of LFEs have templates with diffuse seismic phases, indicating that their measured position is loosely constrained. We therefore remove those families of events from our analysis.

We establish a set of criteria to automatically select LFEs and waveforms suitable for measuring the source parameters. We select 5 s long signals starting 2 s before the cataloged detection time of LFE S-waves.

This allows us to capture the most impulsive part of the S-wave signal, and allows for slight location and detection time imprecisions. This is done so as to trade-off between keeping a maximum of information and not allowing for too much diffuse seismic energy in the form of other seismic phases and events. By isolating the most impulsive part of the wave train, we best constrain the characteristic frequencies and seismic moment of the event.

In a next step, we reject signals with low amplitude and impulsivity to exclude false detections and contaminations with surface waves. First we define two frequency ranges in which we determine both an impulsivity $I$ and an amplitude $A$: they will allow the criterion to be more flexible as for which events are taken into consideration. The impulsivity is calculated as the ratio of the seismic energy in the waveform time window $W$ to the seismic energy contained in a three times larger background time window $B$ around the detection, containing it:

$$I = \frac{\int_{W} s(t)^2 dt}{\int_{B} s(t)^2 dt} \qquad (2.1)$$

where $s(t)$ is the velocity signal. $I$ varies between 0, for least impulsive, and 1, for most impulsive. It is used as a proxy of the detection SNR, but it also addresses the difficulty of defining what is noise within the tremor, where the LFE rate is very high. The amplitude is calculated as the maximum of the envelope of the filtered signal within the 5 s time window. For this purpose,
the waveforms are converted into complex analytical signals $s^A(t)$:

$$ 
\begin{aligned}
s^A(t) &= E(t) \cdot e^{i2\pi f^i(t)} \\
A &= \max_W(E(t))
\end{aligned} 
\qquad (2.2) 
$$

$E(t)$ is the signal envelope and $f^i(t)$ is its instantaneous frequency. The latter will be used later on to determine the corner frequency of the selected LFEs.

Amplitudes and impulsivities are computed for all events on each channel of each station in the 1–2 Hz frequency band. This frequency range is characteristic of Mexican tremor signal in velocity, and has been shown to yield the best SNR for LFEs in this region (Frank et al., 2014; Payero et al., 2008).

So as to robustly ensure that the detection is a correctly located LFE we also base the selection criterion on the correlation coefficient used in the making of the LFE catalog. It is a measure of the correlation of an event waveforms on all 15 detection channels to the stacked-waveforms template of its family (Frank et al., 2014). We consider the highest absolute values of the correlation coefficient to represent a well-constrained location and an impulsive event. As a sum of 15 normalized correlation coefficients, it is normalized by the number of channels used in the detection of the event and thus varies between 0 and 1.

The selected waveforms have to fulfill the following criterion, where $I$ is the waveform impulsivity and $A$ its amplitude:

* The LFE should be impulsive, but waveforms with very high impulsivity are instrumental noise or earthquake wave arrivals:

$$ 0.5 < I < 0.92 \qquad (2.3) $$

* Higher-amplitude LFEs are generally detected in distant earthquakes waveforms. Therefore, we empirically limit the amplitude of a selected detection to $2 \times 10^{-6}$ m.s$^{-1}$. A lower threshold for selection $T_{low}(st, ch)$ is fixed for each station and component (approximately $5 \times 10^{-9}$ m.s$^{-1}$). Both bounding values are manually determined using the impulsivity-amplitude distribution on each channel of each station (see Figure 2.2).

$$ T_{low}(st, ch) < A < 2 \times 10^{-6} \text{m.s}^{-1} \qquad (2.4) $$

For our source analysis, we select LFEs that contain at least 10 different channels satisfying the described criteria and have a network correlation coefficient superior to 0.4. This value has been chosen empirically, as the median of the correlation coefficients in the catalog.

A benefit of the network template matching method is that it allows to detect earthquakes with very close origin times, and thus with overlapping
<figure><img src="image_20.png" /></figure>

**Figure 2.2** — Selection of LFEs. (a.) Schematic recapitulation of the selection procedure. (b.) Amplitude-impulsivity distribution of all detections recorded on the East component of the station TONA. Each dot is an LFE, its color corresponds to the number of channels on which it is recorded that pass the criteria. Selected LFEs are white dots: they are recorded on 10 or more channels passing the criterion and fulfill the recurrence and correlation requirements. The amplitude and first impulsivity criteria (equations 2.4 and 2.3) are shown with dashed lines. LFEs in the criterion domain pass the criterion on TONA.HHE, but 9 other records passing the criterion on other channels of the network are needed for an LFE to be selected.

signals at some stations. Overlapping signals are not suitable for the source parameter analysis. Therefore, if several LFEs are detected within a 10 s time window, only the LFE with the highest correlation coefficient, and thus most reliable waveforms, is kept. Figure 2.2 sums up the selection process.

The selection yields 3498 LFEs representing 822 families out of a total of 1120. Each event is recorded on 10 selected channels over the network and shown with blue circles in Figure 2.1. The most impulsive waveforms of two selected LFEs are displayed in Figure 2.3. The selected LFEs exhibit a similar activity behavior in space and time as Frank et al. (2014) observed for the whole catalog. First, they occur mainly in bursts within the sweet spot region. Second, during the 2006 SSE, bursts are no longer limited to the sweet spot, but extend updip in the transient zone. The activity in those two regions resembles the description Frank et al. (2014), the sweet spot has a nearly con-
<figure><img src="image_21.png" /><figcaption>**Figure 2.3** — Selected LFEs' waveforms. The 4 most impulsive waveforms of LFE 379 and 687 are displayed here, 1–2 Hz band-pass filtered ((a.) and (c.)), and between 1 and 5 Hz ((b.) and (d.)). The blue section of the waveform shows the time window used to define the detected signal, the red patch shows which channel's waveforms are used as an example for source parameters characterization in Figure 2.4. Impulsivity and amplitude are indicated for each channel.</figcaption></figure>

stant activity over the 2-year span, whereas the transient zone is active mainly during the 2006 SSE and more sparse in the inter-SSE period.

## 2.4 Measuring source parameters

### 2.4.1 Seismic moment

We estimate seismic moment $M_0$ from the displacement spectra of S-waves $\tilde{u}(f)$ that converges to a constant value at low frequencies (Brune, 1970). This constant is proportional to $M_0$ (Aki & Richards, 2002). For every selected waveform, we can determine the low-frequency spectral asymptotic value $\Omega_0^{i,j}$ as:

$$ \Omega_0^{i,j} = \lim_{f \to 0} \tilde{u}^{i,j}(f) = C^{i,j} \cdot M_0^{i,j}, \quad (2.5) $$
where *i* corresponds to the LFE number and *j* is an index given to each of the corresponding selected waveforms. The $C^{i,j}$ factor accounts for the source radiation pattern and the seismic wave propagation.

To measure these low-frequency asymptotic values, the 5 s long S-wave velocity waveforms selected in the previous section are first detrended and high-pass filtered above 1 Hz, to remove any influence of spectral leakage of microseismic noise, and then integrated into displacement. The power spectral density (PSD) of the displacement waveforms is computed with the multi-taper method (Thomson, 1982). It gives more robust spectral estimations than a simple discrete Fourier transform, especially for short signals. Taking the square root of this spectral estimation yields the displacement spectral amplitude. The low-frequency value of the displacement spectral amplitude $\Omega_0^{i,j}$ is estimated as the geometric mean of the 1–2 Hz portion of the spectral amplitude of displacement, filtered between 1 and 8 Hz (Equation 2.5, and Figure 2.4).

Then, seismic moment estimates from every waveform $M_0^{i,j}$ can be retrieved by computing the propagation effects factor $C^{i,j}$. For this purpose, we compute synthetic waveforms for each source-station couple, using LFE family locations from Frank et al. (2014) as hypocenters, an average regional seismic velocity model (Iglesias et al., 2010) and an average LFE focal mechanism, determined by Frank et al. (2013). We used the AXITRA software (Coutant, 2008) based on the discrete wavenumber method (Bouchon, 2003).

After the synthetic seismograms are computed, we select a 5 s long time window around the S-waves, process them in the same way as the LFEs' seismogram and compute their displacement spectrum without attenuation $\tilde{u}^s(f)$, where *f* is the frequency. The seismic wave attenuation effects are modeled in the spectral domain by introducing $a(f)$. A corrected synthetic spectral displacement $\tilde{u}^s(f)$ becomes:

$$ \begin{aligned} \tilde{u}^s(f) &= a(f) \cdot \tilde{u}^s(f) \\ a^{i,j}(f) &= \exp\left(\frac{-2\pi f X^{i,j}}{2V_S Q(f)}\right), \text{ where } Q(f) = 273 f^{0.66} \end{aligned} \qquad (2.6) $$

where the frequency-dependent quality factor $Q(f)$ for the Guerrero region has determined by Ordaz and Singh (1992), the average shear wave velocity $V_S$ is chosen to be 3.5 km.s⁻¹ and $X^{i,j}$ is the source-receiver distance.

We determine the low-frequency asymptote of the synthetic spectrum $\tilde{\Omega}_0^{i,j}$ as the maximum of the synthetic spectrum, filtered between 1 and 8 Hz. Knowing the synthetic source moment $\tilde{M}_0$ and following equation (2.5), $C^{i,j}$ can be eliminated:

$$ \frac{\tilde{\Omega}_0^{i,j}}{\tilde{\Omega}_0^{i,j}} = \frac{M_0^{i,j}}{M_0}, \qquad (2.7) $$
<figure><img src="image_22.png" /><figcaption>**Figure 2.4** — Determination of the seismic moment $M_0$ and corner frequency $f_c$, for an example low-frequency earthquake. Example LFE: nb. 379, channel HHE of station TONA. (a.) $\Omega_0$ is determined as the red plateau of the displacement spectrum. (b.) $^s\Omega_0$ is determined as the maximum of the synthetic displacement spectrum, taking into account attenuation. (c.) The corner frequency $f_c$ is measured as the mean instantaneous frequency (MIF) over the waveform time window (see Figure 2.3 for the LFE waveforms), (d.) and as the maximum of the velocity spectrum (MVS) of the waveform. In gray, the un-filtered velocity spectrum, allowing to assess the effects of the 1–8 Hz filter on the spectrum.</figcaption></figure>

and a seismic moment estimated for each waveform of a selected LFE becomes:

$$M_0^{i,j} = {}^s M_0 \frac{\Omega_0^{i,j}}{s \Omega_0^{i,j}}. \qquad (2.8)$$

### 2.4.2 Corner frequency

We use two different techniques to estimate the corner frequency from all selected waveforms for every LFE, $f_c^{i,j}$. For both methods, the waveforms are detrended, corrected for attenuation (see equation 2.6) and band-pass filtered between 1 and 8 Hz to remove any influence of micro-seismic signal and higher frequency noise from the frequency band that correspond to LFEs frequencies. Further discussion on the choice of this frequency band can be found in the appendix, section 2.A.
The first approach consists in automatically picking $f_c^{i,j}$ as the frequency corresponding to the maximum in the velocity amplitude spectrum of every waveform. This method is labeled MVS (maximum of velocity spectrum) in the illustrations. The bandpass filter creates a slight amplitude reduction up to 1.5 Hz, which can be seen when comparing the filtered and unfiltered spectra in Figure 2.4. This could shape the spectrum into a peak slightly above 1 Hz, if the actual peak frequency between 1 and 8 Hz is close to 1 Hz. However, we find that the 1–8 Hz filter band is the best compromise between allowing microseismic noise to bias the corner frequency measurement and slightly constraining the shape of the spectrum. The effect on measured corner frequencies of the low-cutoff of the filter is discussed more extensively in section 2.A and can be assessed in figure S2 of the supporting information file.

In the second approach, we estimate the corner frequency as the dominant instantaneous frequency of the S-waves. The waveforms are converted into complex analytical signals (see equation 2.2) and $f_c^{i,j}$ is estimated as the mean of the smoothed instantaneous frequency over the 5 s time window containing the S-waves. This method is labeled MIF (mean of instantaneous frequency) in the illustrations. For this method too, the effects of the filter band can be assessed in figure S2 of the supporting information file.

For each LFE, we obtain a final estimate of their moment $M_0^i$ and two estimates of their corner frequency $f_c^i$, all computed as medians of those measured from individual waveforms selected for this LFE. For each event, the uncertainty on its source parameters is estimated as the standard deviation of the values of the given source parameter measured on each record of the event.

# 2.5 Results

## 2.5.1 Seismic moment and corner frequency of LFEs

Seismic moments and corner frequencies estimated for the selected LFEs are shown in Figure 2.5 along with the schematic imprint of the $M_0-f_c$ distribution for classical earthquake, which follows $f_c \propto M_0^{-1/3}$, for stress drops ranging from 0.1 to 100 MPa (Allmann & Shearer, 2009), and values of $M_0$ and $f_c$ determined for LFEs in previous studies. Seismic moments and magnitudes estimated for Mexican LFEs range between $M_0 = 4.9 \times 10^{11}$ N.m ($M_w \sim 1.7$) and $M_0 = 2.3 \times 10^{14}$ N.m ($M_w \sim 3.5$). These values are higher than previous estimations in other regions (Bostock et al., 2015; Ide, Beroza, et al., 2007). The median uncertainty on moment magnitude of $\delta M_w \approx 0.23$. It is mainly due to errors in modeling the velocity structure of the crust, assuming a generic
<figure><img src="image_23.png" /></figure>

**Figure 2.5** — Distribution of the measured source parameters. The distribution is displayed and colored according to the point density for both methods of determining $f_c$: (a.) with mean instantaneous frequency (MIF), (b.) as the frequency of the maximum of the velocity spectrum (MVS). Weighted least-square fits to the binned moments (see Figure 2.6) are shown as dotted lines, in each respective plot. Generic values for $M_0$ and $f_c$ measured in other regions are displayed along the distribution (1. Thomas et al. (2016), 2. Supino et al. (2020), 3. Ide, Beroza, et al. (2007), 4. Bostock et al. (2015)) The self-similarity scaling law for classic earthquakes is represented by the purple patch. The higher bounding line corresponds to a stress drop of 100 MPa and the lower to a stress drop of 0.1 MPa (Allmann & Shearer, 2009; Eshelby, 1957; Madariaga, 1976).
radiation pattern for all events and large uncertainties on the hypocentral locations (Frank et al., 2014, supplementary figures 3 to 17). The distribution of seismic moments above a threshold magnitude of $M_{wthr} \sim 2.32$ follows a power-law, with an exponent $\beta = 1.62$ (corresponding b-value $b = 2.43$, see Appendix, Figure 2.7).

The measured corner frequencies mean and extrema for each method are displayed in Table 2.1. They are approximately centered around $f_c \sim 3.0$ Hz and fall within the accepted range of characteristic frequencies for LFEs of 1–6 Hz (Bostock et al., 2015; Ide, Beroza, et al., 2007; Shelly et al., 2007b; Thomas et al., 2016). The median uncertainty on corner frequency is $\delta f_c \approx 1.6$ for the maximum of velocity spectrum method, and $\delta f_c \approx 0.8$ for the mean instantaneous frequency method.

<table><caption>**Table 2.1** — Measured LFE source parameters statistics analysis values</caption><thead><tr><th rowspan="2">Characterization method</th><th colspan="3">Corner frequency f<sub>c</sub> (Hz)</th></tr><tr><th>Min</th><th>Median</th><th>Max</th></tr></thead><tbody><tr><td>Maximum of instantaneous frequency (MIF)</td><td>2.0</td><td>3.2</td><td>5.5</td></tr><tr><td>Maximum of velocity spectrum (MVS)</td><td>1.1</td><td>2.6</td><td>6.5</td></tr></tbody><thead><tr><th rowspan="2">Characterization method</th><th colspan="3">Seismic moment M<sub>0</sub> (N.m)</th></tr><tr><th>Min</th><th>Median</th><th>Max</th></tr></thead><tbody><tr><td rowspan="2">Low-frequency spectral asymptote</td><td>4.9 &times; 10<sup>11</sup></td><td>2.8 &times; 10<sup>12</sup></td><td>2.3 &times; 10<sup>14</sup></td></tr><tr><td>M<sub>w</sub>1.7</td><td>M<sub>w</sub>2.2</td><td>M<sub>w</sub>3.9</td></tr></tbody></table>

## 2.5.2 Moment-duration scaling

The $M_0-f_c$ distribution in Figures 2.5 and 2.6 shows that the corner frequencies seem virtually independent of seismic moment. In order to estimate the scaling law between the source parameters, $M_0-f_c$ points are binned into 35 bins of ascending moment magnitude, containing 100 points each. Moment magnitude of the bin is considered to be the mean $M_w$ of the events in the bin. Corner frequency of the bin $f_c^{bin}$ is the weighted mean of the measured corner frequency of the binned events. The weight associated to the corner frequency of event $k$ is chosen to be the inverse of the variance of corner frequencies measured on each record for event $k$, $w_k = 1/\sigma_k^2$. Figure 2.6 shows
the bin distribution and the associated estimation of variability of corner frequency within bin, calculated as a weighted standard deviation:

$$ \sigma^{bin} = \sqrt{\frac{\sum_k (f_c^{bin} - f_c^k)^2 \times w_k}{\sum_k w_k}} \qquad (2.9) $$

with $f_c^k$ the measured corner frequency for event of index $k$ in the considered bin, $w_k = 1/\sigma_k^2$ its associated weight (see above), and $f_c^{bin}$ as the bin's corner frequency.

Using the bins $M_0$ and $f_c$ and their associated weight $w^{bin} = 1/\sigma^{bin^2}$, a weighted least-square fit to the following power law is performed, in log-log space:

$$ f_c \propto M_0^{-\alpha} \Leftrightarrow \log(f_c) = A - \alpha \log(M_0). \qquad (2.10) $$

The resulting values of the slopes confirm that corner frequencies are weakly dependent on seismic moment: for the mean instantaneous frequency method (MIF) used to obtain $f_c$: $\alpha = 0.052$, for the maximum of velocity spectrum method (MVS): $\alpha = 0.115$, (see Figures 2.5 and 2.6). This result is similar to the moment-duration scaling observed for LFEs in the Cascadia subduction zone by Bostock et al. (2015), where $\alpha \sim 1/10$.

We proceed to estimate the probability density function (PDF) of the exponent of the scaling using a bootstrapping method. For each bin, we keep its moment but pick a new corner frequency from a normal distribution parametrized by the bin's corner frequency $\mu = f_c^{bin}$ and by the bin's weighted standard deviation $\sigma = \sigma^{bin}$ (see equation 2.9). An unweighted, least-square fit to a power law (equation 2.10) is then performed on the bin's seismic moment and re-sampled corner frequency. After 200,000 iterations of the previous steps, we obtain an estimate of the probability distribution function for the scaling parameter $\alpha$. The mean and standard deviation of the obtained distribution give a new estimation of the slope and its associated uncertainty. For the MIF method, $\alpha = 0.054 \pm 0.034$, for the MVS method, $\alpha = 0.12 \pm 0.11$. Assuming that $\alpha$ is normally distributed, the probability that it is higher than 0.25 (thus closer to 1/3) for our $f_c - M_0$ distribution is around 0.12 for the MVS method ($\alpha = 1/3$ is $\sim 1.2\sigma$ away from the mean $\alpha$), and lower than $10^{-8}$ for the MIF method ($\alpha = 1/3$ is $\sim 5.7\sigma$ away from the mean $\alpha$).

The scaling exponent $\alpha$ we measure is thus significantly lower than the earthquake power law exponent $\alpha = 1/3$, for the MIF method at least. The MVS method gives more ambiguous results with regard to the value of the exponent $\alpha$, due to the wider distribution on $f_c$.
<figure><img src="image_24.png" /></figure>

**Figure 2.6** — Binned distribution and uncertainties estimation. The measured $M_0-f_c$ is displayed and colored according to point density, for both the mean instantaneous frequency (MIF) and maximum of velocity spectrum (MVS). Seismic moment and corner frequency are binned, and the weighted standard deviation is shown with the error bars (see text for details). Weighted least-square fits to the binned data are displayed, the fit resulting in the purple scaling has $\alpha = 1/3$ as a constraint.
### 2.5.3 Stress-drop variations

The stress-drop for an event of seismic moment $M_0$ modeled by a circular dislocation of radius $r$ is expressed as follows (Madariaga, 1976):

$$ \Delta\sigma = \frac{7}{16} \frac{M_0}{r^3}. \qquad (2.11) $$

Assuming a rupture velocity proportional to the shear wave velocity, $v_r \propto V_S$, it follows that the source radius is inversely proportional to the corner frequency (Madariaga, 1976):

$$ r \propto \frac{V_S}{f_c}. \qquad (2.12) $$

From equations 2.11 and 2.12, and keeping the assumptions that our events originate from regions of constant $V_S$ and share the average corner frequency of our distribution $f_c \sim 3$ Hz, we can express the stress-drop ratio for two events of seismic moment $M_{01}$ and $M_{02}$:

$$ \begin{gathered} \Delta\sigma \propto M_0 \\ \Rightarrow \quad \frac{\Delta\sigma_1}{\Delta\sigma_2} = \frac{M_{01}}{M_{02}} \end{gathered} \qquad (2.13) $$

Using the extrema of the seismic moment distribution, and within the limit of our assumptions, we observe a stress-drop variation of a factor 300 across the range of measured seismic moments. For illustration purposes, we calculate what values would the stress-drop of our events take under the assumption that they have the same rupture speed as classical earthquakes, $v_r = 0.9V_S$. It yields $r = 0.21V_S/f_c$ (Madariaga, 1976), and using equation 2.11, the selected events have stress-drop ranging from 3 kPa to 1 MPa, with most events around 17 kPa ($M_0 \sim 2 \times 10^{12}$ N.m). In the following section, we discuss a possible mechanism for the observed range of stress-drops.

### 2.5.4 Moment-frequency distribution of LFEs

In Figure 2.7, we represent the complementary cumulative distribution function (abbreviated as CCDF, also called survival function) of seismic moments, $P(M_0 > M_0')$. It is computed with $M_0'$ regularly spaced in the magnitude range.

The higher-moments tail of the distribution can be described as a power-law distribution, with the following probability density function (PDF) and
<figure><img src="image_25.png" /><figcaption>**Figure 2.7** — Moment-frequency distribution of the selected LFEs. The seismic moment complementary cumulative distribution function (CCDF) of our LFEs is modeled by a Gaussian below $M_{wthr}$ 2.32, with mean magnitude $M_w$ 2.4, and by a power law above $M_{wthr}$, with an exponent corresponding to a b-value of 2.4.</figcaption></figure>

CCDF:

$$ 
\begin{align} 
p(M_0) &= \frac{\beta}{M_{0thr}} \left( \frac{M_0}{M_{0thr}} \right)^{-\beta-1} \tag{2.14} \\ 
P(M_0 > M_0') &= \left( \frac{M_0'}{M_{0thr}} \right)^{-\beta} \nonumber 
\end{align} 
$$

where $M_{0thr}$ is the lower bound for which the power-law distribution is valid. The power law exponent $\beta$ relates to the b-value of the Gutenberg-Richter law following: $b = 3/2\beta$. For moments above a threshold $M_{0thr}$, we perform an un-weighted, least-square fit of the logarithm of the moments CCDF to a power-law. We find that $\beta = 1.62$, and a corresponding b-value of $b = 2.43$. For comparison, b-values of classical seismicity are usually around 1. The threshold magnitude $M_{wthr} = 2.32$ (corresponding to $M_{0thr}$) is selected so as to minimize the Kolmogorov-Smirnov statistic for the moments' distribution and the fitted power law.

The deviation from a power law at lower moments is often attributed to observational limitations: below a completeness magnitude, events are too weak to all be detected, and it manifests as a deviation from the otherwise assumed power law distribution of moments. Here, the completeness mag-
nitude would thus be $M_{wthr} = 2.32$. The departure from a power-law at lower
moments can also be the manifestation of a characteristic scale of magnitude
of the events. Several studies report observations of LFEs size distribution
being best modeled by exponential laws and thus of a characteristic scale con-
trolling the LFE size distribution (Chamberlain et al., 2014; Chestler & Crea-
ger, 2017a). In our case, a normal distribution can be fitted to the distribution
of moments below $M_{wthr}$, parametrized by a mean magnitude $\bar{M}_w = \mu = 2.23$
and a standard deviation $\sigma = 0.19$. The mean magnitude $\bar{M}_w = 2.23$ defines
a characteristic scale of the lower-bound of the moment distribution. It is ei-
ther due to a scale limitation of LFEs in Guerrero, or a manifestation of the
limitation of events' detection.

## 2.6 Discussion and interpretation of the observed scaling law

2.6.1 Scaling laws of earthquakes, slow earthquakes and LFEs

The results of our analysis show that the selected LFEs follow a moment-
duration scaling $f_c \propto M_0^{-\alpha}$ with an exponent significantly lower than the
value for regular earthquakes $\alpha = 1/3$ (Allmann & Shearer, 2009; Houston,
2001; Prieto et al., 2004, among others), around $\alpha \sim 0.1$ and $\alpha \sim 0.05$. The
scaling found for selected LFEs in Guerrero concurs with the one found for
LFEs in Cascadia by Bostock et al. (2015), who found a power law exponent
of $\alpha = 1/10$, quite similar to the scaling described in this paper. It is note-
worthy that the detection method used by Bostock et al. (2015) relies on a
matched-filter search, similarly to the catalog used in our study. However,
we argue that performing a matched-filter search in the 1–2 Hz band sets a
lower bound, but no upper bound for the characteristic frequency, for a given
seismic moment. Therefore, the values of corner-frequencies should not be
the product of the detection method, but reflect the real moment-corner fre-
quency distribution of events, within the measured magnitude range. This
argument is developed in appendix 2.A. Interestingly enough, Supino et al.
(2020) found the moment-duration scaling for LFEs in Nankai tremor areas to
be similar to regular earthquakes, $\alpha = 1/3$. This discrepancy will not be thor-
oughly discussed here, but brings out intriguing questions, about the possibly
different nature of the processes recorded in Guerrero, Cascadia and Nankai,
and on the other hand, the potential detection and method biases that can
lead to measuring different scalings.

Using a stochastic model describing a slow earthquake rupture (the Brow-
nian slow earthquake model), Ide (2008); Ide and Maury (2018) showed that
large-scale slow earthquakes as slow slip events' duration and seismic mo-
ment should follow a $T \propto M_0$ scaling, whereas on the scale of LFEs, slow
earthquake events should follow a $T^2 \propto M_0$ scaling, corresponding a value
of $\alpha = 1/2$. In Ide (2008), it is noted that, due to the fact that detection is only
possible above noise levels, the apparent duration of LFEs should be shorter
than expected, bringing them closer to a moment-duration scaling with $\alpha = 1$.
The LFEs we analyzed in Guerrero and the LFEs of Cascadia (Bostock et al.,
2015) cannot be described by either scaling laws.

There is evidence that the shortest slow slip transients measured to this
day are self-similar, exhibiting a $T^3 \propto M_0$ scaling law (Frank & Brodsky,
2019). The exact proportionality measured for larger scales slow earthquakes
(Ide, Beroza, et al., 2007) could be due to the fact that the largest ruptures
are bounded in their growth by the rheological properties of the fault zone,
whereas smaller events growth is un-bounded, and thus self-similar (Gomberg
et al., 2016). This geometric argument is valid both for slow and regular
earthquakes, and further isolates LFEs moment-duration characteristic val-
ues and scaling from the ones of slow and regular slip. For LFEs in Cascadia
and Guerrero, a satisfactory physical description that could reproduce the ob-
servation of the constant duration and relatively low-frequency content still
needs to be devised.

**2.6.2 Models for magnitude-independent duration of LFEs**

2.6.2.1 A source of constant size, but variable slip

The scaling with $\alpha = 1/3$ for regular earthquake is explained by the self-
similarity of the seismic rupture. It develops with a constant ratio between the
characteristic length of the rupture and fault slip, and produces a magnitude-
independent stress drop (Allmann & Shearer, 2009; Kanamori & Anderson,
1975; Prieto et al., 2004). Events with a wide range of magnitudes but with
nearly constant durations could thus be generated by variable amounts of
slip on asperities of nearly constant sizes, producing a very limited range of
recorded signal duration ($f_c$) but variable seismic moments (Bostock et al.,
2015; Lengliné et al., 2014; Lin et al., 2016). Independently from the current
study, Frank and Brodsky (2019) found that, in Guerrero, the seismic mo-
ment rate, calculated from LFEs median amplitude during rapid slow slip
transients scales with the moment rate calculated using GPS displacements
during the same transient events. This reflects that the higher the seismic
moment of LFEs, the higher the slip would be around the LFEs source. Com-
bining those two observations leads to thinking of LFE sources in Guerrero as
patches of fault of nearly constant area, that are forced to rupture by the sur-
rounding slow, aseismic slip, with variable amplitudes of slip directly related
to the surrounding fault slip. A similar model has been designed to explain the observation of a size-limitation for LFEs in Cascadia, describing the LFE family patch as a collection of subpatches randomly activated, totaling to an apparent source radius of about 300 m (Chestler & Creager, 2017a).

### 2.6.2.2 Potential effects of fluid circulation and high fluid pressures

Early work on tectonic tremor and LFEs made the hypothesis that they are generated by fluids circulating in the vicinity of the fault zone (Kao et al., 2005; Katsumata & Kamaya, 2003; Obara, 2002). More recent studies have strengthened the link between non-stationnary, metamorphic fluid flow and migrations of LFEs and slow slip activity, based on numerical models of pore-pressure diffusion (Cruz-Atienza et al., 2018) and on analysis of the large-scale, spatio-temporal behavior of LFE activity (Frank, Shapiro, et al., 2015).

A study by Lengliné et al. (2014) reported that microseismic events occurring during water circulation tests display the same independence of corner frequency on seismic moment that we observe for our LFEs. As for the LFEs in this study, they find that stress drop varies on more than two orders of magnitude across the range of moments of their events. Likewise, Lin et al. (2016) reports that microseismic events detected in Taiwan exhibit the same characteristics. Both studies suggest that these observations can be interpreted as a systematic control of the events’ stress-drop and duration by dynamic fluid injection. Indeed, fluid injection can reduce the effective normal stress on asperities of the fault where the events occur. As fluid is injected or when aseismic slip is active, rapid, strong, heterogeneous increase of pore pressure in the fault can dramatically decrease the effective normal stress and thus the strength of the fault, triggering seismic slip on localized asperities. Under these conditions, Lengliné et al. (2014) and Lin et al. (2016) argue that the heterogeneity of fluid injection is likely to produce the wide range of stress-drops they observe. As noted by the previous papers, reduced effective normal stress have also been suggested as a cause of the low values of stress drop of VLFEs in Japan (Ito & Obara, 2006). In the light of the similarity of our observations and the ones presented in the discussed studies, dynamic triggering of LFEs by fluid injection could be a relevant mechanism to explain the LFEs’ stress drop values.

In order to account for the observation of events’ constant duration, Lin et al. (2016) proposes that slip could only be maintained while the pore pressure is high enough to allow for seismic rupture, that is until the injection stops and the high pore pressure transient diffuses. The duration of the recorded seismic event would thus be controlled by the fluid transport properties in the fault zone. Those conditions of rapid, localized fluid injection are plausible within large active fault zones (N. M. Shapiro et al., 2018). In subduction
zones, fluids released by the dehydration of hydrous minerals in the down-
going slab circulate within a permeable channel formed along or within the
fault zone (Angiboust et al., 2012; van Keken et al., 2011). If the permeability
or width of the channel has strong lateral heterogeneities, N. M. Shapiro et al.
(2018) suggests that strong, localized pore pressure transients can be formed.
The strong, transient gradients of pore pressure can act as a source of seismic
waves. The duration of such a source would thus be the characteristic time
taken by a pore pressure transient to diffuse. N. M. Shapiro et al. (2018) shows
that as with such a mechanism, this duration is only governed by the trans-
port properties of the fault zone and independent of the source magnitude.
Although more precise observations of the source radiation pattern of LFEs
and a substantial modeling effort are needed to confirm this hypothesis, it is
consistent with our observations.

Another lead to account for the apparent magnitude-independent, low-
frequency content of LFEs is to interpret it as the consequence of specific,
near-source propagation effects (Bostock et al., 2017). In this study the au-
thors assume the existence of a narrow region with a very strong shear wave
attenuation around LFEs source and suggest that the observed depletion in
high frequencies can be related to this strong near-source attenuation of seis-
mic waves. Based on accounts of high compressional to shear wave velocity
ratio V<sub>P</sub>/V<sub>S</sub> in seismic cross-section of subduction zones (Bostock et al., 2012;
Shelly et al., 2006), tidal triggering of LFE activity (Beeler et al., 2013; Royer
et al., 2015; Rubinstein et al., 2008, e.g.) and dynamic triggering by teleseis-
mic waves (Rubinstein et al., 2009), the source region of LFEs is thought to
be under nearly lithostatic fluid pressures. Laboratory studies indicate that
such high fluid content can also be the cause of strong shear wave attenuation
(Tompkins & Christensen, 2001).

2.7 Conclusions

We analyzed a subset of carefully selected LFEs from the cataloged events
of Guerrero, Mexico (Frank et al., 2014). The subset of LFEs samples a rela-
tively broad range of seismic moments, between M<sub>w</sub>1.7 and M<sub>w</sub>3, with corner-
frequencies clustering around 3 Hz, between 1.5 and 6 Hz. Overall, we find
corner frequencies to be very weakly dependent on seismic moments for LFEs
in Guerrero. When trying to describe this relationship with the power law
f<sub>c</sub> &prop; M<sub>0</sub><sup>-&alpha;</sup>, we retrieve &alpha; = 0.12 &pm; 0.11 and &alpha; = 0.054 &pm; 0.034 for the two sets
of corner frequency measurements. In both cases those scalings are signifi-
cantly different from a self-similar rupture, for which &alpha; = 1/3.

Recent studies show that the clear departure from self-similarity mea-
sured for LFEs in Guerrero and Cascadia is not a necessary feature of LFEs in
all tremor regions across the world. Exploring this regional specificity might
bring significant insight on LFEs source physics and lead us to understand
more about how deep fault environment affects slow slip processes.
# APPENDICES TO CHAPTER 2

## Contents

<table><tr><td>2.A</td><td>Effects of bandpass filtering on events corner frequencies</td><td>47</td></tr><tr><td>2.1.1</td><td>Analysis bandpass filter</td><td>47</td></tr><tr><td>2.1.2</td><td>Detection bandpass filter</td><td>51</td></tr><tr><td>2.B</td><td>Effects of attenuation on corner frequency f<sub>c</sub> measurements</td><td>52</td></tr><tr><td>2.2.1</td><td>Seismic attenuation in the crust of Guerrero</td><td>53</td></tr><tr><td>2.2.2</td><td>Effects of a highly-attenuating fault zone on LFE spectra</td><td>56</td></tr></table>

## 2.A Effects of bandpass filtering on events corner frequencies

### 2.1.1 Analysis bandpass filter

In Figure 2.A.1, we see that the tremor seismic energy seems really bounded between 0.8 Hz and 8 Hz. Above this band, high-frequency noise, possibly anthropogenic, is active, below this band, the seismic time-series are dominated by low-frequency noise leaking from microseismic modes. The choice of the 1–8 Hz band-pass filter used to measure corner-frequencies is motivated by this observation.

It is possible that the 1–8 Hz band-pass filter we use biases our measure of the corner frequencies to constrain them in a narrow band above 1 Hz. We tested the influence of the lower cutoff frequency of the band-pass filter on the measured corner frequency. Results of this test are shown in Figure 2.A.2. The distribution of corner frequencies is not drastically modified down to a lower cutoff of 0.2 Hz, where corner frequencies suddenly regroup around 0.2 Hz. This sudden shift is due to the filter band not blocking pollution of micro-seismic noise anymore. In addition to it, corner frequencies are slightly
dragged to lower values, as the lower cutoff decreases from 1 Hz to 0.5 Hz. This should be expected: more and more pollution from micro-seismic noise is allowed in records of events on stations that witness higher amplitudes of micro-seismic noise, due to proximity to the ocean or site effects. On those records, the characteristic frequency will be picked lower, due to the influence of micro-seismic noise. Figure 2.4d is a good illustration of this: the peak of the un-filtered spectrum is located at 1 Hz due to micro-seismic noise influence.

The 1–8 Hz band is thus chosen as a compromise. The lower cutoff is taken far enough from micro-seismic frequencies so that the event signal is not polluted, even at stations where micro-seismic noise has the highest amplitude. Ultimately, it does constrain the corner frequency between 1 and 8 Hz: most events have corner frequencies between 2 and 6.5 Hz, well within the filtering band.
<figure><img src="image_26.png" /></figure>

**Figure 2.A.1** — The tremor signal can be observed in the 1–8 Hz frequency, as signal amplitude variations correlated to the hourly rate of event detection. It is clear that there is a nearly constant background noise below 0.8 Hz and above 8 Hz, with an amplitude on the same order as the tremor signal.
<figure><img src="image_27.png" /></figure>

**Figure 2.A.2** — Effects of varying the lower cutoff of the filter are presented for both methods of determination of the corner frequency. In the left column, for the mean instantaneous frequency method, In the right column, for the maximum of velocity spectrum method. The filter band used in the article is 1–8 Hz. For comparison, location of the peak of the current distribution is shown as the black dashed line. Decreasing the lower cutoff of the filter slightly brings the distribution of corner frequencies to lower values, down to 0.5 Hz. Below 0.5, as seen for 0.2 Hz, micro-seismic noise is picked up as characteristic frequency of the signal and pulls the whole distribution to low values.
### 2.1.2 Detection bandpass filter

The LFEs of the catalog used in this study were detected using a matched-filter search (or MFS) approach, using template velocity waveforms, band-pass filtered between 1–2 Hz (Frank et al., 2014). A potential concern is that a MFS using a narrow band-pass would only identify events with a characteristic frequency between 1 and 2 Hz. Assuming it is the case, it would bias our analysis towards the observed magnitude-independent corner frequencies.

However, the use of an MFS approach does not constrain the characteristic frequencies of the event between 1–2 Hz, it simply limits detections to signals with energy within this frequency band. Several studies using MFS to detect classic earthquakes are able to detect events with a range of magnitudes, that, through the earthquake moment-duration scaling law, corresponds to a range of characteristic frequency wider than the chosen band-pass (Warren-Smith et al., 2017, Figure 7.(b)), (Frank & Abercrombie, 2018, Figure 5). In our case most events indeed have corner frequencies above 2 Hz and up to 6.5 Hz (see Table 2.1, Figure 2.5).

A compact form of the Brune spectrum for ground motion velocity can be used to compare the detectability of two events: $\tilde{v}(f)$:

$$ \tilde{v}(f) = C \cdot M_0 \frac{2\pi f}{1 + (f/f_c)^2}, \qquad (2.15) $$

where C is a constant relating to the velocity structure of the medium and the source focal mechanism, $M_0$ the seismic moment of the event and $f_c$ its corner frequency. For any two events with $f_{c1} < f_{c2}$ but identical moment, for any given frequency $f$, the corresponding velocity amplitudes will follow $\tilde{v}_1(f) < \tilde{v}_2(f)$. In Figure 2.A.3, we represent the velocity spectrum $\tilde{v}(f)$ of 5 events with various corner frequencies and the same seismic moments. In Figure 2.A.3b, we show the spectral velocity energy in the 1–2 Hz band for events with various corner frequency, but the same seismic moment. This figure shows that in this case, the higher the corner frequency of the event, the higher the amplitude of velocity between 1–2 Hz, and therefore the more detectable the event by the MFS. Events with a higher characteristic frequency than events we detected thus do not have a lower probability to be detected because of the band-pass filter. However, events with a lower corner frequency than events we detect have a lower probability of detection, because the energy of their velocity waveform is lower in the 1–2 Hz band (Fig. 2.A.3). This effect is especially strong for events with $f_c < 1$ Hz. The matched-filter search band-pass thus imposes a lower bound of about 1 Hz on detections’ corner frequencies in the range of moments we measure, but no higher bound, except the one imposed by the filter used in the analysis, at 8 Hz.
<figure><img src="image_28.png" /><figcaption>Effects of a 1–2 Hz detection on the detectability of LFEs of various $f_c$</figcaption></figure>

**Figure 2.A.3** — (a) Seismic velocity spectra of 5 events with different corner frequencies but the same seismic moment. Corner frequencies are shown with a cross and dashed line. (b) Spectral velocity energy in the 1–2 Hz band for events with various corner frequencies, and same seismic moment. The “energy” of the velocity waveform is defined as the integral of spectral velocity squared, over the 1–2 Hz band. The higher the energy, the more the event should be detectable. In both panels, the 1–2 Hz frequency band used to detect LFEs is shaded in gray. All events used here are described by a Brune spectrum (Brune, 1970).

The 1–8 Hz band-limitation that the MFS should impose should allow to detect significant variation of $f_c$ for the range of seismic moments we measure ($8 \times 10^{11}$–1 × $10^{13}$ N.m) — a 1–10 Hz band limitation was enough for Supino et al. (2020) to measure strong variations of $f_c$, in a similar range of moments. Therefore, we conclude that although the combination of the lower bound on detected $f_c$ inherited from the 1–2 Hz MFS and the upper bound imposed by the analysis low-pass at 8 Hz creates a band-limitation for the detected corner frequency, it is not narrow enough to bias the measurements of the $M_0-f_c$ scaling, so that it is an artifact of band-limited detection.

## 2.B Effects of attenuation on corner frequency $f_c$ measurements

The spectral characteristics of a seismic signal are shaped both by the properties of the source — its duration, the amount of radiated energy, the variations in the rate of energy release at the source — and the effects of the medium on the signal that propagates from the source — geometric spreading, radiation pattern, attenuation through scattering and anelastic energy dissipation.
In particular, measures of the source's corner frequency trade-off with how the propagating medium attenuates seismic energy at high frequencies. A strongly attenuating medium — for instance in the highly heterogeneous, uncosolidated subsurface — can strongly attenuate the high-frequency part of a seismic signal and produce an apparently-constant corner frequency, independent of the events seismic moment (Abercrombie, 1995). In this section, we discuss our choice of attenuation model for the Guerrero crust, and how an unaccounted, strongly-attenuating source region could be the cause of an apparently-constant corner frequency for LFEs in Guerrero.

## 2.2.1 Seismic attenuation in the crust of Guerrero

The effects of attenuation on the spectrum of a seismic event can be modeled by an attenuation term $a(f)$ that multiplies the true source spectrum $\tilde{u}(f)$ to produce the measured, attenuated spectrum $\tilde{u}(f)$ (J. G. Anderson & Hough, 1984):

$$ \tilde{u}(f) = \tilde{u}(f) \cdot a(f). \qquad (2.16) $$

The attenuation term $a(f)$ depends on frequency, usually with an exponentially stronger attenuation at higher frequencies:

$$ a(f) = \exp \left( - \frac{\pi f R}{V_S Q(f)} \right), \qquad (2.17) $$

where $R$ is the source-receiver distance, and $Q(f)$ and the seismic quality factor. In order to measure the true corner frequencies $f_c$ of LFEs, we divide the measured spectrum $\tilde{u}(f)$ of selected events by the modeled effects of attenuation $a(f)$, to obtain the corrected “true” spectrum $\tilde{u}(f)$, in which only the effects of source, radiation pattern and geometric spreading should theoretically remain. Corner frequencies $f_c$ are then measured on $\tilde{u}(f)$.

Several studies found that attenuation in the crust in Guerrero is weakly frequency dependent, and can be described by a quality factor in the form $Q(f) = Q_0 f^\eta$, with $\eta$ lower than, but close to 1 (Castro et al., 1990; Mahdyiar, 1984; Ordaz & Singh, 1992; Ottemöller et al., 2002; Singh et al., 1990). When $\eta = 1$, $Q = Q_0 f$, and $a(f)$ becomes frequency independent (see equation 2.17). Humphrey and Anderson (1992) found that attenuation in the Guerrero crust can be described by such a frequency-independent parametrisation, consistently with the previously mentioned studies. Figure 2.B.1 shows examples of frequency-dependent models for the quality factor $Q$ in Guerrero, measured using S-waves for most of the models. Those attenuation models are quite similar, our choice of attenuation model — where $Q = 273 f^{0.66}$ (Ordaz & Singh, 1992) — should therefore adequately describe the attenuation in the continental crust of Guerrero.
<figure><img src="image_29.png" /></figure>

**Figure 2.B.1** — Models of attenuation for the crust of Guerrero. All modeled frequency-dependent quality factors are measured on the S-waves of regional earthquakes, which wave travel through the crust of Guerrero, except for the model by Ottemöller et al. (2002), who measured it on the Lg phase.

As values of $\eta$ close to 1 imply that $a(f)$ is nearly frequency-independent, the measured, attenuated spectrum of LFEs in Guerrero and their source spectrum should have the same spectral shape, only varying in amplitude (see equation 2.16). Therefore small imprecisions in the attenuation modeling should have a negligible effect on the measured corner frequency, assuming that the continental crust is the main attenuating layer between LFEs' sources and the instrument at the surface. Figure 2.B.2 shows the effect of various attenuation models on the spectra and measured $f_c$ of events of different seismic moment $M_0$ and corner frequency $f_c$. For all expressions of the quality factor for the crust in Guerrero, the way attenuation modifies the true spectrum are very similar. It is also clear that such crustal attenuation does not have a strong impact on the spectral shape between 1–10 Hz, and therefore on the apparent corner frequency. Our choice of model for the attenuation of the crust can therefore be excluded as being responsible for the $M_0-f_c$ scaling we obtain. The details of calculations to obtain Fig. 2.B.2 are described in the next section.

Measures of attenuation in other subduction zones show that the values found in Guerrero — $Q \approx 300$ at 1 Hz, $Q = 1000$ at 8 Hz — are consistent with the high values of attenuation found in the continental crust found of tremor areas. Using earthquake-receiver paths that sample the whole continental crust and parts of the subduction fault under Vancouver Island (Cascaedia) Atkinson (1995) and Shikoku (Japan) Takahashi et al. (2014), studies
Source-receiver geometry, and propagation model

a

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

Ordaz and Singh (1992) $Q_s = 273 f^{0.66}$

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

Castro et al. (1990) $Q_s = 278 f^{0.92}$

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

Singh et al. (1990) $Q_s = 100 f$

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

Ottemoller et al. (2002) $Q_{Lg} = 193 f^{0.87}$

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

Mahdiyar et al. (1984) $Q_s = 87 f^{0.78}$

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

**Figure 2.B.2** — Effects of attenuation on the spectra and apparent corner frequency of earthquakes of different magnitudes and corner frequencies. (a) Source-receiver geometry and physical properties of the source and propagating medium, here only composed of a model continental crust. (b)–(f), True and attenuated spectra for different expressions of the quality factor. The true corner frequency of the model event is symbolized with a cross, and the apparent, attenuated corner frequency is represented as a circle. See section 2.2.2 for more details on the modeling.
respectively found $Q(f) = 263 f^{0.49}$ and $Q(f) = 315 f^{0.5}$. Those laws are sensibly similar to the ones measured in Guerrero. Other studies measure the attenuation using tremor waveforms, and find that at $f = 4.5$ Hz, the full-path attenuation can be described by a quality factor of $Q \approx 500–800$ in Shikoku, Jalisco (Mexico) and Cascadia (Baltay & Beroza, 2013; Yabe et al., 2014). Those results are especially relevant to our case, as they measure the attenuation on the typical path travelled from LFE source to receiver. The fact that these measures are consistent between subduction zones tend to show that the crustal attenuation model should be adequate for Guerrero.

### 2.2.2 Effects of a highly-attenuating fault zone on LFE spectra

In studies by Z. Wang and Zhao (2019) and Yabe et al. (2014), the presence of intraslab earthquakes in addition to other seismic event at the plates' interface and within the crust allow the authors to measure the attenuation within the fault zone, where LFEs are located. Both studies report that the fault zone should behave as a few-kilometer-thin, highly-attenuating layer, with $Q \approx 100$ for a thickness of 5 km (Yabe et al., 2014). In the next paragraphs, we investigate the effects of a thin, highly-attenuating layer located just above the LFE sources, and how it could bias corner frequency measurements towards the observed scaling. Such strong attenuation of the high frequency content of a signal can be responsible of an apparent constant corner frequency for events of different sizes, as pointed out by Abercrombie (1995), who showed that the observation of a constant $f_c$ for small events in California by Archuleta et al. (1982) was likely due to a highly attenuating subsurface ($Q \approx 5–30$ in the shallowest few hundreds of meters). A similar reasoning could be made for small earthquakes in different parts of the world B. Chouet et al. (1978). We will use this reasoning here to infer what value of the quality factor in the fault zone is necessary to obtain an apparent $M_0-f_c$ scaling similar to the one we measure for LFEs in Guerrero. This approach was developed by J. G. Anderson et al. (1986).

The hypothesis we want to test is that a population of small, self-similar earthquakes which waves propagate in a thin, strongly attenuating fault-zone could exhibit an apparently-constant corner frequency $f_c$, independent of seismic moment $M_0$. If values of the quality factor necessary for such a scaling to occur are realistic, the scaling we measure for LFEs could be an artifact due to the un-modeled fault-zone attenuation in our study. The population of small earthquakes are located directly below the receiver, at a depth $R$. They are self-similar, and their true, un-attenuated displacement spectrum $\tilde{u}(f)$ is
described by a simple Brune's spectrum (Brune, 1970):

$$0\tilde{u}(f) = \frac{0.84M_0}{4\pi\rho V_S^3 R} \times \frac{1}{1 + (f/f_c)^2}. \qquad (2.18)$$

For a self-similar rupture, the corner frequency depends on the stress drop $\Delta\sigma$ and seismic moment $M_0$ as

$$f_c = 0.49 \left( \frac{\Delta\sigma}{M_0} \right)^{1/3} V_S. \qquad (2.19)$$

The population of small earthquakes we use is defined by seismic moments $M_0 = 1e8-1e16$ N.m, a stress drop of $\Delta\sigma = 4$ MPa (median stress drop in Allmann & Shearer, 2009), and true corner frequencies $f_c = 1-500$ Hz. Following J. G. Anderson and Hough (1984), we account for attenuation in the crust and fault-zone with a multiplicative term for each, to obtain the measured spectrum $\tilde{u}(f)$:

$$\tilde{u}(f) = 0 \tilde{u}(f) a_C(f) a_F(f). \qquad (2.20)$$

where both the crustal attenuation factor $a_C(f)$ and fault attenuation factor $a_F(f)$ are described by an exponential term (see equation 2.17). The “apparent” corner frequency is defined as the maximum of the measured, attenuated velocity spectrum $\tilde{u}(f)$.

We first test the effect of crustal attenuation without a highly-attenuating fault zone, and set $a_F(f) = 1$. Figure 2.B.2a describes the source receiver geometry and the velocity and attenuation structure. The source is located at $R = h_C = 40$ km depth, below a homogeneous crust with $V_S = 3.5$ km.s$^{-1}$, $\rho = 2700$ kg.m$^3$ and a frequency-dependent quality factor $Q(f) = Q_0 f^\eta$. Panels (b)–(f) show the true and attenuated spectra for events with different $M_0$, $f_c$. In Guerrero, the crustal attenuation does not strongly affect the spectral shape, and should not have a strong effect on the apparent corner frequency in the 1–10 Hz frequency band we are interested in for the study of LFEs.

In a second experiment, we add a $h_F = 1$ km-thick fault zone above the source and below the $h_C = 40$ km-thick continental crust. Both layers share the same density and shear wave velocity. Crustal attenuation is described by our choice of model $Q(f) = 273 f^{0.66}$ (Ordaz & Singh, 1992). Fault attenuation is described by a frequency-independent quality factor $Q_F$. Figure 2.B.3a describes the simulation set up. We test values of the fault-zone quality factor $Q_F = 1-1000$, and for each we obtain the distribution of apparent corner frequency $f_c$ and seismic moment $M_0$. Figure 2.B.3c shows that for a relatively low quality factor in the fault zone $Q_F = 100$ (similar to what Yabe et al., 2014, found), earthquakes of the same magnitude as the LFEs in Guerrero
$(M_0 \approx 5e11-1e14 \text{ N.m})$ have a higher apparent corner frequency than LFEs ($f_c \approx 10-30 \text{ Hz}$), and their $M_0-f_c$ distribution is very close to their true value, where $f_c \propto M_0^{-1/3}$. Only for an extremely low quality factor $Q_F = 3$ — very high attenuation — the effect of attenuation on the spectrum are significant between 1–10 Hz, and therefore could strongly affect the shape of the spectrum, biasing the measured corner frequency (Figure 2.B.3b and c). In this case, if the attenuation of the fault zone is not corrected for in the measurement of $f_c$, the apparent $M_0-f_c$ scaling resembles the one we measure for LFEs in Guerrero. Bostock et al. (2017) found with a similar calculation that the $M_0$-independent corner frequency scaling of LFEs under Vancouver Island could result from the attenuation of a 1 km-thick fault-zone with $Q_S = 3$.

A value of $Q_F = 3$ for S-waves seems unrealistically low, as even the highly-attenuating, sedimentary subsurface rarely goes down to $Q = 5$ (Abercrombie, 1995). More than the absolute value of $Q_F$, its frequency independence is responsible for the apparent scaling, as it produces a steep high-frequency attenuation: $Q_F = Q_0 f^0 \implies a_F(f) = \exp(-Cf)$, where C is a constant. A very highly-permeable ($k \ge 1e-12 \text{ m}^2$), fluid-saturated layer could generate strongly frequency-dependent attenuation, with a low $Q_0 \approx 100$ and $\eta \approx 0$ or even $\eta < 0$ (Müller et al., 2010; Takahashi et al., 2014). If the subduction fault zone can be described as such a highly permeable layer, its attenuation could be responsible for the band-limited nature of tremor and LFEs. However, the extraordinary nature of the attenuation needed to explain this scaling as resulting from propagation effects tends to show that the scaling we measure can be reasonably interpreted as a result of the source properties of LFEs. In the end, more research remains to be done on the attenuation characteristics of the subduction fault-zone, and the consequences of such a high permeability, to assess if propagation effects could realistically be responsible for the apparent $M_0$-independent corner frequency.
<figure><img src="image_36.png" /><figcaption>**Figure 2.B.3** — (a) Source-receiver geometry and physical properties of the source and propagating medium, composed of a model continental crust and fault zone. (b) True and attenuated spectra for events of different sizes and magnitudes, for a quality factor in the fault zone of $Q_F = 3$. The true corner frequencies are symbolized with crosses, and the apparent, attenuated corner frequencies are represented as circles. (c) Relationship between seismic moment $M_0$ and apparent corner frequency $f_c$ for different values of the quality factor in the fault zone $Q_F$. In (b) and (c), the $M_0-f_c$ distribution for LFEs in Guerrero (present study) is represented approximately in purple.</figcaption></figure>
<figure><img src="image_37.png" /></figure>
# Part III

## Spatio-temporal patterns of tremor activity
<figure><img src="image_38.png" /></figure>
# CHAPTER 3

## TREMOR GENERATION THROUGH FAST FLUID PRESSURE TRANSIENTS IN THE PERMEABLE SUBDUCTION INTERFACE

Based on: G. FARGE, C. JAUPART, N. M. SHAPIRO, *Episodicity and Migration of Low Frequency Earthquakes Modeled With Fast Fluid Pressure Transients in the Permeable Subduction Interface*. Journal of Geophysical Research: Solid Earth, 126(9), September 2021. ISSN 2169-9313, 2169-9356. https://doi.org/10.1029/2021JB021894

## Contents

<table><tr><td>3.1</td><td>Introduction</td><td>64</td></tr><tr><td>3.2</td><td>Tremor patterns observed in subduction zones</td><td>68</td></tr><tr><td>3.2.1</td><td>Spatial Segmentation</td><td>69</td></tr><tr><td>3.2.2</td><td>Tremor Migration</td><td>70</td></tr><tr><td>3.2.3</td><td>Temporal Patterns</td><td>70</td></tr><tr><td>3.2.4</td><td>Tremor activity Patterns in Guerrero, Mexico</td><td>71</td></tr><tr><td>3.3</td><td>Model design</td><td>71</td></tr><tr><td>3.3.1</td><td>Diffusion equation for pore pressure variations</td><td>72</td></tr><tr><td>3.3.2</td><td>A mechanism for large changes of permeability</td><td>73</td></tr><tr><td>3.3.3</td><td>A valve mechanism</td><td>75</td></tr><tr><td>3.3.4</td><td>Valves as elementary seismic sources</td><td>77</td></tr><tr><td>3.3.5</td><td>Permeability values, boundary and initial conditions</td><td>77</td></tr><tr><td>3.3.6</td><td>Numerical implementation</td><td>80</td></tr><tr><td>3.4</td><td>Time-dependent valve behaviour</td><td>80</td></tr><tr><td>3.4.1</td><td>A single valve</td><td>80</td></tr><tr><td>3.4.2</td><td>Valve-valve interaction via pore pressure diffusion</td><td>83</td></tr><tr><td>3.5</td><td>The collective behavior of interacting valves</td><td>84</td></tr></table>
<table><tr><td>3.5.1</td><td>Characteristics of valve activity in the permanent regime</td><td>84</td></tr><tr><td>3.5.2</td><td>Cascades and migrations of synthetic activity</td><td>87</td></tr><tr><td>3.5.3</td><td>Large-scale valving behavior</td><td>91</td></tr><tr><td>3.5.4</td><td>Random valve distributions</td><td>91</td></tr><tr><td>3.6</td><td>Discussion</td><td>92</td></tr><tr><td>3.6.1</td><td>Comparison with LFE activity patterns in Guerrero</td><td>92</td></tr><tr><td>3.6.2</td><td>Scaling numerical results to LFE activity at Guerrero</td><td>93</td></tr><tr><td>3.6.3</td><td>Perspectives for diagnosing hydraulic conditions in fault zones</td><td>95</td></tr><tr><td>3.6.4</td><td>Perspectives for exploring model configurations</td><td>96</td></tr><tr><td>3.6.5</td><td>Coupling between hydraulic transients and fault slip</td><td>97</td></tr><tr><td>3.7</td><td>Conclusion</td><td>98</td></tr></table>

## 3.1 Introduction

The extensive records of weak low-frequency seismicity detected in many fault zones over the last two decades provide powerful information on their dynamics. At intermediate depths (30–50 km) in several subduction zones, seismic activity takes the form of tremor (e.g., A. L. Husker et al., 2012; Kao et al., 2005; Obara, 2002; Payero et al., 2008). It now seems clear that tremor is made of swarms of low-frequency earthquakes (LFEs) and many authors use the term “tremor” as a synonym for LFE activity (e.g., Audet & Kim, 2016; Bostock et al., 2012; Frank et al., 2014, 2013; Peacock et al., 2011; Rogers, 2003; Shelly et al., 2006; Song et al., 2009). These LFEs appear to be due to repetitive, impulsive and localized sources (Beroza & Ide, 2011; Frank et al., 2014). High rates of tremor activity are associated with geodetically observed slow deformation episodes called slow-slip events (SSEs) (e.g., Delbridge et al., 2020; Frank, Radiguet, et al., 2015; Kostoglodov et al., 2010; Rousset et al., 2019). The correlation between LFE and SSE activity gave rise to the term Episodic Tremor and Slip (ETS) (Rogers, 2003). In a more general acceptation, transient seismic and geodetic signals that are both generated during inter-seismic phases are often lumped together as slow earthquake phenomena. This new class of seismic and geodetic observations holds great promise for improving our understanding of the seismo-tectonic cycle.

According to the framework that prevails today, slow earthquakes are the “slow” analogs of “regular” earthquakes. Geodetic and seismic events are manifestations of the same process of slow slip at transitional depths in fault zones, in between unstable stick-slipping and stable sliding areas.
The transitional part is made of slow sliding patches that are responsible for geodetically detected SSEs and that contain many small embedded frictionally unstable asperities generating LFEs. This framework explains why SSEs and tremor activity are correlated in time and accounts for the S-wave dominated LFE radiation pattern (e.g., Shelly et al., 2007b). Numerical models of these phenomena (e.g., Barbot, 2019; Ben-Zion, 2012; Ide, 2014; Liu & Rice, 2007; Luo & Liu, 2019; Shibazaki & Iio, 2003) are often based on the rate-and-state friction (RSF) formalism developed for “regular” earthquakes (Dieterich, 1992). This efficient modeling approach has two significant shortcomings, however (van den Ende et al., 2018). First, all the dynamic variability of the interface and its surroundings is collapsed onto a planar boundary. Second, it depends on empirical RSF parameters that are not inferred directly from processes occurring within the fault zone. Thus, it is difficult to specify RSF parameters and account for internal fault zone dynamics in a self-consistent manner. These shortcomings may be limiting in the brittle-to-ductile transition zone of subduction faults because the deformation proceeds over an increasingly wide off-fault region as depth increases (Platt et al., 2018). Source dimensions for tremor and LFEs are typically inferred to be a few hundred meters (e.g., Chestler & Creager, 2017a; Sammis & Bostock, 2021; Supino et al., 2020), close to the thickness of the highly deformed permeable fault zone (e.g., Angiboust et al., 2015). This calls for more complete models dealing with processes in the interior of these zones.

There can be little doubt that pore fluid pressures and their variations play a key role in fault zone dynamics (Saffer & Tobin, 2011). From a frictional perspective, high pore pressures reduce the fault's strength, bringing it closer to unstable sliding conditions. Numerous geophysical and geological observations provide evidence for fluid saturation in the subduction interface (Saffer & Tobin, 2011). It is widely believed that fluids are supplied by metamorphic dehydration reactions in the downgoing plate (Frezzotti & Ferrando, 2015; Hyndman et al., 2015; Manning, 1997; van Keken et al., 2011). High values of the $V_p/V_s$ ratio and electrical conductivity indicate that pore fluids are at near-lithostatic pressures (Audet & Kim, 2016; Peacock et al., 2011; Wannamaker et al., 2014, among others), implying that the fault should be very weak. This explains why seismic activity is highly sensitive to very small dynamic stress changes generated by passing teleseismic waves (Rubinstein et al., 2009) or tidal deformation (Beeler et al., 2013; Royer et al., 2015; Rubinstein et al., 2008; Thomas et al., 2009).

There is mounting evidence for large variations of pore pressure on the time scale of the slow-earthquake cycle. During slow slip events, seismic images, gravimetric studies as well as earthquake source mechanisms indicate that pore fluid pressure around the subduction interface drops markedly
(Gosselin et al., 2020; Tanaka et al., 2018; Warren-Smith et al., 2019). The migration of LFE activity is highly suggestive of pore-pressure diffusion (Cruz-Atienza et al., 2018; Frank, Shapiro, et al., 2015). In exhumed subduction interfaces, the brittle-to-ductile transition seems to be intrinsically coupled to transient pore-pressure cycling and fluid transport (Angiboust et al., 2015; Behr & Bürgmann, 2021; Platt et al., 2018). Episodes of fluid accumulation and drainage have been documented in shallow thrust faults and are recorded by metasomatic reactions in host rocks (Angiboust et al., 2014; Saffer, 2015; Taetz et al., 2018). They have been attributed to changes of permeable pathways due to tectonic deformation, precipitation/dissolution or deposition/erosion processes (Giger et al., 2007; Jäger et al., 2017a; Williams et al., 2019). In recent efforts to model ETS phenomena, the emphasis has shifted from a static view of the impact of fluid pressure on RSF parameters to dynamic coupling between time-dependent fluid pressure and slip (e.g., Bernaudin & Gueydan, 2018; Cruz-Atienza et al., 2018; Zhu et al., 2020).

With the recognition that pore fluid pressures are high and variable in the subduction interface, it may not be appropriate to attribute seismic events to seismic slip only, without considering other processes that are active in the fault zone. Indeed, there are indications that the high fluid pressure transients that seem fundamental in shaping tremor patterns may well generate seismic events on their own. Observations do not demonstrate unambiguously that slow slip events and low-frequency seismic radiation are generated by the same process. It is true that these two types of activity are frequently correlated in time, but the sources of seismic and geodetic events are not systematically co-located (e.g., Hirose et al., 2010; Kostoglodov et al., 2010). In addition, in some subduction zones, LFEs do not conform to the moment-duration scaling of classical shear rupture (Bostock et al., 2015, and Chapter 2) and do not follow the Gutenberg-Richter size distribution characteristic of regular earthquake activity (Sammis & Bostock, 2021). Finally, it should be noted that the predominance of S-waves in the seismic signals does not necessarily establish a shear rupture mechanism. S-wave dominated signals can arise from volumetric sources with non-spherical geometries (e.g., Melnik et al., 2020a) or from single-force sources associated with fast local pressure variations (N. M. Shapiro et al., 2018).

Strong arguments in favor of a fluid-related origin for LFEs comes from laboratory studies and volcano seismology. In volcanic systems, shear deformation and slip do not play a predominant role, and yet tremors and LFEs (labeled as long-period earthquakes) are generated in abundance. These events have been attributed to fluid related processes, such as resonance in fluid-filled cracks (B. A. Chouet, 1996), strong flow in tortuous channels (Benson et al., 2008; Julian, 1994) and rapid pressurization of volatiles (Melnik et al.,
2020a; Wech et al., 2020), among others. In some cases, their sources are dis-
tributed over large depth intervals of a few tens of kilometers and appear
to migrate on time scales of a few months or days (N. M. Shapiro, Droznin,
Droznina, Senuykov, et al., 2017). These migrations have been attributed to
pore-pressure transients associated with hydraulic connections between dif-
ferent parts of a volcanic plumbing system.

A key piece of information on fault zone dynamics is the distribution of LFEs in space and time. Here, we explore the potential of pressure transients in a permeable subduction interface to generate this type of activity. The structure and hydraulics of such a channel have rarely been studied in relation to seismicity (Aochi et al., 2013; Cruz-Atienza et al., 2018; Segall & Rice, 1995; Zhu et al., 2020). We investigate the behavior of plugs of low permeability in a heterogeneous channel. Pore pressure variations can arise in two different ways. They may be generated by a time-dependent fluid supply, but there are no known mechanisms for strongly fluctuating rates of metamorphic dehydration reactions. Alternatively, one can invoke changes of permeability, as in the classical fault-valve model (R. Sibson, 1992). In this model, an increase of pore fluid pressure can generate an earthquake rupture, which opens up or reactivates fluid pathways leading to an increase of permeability and a decrease of pore pressure. Clogging of pores due to precipitation from circulating fluids or to the trapping of particles then acts to decrease permeability, implying an increase of pore pressure and a new cycle of pressure change. The permeability increase of the first phase of the cycle can be achieved in two different ways. One is shear displacement in association with fault dilation (Im et al., 2019a) — although it appears that, in some cases, this may in fact act to decrease permeability in the slip direction (Auradou et al., 2005). In the other mechanism, high pore pressures can induce fracturing and/or unclogging of small pores and pore throats. In the laboratory, impulsive flushing events similar to hydrofracturation have been observed in densely packed permeable beds due to sudden particle detachment and called erosive bursts (Bianchi, Thielmann, et al., 2018; Jäger et al., 2017b). These bursts are associated with sharp pore pressure drops and recur as particles are redeposited. Coupled variations of permeability and pore fluid pressure due to particle clogging and un-clogging have been observed in many laboratory experiments (Bianchi, Thielmann, et al., 2018; Candela et al., 2014) and have been surmised for fault zones (Manga et al., 2012). These phenomena have also been documented in large-scale systems. In some oil fields, for example, long-term records show that the bulk permeability of the reservoir oscillates in relation to cycles of deposition and removal of asphalt aggregates (Sahimi et al., 1999). We therefore argue that these mechanisms are bound to be active in faults, where fluid circulates under high pressure gradi-
ents, carrying precipitates and cataclastic fines that can be deposited and dis-
lodged. In this framework, fast local pore pressure variations (N. M. Shapiro
et al., 2018) and/or rupture of brittle heterogeneities due to increased fluid
pressure (Fagereng & Sibson, 2010; Kotowski & Behr, 2019) may be respon-
sible for LFEs. We expect that the unclogging or rupture of an individual
plug induces pore pressure variations that propagate to neighboring plugs
and generates swarms of similar events.

Through the rapid pressure drop that occurs, the unclogging of an in-
dividual plug is responsible for a local impulsive force on the solid matrix
which generates elastic waves with both compressional and shear polarisa-
tion components. N. M. Shapiro et al. (2018) have shown that the radiation
pattern of such a single-force source is consistent with seismic observations
in subduction zones due to the constrained source-receiver geometry, and
that it can account for the weak dependency of LFEs source duration on their
magnitude, as observed in Cascadia and Guerrero (Bostock et al., 2015, and
Chapter 2). Here, we build on this study to evaluate the seismic patterns that
can be generated by a series of interacting plugs in a permeable channel. The
clogging/unclogging mechanism does not involve shear rupture explicitly
and we will discuss a *posteriori* how it should reflect on slip and its sensitivity
at high fluid pressures.

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

## 3.2 Tremor patterns observed in subduction zones

We summarize the most salient features of tremor activity in subduction zones.
The Guerrero, Mexico, area is particularly appropriate to study LFE activity
because it stretches over a distance of more than 100 km along-dip, owing to
a low dipping angle. A very large catalog of events is available there (Frank
et al., 2014).
<figure><img src="image_39.png" /></figure>

**Figure 3.1** — The Guerrero (Mexico) tremor zone in (a) map view and (b) vertical cross-section. Each dot represents a low-frequency earthquake family from Frank et al. (2014). Thick black lines outline the trench and the subduction interface (data from Hayes (2018)). (c) and (d), LFE activity in a time-location chart from the Frank et al. (2014) catalog. Activity occurs in two patches centered on the subduction interface, a dense downdip one called the *sweet spot* and a sparser updip one called the *transient zone*. Tremor activity migrates along the subduction interface in both the updip and downdip directions, as shown by the orange arrows.

### 3.2.1 Spatial Segmentation

At a local scale, the fact that LFEs occur in families of recurring events with similar waveforms indicates that they are due to the repeated activation of individual sub-kilometric scale sources (Chestler & Creager, 2017a). In subduction zones, activity is often segmented in the along-strike direction and concentrated in patches (Brudzinski & Allen, 2007; Poiata et al., 2021). In Shikoku, Japan, and elsewhere, events delineate pluri-kilometric striations that are aligned with the subduction direction and that have been linked to the subduction of seamount chains (Ide, 2010). Shorter striations have also been observed in the Olympic Peninsula, Cascadia (Ghosh et al., 2010). In Guerrero, Mexico, activity is split into two well-separated patches along the dip of the interface (Frank et al., 2014; A. L. Husker et al., 2012; Kostoglakov et al., 2010; Payero et al., 2008) (Figure 3.1).
### 3.2.2 Tremor Migration

Bursts of tremor activity migrate along the subduction interface on time scales of seconds to weeks, and spatial scales from a few to a few hundreds of kilometers. Short-lived, short-distance migrations are often restricted to an individual tremor patch and proceed at rates in the 10–100 km/hr range (on time scales *T* of minutes to hours over length scales *L* of tens of kilometers) (Cruz-Atienza et al., 2018; Ghosh et al., 2010). In marked contrast, long-lived, long-distance migrations that connect spatially separated patches are slower, with velocities in the ~ 1–10 km/day range (*T* ~ weeks to months, *L* ~ 100 km) (Frank, Shapiro, et al., 2015; Kao et al., 2009). Such systematics between time and space scales are reminiscent of a diffusion process, such that *L*² ∝ *T* and therefore *V* ~ *L*/*T* ∝ 1/*L*.

In Japan and Cascadia, the largest transients of activity (*T* ~ weeks, *L* ~ 100 km) propagate along-strike in a belt lying in a 30 to 60 km depth range. It has been suggested that this depth interval delineates a zone rheologically conducive to slow-slip (Gomberg et al., 2016; Ide, 2012). In Guerrero, Mexico, the active segment of the subducted slab is dipping at a low angle and stretches over as much as 200 km in a 40-50 km depth range (Kim et al., 2010; Pardo & Suárez, 1995). Due to this particularly favorable configuration, it is clear that the largest migrations systematically originate from the deep end of the fault and slowly propagate updip (Frank, Shapiro, et al., 2015).

Short-lived migrations (*T* ~ minutes to hours) often proceed along striations (Ide, 2012), sometimes in rapid sweeps (Ghosh et al., 2010). Remarkably, these rapid migrations sometimes occur within larger scale migrations, propagating in the reverse direction with respect to the large scale propagation. This has been observed in the along-dip direction in Guerrero, Mexico (Cruz-Atienza et al., 2018; Frank et al., 2014) and in the along-strike direction in Cascadia (Houston et al., 2011).

### 3.2.3 Temporal Patterns

Tremor activity occurs during episodic bursts of events, which are composed of shorter bursts. Temporal clustering of activity covers time scales of years for the largest bursts — often correlated with slow-slip events and long-distance migrations — to days or even minutes for smaller bursts (Idehara et al., 2014; Wech & Creager, 2011) — called cascades. The imbrication of small bursts within larger ones makes it difficult to determine a characteristic duration. The shortest bursts proceed over a few minutes, whereas the largest ones can persist for days or even months (e.g., Frank et al., 2014; Ghosh et al., 2010; Shelly, 2015) As a temporal point process, this activity is clustered and does not conform to simple memoryless realizations. It has been analyzed using
non-Poissonian statistics (Beaucé et al., 2019; Frank et al., 2016).

The recurrence interval between tremor bursts decreases with increasing
depth along the subduction interface (Wech & Creager, 2011). In addition, the
time-clustering of activity can change from one tremor patch to its neighbors
(Frank et al., 2016; Idehara et al., 2014), and can evolve as a function of time
(Frank et al., 2016). This has been interpreted as evidence for variations of
fault zone conditions in both space and time (Frank et al., 2016).

**3.2.4 Tremor activity Patterns in Guerrero, Mexico**

Figure 3.1 shows the spatio-temporal patterns of LFE activity in the Guerrero tremor subduction segment, which is concentrated in two separate patches. One, called the *sweet spot*, is located at the downdip end of the segment and has a high rate of activity with a high level of time clustering (Frank et al., 2016). The other patch, called *transient zone* (Frank et al., 2014), is further updip and on the whole more quiescent. It hosts episodic bursts of activity with a recurrence interval of about 3 months (Frank et al., 2014), which coincide with short-term SSEs recurrence (Frank, Radiguet, et al., 2015). One of these episodes is displayed in Figure 3.1c. Activity migrates away from the sweet spot over a hundred kilometers updip. These episodes are built from many shorter bursts that may drift in both updip and downdip directions (Figure 3.1d).

The Guerrero subduction zone is also known for its very large SSEs that recur at intervals of about 4 years (e.g., Radiguet et al., 2012b) and are accompanied by intense tremor activity in both the transient zone and sweet spot. These events persist for several months with slow slip occurring mostly updip of the tremor-generating segment. The unstable slipping segment of the fault encroaches upon the transient zone over a small distance, where short-term SSEs are concentrated. During both short and long-term SSEs, tremor activity in the transient zone switches from episodic to highly clustered (Frank et al., 2016).

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

**3.3 Model design**

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

### 3.3.1 Diffusion equation for pore pressure variations

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

$$q(x, t) = -\frac{k(x)\rho}{\eta} \left( \frac{\partial P}{\partial x}(x, t) + \rho g \sin \alpha \right) \quad (3.1)$$

where $\rho$ is the fluid's density, $\eta$ is the fluid viscosity, $\alpha$ is the dip angle of the subduction interface and $g$ the acceleration of gravity. Mass balance dictates that:

$$-\frac{\partial q}{\partial x}(x, t) = \frac{\partial(\rho\phi)}{\partial t}. \quad (3.2)$$

where $\phi$ is porosity. We consider that pores are saturated with fluid and neglect temperature changes, which is justified by the small thickness of the channel and the small flow rates that are expected. Changes of the local mass of fluid can be due to the fluid compressibility as well as expanding/contracting pores. This is accounted for by an effective rock-fluid compressibility $\beta$, such that:

$$\begin{align}
\beta &= \frac{1}{\rho\phi}\frac{\partial(\rho\phi)}{\partial P} \nonumber \\
\implies -\frac{\partial q}{\partial x}(x, t) &= \rho\phi\beta\frac{\partial P}{\partial t}(x, t). \tag{3.3}
\end{align}$$

For simplicity, we neglect density variations in the body force term in equation (3.1), because they are due to changes of flow-induced dynamic pressure, which account for a fraction of the total pressure, and because they are not expected to play a significant role. We split pressure into a hydrostatic component and a time-dependent dynamic component due to flow noted $p$. Neglecting changes of fluid viscosity, we obtain the following equation for the dynamic pore pressure component:

$$\frac{\partial p}{\partial t} = \frac{1}{\phi\beta\eta}\frac{\partial}{\partial x}\left(k(x)\frac{\partial p}{\partial x}\right). \quad (3.4)$$

We assume that permeability is piecewise constant. Thus, in each segment, the dynamic pressure obeys a diffusion equation:

$$\frac{\partial p}{\partial t} = \frac{k}{\phi\beta\eta}\frac{\partial^2 p}{\partial x^2}, \quad (3.5)$$
<figure><img src="image_40.png" /></figure>

**Figure 3.2** — (a) Schematic representation of a permeable fault zone. Below the tremor source region, metamorphic dehydration reactions release fluid, which is channeled in the permeable subduction interface as it circulates under high fluid pressures. (b) Permeability structure for a 1-D hydraulic model of this system, involving a number of randomly distributed plugs that can open and close under certain conditions. The channel is fed at its base by input fluid flux $q_{in}$. At the top, the channel is connected to a high permeability region such that the fluid pressure is close to the hydrostatic pressure $p_{out}$. (c) Initial conditions correspond to steady-state flow regime with a small input flux and all plugs closed.

with diffusivity $D = k/(\phi\beta\eta)$. Along the modeled channel, $D$ is therefore piecewise constant.

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

### 3.3.2 A mechanism for large changes of permeability

Our model has two main ingredients, which deal with the nature of the fluid phase and its interaction with the network of pores in the permeable medium. Subduction zone metamorphic fluids are complex silica-saturated mixtures (Frezzotti & Ferrando, 2015; Manning, 1997). As they rise, they undergo decompression and, due to the pressure dependence of solubility, are expected to precipitate either amorphous silica (silica gel) or quartz nanoparticles as well as other mineral phases. Depending on the availability of nucleation sites, precipitation may proceed along the pore walls or in the bulk fluid. High temperature decompression experiments on silica-rich fluids demonstrate that the silica nanoparticles that are generated are readily mobilized by
fluid flow (Amagai et al., 2019; Okamoto, 2019). Small clasts generated by grinding motion along the fault are also likely involved (e.g., Han & Hirose, 2012). As in other deep permeable systems, the interstitial fluid phase is thus expected to carry fines and colloids (Manga et al., 2012).

Natural permeable media are made of pores and fractures connected by thinner throats which have a strong effect on permeability because they determine the effectiveness of fluid pathways through a volume. Capture and/or precipitation of even small amounts of solids within throats can decrease permeability by orders of magnitude (McDowell-Boyer et al., 1986). Calculations by Beckingham (2017) show that, depending on whether pore throats are open or clogged, permeability can change by more than a factor of ten. These changes are associated with negligible porosity variations because pore throats occupy a very small fraction of the total interstitial volume. Thus, the transition between high and low permeability states is very abrupt, as plugging thin throats between larger pores involves very small amounts of sediment. Even if they are much smaller than a pore throat, fines and colloids can get stuck through a jamming process at rates that depend on their concentration (Agbangla et al., 2014; Civan, 2016; Delouche et al., 2020). Clogging can be achieved in a matter of minutes (Candela et al., 2014; Delouche et al., 2020; Manga et al., 2012).

Material that gets stuck in pore throats may be remobilized by the flow (Amagai et al., 2019; Bianchi, Wittel, et al., 2018; Civan, 2016; Jäger et al., 2017b; Kudrolli & Clotet, 2016). Flow exerts a shear stress on the encasing solid, noted $\sigma$, that scales with the average pathway diameter $d$ and the pressure gradient $\partial p/\partial x$. If this stress exceeds a certain threshold value noted $\sigma_c$, infillings are eroded and transported by the fluid (Civan, 2016). This may occur in catastrophic fashion in events called *erosive bursts* (Bianchi, Thielmann, et al., 2018). Such bursts may be composites of several smaller unclogging events (Bianchi, Thielmann, et al., 2018). The process has been studied in detail by Kudrolli and Clotet (2016). In these experiments, a pure fluid is fed into a polydisperse particle bed. When the flow rate exceeds a threshold value, such that the flow-induced stress exceeds the threshold value $\sigma_c$, the flow dislodges and mobilizes the smaller particles. Interestingly, the authors observe small precursor events when the flow rate is just below critical and a large avalanche when the flow rate reaches the critical value. This is due to the collective behavior of particles, as some of them are not able to move if they are in contact with stronger ones that resist detachment. Thus, these particles are only dislodged when their neighbors are, which triggers an erosion avalanche. As long as $\sigma > \sigma_c$, particles cannot remain stuck and incipient coatings of pore walls are eroded, ensuring that fluid pathways stay fully open. Clogging thus requires that $\sigma \le \sigma_c$. The alternation between ero-
sion and deposition phases leaves little opportunity for particles and colloids to become permanently cemented to one another and to pore walls. Clogging and unclogging events therefore proceed rapidly, and can be considered as instantaneous at the time-scales of LFE activity.

As shown above, these processes have been observed in laboratory experiments on different types of porous media. Owing to their local character, they have not been documented in the field but, on a large scale, simultaneous variations of permeability and fluid pressure in faults (Manga et al., 2012) and oil reservoirs (Sahimi et al., 1999) have been attributed to erosion/deposition processes.

### 3.3.3 A valve mechanism

We model clogging and unclogging events as follows. Within a long permeable channel, let us consider a clogged plug where fluid pathways are narrow and contorted because they wrap around captured particles (closed state). In this case, the average pathway aperture $d$ and permeability $k$ are both small, with values $d_{lo}$ and $k_{lo}$, respectively. If and when $\delta p$, the pressure difference across the plug, reaches a critical value noted $\delta p_c^{break}$ such that the shear stress on the plug particles $\sigma$ is equal to critical value $\sigma_c$, unclogging occurs. As a consequence, pore throats are no longer obstructed by particles and permeability increases to $k_{hi} \gg k_{lo}$, in association with an increase of $d$ from $d_{lo}$ to $d_{hi}$ (open state). This drives a lowering of pore pressures upstream of the plug, and an increase downstream. This pressure transient is very rapid initially, and then decelerates, in a diffusion-induced flow transient. A key point is that $\sigma$ initially increases. Just before unclogging, it is equal to:

$$ \sigma_{-} = \lambda_{1}d_{lo}\frac{\delta p_{c}^{break}}{w} = \sigma_{c} \qquad (3.6) $$

where $\lambda_1$ is a coefficient and $w$ is the plug length along the channel. As plug particles are set in motion by the hydraulic stress and can no longer adhere to throat walls, the effective pore aperture $d$ increases from $d_{lo}$ to $d_{hi}$. In the permeable medium, variables are defined as spatial averages over a representative volume, and this is true for porosity and pressure in particular. As the solid particles get flushed out of the thin throats that connect larger pores, porosity is barely affected and, by the same token, so is the fluid pressure, because it is the average pressure for a group of pores and pore throats, with the larger pores acting to buffer local pressure variations at the pore throat scale. Thus, the pressure differential across the pore is maintained at $\delta p_c^{break}$
as the plug opens up (Jäger et al., 2018). Therefore, $\sigma$ increases to $\sigma_+$:

$$ \sigma_+ = \lambda_1 d_{hi} \frac{\delta p_c^{break}}{w} \qquad (3.7) $$

which is larger than $\sigma_c$. $\sigma$ then decreases gradually as the pressure distribution adapts to the new permeability structure. The plug remains unobstructed as long as $\sigma$ is larger than $\sigma_c$. Clogging occurs when $\sigma = \sigma_c$ such that $\delta p$ reaches a second threshold value of pressure differential noted $\delta p_c^{clog}$:

$$ \lambda_1 d_{hi} \frac{\delta p_c^{clog}}{w} = \sigma_c. \qquad (3.8) $$

Clogging causes the permeability to drop back to $k_{lo}$, driving an increase in $\delta p$ until unclogging occurs again, thus starting a new cycle. This is a valve behavior akin to the toggle-switch of Miller and Nur (2000).

This framework allows an internally consistent mechanism for a succession of unclogging and clogging events. One requirement is that $\delta p_c^{clog} < \delta p_c^{break}$. From equations (3.6) and (3.8), we deduce that:

$$ \frac{\delta p_c^{break}}{\delta p_c^{clog}} = \frac{d_{hi}}{d_{lo}}. \qquad (3.9) $$

Thus, the fact that $d_{hi} > d_{lo}$ automatically ensures that $\delta p_c^{clog} < \delta p_c^{break}$. One can go one step further using relationships between the average width of fluid pathways $d$ and permeability $k$. A common parameterization is the Kozeny-Carman equation (Civan, 2016):

$$ k = \lambda_2 \frac{\phi d^2}{\tau} \qquad (3.10) $$

where $\lambda_2$ is a coefficient that depends on the nature of the fluid pathways, e.g. tubes or slits, $\phi$ is porosity and $\tau$ is tortuosity. We may neglect variations of porosity, as discussed above, and obtain:

$$ \frac{k_{hi}}{k_{lo}} = \frac{\tau_{closed}}{\tau_{open}} \left( \frac{d_{hi}}{d_{lo}} \right)^2 \qquad (3.11) $$

where $\tau_{open}$ and $\tau_{closed}$ are the tortuosities of the open and closed plugs. Fluid pathways are more irregular and contorted in a closed plug than in an open one, implying that $\tau_{closed} > \tau_{open}$. Using equation (3.9), one deduces that:

$$ \frac{\delta p_c^{break}}{\delta p_c^{clog}} = \frac{d_{hi}}{d_{lo}} = \left( \frac{\tau_{open} k_{hi}}{\tau_{closed} k_{lo}} \right)^{1/2}. \qquad (3.12) $$

For the small porosity media that are involved here, tortuosity ratio $\tau_{open}/\tau_{closed}$
is likely to be in a 1/2 – 1/4 range (Ghanbarian et al., 2013). As stated above, permeability changes induced by clogging are typically larger than one order of magnitude (Beckingham, 2017; McDowell-Boyer et al., 1986), which guarantees that $\delta p_c^{break} > \delta p_c^{clog}$. For the sake of further discussion in section 3.4.1, it should be noted that such values of the tortuosity ratio imply that $k_{hi}/k_{lo} > \delta p_c^{break}/\delta p_c^{clog}$.

Figure 3.3 summarizes how these processes are implemented in our model. In valve segments of the channel, permeability responds to the local fluid pressure gradient according to the hysteretical cycle shown in Figure 3.3b.

### 3.3.4 Valves as elementary seismic sources

At each valve-opening, the fluid pressure gradient that has built up is suddenly dissipated in a rapid, localized diffusion transient. Variations of fluid pressure on each sides of the valve are responsible for temporal changes of the force applied to the solid walls of the permeable medium, which are able to generate seismic waves, as shown by N. M. Shapiro et al. (2018). Other seismic processes may also be triggered. For example, high fluid pressures may induce the seismic rupture of critically-stressed brittle heterogeneities (e.g., Taetz et al., 2018), and an impulsive turbulent flow transient may generate resonance in tortuous pathways (Benson et al., 2008). We thus consider that as they open, valves act as elementary seismic sources and generate LFE-like events.

### 3.3.5 Permeability values, boundary and initial conditions

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

Metamorphic reactions occur at different depths depending on the dehydrating material, which includes sediments, serpentinized mantle and hydrated oceanic crust. Focussing on the oceanic crust, dehydration proceeds at depths exceeding the tremor source region (van Keken et al., 2011) and most of this fluid gets trapped in the high permeability subduction interface, sealed from above by the low permeability mantle wedge (Hyndman et al., 2015). Thus, the channel is fed at its base by an input fluid flux $q_{in}$. For a first instalment of the model, it is instructive to maintain this flux at a constant value. At its updip end, the channel is connected to a high permeability
<figure><img src="image_41.png" /><figcaption>**Figure 3.3** — (a) The 4 stages of the valve mechanism. 1, the valve has been closed and permeability is $k_{lo}$. $\delta p$, the pressure difference across the valve, increases. 2, when $\delta p$ reaches the threshold value $\delta p_c^{break}$, an impulsive unclogging event (an erosive burst) occurs. Particles get flushed out, the valve opens and permeability rises to $k_{hi}$. 3, valve opening generates a rapid pressure transient that is associated with a low frequency earthquake. 4, when $\delta p$ decreases to $\delta p_c^{clog}$, the deposition of fines and colloids resumes and permeability drops to $k_{lo}$. The valve closes and a new cycle begins. (b) Hysteretical cycle of pressure difference $\delta p$ and valve permeability $k$. (c). Evolution of $\delta p$ through 3 successive cycles of valve opening and closure for an isolated valve.</figcaption></figure>
region which may be either continental crust (Hyndman et al., 2015) or an accretionary prism. In this case, pore pressures are close to hydrostatic values, such that the dynamic pressure is zero ($p_{out} = 0$).

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

<table><caption>**Table 3.1** — Characteristic dimensions used to scale physical variables</caption><thead><tr><th>Variable</th><th>Characteristic scale</th></tr></thead><tbody><tr><td><i>x</i> Distance along channel</td><td><i>X</i> = <i>L</i></td></tr><tr><td><i>z</i> Depth</td><td><i>Z</i> = <i>X</i> sin(<i>&alpha;</i>)</td></tr><tr><td><i>t</i> Time</td><td><i>T</i> = <i>X</i><sup>2</sup>/<i>D</i></td></tr><tr><td><i>p</i> Pore pressure</td><td><i>P</i> = (<i>&rho;</i><sub>r</sub> &minus; <i>&rho;</i>)g<i>Z</i></td></tr><tr><td><i>q</i> Massic flow rate<sup>b</sup></td><td><i>Q</i> = <i>k</i><sub>hi</sub><i>&rho;</i>/<i>&eta;</i> &times; <i>P</i>/<i>X</i></td></tr><tr><td><i>k</i> Permeability</td><td><i>K</i> = <i>k</i><sub>hi</sub></td></tr></tbody></table>

<table><caption>**Table 3.2** — Valve parameters for all simulations presented in this chapter</caption><thead><tr><th>Parameter</th><th>Value (scaled)</th></tr></thead><tbody><tr><td><i>w</i> Valve width</td><td>0.02</td></tr><tr><td><i>k</i><sub>hi</sub> Open valve permeability</td><td>1</td></tr><tr><td><i>k</i><sub>lo</sub> Closed valve permeability</td><td>0.05</td></tr><tr><td>δ<i>p</i><sub>c</sub><sup>break</sup> Threshold <i>p</i> difference for opening</td><td>0.04</td></tr><tr><td>δ<i>p</i><sub>c</sub><sup>clog</sup> Threshold <i>p</i> difference for closing</td><td>0.021</td></tr><tr><td><i>q</i><sub>c</sub><sup>break</sup> Flux above which a closed valve can open</td><td>0.1</td></tr><tr><td><i>q</i><sub>c</sub><sup>clog</sup> Flux below which an open valve can close</td><td>1.05</td></tr></tbody></table>

<footnote><sup>a</sup> along a fully open channel, $D = k_{hi}/(\rho\eta\beta)$</footnote>

<footnote><sup>b</sup> per unit channel cross section area</footnote>