# Où et quand se produisent les grands séismes ?
Comment les anticiper ?

Quentin Bletery

► To cite this version:

Quentin Bletery. Où et quand se produisent les grands séismes ? Comment les anticiper ?. Géophysique [physics.geo-ph]. Université Côte d'Azur; ED n°364 : Sciences fondamentales et appliquées, 2022. tel-03869413

HAL Id: tel-03869413

https://hal.science/tel-03869413

Submitted on 24 Nov 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
# HABILITATION À DIRIGER DES RECHERCHES

DE L'UNIVERSITÉ CÔTE D'AZUR

ÉCOLE DOCTORALE SCIENCES FONDAMENTALES ET APPLIQUÉES

Présentée par

Quentin Bletery

## Où et quand se produisent les grands séismes ?

## Comment les anticiper ?

Soutenue le 4 Novembre 2022
Devant la Commission d'Examen composée de :

## JURY

Isabelle Manighetti

Cécile Lasserre

Jean-Philippe Avouac

Jean-Paul Montagner

Cécile Cornou

Michel Campillo

Physicienne à l'Observatoire de la Côte d'Azur

Directrice de Recherche au CNRS

Professeur au California Institute of Technology

Professeur à l'Institut de Physique du Globe de Paris

Directrice de Recherche à l'IRD

Professeur à l'Université Grenoble Alpes

Présidente du Jury

Rapportrice

Rapporteur

Rapporteur

Examinatrice

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

Je remercie les membres du Jury, Isabelle Manighetti, Cécile Lasserre, Cécile Cornou, Jean-
Philippe Avouac, Jean-Paul Montagner et Michel Campillo d'avoir considéré mon travail avec
intérêt et pris le temps de le juger.

Je remercie les personnes qui ont participé à la rédaction des publications qui composent, en
partie ce manuscrit, Amanda M. Thomas, Alan W. Rempel, Leif Karlstrom, Anthony Sladen,
Louis de Barros, Jeanne L. Hardebeck, Olivier Cavalié, Jean-Mathieu Nocquet, Théa Ragon,
Jessica C. Hawthorne, Robert M. Skarbek, Randy D. Krogstad, Andrea Licciardi, Bertrand
Rouet-Leduc, Jean-Paul Ampuero et Kévin Juhel.

Je remercie plus personnellement les chercheurs qui m'ont formé au métier de la recherche,
Luis Rivera, Mark Simons, Anthony Sladen, Bertrand Delouis, Amanda M. Thomas et Alan W.
Rempel ainsi que les doctorants et postdoctorants avec qui j'ai la chance de travailler, Gabriela
Arias, Pablo Lara, Andrea Licciardi et Kévin Juhel.

Merci au Seismological Laboratory du California Institute of Technology, au département
de Sciences de la Terre de l'Université d'Oregon et à Géoazur, ainsi qu'à tous les chercheurs
qui y travaillent, de m'avoir accueilli aux différentes étapes de ma carrière.

Merci à ma famille, à mes amis, ma compagne, pour avoir traverser des moments difficiles
avec moi, Aline, Anne-Laure, Nicolas, Joris, Niels, François-Xavier, Raphaël, Agnès, Jacques
et bien sûr Laure. Merci infiniment.

« Le Vrai est le Tout »

Georg Wilhelm Friedrich Hegel
À mon père,
Table des matières

Introduction 1

I Anticiper les séismes, pourquoi est-ce si difficile? 3

1 Qu'est-ce qu'un séisme? 3

1.1 Observations en surface 3

1.2 Mécanisme. 5

1.3 Que se passe-t'il en profondeur? 7

2 Que font les failles entre deux séismes? 9

2.1 Glissement post-sismique. 10

2.2 Glissement inter-sismique 11

2.3 Glissement transitoire 12

3 Est-il possible de prédire les séismes? 14

3.1 Approche physique 14

3.2 Approche statistique 16

3.3 Approche observationnelle 18

4 Anticiper les séismes autrement : l'alerte précoce. 19

4.1 L'alerte sismique 19

4.2 L'alerte tsunami. 21

4.3 La découverte prometteuse des PEGS. 22

5 Possibilités offertes par l'émergence de l'IA 24

5.1 Principes simplifiés de l'IA à l'usage des non-spécialistes 24

5.2 Accomplissements récents en sismologie. 26
Introduction

Pourquoi la Terre tremble ? Cette question, que se posent les hommes – probablement – depuis toujours, a trouvé, à la fin des années 1960, une réponse partielle avec la théorie de la tectonique des plaques : des plaques rigides, entraînées par un manteau asthénosphérique en mouvement, glissent les unes par rapport aux autres. Leur mouvement relatif est accommodé lors d'évènements rares mais violents : les séismes. Si ce cadre conceptuel simple a permis de comprendre l'origine des tremblements de terre, il ne permet pas de répondre aux deux questions fondamentales qui permettraient de limiter le coût humain de ces catastrophes : Où et Quand se produisent les grands séismes ? Ce mémoire synthétise mes travaux de recherche réalisés depuis la soutenance de mon doctorat qui, par divers angles, ont tenté de répondre à ces questions.

Le chapitre 1 présente un bref état de l'art autour des concepts de rupture sismique, cycle sismique, séismes lents, des observations qui les contraignent ainsi que des outils et méthodes utilisés pour étudier le comportement des failles. Je mettrai en relief les difficultés à anticiper les ruptures sismiques, un phénomène impossible à prédire à partir des observations et outils que l'on utilise à l'heure actuelle. Je présenterai brièvement les méthodes d'alerte sismique et tsunami, une alternative à la prédiction pour réduire le coût sociétal des tremblements de terre. Enfin, je présenterai un outil émergent ayant le potentiel pour révolutionner la recherche autour de l'alerte : l'Intelligence Artificielle (IA).

Le chapitre 2 s'intéresse à la question du « Où ? ». Où les grands séismes sont-ils le plus susceptible de se produire ? Je présenterai dans ce chapitre trois études tentant de répondre à cette question par trois angles différents. La première montre que les grandes ruptures sismiques se produisent sur des portions de failles particulièrement planes. La seconde décrit un cadre théorique permettant d'évaluer la distribution spatiale du seuil de rupture en contraintes le long
des failles. Appliquée au Japon cette approche révèle que les grands séismes se produisent sur
des portions de faille présentant un seuil de rupture particulièrement homogène. La troisième
approche est une inversion Bayésienne de la distribution spatiale du couplage inter-sismique le
long des failles Nord- et Est- Anatolienne, en Turquie. Cette approche permet d'identifier où
les contraintes s'accumulent le long de ces grandes failles, et donc où celles-ci sont susceptibles
d'être relâchées au cours d'un évènement sismique à venir.

Le chapitre 3 s'intéresse à la question du « Quand ? ». Répondre à cette question est un échec
notoire de la sismologie contemporaine. Mais est-ce réellement impossible ? Serait-il possible de
connaître la dynamique des failles au point de pouvoir « prédire » quand de grands évènements
pourraient se produire. Comment étudier cette dynamique ? L'identification d'évènements de
glissement transitoires – les séismes lents – ouvre une fenêtre pour monitorer les failles en
ce sens. Je présenterai dans ce chapitre une étude sur la dynamique interne des séismes lents,
laquelle a mis en évidence que ces derniers sont constitués de sous-épisodes de glissement. Puis je
présenterai une étude dans laquelle nous montrons que la coalescence de séismes lents distincts
génère un « emballement » du glissement, un mécanisme observable sur les failles naturelles qui
pourrait être à l'origine de la nucléation de grands séismes.

Si le chapitre 3 suggère la possibilité de la prédiction des séismes, celle-ci reste pour l'heure
très hypothétique. Le chapitre 4 vise à dépasser la question de la prédicibilité pratique ou
théorique des séismes à travers une approche concrète pour anticiper les tremblements de terre
et les tsunamis. Je présenterai PEGSnet, un algorithme d'IA que nous avons développé, qui uti-
lise un signal, identifié récemment, d'amplitude très faible mais précédant les premières ondes
sismiques. Ce signal d'origine gravitationnelle a la propriété de se propager à la vitesse de la lu-
mière. Grâce a cette propriété, PEGSnet est capable d'évaluer quasi-instantanément l'évolution
de la magnitude des grands séismes et ainsi apporter des informations crucales dans la décision
de lancer une alerte sismique/tsunami avant les premières secousses et les raz de marée. Le cha-
pitre se poursuivra par quelques perspectives sur lesquelles nous travaillons actuellement avec
les doctorants et postdoctorants que j'encadre. Je présenterai quelques uns des développements
en cours autour de PEGSnet, l'un de nos grands objectifs étant de développer un algorithme
d'alerte global capable de suivre l'évolution de la magnitude des grands séismes en tout point
de la Terre à partir de sismomètres installés partout sur Terre. Je présenterai ensuite un second
algorithme d'IA qui vise à combler le déficit de performance de PEGSnet pour les séismes de
magnitude inférieure à 8.3. Cet algorithme n'utilise que 3 secondes d'enregistrement sur une
unique station, procurant ainsi des informations pertinentes pour l'alerte en des temps records.

Enfin, ce mémoire s'achèvera par quelques réflexions sur les possibilités apportées par l'IA
en sismologie : la prédiction est-elle réellement impossible ?
# Chapitre I

## Anticiper les séismes, pourquoi est-ce si difficile ?

### 1 Qu'est-ce qu'un séisme ?

Un séisme est un glissement rapide le long d'une faille, interface entre deux « blocs de Terre » en mouvement relatif l'un par rapport à l'autre. Il existe des séismes de toutes tailles. Les plus petits se produisent de manière quasi-continue sans qu'on les ressente. Les plus grands génèrent des secousses capables de détruire des villes. Certains peuvent générer des tsunamis qui font encore plus de dégâts. Le bilan humain de ces catastrophes naturelles s'élève à près d'un million de morts sur les trois dernières décennies. Ces phénomènes apparaissent comme des évènements chaotiques que la Science moderne n'arrive pas à prédire. Nous disposons pourtant de nombreuses observations à partir desquelles nous essayons de comprendre ce qu'est un séisme, dans l'espoir, un jour, de savoir où et quand la Terre tremblera.

#### 1.1 Observations en surface

Les séismes se produisent sur des failles à plusieurs kilomètres de profondeur. Les observations de ces évènements dont l'on dispose sont malheureusement toutes obtenues à partir de capteurs à la surface de la Terre. Nous ne disposons donc que d'observations indirectes pour étudier ces phénomènes. Les deux grands types d'observations que l'on utilise à cette fin sont les observations sismologiques et les observations géodésiques.

#### Observations sismologiques

Les ruptures sismiques génèrent un champ d'onde responsable des « tremblements » de la Terre ressentis aux alentours de la rupture. Ce champ d'onde est enregistré localement, et globalement sur la planète, par des appareils de mesure appelés sismomètres. Il existe deux grands types de sismomètres : les vélocimètres qui enregistrent la vitesse des mouvements du
sol et les accéléromètres qui enregistrent leur accélération.

Passer d'un sismogramme en vitesse à un sismogramme en accélération est très simple et ne nécessite que de dériver le signal. La réelle différence entre les deux types d'instruments est leur sensibilité à différentes fréquences. Les vélocimètres sont également appelés sismomètres large-bande car ils sont capables d'enregistrer des signaux sur une gamme de fréquences très large (~ 0.01 - 20 Hz). Ce sont les instruments les plus employés pour étudier les séismes. Ils contraignent par exemple les signatures basse fréquence des séismes qui permettent de déterminer leur magnitude [Hanks and Kanamori, 1979].

Les accéléromètres sont beaucoup plus bruités à basse fréquence (< 1 Hz) mais sont utiles pour enregistrer les mouvements forts du sol. En effet, proches d'une forte source sismique ils saturent moins que les vélocimètres. Ils sont ainsi utiles pour caractériser rapidement la source sismique dans une optique d'alerte [Wu and Kanamori, 2008].

## Observations géodésiques

Les observations géodésiques sont des mesures de déplacement du sol au cours du temps. Elles sont très utiles pour contraindre la distribution spatiale des ruptures sismiques. Les deux techniques géodésiques les plus utilisées pour étudier la source sismique sont le *Global Navigation Satellite System* (GNSS) et l'interférométrie radar.

Le système GNSS est constitué d'une constellation de satellites permettant de localiser précisément des récepteurs à la surface de la Terre. Le déploiement d'un grand nombre de récepteurs – appelés stations GNSS (ou abusivement stations GPS) – permet de déterminer de manière précise l'orbite des satellites dans un référentiel terrestre régulièrement mis à jour, l'International Terrestrial Reference Frame (ITRF) [Altamimi et al., 2002, 2007, 2011, 2016]. Après un post-traitement élaboré [King, 1995], les mesures GNSS peuvent atteindre une précision millimétrique. Lors d'un séisme, le GNSS enregistre le déplacement du sol à l'emplacement des stations sur les trois composantes de l'espace. Ces déplacements peuvent être de plusieurs mètres pour de grands évènements [e.g., Ozawa et al., 2011]. Historiquement utilisé pour des mesures journalières – servant en particulier à la caractérisation du champ de déplacement long-terme [e.g., DeMets et al., 2000; Nocquet and Calais, 2004; Calais et al., 2006] – le GNSS permet aujourd'hui des mesures haute fréquence (jusqu'à 10 Hz) [Melgar et al., 2020]. Toutefois, la précision dépend d'un « moyennage » des erreurs sur un grand nombre de mesures. Elle a donc tendance à décroitre lorsque l'on échantillonne à plus haute fréquence.

L'autre grande technique de géodésie spatiale utilisée pour étudier les tremblements de Terre est l'interférométrie radar, ou InSAR pour Interferometric Synthetic Aperture Radar. Le
principe de cette méthode est d'utiliser les différences de phase entre deux images radar (SAR). Un interférogramme entre une image SAR acquise après un séisme et une image acquise avant donne le champ de déplacement généré par le séisme en surface dans la direction de visée du satellite. Là où le GNSS permet des mesures ponctuelles, l'InSAR permet des mesures sur de vastes étendues spatiales. Obtenir un interférogramme nécessite toutefois un certain nombre d'étapes de traitement, l'étape la plus cruciale étant la correction atmosphérique [Doin et al., 2009; Jolivet et al., 2011]. La résolution d'un interférogramme bien corrigé est de l'ordre du centimètre mais dépend de la fréquence du radar, qui varie d'un satellite à l'autre. La résolution temporelle est déterminée par la fréquence de passage des satellites. Historiquement de l'ordre du mois avec les satellites d'ancienne génération ERS-1, ERS-2, ALOS, ENVISAT, elle est depuis le lancement du programme SENTINEL de moins d'une semaine. Comme le GNSS, l'InSAR est utilisé pour mesurer la déformation lors d'un séisme [e.g., Massonnet et al., 1993; Lasserre et al., 2005], mais également la déformation long-terme [e.g., Cavalié et al., 2008].

Autres types d'observations

La manifestation en surface la plus impressionnante des séismes est certainement le tsunami qui suit les plus grands d'entre eux. Des bouées dites DART (pour Deep-ocean Assessment and Reporting of Tsunamis) ou encore des capteurs de pression sous-marin [Hino et al., 2001; Maeda et al., 2011] permettent de mesurer l'élévation du niveau de l'eau et donc le passage éventuel d'un tsunami avant qu'il ne frappe les côtes. Des mesures de la perturbation des signaux GNSS dans l'ionosphère permettent également une mesure indirecte du déplacement co-sismique [Calais and Minster, 1995] et possiblement du tsunami [Heki et al., 2006; Rolland et al., 2011; Bagiya et al., 2020; Sunil et al., 2021].

La technologie la plus prometteuse en développement est certainement celle du *Distributed Acoustic Sensing* (DAS) appliquée à des câbles de télécommunication en fibre optique. Cette technologie permet virtuellement d'avoir accès à des milliers de sismogrammes le long des câbles, une densité instrumentale sans précédent à faible coût, qui pourrait révolutionner l'observation sismique en particulier en mer [Sladen et al., 2019].

1.2 Mécanisme

Expliquer les diverses manifestations de la source sismique observées à la surface de la Terre nécessite de résoudre différents problèmes physiques, selon les type de données. Les équations décrivant la résolution de ces divers problèmes sont détaillées dans ma thèse [Bletery, 2015]. Je me contenterai ici d'exposer succinctement les moyens couramment employées pour modéliser
les ondes sismiques, le déplacement du sol et le tsunami.

## Les ondes sismiques

Les ondes sismiques enregistrées sur les sismomètres sont bien expliquées par le théorème de représentation, qui donne le champ d'onde généré par une dislocation cisaillante ponctuelle dans un milieu isotrope, homogène, infini [Aki and Richards, 2002]. Certains termes peuvent être négligés en champs lointain (si les stations sont loin de la source), comme dans le code développé par Nabelek [1984] dans une Terre stratifiée (1D). Dans le cas de stations proches, tous les termes de radiations doivent être pris en compte, ce que permet la méthode des nombres d'onde discrets [e.g., Bouchon, 1981], également dans une Terre stratifiée. Des méthodes de résolution numérique en 3D existent, comme SPECFEM ou sa version simplifiée (plus rapide) AXISEM [Nissen-Meyer et al., 2014], mais sont limitées par leur coût en terme de temps de calcul ainsi que par notre connaissance limitée de la complexité 3D de la Terre, limitant l'intérêt d'une complexification numérique top poussée pour l'étude de la source.

## Le déplacement du sol

Le déplacement du sol enregistré par les données géodésiques est donné par une solution analytique pour une dislocation « rectangulaire » dans un milieu élastique homogène semi-infini. Cette solution analytique connue sous le nom de solution d'Okada a en fait été résolue par Mansinha and Smylie [1971] à partir d'une intégrale proposée par Steketee [1958]. Okada [1985] a, des années plus tard, implémenté cette solution analytique dans un script utilisé par toute la communauté depuis. Des codes intégrants des degrés de complexité supplémentaires ont depuis été développés, comme celui de Zhu and Rivera [2002] qui considère une Terre stratifiée par la technique des matrices propagatives de Thompson-Haskell ou celui de Pollitz [1996] qui permet de considérer une Terre sphérique stratifiée à partir de calculs de modes propres. Aujourd'hui, des codes comme PyLith permettent de prendre en compte la complexité 3D de la Terre [e.g., Romano et al., 2014], mais comme pour le champ d'onde sismique, notre connaissance limitée de cette complexité rend sa prise en compte presque superflue pour l'étude de la source.

## Le tsunami

Le mécanisme à l'origine des tsunamis est assez simple : un séisme entraine une élévation du fond de l'océan qui peut, pour les plus grands évènements, atteindre plusieurs mètres et s'étendre sur des centaines de kilomètres. La colonne d'eau située au dessus de la partie du fond de l'océan sur-élevée par le séisme est automatiquement translatée verticalement. Un immense volume
d'eau se retrouve alors au dessus du niveau de l'équipotentiel. Le ré-équilibrage du niveau de
la mer se fait par propagation d'une onde de gravité : un tsunami. Des codes de simulation
numérique comme COMCOT ou NEOWAVE [Yamazaki et al., 2009, 2011] permettent calculer
le tsunami partout sur Terre à partir des paramètres du séisme. Ces codes négligent toutefois la
contribution du déplacement horizontal du fond de l'océan [Tanioka and Satake, 1996; Bletery
et al., 2015] ainsi que des phénomènes complexes de dispersion induits par la réponse élastique
de la Terre et l'hétérogénéité de la densité des couches d'eau [Tsai et al., 2013; Watada, 2013;
Allgeyer and Cummins, 2014].

1.3 Que se passe-t'il en profondeur ?

Les observations sismologiques permettent depuis de nombreuses années de localiser la
source des séismes en profondeur [e.g., Kennett and Engdahl, 1991]. On dispose ainsi de vo-
lumineux catalogues de sismicité nous renseignant sur où les séismes se produisent à travers
le monde. Ces catalogues ont permis de modéliser la géométrie des grandes failles (Figure I.1)
[Hayes et al., 2012, 2018].

<figure><img src="image_3.png" /><figcaption>**Figure I.1** – Modèle slab2.0 décrivant la géométrie des failles de subdcution [Hayes et al., 2018].</figcaption></figure>
<figure><img src="image_4.png" /></figure>

**Figure I.2** – Glissement sur la faille survenu lors du séisme de Tohoku-Oki (M<sub>w</sub> 9.0, 2011). Les formes d'ondes autour de la carte sont des exemples de données sismologiques, géodésiques et tsunami (couleur selon le type de données) expliquées par le modèle (courbes noires). La figure en insère donne le moment relâché en fonction du temps [Bletery et al., 2014].
Le long de ces grandes failles se produisent les plus grands séismes. Une approche populaire pour étudier la source sismique est d'imager le glissement le long des failles pendant les plus grands séismes. Cette approche consiste à chercher la distribution de glissement qui explique le mieux les différentes observations du séismes, i.e. les données sismologiques, géodésiques et éventuellement les mesures du tsunami. De nombreuses techniques, dites d'inversion, permettent ainsi d'obtenir un modèle de rupture en accord avec les données. Ces différentes techniques sont présentées en détails dans ma thèse [Bletery, 2015]. Leur principe de base est de discrétiser la faille rompue en sous-failles et chercher quelle amplitude de glissement sur chaque sous-faille à chaque instant explique le mieux les données. Cela permet d'obtenir une « image » de la rupture sur la faille (Figure I.2), et même un film [e.g., Bletery et al., 2014].

Les résultats donnés par les techniques d'inversion présentent toutefois une faiblesse fondamentale : la solution donnée par l'inversion n'est jamais unique. Cette faiblesse limite l'interprétation que l'on peut faire des modèles car différents groupes d'auteurs peuvent aboutir à des modèles très différents pour un même séisme [e.g., Lay, 2018]. Pour résoudre ce problème fondamental, des régularisations sont appliquées pour « guider » l'inversion vers la solution la plus lisse. Ce biais dans les modèles a été introduit car l'on pense qu'un modèle lisse (avec un glissement relativement homogène) est plus physique qu'un modèle rugueux (avec un glissement plus hétérogène), mais cet a priori n'est pas réellement fondé et les modèles les mieux contraints sont plutôt hétérogènes [e.g., Lay, 2018]. Une approche plus sophistiquée, permettant d'obtenir l'ensemble des solutions du problème et ainsi d'évaluer les incertitudes sur les modèles, a émergé récemment [Minson et al., 2013]. Le cadre théorique de cette approche permet même de prendre en compte les erreurs épistémologiques afin d'obtenir l'ensemble des solutions expliquant de manière satisfaisante les observations en fonction des incertitudes sur les données et sur les paramètres du modèle, comme les paramètres élastiques de la Terre ou la géométrie de la faille [Duputel et al., 2014; Bletery et al., 2016a; Ragon et al., 2018, 2019a, 2021].

## 2 Que font les failles entre deux séismes ?

Que se passe-t-il sur une faille lorsqu'il ne se passe rien ? Cette question est fondamentale pour comprendre ce que sont les séismes, comment, où et quand ils se produisent. La technique d'inversion décrite au paragraphe précédant peut être adaptée pour imager le glissement le long des failles entre 2 séismes. Celle-ci a permis de mettre en évidence différentes phases du cycle sismique : la phase post-sismique suivant les grands évènements, la phase inter-sismique et des épisodes de glissement transitoires.
## 2.1 Glissement post-sismique

Les inversions de données géodésiques enregistrées après les grands séismes ont révélé qu'après une grande rupture, le glissement sur la faille ne s'arrête pas soudainement, il décelère progressivement (suivant une décroissance logarithmique, Figure I.3). Les plus grandes amplitudes de glissement post-sismique semblent, par ailleurs, se concentrer en bordure des zones rompues co-sismiquement (Figure I.3).

<figure><img src="image_5.png" /><figcaption>**Figure I.3** – Distribution spatiale (a) et temporelle (b - e) du glissement post-sismique suivant les séismes de Nias (M<sub>w</sub> 8.6, 2005) et Mentawai (M<sub>w</sub> 7.9, 2007) [Avouac, 2015].</figcaption></figure>

La transition entre glissement co-sismique et post-sismique n'est pas évidente et une grande
partie du signal que l'on attribue généralement au glissement co-sismique pourrait bien être a-
sismique [Twardzik et al., 2019; Ragon et al., 2019b]. Pour les très grands séismes, le glissement
post-sismique peut durer des années. Au bout de quelques mois/années, le signal dans les
données géodésiques devient toutefois dominé par la relaxation visco-élastique du milieu et non
plus par le glissement sur la faille [e.g., Sun and Wang, 2015].

## 2.2 Glissement inter-sismique

Une fois la phase post-sismique terminée, le glissement sur les failles reprend progressivement
son comportement long terme. Ce glissement inter-sismique n’est généralement pas homogène,
certaines portions de failles « glissent » à la vitesse relative des plaques tectoniques, d’autres
portions restent « bloquées », ou couplées.

<figure><img src="image_6.png" /><figcaption>**Figure I.4** – Distribution du couplage intersismique le long de la faille de subduction du Mexique [Rousset *et al.*, 2016].</figcaption></figure>

Le degrés de couplage – variant entre 0, si la faille glisse à la vitesse des plaques, et 1 si
elle est complètement couplée – peut être imagé le long des grandes failles par inversion de données géodésiques long terme [e.g., Loveless and Meade, 2011; Stevens and Avouac, 2015; Jolivet et al., 2015; Metois et al., 2016; Rousset et al., 2016; Villegas-Lanza et al., 2016; Jolivet et al., 2020]. Imager la distribution du couplage le long des grandes failles permet d'identifier les zones où des déficits de glissement s'accumulent (Figure I.4). Les zones à fort couplage sont ainsi, dans le schéma conceptuel classique, considérées comme les zones les plus à risque de produire des grands séismes [e.g., Avouac, 2015].

## 2.3 Glissement transitoire

Le début des années 2000 a vu la découverte de signaux transitoires dans les données GNSS de la zone de subduction des Cascades [Dragert et al., 2001; Miller et al., 2002]. Parallèlement, des signaux micro-sismiques, du nom de trémors, ont été reliés à l'activité tectonique au Japon [Obara, 2002]. Rogers and Dragert [2003] ont alors observé une forte corrélation temporelle entre ces deux types de signaux dans la zone des Cascades (Figure I.5).

<figure><img src="image_7.png" /><figcaption>**Figure I.5** – Signal journalier enregistré sur une station GNSS aux Cascades (cercles bleus) et activité des trémors (courbe noire) [Rogers and Dragert, 2003].</figcaption></figure>

Au delà de la corrélation temporelle frappante (Figure I.5), la périodicité du signal a énor- mément intéressé la communauté [Schwartz and Rokosky, 2007]. L'inversion des signaux GNSS a mis en évidence que ces phénomènes sont des épisodes de glissement transitoire sur la faille [e.g., Bartlow et al., 2011], accompagnés de millions de micro-séismes dont la superposition gé- nère les trémors enregistrés par les sismomètres [Shelly et al., 2007]. Les épisodes de glissement
transitoire ont été nommés *Slow Slip Events* (SSEs), ou séismes lents, les micro-séismes qui les accompagnent *Low-Frequency Earthquakes* (LFEs) car ces évènements apparaissent comme dépourvus en hautes fréquences [Shelly et al., 2006]. Une propriété intéressante des LFEs est leur sensibilité à des perturbations de contraintes extrêmement faibles [Thomas et al., 2009]. Leur activité représente, par ailleurs, une fenêtre d'observation des conditions physiques en profondeur [Thomas et al., 2018]. D'autres évènements, encore plus basses fréquences, nommés *Very-Low-Frequency Earthquakes* (VLFEs) ont été identifiés [Ito et al., 2007]. Le cadre conceptuel classique est que les SSEs, des glissements lents a-sismiques, déstabilisent des micro-aspérités, générant ainsi des LFEs et des VLFEs que l'on observe dans les données sismologiques sous forme de trémors. Toutefois la relation entre SSEs et LFEs n'est pas totalement établie, les LFEs n'étant pas toujours observés au centre des zones de glissement lent, mais souvent en bordure. La relation entre tous ces phénomènes, la distribution du couplage [Jolivet and Frank, 2020] et les grands séismes (Figure I.4) comporte peut-être la clé de la compréhension et de l'anticipation des séismes [Obara and Kato, 2016].

<figure><img src="image_8.png" /><figcaption>**Figure I.6** – Loi d'échelle proposée par Ide et al. [2007] liant SSEs, LFEs et VLFEs, séparée de la loi d'échelle des séismes par un gap d'observations.</figcaption></figure>
Des techniques de plus en plus fines ont permis de détecter des épisodes de glissement transitoire de plus en plus petits [Rousset et al., 2017, 2019] et de partout sur la planète [Obara and Kato, 2016, et références citées]. S'est alors posée la question d'un continuum des phénomènes de glissement lent : existe-t-il des épisodes de glissement lent de toute magnitude, de toute durée [Frank et al., 2018; Jolivet and Frank, 2020] ? Ce continuum inclut-il les évènements sismiques ou les glissements lents sont-ils un mode de glissement distinct des séismes ? Sur la base d'une loi d'échelle entre moment et durée (Figure I.6), Ide et al. [2007] ont proposé que SSEs, LFEs et VLFEs sont différentes manifestations d'un même phénomène (à différentes échelles), distinct des séismes. Des études postérieures ont remis en question la pente (moment proportionnel à la durée) de la Figure I.6 pour les SSEs, la rapprochant plus ou moins de la pente observée pour les séismes (moment proportionnel au cube de la durée) [Gomberg et al., 2016; Michel et al., 2019]. Mais reste entière la question du gap d'observations entre « séismes lents » et « séismes rapides » (Figure I.6). Dans une gamme de magnitudes et de durées en théorie observables par l'instrumentation sismique et géodésique, pourquoi n'observe-t-on rien ?

# 3 Est-il possible de prédire les séismes ?

## 3.1 Approche physique

Une approche naïve pour prédire les séismes serait de modéliser l'ensemble du cycle sismique (i.e. les phases post-sismique, inter-simique, co-sismique ainsi que les épisodes transitoires) et de « calibrer » les différents paramètres physiques contrôlant le comportement des failles selon les observations enregistrées sur chaque faille. Une modélisation calibrée finement selon les propriétés d'une région donnée pourrait, en théorie, décrire le comportement passé et futur d'une faille. Si aucune prédiction n'a encore jamais pu être réalisée ainsi, des modèles de cycle sismique de ce type existent. Leur formalisme, dit *rate-and-state*, est basé sur des expériences réalisées en laboratoire [Marone, 1998]. Plusieurs lois ont été proposées [Dieterich, 1979; Ruina, 1983; Marone, 1998]. La plus populaire s'écrit

$$ \mu = \mu_0 + a \ln \frac{V}{V_0} + b \ln \frac{\theta}{\theta_0}, \qquad (I-1) $$

où $\mu$ est le coefficient de friction à l'équilibre (*steady-state*), $V$ est la vitesse de glissement à l'équilibre, $a$ et $b$ sont des scalaires, l'indice 0 indique un état de référence et $\theta$ est une variable d'état suivant la loi, dite de vieillissement (*aging law*) – une autre loi très utilisée est dite de
glissement (slip law) – suivante :

$$ \frac{d\theta}{dt} = 1 - \frac{V\theta}{d_c}, \qquad (I-2) $$

$d_c$ étant une distance de glissement caractéristique.

À l'état d'équilibre, la loi ne dépend plus de la variable d'état, mais uniquement de la vitesse de glissement $V$ :

$$ \mu_{eq} = \mu_0 + (a - b) \ln \frac{V}{V_0}. \qquad (I-3) $$

Une analyse de stabilité de cette loi montre que le glissement ne peut être stable que si $a-b > 0$ et qu'à l'inverse il ne peut être instable que si $a-b < 0$ [Scholz, 2019]. Un glissement instable se produit si $a-b$ est inférieur à une valeur critique négative, ce qui définit un domaine intermédiaire « conditionnellement stable ». Dans le cas d'un crack de taille $L$ dans un milieu élastique avec un module de cisaillement $G$, cette condition d'instabilité s'écrit

$$ b - a > \frac{\lambda G d_c}{L(\sigma_n - p)} \qquad (I-4) $$

où $\lambda$ est de l'ordre de 1, $\sigma_n$ est la contrainte normale et $p$ la pression de pore [Avouac, 2015].

Ce formalisme dérivé d'expériences en laboratoire a permis le développement de modèles numériques du cycle sismique. Ces modèles reproduisent un certain nombre d'aspects caractérisants le comportement des failles [Tse and Rice, 1986; Dieterich, 1987, 1994], en particulier la décroissance logarithmique du glissement post-sismique avec le temps [Marone et al., 1991]. De manière intéressante, ils prédissent également que les séismes s'initient avec une phase de nucléation, puisque le crack doit attendre une taille critique avant que le glissement ne devienne instable [Dieterich, 1987; Rubin and Ampuero, 2005; Ampuero and Rubin, 2008].

En considérant des distributions de propriétés frictionnelles hétérogènes le long de la faille (i.e. des zones où $a-b > 0$ et des zones où $a-b < 0$) (Figure I.7), les modèles récents sont capable de reproduire des distributions de glissement très complexes [e.g., Kaneko et al., 2010; Barbot et al., 2012; Ader et al., 2014; Jiang and Lapusta, 2016] cohérents avec les modèles de glissement co-, post- et inter-sismique obtenus par inversion des observations en surface.

Dans ce cadre conceptuel, une connaissance fine de la distribution des propriétés frictionnelles le long des failles pourrait ainsi, en théorie, permettre aux modèles de prédire le comportement passé, présent et futur des failles, et ainsi de prédire les grandes ruptures. Une telle connaissance est toutefois tout-à-fait hors d'atteinte à partir des observations dont l'on dispose actuellement et il est peu probable que l'on ait un jour accès au degré d'information nécessaire à la prédiction à partir d'observations de surface. L'approche physique s'avère ainsi, en pratique,
<figure><img src="image_9.png" /></figure>

**Figure I.7** – Modèle conceptuel de comment la distribution des propriétés frictionnelles – en zones dites *rate-weakening* (i.e ou $a - b < 0$) et en zones dites *rate-strengthening* (i.e ou $a - b > 0$) – le long d'une faille contrôle le comportement de la faille lors des phases co-, post- et inter-sismique [Avouac, 2015].

incapable de prédire les séismes.

## 3.2 Approche statistique

Une autre approche naïve pour prévoir les séismes serait de décrire statistiquement la sismicité afin d'extrapoler l'occurrence de futurs évènements à partir de séquences sismiques passées. Une caractéristique évidente lorsque l'on observe la sismicité, est qu'après un grand évènement, un certain nombre de répliques se produisent. À partir d'observations suivant le séisme de Nobi (1891), *Omori* [1894] a décrit ce phénomène de manière empirique et conclu que le nombre de répliques décroit en $\frac{K}{t+c}$, où K et c sont des constantes [Utsu et al., 1995]. Utsu [1961] a modifié cette célèbre loi d'Omori, observant que les séquences de répliques décroissent avec différents exposants. D'autres lois alternatives ont été proposées [Mignan, 2015, 2016]. Aujourd'hui, l'activité des répliques est généralement décrite comme un processus en cascade et le modèle ETAS (Epidemic-Type Aftershock Sequence) [Ogata, 1988, 1998] donne la meilleure description de l'ac-
tivité sismique [Marzocchi et al., 2017]. Même si ces modèles sont purement empiriques, des explications physiques, souvent dérivées d’expériences en laboratoire [e.g., Toda et al., 2005], ont été proposées pour expliquer des augmentations du taux de sismicité pouvant aller jusqu’à un facteur 1000 dans certaines régions ou même engendrer des diminution dans d’autres [Wyss and Wiemer, 2000; Gerstenberger et al., 2005].

Une autre loi empirique célèbre, la loi de Gutenberg-Richter [Gutenberg and Richter, 1944] donne le nombre N de séismes de magnitude supérieure ou égale à M par la formule

$$ \log N = a - bM, \qquad (I-5) $$

où a et b n'ont rien à voir avec les constantes du formalisme *rate-and-state* décrit dans le précédent paragraphe mais où a décrit le volume global de la sismicité d'un catalogue donné (i.e. le nombre d'évènements recensé dans une région donnée au cours d'une période donnée) et b décrit la proportion entre le nombre de grands séismes et de petits séismes. La valeur de b – connue sous le terme de *b-value* – est généralement proche de 1, ce qui veut dire que pour 1 séisme de magnitude 9, il y a 10 séismes de magnitude 8, 100 de magnitude 7 (etc.). Une valeur légèrement inférieure à 1 signifie qu'il y a une plus grande proportion de séismes de forte magnitude, une valeur légèrement supérieure à 1 signifie l'inverse. Des expériences en laboratoires ont montré un lien entre la *b-value* et l'état de contraintes; plus les contraintes différentielles appliquées ($\sigma_1 - \sigma_3$) sont élevées plus la *b-value* est faible [Scholz, 1968; Amitrano, 2003; Goebel et al., 2013]. Des observations faites dans divers contextes tectoniques semblent cohérentes avec un contrôle de l'état de contraintes sur b [Schorlemmer et al., 2005], et suggèrent également des liens entre la valeur de b et le mécanisme au foyer [Gulia and Wiemer, 2010], la profondeur [Spada et al., 2013] et la pression fluide [Bachmann et al., 2012].

Les variations spatiales et temporelles de b ont suscité un grand intérêt dans la communauté car si celles-ci sont reliées à un changement d'état de contraintes, elles pourraient nous renseigner sur l'imminence d'un grand évènement. Une augmentation de la *b-value* a ainsi été observée après quelques grands séismes [Wiemer and Katsumata, 1999; Wiemer et al., 2002; Tormann et al., 2012; Ogata and Katsura, 2014; Tormann et al., 2014, 2015; Tamaribuchi et al., 2018], puis de manière plus systématique [Gulia et al., 2018]. De manière intéressante, les quelques exceptions – i.e. des séismes suivis d'une diminution de la *b-value* – ont été suivis d'évènements de magnitude plus forte [Gulia and Wiemer, 2019]. Cette observation suggère qu'un monitoring en temps réel de la *b-value* pourrait aboutir à un système de prévention des grands évènements (dont le cadre d'application resterait à inventer). Pour être validé statistiquement, ce système devra toutefois être confronté à un grand nombre de données, lesquelles prendront des an-
nées à être acquises [Gulia and Wiemer, 2019]. Malgré des éléments encourageants, l'approche statistique s'avère ainsi, à l'heure actuelle, incapable de prévoir les séismes.

## 3.3 Approche observationnelle

Devant l'échec des approches physique et statistique à prédire les grands séismes, l'approche observationnelle apparaît comme une solution simple et peut-être plus prometteuse. Si l'on arrivait à observer un signal, quel qu'il soit, systématiquement avant un grand séisme, on pourrait se baser sur l'occurrence de ce signal pour prévoir l'activité sismique à venir. Cette idée simple est à l'origine de la quête des signaux précurseurs. L'existence de tels signaux est hypothétique et dépend de l'existence d'une phase précurserice. Il est possible que les séismes soient des phénomènes de nature chaotique [e.g., Geller, 1997; Kagan, 1997]. Dans ce cas, aucune phase précurserice n'existe et aucun signal précurseur non plus. Si, à l'inverse, les séismes sont des phénomènes de nature déterministe [e.g., Dieterich, 1987; Rubin and Ampuero, 2005; Ampuero and Rubin, 2008], alors les signaux précurseurs existent. Dans la deuxième hypothèse, tout l'enjeu est de les détecter.

Si aucune détection systématique de signaux précurseurs n'a jamais pu être démontrée, de nombreuses observations sur des séismes ponctuels ont été proposées [Bouchon et al., 2011, 2013; Kato et al., 2012; Ruiz et al., 2014; Brodsky and Lay, 2014; Ellsworth and Bulut, 2018]. Le cadre conceptuel sous-jacent est que les grands séismes s'initient avec une phase de glissement a-sismique lent sur la faille, associée à une augmentation de la micro-sismicité dans la zone qui en viendra à produire la rupture [Bouchon et al., 2011, 2013; Kato et al., 2012; Ruiz et al., 2014]. Ces études suggèrent que les précurseurs sismiques existent et que les grands séismes pourraient être anticipés des heures [Bouchon et al., 2011], des jours [Kato et al., 2012], des semaines [Ruiz et al., 2014] ou des mois [Bouchon et al., 2013] avant qu'ils ne se produisent.

Si ces observations sont intéressantes, elles se limitent à un petit nombre d'évènements. L'absence de généralisation à l'ensemble – ou à défaut à un grand nombre – des évènements connus rend impossible la mise en place d'un système de prévision, d'autant que des signaux similaires associés à des glissements a-sismiques transitoires se produisent couramment sans provoquer de grands séismes [e.g., Schwartz and Rokosky, 2007; Gomberg et al., 2010; Obara and Kato, 2016; Wallace, 2020].
# 4 Anticiper les séismes autrement : l’alerte précoce

Détecter des signaux précurseurs, i.e. des signaux précédant l'origine des séismes, est un objectif très (trop ?) ambition, qui n'a jamais mené à un système de prévention opérationnel. Des signaux précoces suivant directement l'origine des séismes – mais ne les précédant pas – peuvent en revanche être utilisés pour alerter les populations de manière concrète. L'exploitation de ces signaux précoces non-précurseurs permet en effet d'alimenter des systèmes d'alerte sismique. Les algorithmes dits d'alerte estiment la magnitude et la localisation des évènements dès les premiers signaux recueillis sur les capteurs les plus proches de la source afin de lancer (ou non) l'alerte dans des régions habitées, plus éloignées, avant que les premières secousses n'y soient ressenties et avant qu'un éventuel tsunami ne frappe les côtes. L'alerte précoce est ainsi moins ambitieuse que la prédiction car elle ne cherche pas à connaître le futur, seulement à déduire très vite des signaux disponibles une information utile pour alerter les populations.

## 4.1 L'alerte sismique

Si énormément d'espoir a été mis dans la prédiction dans le début des années 1970 [Scholz et al., 1973], sa faisabilité a été sérieusement remise en question par des études théoriques sug-
gérant que les séismes sont un processus non-linéaire sensibles à des variations non mesurables
des conditions physiques en profondeur [Geller, 1997; Kagan, 1997]. Les approches probabilistes
(statistiques) développées par la suite n'ont jamais apporté suffisamment d'information pour
être implémentées dans un système d'alerte opérationnel [Jordan et al., 2011; Zechar et al.,
2016]. Alors que l'espoir de la prédiction des tremblements de terre s'est progressivement es-
tompé, le domaine de l'alerte sismique précoce a pris de plus en plus d'importance. Il a donné
vie à une variété de techniques aujourd'hui implémentées dans des systèmes d'alerte opération-
nels dans de nombreux pays exposés au risque sismique et tsunami. L'alerte sismique repose
sur le fait que les ondes sismiques se propagent à une vitesse finie, très rapide (quelques km/s),
mais beaucoup plus lente que les messages d'alerte, les systèmes de télécommunication trans-
mettant à la vitesse de la lumière. On peut donc utiliser des signaux sismiques enregistrés par
des stations sismiques situées proche de l'épicentre pour alerter des régions habitées situées plus
loin de la source avant que les ondes sismiques n'arrivent. De plus, les plus forts mouvements
du sol sont produits par les ondes S qui se propagent plus lentement que les premières ondes
enregistrées, les ondes P, ce qui donne un peu plus de temps pour l'alerte (Figure I.8). On ne
parle pour autant que de quelques secondes – dans le meilleur des cas quelques dizaines de
secondes – de temps d'alerte avant les secousses.
<figure><img src="image_10.png" /><figcaption>**Figure I.8** – Illustration des principes de l’alerte sismique (ETH Zurich)</figcaption></figure>

L’alerte sismique peut se décomposer en 4 catégories d’approches.

(1) Les approches in situ utilisent un capteur à un endroit donné pour alerter directement la région alentour [Böse et al., 2009; Meier et al., 2015]. Ces approches peuvent s’avérer efficaces sur une région étendue si de nombreux capteurs sont disponibles, c’est l’ambition de l’application MyShake qui utilise les smartphones comme capteurs sismiques [Kong et al., 2016a,b].

(2) Les algorithmes « point-source », tels que PRESTo [Satriano et al., 2011] ou ElarmS [Serdar Kuyuk et al., 2014; Sheen et al., 2017; Chung et al., 2019], procurent des estimations plus rapides des caractéristiques du séisme (magnitude, localisation) en utilisant les toutes premières secondes d’ondes P enregistrées.

(3) Les algorithmes de modélisation des mouvements du sol, tels que FinDer [Böse et al., 2012, 2015, 2018] ou PLUM [Hoshiba and Aoki, 2015], prédisent les mouvements du sol à venir à partir des mouvements du sol observés. Leurs performances sont meilleures que celles des algorithmes « point-source ».

Ces 3 familles d’algorithmes fonctionnent relativement bien dans des systèmes régionaux, et pour des séismes de magnitudes modérées, car ils sont assez rapides. Ils ont, en revanche, tendance à systématiquement sous-estimer la magnitude des grands séismes (magnitude $M_w$ > 7). En d’autres termes, l’estimation de la magnitude sature à $M_w$ = 7 : pour des séismes de magnitude supérieure l’estimation sera $M_w$ = 7. Ce problème a pour origine principale la saturation des instruments ; l’amplitude des ondes P enregistrées proche de la source d’un séisme
de magnitude supérieure à 7 est supérieure à ce que les accéléromètres (instruments les plus performants pour enregistrer les mouvements forts du sol) sont capables de mesurer. S’ajoute à ce problème instrumental, un problème de modélisation : l’approximation « point-source » devient, en effet, largement fausse proche de la source pour des séismes d’une telle magnitude.

(4) Les algorithmes « faille-finie », comme G-larmS [Grapenthin et al., 2014a,b], G-FAST [Crowell et al., 2016, 2018], BEFOREES [Minson et al., 2014] ou REGARD [Kawamoto et al., 2016, 2017], inversent le glissement sur la faille à partir des données GNSS. Ces algorithmes règlent en partie le problème de saturation de l’estimation de la magnitude mais sont significativement plus lents que les 3 autres types d’algorithmes [Ruhl et al., 2019; Allen and Melgar, 2019]. Il y a ainsi un compromis à trouver entre la rapidité de l’estimation et sa fiabilité [Meier, 2017].

## 4.2 L’alerte tsunami

Le problème de saturation de l’estimation de la magnitude est particulièrement problématique, dans le cas des grands séismes de subduction, pour évaluer l’amplitude du tsunami, lequel peut être beaucoup plus destructeur que les secousses dans le cas de séismes de magnitude supérieure à 9, ou inexistant pour des séismes de magnitude inférieure à 8. L’estimation de la magnitude des séismes est ainsi cruciale pour prévoir l’ampleur du tsunami, mais les algorithmes d’alerte sismique ont des difficultés à estimer rapidement et de manière fiable la magnitude des grands évènements. Les tsunamis atteignent les côtes entre 30 min et 2 h après le séisme, ce qui laisse plus de temps pour l’alerte que pour les secousses. Cependant, réduire le temps d’alerte est essentiel pour arrêter les infrastructures sensibles, telles que les centrales nucléaires, et évacuer les régions côtières densément peuplées.

Le calcul de la magnitude d’un séisme est normalement réalisé à partir des ondes de surface enregistrées partout sur Terre [Kanamori and Given, 1981], mais ces ondes se propagent beaucoup plus lentement que les ondes P. Le calcul « classique » de la magnitude ne peut donc pas être effectué avant l’arrivée du tsunami. Les systèmes d’alerte tsunami modernes utilisent ainsi un autre type d’onde, de faible amplitude et de très longue période, appelé phase W [Kanamori and Rivera, 2008], qui peut être enregistré plus tôt que les ondes de surface (Figure I.9) et procure de meilleurs estimations de la magnitude que les systèmes d’alerte basés sur les ondes P [Duputel et al., 2011, 2012]. Les algorithmes basés sur la phase W permettent ainsi une estimation plus fiable de la taille des tsunamis mais sont beaucoup plus lents (≈ 20 min) que les algorithmes basés sur les ondes P [Duputel et al., 2011].
<figure><img src="image_11.png" /><figcaption>**Figure I.9** – Sismogramme du séisme péruvien de 2001 (M<sub>w</sub> 8.4) (black line) indiquant les temps d'arrivée des ondes P, S et des ondes de surface, enregistré à la station HRV, située à une distance épicentrale de 58.5°. La courbe rouge montre une forme d'onde synthétique de la phase W [Kanamori and Rivera, 2008].</figcaption></figure>

## 4.3 La découverte prometteuse des PEGS

Les déformations statiques et transitoires générées par les séismes entrainent des perturbations de densité qui, à leur tour, provoquent des perturbations à grande échelle du champ de gravité terrestre de manière instantanée [Harms et al., 2015; Harms, 2016]. Ainsi, les signaux les plus précoces générés par un séisme ne sont pas les ondes P mais les perturbations du champ de gravité terrestre – i.e. des variations infimes de l'accélération de la pesanteur *g* – qui se propagent à la vitesse de la lumière. Cependant, ces perturbations sont trop faibles pour être observées systématiquement sur tous les séismes, ce qui a, jusqu'ici, considérablement limité leur utilité pour l'alerte précoce.

La première observation d'un signal dans le champ de gravité terrestre induit par un séisme a été faite par un gravimètre superconducteur lors du séisme de Tohoku-Oki (M<sub>w</sub> 9.0) survenu en 2011 au large du Japon [Montagner et al., 2016]. Cette découverte a été suivie de multiples observations sur des sismogrammes enregistrés lors du même séisme (Figure I.10) – la gravité étant une accélération, les variations de *g* peuvent être mesurées sur un accélérogramme ou un vélocigramme dérivé (plus sensible aux basses fréquences). Dans le détail, les signaux observés sur ces sismogrammes sont composés d'un terme direct, i.e. la perturbation de *g*, et d'un terme induit par la réponse élastique de la Terre à la perturbation gravitationnelle. Ce terme induit est de signe opposé et d'amplitude similaire au terme direct, « annulant » une grande partie du signal [Vallée et al., 2017]. La résultante des 2 termes a été nommée *Prompt Elasto-Gravity Signal* (PEGS). Découpler les termes direct et induit rendrait le signal plus facilement
observable mais la technologie nécessaire à ce découplage n'est malheureusement pas encore disponible.

<figure><img src="image_12.png" /><figcaption>**Figure I.10** – PEGS observés (rouge) et prédits (noir) pour le séisme de Tohoku-Oki (étoile rouge) à différentes stations (triangles noirs sur la carte). Les formes d'ondes sont coupées à l'arrivée des ondes P. Les courbes pointillées bleues montrent le signal prédit pour un séisme de magnitude M<sub>w</sub> 8.5 [Vallée et al., 2017].</figcaption></figure>

Les PEGS montrent toutefois un fort potentiel pour l’alerte précoce puisque ces signaux présentent des caractéristiques exceptionnelles : ils sont très sensibles à la magnitude, ils ne saturent pas (au contraire des ondes P dans les systèmes « classiques »), et surtout ils se propagent à la vitesse de la lumière. Cela veut dire que des signaux très précoces, de faible amplitude, contenant l’information de la magnitude des grands séismes – ce que les premières secondes des ondes P ne contiennent probablement pas [Meier et al., 2017] – existent dans les formes d’ondes sismiques, cachés dans le bruit. Extraire ces signaux du bruit sismique en temps réel est cependant un réel challenge. L’observation directe des PEGS lors du séisme de Tohoku n’a été rendue possible que par la magnitude exceptionnelle de l’évènement (M<sub>w</sub> 9.0) conjuguée à la qualité des instruments de mesure installés au Japon. Les séismes de magnitude 9 sont rares et ne se produisent pas toujours dans des régions aussi bien instrumentées. Un séisme de magnitude 8.5 aurait été difficilement observable par le même réseau sismique (Figure I.10). Une recherche systématique des PEGS générés par les grands séismes récents a été réalisée et n’a été concluante que pour 5 évènements [Vallée and Juhel, 2019].
La faisabilité de la détection rapide et systématique des PEGS n’est ainsi pas démontrée. Une piste prometteuse est le développement d’une nouvelle génération d’instruments spécifiquement conçus pour n’enregistrer que le terme gravitationnel direct des PEGS, découplé du terme induit [Juhel et al., 2018]. Nous explorerons dans ce manuscrit une piste alternative : l’utilisation de la technologie émergente de l’Intelligence Artificielle (IA) pour rapidement détecter les PEGS et estimer la magnitude et la localisation des grands séismes. Cette technologie s’avère, en effet, être extrêmement efficace pour extraire du bruit des signaux de faible amplitude.

## 5 Possibilités offertes par l’émergence de l’IA

L’IA – ou plus exactement sa sous-branche appelée Apprentissage Machine (Machine Learning) – offre une variété d’outils pour extraire de l’information de n’importe quel type de données [Jordan and Mitchell, 2015]. Les algorithmes de Machine Learning – que l’on appellera improprement algorithmes d’IA par simplicité – sont conçus pour apprendre à reconnaître des structures dans les données à partir d’un grand nombre d’exemples. Ces algorithmes suivent une approche différente des techniques classiques d’analyse de données, une approche qui est robuste, rapide et permet l’exploration d’espaces de très grandes dimensions. Ces caractéristiques rendent l’IA particulièrement adaptée à la résolution de nombreux problèmes en Sciences de la Terre [Bergen et al., 2019]. Les algorithmes d’IA sont utilisés en sismologie depuis plus de 30 ans [Dowla et al., 1990; Dysart and Pulli, 1990; Dai and MacBeth, 1995] mais ont été accueillis avec peu d’enthousiasme par la communauté dans les années 1990. Le développement exponentiel des capacités de calcul a depuis révélé le potentiel l’IA.

### 5.1 Principes simplifiés de l'IA à l'usage des non-spécialistes

Si la théorie et le développement de schémas algorithmiques nouveaux sont des domaines en-
core réservés aux spécialistes, l'utilisation d'algorithmes d'IA « pré-implémentés » s'est considé-
rablement démocratisée ces dernières années grâce au développement de librairies d'algorithmes
en libre d'accès (e.g., Tensorflow, Keras, PyTorch). Écrire un algorithme d'IA est aujourd'hui
largement accessible à un non-spécialiste. Cette démocratisation ouvre la voie à des découvertes
majeures dans tous les champs disciplinaires. Elle procure un nouvel outil que des spécialistes
de tous les domaines peuvent s'approprier pour résoudre des problèmes qui n'ont, pour l'heure,
pas pu être résolus par les techniques « standards ».

Si l'apprentissage non-supervisé est un domaine en développement, relativement abstrait, et
encore peu accessible aux non-spécialistes, le principe des algorithmes d'apprentissage supervisé
(Supervised) Machine Learning

<figure><img src="image_13.png" /><figcaption>**Figure I.11** – Principes de l'apprentissage machine supervisé (Supervised Machine Learning).</figcaption></figure>

– de loin les plus utilisés aujourd’hui – peut être résumé de manière assez simple : il permet de prédire n’importe quel type de sortie (scalaire, vecteur, matrice, image, son, vidéo...) à partir de n’importe quel type d’entrée, en apprenant les liens entre entrée et sortie sur un grand nombre d’exemples. Lors de la phase dite d’entrainement, l’algorithme va chercher dans un espace de grandes dimensions des « liens » entre les jeux de données en entrée et en sortie. Si l’apprentissage est concluant, l’algorithme sera capable lors de la phase de prédiction, de prédire la sortie à partir d’une entrée qu’il n’a jamais vue (Figure I.11). Les 2 seules conditions pour que l’apprentissage fonctionne sont (1) d’avoir suffisamment d’exemples (au moins plusieurs dizaines de milliers) non-redondants et (2) que les jeux de données en entrée et en sortie ne soient pas indépendants.

L’exemple le plus utilisé pour illustrer l’apprentissage supervisé est celui de la reconnaissance d’image. Les données en entrée sont dans l’exemple que l’on prendra des photographies de chiens et de chats. Les données en sortie associent à chaque photographie un scalaire, 0 pour les photographies de chiens, 1 pour celles de chats. Une fois entrainé sur des millions d’exemples, l’algorithme s’avère capable de « reconnaître » des chiens et des chats qu’il n’a jamais vus, i.e. il prédira 0 ou 1 de manière correcte à partir d’une donnée d’entrée absente du jeu d’entrainement. Cet exemple laisse entrevoir les possibilités immenses de l’IA. Celle-ci permet en effet de
répondre à n'importe quelle question formulée en des termes d'entrée et sortie si l'on dispose de suffisamment de données.

## 5.2 Accomplissements récents en sismologie

Les succès de l'IA dans de nombreuses applications couvrant des disciplines très différentes (reconnaissance d'image, traduction automatique, etc.) ont entraîné un regain d'intérêt dans la communauté sismologique. Bergen et al. [2019] ont classé les applications récentes de l'IA en Sciences de la Terre solide en 3 catégories.

(1) L'automatisation. L'IA peut, par exemple, être utilisée pour automatiser des taches comme la détection de séismes [Zhang et al., 2014; Yoon et al., 2015; Tibi et al., 2017; Bergen and Beroza, 2018; Rong et al., 2018; Perol et al., 2018; Ross et al., 2018a], le pointé des temps d'arrivée des ondes P et S – une tache toujours largement réalisée manuellement aujourd'hui – [Ross et al., 2018b; Zhu and Beroza, 2019], l'association des phases sismiques [Ross et al., 2019a], et a même été utilisée dans des algorithmes d'alerte précoce [Reddy and Nair, 2013; Kong et al., 2016b; Ochoa et al., 2018; Yin et al., 2018; Lin et al., 2021]. L'automatisation par IA donne des résultats plus rapides, plus robustes et moins biaisés que les outils standards, ce qui permet la constitution de catalogues de sismicité beaucoup plus complets [Ross et al., 2019b].

(2) La modélisation. La modélisation par IA inclue des applications comme la prédiction des mouvements du sol à partir des paramètres de la source sismique [Derras et al., 2014; Trugman and Shearer, 2018], l'accélération des calculs de simulations numériques [DeVries et al., 2017; Valera et al., 2018] ou encore des problèmes d'inversion [Käufl et al., 2016; McCann et al., 2017].

(3) La découverte. Si l'apprentissage machine procure des outils plus performants que les techniques standards d'automatisation et de modélisation, le plus gros potentiel de l'IA se situe probablement dans la data-driven discovery, i.e. la découverte orientée par la donnée. L'IA procure en effet un moyen unique d'extraire de l'information dans de gros volumes de données qu'aucune autre technique ne permet. L'IA peut ainsi permettre de découvrir des signaux ou des processus physiques inconnus [Holtzman et al., 2018; Rouet-Leduc et al., 2018, 2019; Hulbert et al., 2019] et même prédire des séismes en laboratoire [Rouet-Leduc et al., 2017]...
# Chapitre II

## Où se produisent les grands séismes ?

### 1 Liens entre failles et séismes

La magnitude d'un séisme dépend principalement de la surface de la faille rompue. Pour faire un grand séisme il faut donc une grande faille. La première réponse à la question « Où se produisent les grands séismes » est donc : sur les grandes failles (Figure II.1).

<figure><img src="image_14.png" /><figcaption>**Figure II.1** – Distribution des grandes failles et des grands séismes à travers le monde [Hicks, 2015].</figcaption></figure>

Toutefois, la relation entre failles et séismes est plus profonde et complexe que cela. Un séisme n'est pas un objet physique hors sol mais l'un des modes de glissement entre deux blocs de Terre le long d'une faille. Les propriétés des failles conditionnent ainsi les propriétés fondamentales des séismes, comme la chute de contraintes [Cao and Aki, 1986; Choy and Kirby, 2004; Manighetti
et al., 2007; Hecker et al., 2010], la longueur de rupture et l'amplitude du glissement [Manighetti et al., 2007], le nombre de segments rompus [Wesnousky, 1988; Manighetti et al., 2007, 2009], la vitesse de la rupture [Perrin et al., 2016] et même les mouvements forts du sol [Radiguet et al., 2009] ainsi que la b-value [Stirling et al., 1996].

Peut-on alors déduire des propriétés des failles où les prochains grands séismes sont le plus susceptibles de se produire ? Quelles propriétés des failles ont-elles un contrôle significatif sur la magnitude des grands évènements ? Schellart and Rawlinson [2013] ont passé en revue 24 paramètres physiques caractérisant les plus grandes failles de la planète – les failles de subduction – découpées pour l'occasion en 44 segments. Il n'ont trouvé que de faibles corrélations entre ces paramètres et la magnitude des séismes historiques, et pour seulement 6 des paramètres testés.

## 2 Les grands séismes se produisent sur les grandes failles planes

Toutes les expériences en laboratoire – à partir desquelles ont été développées les équations utilisées pour modéliser la rupture et le cycle sismique (voir Chapitre I.3.1) – ont été réalisées sur des échantillons de roche de quelques centimètres et présentant une faille pré-existante plane. Or les failles réelles – en particulier les plus grandes d'entre elles, les failles de subdcution – sont longues de milliers de kilomètres et ne sont pas planes. Imaginez que l'on ajoute une petite aspérité géométrique dans la topographie d'une faille centimétrique en laboratoire. On intuitue facilement que le comportement d'une telle micro-faille va être différent du comportement de la même faille plane. Comment cette intuition pourrait-elle se traduire sur des failles naturelles, longues de milliers de kilomètres ? J'ai voulu testé cette idée en calculant le gradient du pendage – une mesure de la topographie – le long des grandes failles de la planète. Il est apparu que les grands séismes historiques se sont produits sur des portions de failles particulièrement planes (à faible gradient de pendage). Notre interprétation de cette observation est que les différentes portions d'une faille ont plus de chances de rompre de manière synchrone, si celle-ci est plane, rendant ainsi le système plus homogène.

Cette étude a été réalisée à l'Université d'Oregon en collaboration avec Amanda M. Thomas, Alan W. Rempel, Leif Karlstrom, Anthony Sladen et Louis De Barros. Elle a été publiée en 2016 dans le journal Science [Bletery et al., 2016b].
<figure><img src="image_15.png" /></figure>

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

Science

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

AAAS

**Mega-earthquakes rupture flat megathrusts**

Quentin Bletery, Amanda M. Thomas, Alan W. Rempel, Leif Karlstrom, Anthony Sladen and Louis De Barros (November 24, 2016)

*Science* **354** (6315), 1027-1031. [doi: 10.1126/science.aag0482]

**Editor's Summary**

**Mega-earthquakes go the flat way**

Megathrust faults in subduction zones cause large and damaging earthquakes. Bletery *et al*. argue that certain geometric features of the subduction zones relate to earthquake size. The key parameter is the curvature of the megathrust. Larger earthquakes occur where the subducting slab is flatter, providing a rough metric for estimating where mega-earthquakes may occur in the future.

*Science*, this issue p. 1027

This copy is for your personal, non-commercial use only.

**Article Tools**
Visit the online version of this article to access the personalization and article tools:
http://science.sciencemag.org/content/354/6315/1027

**Permissions**
Obtain information about reproducing this article:
http://www.sciencemag.org/about/permissions.dtl

*Science* (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright 2016 by the American Association for the Advancement of Science; all rights reserved. The title *Science* is a registered trademark of AAAS.
in combination with periodic reservoir reloading from a cold atom source (such as a MOT), could be used to maintain arrays indefinitely.

Atom-by-atom assembly of defect-free arrays forms a scalable platform with unique possibilities. It combines features that are typically associated with ion-trapping experiments, such as single-qubit addressability (32, 33) and fast cycling times, with the flexible optical trapping of neutral atoms in a scalable fashion. Furthermore, in contrast to solid-state platforms, such atomic arrays are highly homogeneous (31) and mostly decoupled from their environment. The homogeneity of our array should also allow for cooling of the atomic motion via simultaneous sideband cooling in all tweezers at once (34, 35).

These features provide an excellent starting point for multiqubit experiments, for studies of quantum many-body effects, and for exploring future applications. The required interactions between the atoms can be engineered using several approaches. Even without sideband cooling, exciting the atoms into high-lying Rydberg states would introduce strong dipole interactions that can be used for fast entangling gates (24, 25, 27). The parallelism afforded by our flexible atom rearrangement enables efficient diagnostics of such Rydberg-mediated entanglement. These interactions may also enable approaches to quantum simulations that involve both coherent coupling and engineered dissipation (26, 27), as well as large-scale entangled quantum states for applications in precision measurements (36).

An alternative approach to engineering interactions involves the integration of atom arrays with nanophotonic platforms as demonstrated previously (28, 29). These enable photon-mediated interactions that can be employed to couple the atoms within a local multiqubit register or for efficient communication between the registers using a modular quantum network architecture (3).

Finally, our platform could enable new bottom-up approaches to studying quantum many-body physics in Hubbard models (15, 16, 30), where atomic Mott insulators with fixed atom number and complex spin patterns could be directly assembled. This requires atom temperatures close to the ground state, coherent tunneling between the traps, and sizable on-site interactions. With side-band cooling, ground-state fractions in excess of 90% have already been demonstrated (34, 35) and can likely be improved via additional optical trapping along the longitudinal tweezer axes, which would also increase on-site interaction strengths. Coherent tunneling of Rb atoms between similarly sized tweezers has been observed before by reducing the tweezer distance (15, 16). The parametric heating, currently limiting the minimal distance between our traps, could be reduced by working with shallower traps, as needed for tunneling, and by employing fewer traps to increase the frequency separation between neighboring traps. Eventually, this approach could be applied to create ultracold quantum matter composed of exotic atomic species or complex molecules (37, 38) that are difficult to cool evaporatively.

## REFERENCES AND NOTES

1. S. Haroche, *Ann. Phys.* **525**, 753–776 (2013).

2. D. J. Wineland, *Rev. Mod. Phys.* **85**, 1103–1114 (2013).

3. C. Monroe, J. Kim, *Science* **339**, 1164–1169 (2013).

4. M. H. Devoret, R. J. Schoelkopf, *Science* **339**, 1169–1174 (2013).

5. W. S. Bakr et al., *Science* **329**, 547–550 (2010).

6. J. F. Sherson et al., *Nature* **467**, 68–72 (2010).

7. C. Weitenberg et al., *Nature* **471**, 319–324 (2011).

8. D. S. Weiss et al., *Phys. Rev. A* **70**, 040302 (2004).

9. J. Vala et al., *Phys. Rev. A* **71**, 032324 (2005).

10. N. Schlosser, G. Reymond, I. Protsenko, P. Grangier, *Nature* **411**, 1024–1027 (2001).

11. M. Weber, J. Volz, K. Saucke, C. Kurtsiefer, H. Weinfurter, *Phys. Rev. A* **73**, 043406 (2006).

12. K. D. Nelson, X. Li, D. S. Weiss, *Nat. Phys.* **3**, 556–560 (2007).

13. M. J. Piotrowicz et al., *Phys. Rev. A* **88**, 013420 (2013).

14. F. Nogrette et al., *Phys. Rev. X* **4**, 021034 (2014).

15. A. M. Kaufman et al., *Science* **345**, 306–309 (2014).

16. A. M. Kaufman et al., *Nature* **527**, 208–211 (2015).

17. Y. Miroshnychenko et al., *Nature* **442**, 151 (2006).

18. J. Beugnon et al., *Nat. Phys.* **3**, 696–699 (2007).

19. M. Schlosser et al., *New J. Phys.* **14**, 123034 (2012).

20. H. Kim, et al., *Nat. Commun.* **7**, 13317 (2016).

21. T. Grünzweig, A. Hilliard, M. McGovern, M. F. Andersen, *Nat. Phys.* **6**, 951–954 (2010).

22. B. J. Lester, N. Luick, A. M. Kaufman, C. M. Reynolds, C. A. Regal, *Phys. Rev. Lett.* **115**, 073003 (2015).

23. Y. H. Fung, M. F. Andersen, *New J. Phys.* **17**, 073011 (2015).

24. D. Jaksch et al., *Phys. Rev. Lett.* **85**, 2208–2211 (2000).

25. M. Saffman, T. G. Walker, K. Mølmer, *Rev. Mod. Phys.* **82**, 2313–2363 (2010).

26. H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, H. P. Büchler, *Nat. Phys.* **6**, 382–388 (2010).

27. A. Browaeys, D. Barredo, T. Lahaye, *J. Phys. B* **49**, 152001 (2016).

28. J. D. Thompson et al., *Science* **340**, 1202–1205 (2013).

29. A. Goban et al., *Nat. Commun.* **5**, 3808 (2014).

30. S. Murmann et al., *Phys. Rev. Lett.* **115**, 215301 (2015).

31. See supplementary materials on *Science Online*.

32. T. Xia et al., *Phys. Rev. Lett.* **114**, 100503 (2015).

33. Y. Wang, X. Zhang, T. A. Corcovilos, A. Kumar, D. S. Weiss, *Phys. Rev. Lett.* **115**, 043003 (2015).

34. A. M. Kaufman, B. J. Lester, C. A. Regal, *Phys. Rev. X* **2**, 041014 (2012).

35. J. D. Thompson, T. G. Tiecke, A. S. Zibrov, V. Vuletić, M. D. Lukin, *Phys. Rev. Lett.* **110**, 133001 (2013).

36. P. Kömär et al., *Phys. Rev. Lett.* **117**, 060506 (2016).

37. J. F. Barry, D. J. McCarron, E. B. Norrgard, M. H. Steinecker, *De Mille, Nature* **512**, 286–289 (2014).

38. N. R. Hutzler, L. R. Liu, Y. Yu, K.-K. Ni, https://arxiv.org/abs/1605.09422 (2016).

39. D. Barredo, S. de Léséleuc, V. Lienhard, T. Lahaye, A. Browaeys, *Science* **354**, 1021–1023 (2016).

## ACKNOWLEDGMENTS

We thank K.-K. Ni, N. Hutzler, A. Mazurenko, and A. Kaufman for insightful discussion. This work was supported by NSF, Center for Ultracold Atoms, National Security Science and Engineering Faculty Fellowship, and Harvard Quantum Optics Center. H.B. acknowledges support by a Rubicon Grant of the Netherlands Organization for Scientific Research (NWO). During the completion of this work, we became aware of a related approach (39).

## SUPPLEMENTARY MATERIALS

www.sciencemag.org/content/354/6315/1024/suppl/DC1

Materials and Methods

Figs. S1 to S3

Movies S1 to S3

References (40, 41)

17 June 2016; accepted 17 October 2016

Published online 3 November 2016

10.1126/science.aah3752

## GEOPHYSICS

# Mega-earthquakes rupture flat megathrusts

Quentin Bletery,<fnref n="1" />* Amanda M. Thomas,<fnref n="1" /> Alan W. Rempel,<fnref n="1" /> Leif Karlstrom,<fnref n="1" /> Anthony Sladen,² Louis De Barros²

The 2004 Sumatra-Andaman and 2011 Tohoku-Oki earthquakes highlighted gaps in our understanding of mega-earthquake rupture processes and the factors controlling their global distribution: A fast convergence rate and young buoyant lithosphere are not required to produce mega-earthquakes. We calculated the curvature along the major subduction zones of the world, showing that mega-earthquakes preferentially rupture flat (low-curvature) interfaces. A simplified analytic model demonstrates that heterogeneity in shear strength increases with curvature. Shear strength on flat megathrusts is more homogeneous, and hence more likely to be exceeded simultaneously over large areas, than on highly curved faults.

P
ast mega-earthquakes, such as the magnitude (M) 9.6 Chile earthquake in 1960 and the M 9.3 Alaska earthquake in 1964, occurred in areas where the subducting lithosphere was relatively young (and buoyant) and the plate convergence rate was relatively high (I). These observations led some authors to hypothesize that maximum earthquake

size is controlled by these two geological parameters (2, 3). The development of space-based geodesy enabled refined measurements of plate motion that challenged the role of convergence rate (4–6). Additionally, the moment magnitude (Mw) 9.0 Tohoku-Oki earthquake (7) ruptured lithosphere that is over 120 million years old (8), ruling out lithospheric age as the dominant control. Weak correlations appear in recent data sets among a variety of parameters, including forearc structure (9, 10); age, density, and buoyancy of the slab (6); upper plate motion (II); upper plate strain (I2); long-term trench migration (II); trench sediment thickness (I2); and width of the seismogenic

<footnote n="1">¹Department of Earth Sciences, University of Oregon, 1272 University of Oregon, Eugene, OR 97403, USA. ²Université Côte d'Azur, CNRS, OCA, IRD, Géoazur, 250 rue Albert Einstein, Sophia Antipolis, 06560 Valbonne, France.</footnote>

<footnote>*Corresponding author. Email: qbletery@uoregon.edu</footnote>
<figure><img src="image_18.png" /><figcaption>**Fig. 1.** Variations in megathrust curvature and distribution of historical mega-earthquakes. Along-dip curvature $K_s = d\theta/ds$ on the megathrust is shown for the 13 main subduction zones of the world, overlaid with the estimated slip contours of known historical mega-earthquakes ($M \ge 8.5$) (see table S1 for a list of events). Tick marks in each panel are spaced at 5° increments.</figcaption></figure>

zone (11, 13, 14). However, aside from dip angle and seismogenic width [coefficient of correlation $r = 0.61$ (14)], none of these parameters correlate strongly ($|r| < 0.5$) with the maximum earthquake magnitude recorded in the different zones (11). Stronger correlations are found with fault properties such as thrust interface dip angle (II) and apparent friction derived from heat flow measurements (15).

Today, forecasts of potential large earthquakes focus mainly on imaging the slip deficit with respect to the relative plate motion, because such deficits are thought to result in stress loading that is ultimately released in earthquakes [e.g., (16)]. We followed a complementary approach by analyzing large-scale geometrical features of subduction faults and assessing their possible influence on the physical conditions favoring large earthquake ruptures. Seismic moment $M_0 \propto \Delta\sigma S^{3/2}$ (17) depends mainly on the surface area $S$ over which an earthquake ruptures, because the stress drop $\Delta\sigma$ (the difference in stress before and after an earthquake) is roughly constant among earthquakes over a large range of magnitudes (18). Here we demonstrate that slab interface curvature

exerts a leading-order control on the spatial extent of potential ruptures in subduction zones and hence on the magnitude of the largest earthquakes.

Data sets from active- and passive-source seismology have been combined to constrain the slab1.0 model for the geometry of the world's major subduction zones (14) [details and limitations of the slab1.0 model are discussed in the supplementary materials (19)]. We computed maps of the along-dip interface curvature $K_s = d\theta/ds$, hereafter referred to as curvature [where $\theta$ is the dip angle and $s$ is the tangent to the interface pointing in the down-dip direction (19)], using the slab1.0 model (14, 20) (Fig. 1). Because $K_s$ is almost always positive (Fig. 1), we hereafter interchangeably refer to low-$K_s$ fault regions as flat, low-curvature, or planar. Similarly, we computed the along-strike gradient of the dip angle $K_t = d\theta/dt$ (fig. S1). Comparison with a catalog of historical events (19) reveals a tendency for mega-earthquakes to occur on relatively flat (low-$K_s$) megathrusts (Fig. 1). The curvature is particularly small in the Japan-Kuriles-Kamchatka, Alaska-Aleutians, Sumatra-Java, South America, and Cascadia subduction zones, which are known

to produce $M \ge 9.0$ earthquakes. Subduction zones with large curvatures, such as the Philippines, Solomon Islands, Izu-Bonin, Santa Cruz-Vanuatu-Loyalty, and Tonga-Kermadec zones, lack historic mega-earthquakes. At a smaller scale, we also observe that the down-dip limit of mega-earthquakes offshore of Sumatra—and to a lesser extent, in the Aleutians, Alaska, Nankai, and Kamchatka—coincides with an abrupt change in the slope, suggesting that the conditions required to generate large earthquakes are directly related to the local curvature along the megathrust interface.

To further explore the role of subduction geometry, we calculated the average dip angle $\bar{\theta}$ (Fig. 2A), the average curvature $\bar{K}_s$ (Fig. 2B), and the average along-strike gradient of dip angle $|\bar{K}_t|$ (fig. S2) from the trench to 60 km depth, and we compared these values with the magnitude of the largest megathrust earthquake $M_{max}$ recorded in each subduction zone. We found that $M_{max}$ anticorrelates ($r = -0.72$) with $\bar{\theta}$ (Fig. 2A) and anticorrelates even more strongly ($r = -0.80$) with $\bar{K}_s$ (Fig. 2B). The anticorrelation with $|\bar{K}_t|$ is weaker ($r = -0.64$) but still significant (fig.
<figure><img src="image_19.png" /><figcaption>**Fig. 2.** Correlation of average dip angle and curvature with maximum earthquake magnitude. (A) Average dip angle $\bar{\theta}$ of the 13 main subduction zones of the world as a function of the maximum event magnitude $M_{max}$ recorded on the megathrust (see table S2 for a list of events). The coefficient of correlation *r* is -0.72. (B) Same as (A), but for the average along-dip curvature $\bar{K}_s$ (*r* = -0.80). Light blue bars show the standard deviation for the distribution of $\bar{\theta}$ (A) and $\bar{K}_s$ (B) in each subduction zone. These bars should</figcaption></figure>

not be considered as error bars because we are reasoning with averages; they rather indicate how $\bar{\theta}$ and $\bar{K}_s$ are distributed along the different subduction zones. Blue dashed lines are the linear regressions. Ala, Alaska-Aleutians; CAm, Central America; Cas, Cascadia; Izu, Izu-Bonin; Jap, Japan-Kuriles-Kamchatka; Phi, Philippines; Ryu, Ryukyu-Nankai; SAm, South America; Sco, Scotia arc; Sol, Solomon Islands; Sum, Sumatra-Java; Ton, Tonga-Kermadec; Van, Santa Cruz-Vanuatu-Loyalty.

S2). We suggest, on the basis of our observations, that planar slab interfaces are prone to hosting larger earthquakes. Multiple factors, including the subducting (21) or overriding (22) plate thickness and viscosity contrast (23), have been invoked as controls on subduction curvature; these studies (21–23) provide a useful introduction to this field of research.

The average curvature within mega-earthquake slip contours is consistent with the linear regression in Fig. 2B, with observed values always smaller than $3.76 \times 10^{-6} \text{ m}^{-1}$ (fig. S5B). Given the distribution of $K_s$ across the different megathrusts of the world, the likelihood that all known mega-earthquakes occurred by chance on $K_s \le 3.76 \times 10^{-6} \text{ m}^{-1}$ areas is 0.8% (19), allowing us to affirm with >99% confidence that earthquake magnitude is related to megathrust curvature. The distribution of the average dip angle within mega-rupture contours is much broader (fig. S5A), confirming the stronger relationship of earthquake size with curvature than with dip angle (19). Nevertheless, known mega-earthquakes do not exhibit a strong tendency to rupture the lower-

curvature portions of their host subduction zones (19). This suggests that over the long term, mega-earthquakes might rupture any portion of a low-curvature subduction thrust, such as those of the Cascadia, South America, Alaska-Aleutians, Sumatra-Java, Japan-Kuriles-Kamchatka, Ryukyu-Nankai, and Central America subduction zones. Large earthquakes are preferentially hosted by megathrusts with low dip angles in part as a consequence of larger down-dip seismogenic events making ruptures possible over wider fault areas (II, I3, I4). However, the stronger relationship of earthquake size with curvature than with dip angle noted above (Fig. 2B) warrants a more thorough description of fault loading.

The present understanding of seismic ruptures can be framed in terms of the asperity model: Earthquakes are caused by the sudden failure of locked macroscopic asperities that accumulate stress during interseismic periods (24). In this context, the extent of an earthquake has a leading-order dependence on the size of the rupturing asperity and the neighboring region that is pushed to failure by coseismic stress change.

Recent efforts aimed at characterizing the nature of asperities and their immediate surroundings have mostly focused on potential variations in friction along faults [e.g., (25–28)]. The capacity of a fault segment to accumulate stress is bounded by the Coulomb failure criterion

$$|\tau^c| = \mu(\sigma_n^c - p) \quad (1)$$

where $\tau^c$ is the critical shear stress required to initiate a rupture (hereafter referred to as shear strength), $\mu$ is the coefficient of friction, $p$ is the pore pressure, $\sigma_n^c$ is the normal stress when the rupture initiates, and cohesion is neglected. The Coulomb failure criterion is met when the accumulated shear stress reaches the frictional shear strength and an earthquake initiates. From Eq. 1, one can define an asperity as a fault area characterized by a particularly high friction coefficient $\mu$. Alternatively, one can also explain the time-dependent behavior of asperities by appealing to variations in pore pressure $p$ that accompany fluid migration [e.g., (29)]. We explore here a different hypothesis, which is that the locations of
asperities can be determined by analyzing variations in normal stress $\sigma_n$ that are governed by large-scale geometrical characteristics of megathrust interfaces. Assuming constant friction $\mu$ and hydrostatic pore pressure $p = \rho_w gh$ (with $\rho_w$ being the water density, $g$ the acceleration of gravity, and $h$ the local interface depth), we illustrate the effect of variations in normal stress $\sigma_n$ on variations in shear strength [this method is reminiscent of the consideration of normal force in (16) but invokes different causative factors, i.e., curvature instead of plate tectonic forces].

We considered a simple two-dimensional (2D) fault model: a fault interface subjected to the pressure of the upper plate mass and an unknown horizontal tectonic stress. We also treated the convergence direction as aligned with the horizontal projection of the local along-dip direction (i.e., pure dip-slip convergence). This latter simplification is justified by the dominantly dip-slip mechanisms of most megathrust earthquakes [e.g., figures 14, 15, and 16 of (30)]. Even though subduction zones present, to varying degrees, oblique convergences, this obliquity is usually accommodated for the most part by large strike-slip faults in the back arc. One notable exception is the Solomon Islands subduction zone, which presents a very strong oblique convergence and hosts megathrust earthquakes that typically exhibit large strike-slip components. Our idealized description of the convergence direction allows us to treat the 3D problem with a 2D model (fig. S6). Within this framework, the shear strength—defined as the critical shear stress required to initiate slip at a given location—can be expressed as a function of the crustal density $\rho$, the friction $\mu$, the depth $h$, the dip angle $\theta$, and the angle $\psi$ between the maximum principal component of stress and the horizontal

$$ \tau^c = \frac{gh\mu(\rho-\rho_w)(\sin 2\theta + \tan 2\psi \cos 2\theta)}{\sin 2\theta - \mu(1-\cos 2\theta) + \tan 2\psi(\cos 2\theta - \mu \sin 2\theta)} \quad (2) $$

$\tau^c$ increases with $\theta$ ($\partial\tau^c/\partial\theta \ge 0$) (19) meaning that large dip angles imply greater shear strength. If the dominant control of earthquake magnitude were the amplitude of shear strength, we might expect a positive correlation between $M_{\max}$ and $\theta$. The negative correlation that we observe (Fig. 2A) points to other factors.

Because the magnitude of an earthquake is primarily controlled by the area of rupture, homogeneous distributions of shear strength over large fault areas may favor the occurrence of mega-earthquakes. In this framework, the shear-strength gradient $d\tau^c/ds$ is a critical parameter. Treating $\mu$, $\rho$, $\rho_w$, and $g$ as constants, spatial variations of $\tau^c$ obtained by differentiation of Eq. 2 satisfy

$$ \frac{d\tau^c}{ds} = g\mu(\rho-\rho_w) \left( A_\mu(\theta, \psi)\sin\theta + B_\mu(\theta, \psi)hK_s + C_\mu(\theta, \psi)\frac{d\psi}{ds} \right) \quad (3) $$

where $\rho > \rho_w$ and the functions $A_\mu(\theta, \psi)$, $B_\mu(\theta, \psi)$, and $C_\mu(\theta, \psi)$ are positive for $\mu = 0.6$ and values of $\theta$ and $\psi$ typically encountered in subduction zones (19). Therefore, $|d\tau^c/ds|$ increases with $K_s$ [a similar argument can be made to show that $|d\tau^c/dt|$ increases, on average, with $|K_t|$ (19)]. Shear-strength heterogeneity increases with curvature: The flatter the megathrust interface (lower $K_s$), the more homogeneous the shear-strength distribution. Stress accumulation along faults is complex, and heterogeneities in driving stress can result from the effect of stress concentrations inherited from past events or changes in coupling associated with variations in the friction coefficient or pore fluid pressure. However, dynamic stress perturbations induced by earthquake propagation are more likely to overstep the resistance to failure over broad areas if shear strength is homogeneously distributed along a large fault surface. In contrast, barriers to earthquake propagation are more likely for heterogeneous shear-strength distributions, because the dynamic stress perturbations required to reach local values of shear strength will tend to be larger in some regions. The idea that geometrical heterogeneity may limit earthquake propagation has been previously invoked (15, 31) and is consistent with the observation that many large earthquake ruptures terminate near subducting seamounts or other structural heterogeneities (31, 32).

During an earthquake, fault slip will propagate as long as the dynamic stress perturbation induced by the earlier phase of the rupture, $\delta t$, is larger than the difference between shear strength and shear stress in the surrounding area, $\Delta t$. Hence, one way to generate mega-earthquakes is to initiate a large dynamic stress perturbation by the failure of a strongly locked asperity so that $\delta t > \Delta t$ over a large surrounding fault area. In this scenario, mega-earthquakes are characterized by the rupture of one or several asperities where shear strength is particularly high. Our finding that large earthquakes rupture flat megathrusts suggests that another way to generate mega-earthquakes is if $\Delta t$ is small in a large area surrounding rupture nucleation, and thus even relatively small dynamic stress perturbations $\delta t$ can continue to propagate rupture. The 2004 Sumatra-Andaman earthquake, a unilaterally propagating 1600-km-long rupture without evidence of a localized high stress drop (33), exemplifies this second scenario, because it is unlikely that the energy released during the early phases of rupture could have initiated slip on fault portions as far as 1600 km away if the entire ruptured area was not already close to failure. In this situation, mega-earthquakes are not characterized by the rupture of asperities in regions where strength is much greater than in the surrounding area but, on the contrary, by the absence of strong variations in strength. The limiting case of homogeneous null shear strength results in constant creep and the absence of earthquakes of any magnitude. Hence, although the amplitude of $\tau^c$ is an important parameter, the observation that mega-earthquakes preferentially occur on flat subduction zones (Fig. 2) seems to indicate that the homogeneity of the

shear-strength distribution is among the most critical factors in enabling the generation of mega-earthquakes.

It has recently been proposed that any subduction zone may produce a $M \ge 9$ earthquake (34). Our results suggest that earthquake size is limited by curvature and, if our interpretation is correct, mega-earthquakes might be physically incapable of rupturing highly curved subduction zones with large shear-strength gradients, such as the Philippines, the Solomon Islands, Scotia arc, and Santa Cruz-Vanuatu-Loyalty; heterogeneous shear strength creates natural barriers that stop earthquake propagation. It is possible that heterogeneous shear stress could build to match a heterogeneous shear-strength distribution and allow rare large events, but such scenarios are expected to be less likely, and therefore more infrequent, than on megathrusts with more homogeneous shear strength. As indicated by Fig. 2B, some subduction zones, such as Izu-Bonin or Central America, may yet host earthquakes larger than previously recorded. Cascadia is noteworthy because the last mega-earthquake that it hosted is thought to have ruptured the entire subduction fault and therefore reached the maximum possible magnitude for this subduction zone (35). At a smaller scale, some areas, such as Peru, Java, or the large low-to-negative $K_s$ fault region in Central America (Fig. 1), which includes the Guerrero gap but extends much farther, show favorable features for a possible very large rupture—implying in the Guerrero case the potential for an even larger earthquake than that inferred from the extent of the gap. Moreover, large flat portions of subduction faults may sometimes rupture as one mega-earthquake and sometimes as several smaller earthquakes. Such behavior is documented, for instance, in Nankai (36), where two nearby asperities sometimes break as two separate consecutive earthquakes (as in 1854 and 1944–1946) and sometimes as one larger earthquake (1707; blue contour in the Ryukyu-Nankai box). This suggests that high-curvature regions of generally flat subduction zones may act as barriers to rupture most of the time and still sometimes be overcome by coseismic stress changes. Planar fault areas may thus host moderate-sized earthquakes for a long time and still have the potential to generate mega-earthquakes in the future.

It is possible that frictional properties and geometry are related. Such mechanical relationships are well described in accretionary wedges (37) and are certainly expected, though more difficult to constrain, along the deeper portions of megathrusts. Slab curvature has also been proposed to promote higher permeability (38–41) and thus to have an impact on the pore fluid pressure as well. Nevertheless, we have shown that spatial variations in friction and pore fluid pressure are not required to explain the gross distribution of mega-earthquakes. The observation that mega-earthquakes preferentially rupture flat megathrusts (Figs. 1 and 2B) is consistent with the inference that shear strength tends to be more homogeneously distributed along such subduction interfaces (Eq. 3), facilitating synchronized
failure over large areas. This implies that the critical feature at play in the generation of mega-earthquakes is not the amplitude of shear strength but its spatial variations. Thus, the absence of asperities on large faults may counterintuitively be a source of higher hazard. Though our study focused on subduction earthquakes, flatness may favor large earthquakes on long strike-slip faults as well.

### REFERENCES AND NOTES

1. H. Kanamori, *Philos. Trans. R. Soc. A* **364**, 1927-1945 (2006).

2. L. Ruff, H. Kanamori, *Phys. Earth Planet. Inter.* **23**, 240-252 (1980).

3. L. Ruff, H. Kanamori, *Tectonophysics* **99**, 99-117 (1983).

4. R. McCaffrey, *Geophys. Res. Lett.* **21**, 2327-2330 (1994).

5. S. Stein, E. A. Okal, *Bull. Seismol. Soc. Am.* **97**, S279-S295 (2007).

6. T. Nishikawa, S. Ide, *Nat. Geosci.* **7**, 904-908 (2014).

7. Q. Bletery et al., *J. Geophys. Res. Solid Earth* **119**, 7636-7653 (2014).

8. R. D. Muller, M. Sdrolias, C. Gaina, W. R. Roest, *Geochem. Geophys. Geosyst.* **9**, Q04006 (2008).

9. T.-R. A. Song, M. Simons, *Science* **301**, 630-633 (2003).

10. R. E. Wells, R. J. Blakely, Y. Sugiyama, D. W. Scholl, P. A. Dinterman, *J. Geophys. Res. Solid Earth* **108**, 2507 (2003).

11. W. Schellart, N. Rawlinson, *Phys. Earth Planet. Inter.* **225**, 41-67 (2013).

12. A. Heuret, C. Conrad, F. Funiciello, S. Lallemand, L. Sandri, *Geophys. Res. Lett.* **39**, L05304 (2012).

13. J. F. Pacheco, L. R. Sykes, *Bull. Seismol. Soc. Am.* **82**, 1306 (1992).

14. G. P. Hayes, D. J. Wald, R. L. Johnson, *J. Geophys. Res. Solid Earth* **117**, B01302 (2012).

15. X. Gao, K. Wang, *Science* **345**, 1038-1041 (2014).

16. C. H. Scholz, J. Campos, *J. Geophys. Res. Solid Earth* **117**, B05310 (2012).

17. R. Madariaga, *J. Geophys. Res.* **84**, 2243 (1979).

18. H. Kanamori, D. L. Anderson, *J. Geophys. Res.* **80**, 1075-1078 (1975).

19. Materials and methods are available as supplementary materials on *Science Online*.

20. P. A. McCrory, J. L. Blair, F. Waldhauser, D. H. Oppenheimer, *J. Geophys. Res. Solid Earth* **117**, n/a (2012).

21. W. P. Schellart, J. Freeman, D. R. Stegman, L. Moresi, D. May, *Nature* **446**, 308-311 (2007).

22. A. F. Holt, B. A. Buffett, T. W. Becker, *Geophys. Res. Lett.* **42**, 3802-3810 (2015).

23. W. Schellart, *J. Geophys. Res. Solid Earth* **115**, B11406 (2010).

24. T. Lay, H. Kanamori, L. Ruff, *Earthq. Predict. Res.* **1**, 3 (1982).

25. S. T. Tse, J. R. Rice, *J. Geophys. Res.* **91**, 9452 (1986).

26. J. Dieterich, *J. Geophys. Res. Solid Earth* **99**, 2601-2618 (1994).

27. C. H. Scholz, *Nature* **391**, 37-42 (1998).

28. C. Marone, *Nature* **391**, 69-72 (1998).

29. Y. Guglielmi, F. Cappa, J.-P. Avouac, P. Henry, D. Elsworth, *Science* **348**, 1224-1226 (2015).

30. B. P. Allmann, P. M. Shearer, *J. Geophys. Res. Solid Earth* **114**, B01310 (2009).

31. K. Wang, S. L. Bilek, *Tectonophysics* **610**, 1-24 (2014).

32. K. Wang, S. L. Bilek, *Geology* **39**, 819-822 (2011).

33. Q. Bletery, A. Sladen, J. Jiang, M. Simons, *J. Geophys. Res. Solid Earth* **121**, S116-S135 (2016).

34. R. McCaffrey, *Geology* **36**, 263 (2008).

35. K. Satake, K. Wang, B. F. Atwater, *J. Geophys. Res. Solid Earth* **108**, 2535 (2003).

36. T. Furumura, K. Imai, T. Maeda, *J. Geophys. Res. Solid Earth* **116**, B02308 (2011).

37. F. Dahlen, *Annu. Rev. Earth Planet. Sci.* **18**, 55-99 (1990).

38. R. J. Lisle, J. M. Robinson, *J. Struct. Geol.* **17**, 739-750 (1995).

39. C. R. Ranero, J. Phipps Morgan, K. McIntosh, C. Reichert, *Nature* **425**, 367-373 (2003).

40. M. Faccenda, T. V. Gerya, L. Burlini, *Nat. Geosci.* **2**, 790-793 (2009).

41. T. Nishikawa, S. Ide, *Geophys. Res. Lett.* **42**, 7081-7089 (2015).

### ACKNOWLEDGMENTS

We thank R. Bürgmann for his valuable comments on the manuscript. The slab1.0 model is available online at earthquake.usgs.gov/data/slab/. The U.S. Geological Survey catalog of historical earthquakes (M ≥ 8.0 since 1900) that we used in this

study can be found at http://earthquake.usgs.gov/earthquakes/world/historical.php. For the subduction zones that have not experienced any M ≥ 8.0 earthquakes since 1900, this catalog was complemented by the Global Centroid Moment Tensor catalog (www.globalcmt.org/). We used Generic Mapping Tools to compute the distributions of dip-angle gradients (gmt.soest. hawaii.edu/). This work was supported by NSF grant EAR-1520238, ANR project TO-EOS, and the French Ministry of Research and Education.

### SUPPLEMENTARY MATERIALS

www.sciencemag.org/content/354/6315/1027/suppl/DC1
Materials and Methods
Figs. S1 to S8
Tables S1 to S2
References (42-53)

6 May 2016; accepted 19 October 2016
10.1126/science.aag0482

## CATALYSIS

# Direct and continuous strain control of catalysts with tunable battery electrode materials

Haotian Wang,<fnref n="1" /> Shicheng Xu,² Charlie Tsai,<fnref n="3" />,⁴ Yuzhang Li,<fnref n="5" /> Chong Liu,<fnref n="5" /> Jie Zhao,<fnref n="5" /> Yayuan Liu,<fnref n="5" /> Hongyuan Yuan,⁶ Frank Abild-Pedersen,⁴ Fritz B. Prinz,²,<fnref n="5" /> Jens K. Nørskov,<fnref n="3" />,⁴ Yi Cui<fnref n="5" />,<fnref n="7" />*

We report a method for using battery electrode materials to directly and continuously control the lattice strain of platinum (Pt) catalyst and thus tune its catalytic activity for the oxygen reduction reaction (ORR). Whereas the common approach of using metal overlayers introduces ligand effects in addition to strain, by electrochemically switching between the charging and discharging status of battery electrodes the change in volume can be precisely controlled to induce either compressive or tensile strain on supported catalysts. Lattice compression and tension induced by the lithium cobalt oxide substrate of ~5% were directly observed in individual Pt nanoparticles with aberration-corrected transmission electron microscopy. We observed 90% enhancement or 40% suppression in Pt ORR activity under compression or tension, respectively, which is consistent with theoretical predictions.

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

Highly efficient electrocatalysts for renewable energy conversion processes, such as in H₂ fuel cells and water-splitting electro-catalysis, is becoming increasingly important (1-3). One strategy for systematically improving the activities of known catalysts is to modify their electronic structure (4-6). Numerous examples have been demonstrated in H₂O-O₂-H₂ electrocatalysis, such as the changing of d band filling in perovskite oxides for oxygen evolution (7), transition-metal alloying for the oxygen reduction reaction (ORR) (4, 8-11), and our recent studies of using lithium (Li)-ion intercalation and extraction in layered materials for water-splitting (12, 13).

Lattice strain, either compressive or tensile, can alter the surface electronic structure by mod-

ifying the distances between surface atoms and in turn catalytic activity (14-17). For platinum (Pt), previous studies have suggested that the 5d-band center of Pt can be shifted by ~0.1 eV with only 1% lattice strain (18), which can appreciably strengthen or weaken bonding of reaction intermediates to the surface (14, 18). Lattice-mismatch between metals can be generated by directly synthesizing core-shell structures (19-23) or by selectively removing atoms from an alloy (for example, stripping away Cu from a Pt-Cu alloy) (8, 14, 24-26). However, because of the larger lattice of Pt as compared with that of most metal cores, this method is typically restricted to compressive strain (14, 27). Additionally, both electronic charge transfer between the different metal atoms (ligand effects) and changes in the surface stability—and thus surface area—are present, making it difficult to identify and control the effects of strain alone (14, 25). Another strategy is to deposit catalysts onto flat substrates that undergo physical transformations as external forces are applied or the temperatures varied (28, 29). Those flat and tunable substrates present great importance to fundamental analysis, but only a few of them have been successfully demonstrated effective in electrocatalysis (28). Thus, new methods that can flexibly and effectively control both tensile and compressive lattice strain in catalysts without introducing additional effects are needed.

<footnote n="1">¹Department of Applied Physics, Stanford University, Stanford, CA 93205, USA. ²Department of Mechanical Engineering, Stanford University, Stanford, CA 93205, USA.</footnote>

<footnote n="3">³SUNCAT Center for Interface Science and Catalysis, Department of Chemical Engineering, Stanford University, Stanford, CA 94305, USA. ⁴SUNCAT Center for Interface Science and Catalysis, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA.</footnote>

<footnote n="5">⁵Department of Material Science and Engineering, Stanford University, Stanford, CA 94305, USA. ⁶Department of Physics, Stanford University, Stanford, CA 94305, USA.</footnote>

<footnote n="7">⁷Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA.</footnote>

<footnote>*Corresponding author. Email: yicui@stanford.edu</footnote>
3 Imager le seuil de rupture le long des grandes failles

Nous avons proposé dans l'étude précédente que les grandes ruptures sismiques sont plus probables le long d'une faille aux conditions physiques homogènes, rendant le comportement du « système faille » plus propice à des ruptures synchrones. Cette idée peut se traduire en des termes plus physiques. Une façon de modéliser le comportement des failles est de considérer que chaque point de la faille a une valeur de contraintes cisaillantes τ. Dans ce cadre conceptuel, chaque point de la faille est susceptible de rompre lorsque τ attend une valeur critique τ^c. Nous avons, dans l'étude qui suit, trouvé un moyen d'estimer la distribution spatiale de τ^c le long des grandes failles. Après avoir décrit le cadre théorique et dérivé l'équation donnant τ^c en fonction d'un certain nombre de paramètres physiques – incluant notamment des paramètres décrivant la géométrie de la faille mais aussi l'orientation des contraintes – nous avons « imagé » les variations spatiales de τ^c le long des failles de subduction japonaises et observé que les deux plus grands séismes historiques sur ces failles s'étaient produits sur des zones présentant une distribution spatiale de τ^c particulièrement homogène. Cette observation renforce notre intuition que plus le seuil de rupture est similaire d'un bout à l'autre d'une faille, plus il est probable qu'une grande portion de cette faille puisse rompre de manière synchrone, générant ainsi un très grand séisme.

Cette étude a été réalisée à l'Université d'Oregon en collaboration avec Amanda M. Thomas, Alan W. Rempel et Jeanne L. Hardebeck. Elle a été publiée en 2017 dans le journal Geophysical Research Letters [Bletery et al., 2017a].
<figure><img src="image_21.png" /></figure>

# Geophysical Research Letters

## RESEARCH LETTER

10.1002/2017GL075501

**Key Points:**

* We relate earthquake rupture threshold to physical parameters

* We estimate the variation of the rupture threshold along an actual subduction fault

* We relate estimates of rupture threshold variations to earthquake magnitude

**Supporting Information:**

* Supporting Information S1

**Correspondence to:**

Q. Bletery, qbletery@uoregon.edu

**Citation:**

Bletery, Q., Thomas, A. M., Rempel, A. W., & Hardebeck, J. L. (2017), Imaging shear strength along subduction faults. Geophysical Research Letters, 44. https://doi.org/10.1002/2017GL075501

Received 1 SEP 2017
Accepted 22 OCT 2017
Accepted article online 26 OCT 2017

# Imaging Shear Strength Along Subduction Faults

Quentin Bletery<fnref n="1" /> ID, Amanda M. Thomas<fnref n="1" /> ID, Alan W. Rempel<fnref n="1" /> ID, and Jeanne L. Hardebeck² ID

**Abstract** Subduction faults accumulate stress during long periods of time and release this stress suddenly, during earthquakes, when it reaches a threshold. This threshold, the shear strength, controls the occurrence and magnitude of earthquakes. We consider a 3-D model to derive an analytical expression for how the shear strength depends on the fault geometry, the convergence obliquity, frictional properties, and the stress field orientation. We then use estimates of these different parameters in Japan to infer the distribution of shear strength along a subduction fault. We show that the 2011 Mw 9.0 Tohoku earthquake ruptured a fault portion characterized by unusually small variations in static shear strength. This observation is consistent with the hypothesis that large earthquakes preferentially rupture regions with relatively homogeneous shear strength. With increasing constraints on the different parameters at play, our approach could, in the future, help identify favorable locations for large earthquakes.

**Plain Language Summary** Subduction faults accumulate stress during hundreds of years until this stress reaches a threshold. When this threshold is reached, an earthquake occurs. We derive an expression of the rupture threshold depending on different physical parameters, which vary along subduction faults. We then use estimates of the spacial variations of these different parameters in Japan to infer the distribution of the rupture threshold along the Japanese subduction faults. Our results suggest that very large earthquakes, such as the 2011 Mw 9.0 Tohoku earthquake, occur on fault portions of particularly homogeneous rupture thresholds. Our interpretation is that the rupture threshold is more likely to be reached simultaneously on broad fault portions if it does not vary much from one point of the fault to another. Imaging spatial variations of the rupture threshold along large faults may then reveal the location of possible future large earthquakes.

## 1. Introduction

Earthquakes are the result of relative motion between elastic blocks that is accommodated by slip on active faults. On convergent margins, during prolonged time periods, the contacts between crustal rocks are strong enough to prevent significant slip in the first 40–70 km beneath Earth’s surface. As a result, stress gradually accumulates along these “locked” portions of subduction faults. The accumulated stress is ultimately released when it reaches a threshold that overcomes the shear strength along a fault segment of sufficient extent for an earthquake to nucleate. Once initiated, an earthquake can propagate as long as the dynamic stress perturbation $\delta\tau$ inherited from earlier phases of the rupture is larger than the difference between the local shear strength $\tau^c$—defined as the critical shear stress required to slip at a given location—and the background shear stress $\tau$, so that

$$ \delta\tau \geq \tau^c - \tau. \qquad (1) $$

We see from equation (1) that rupture propagation past any given point is favored by large dynamic stress perturbations $\delta\tau$ and/or small values of $\tau^c - \tau$. Large $\delta\tau$ is expected from the sudden release of large stress concentrations that were able to build up in regions where the local shear strength $\tau^c$ was higher than that of surrounding fault segments. On the other hand, high $\tau^c$ fault segments also produce barriers to rupture propagation by making $\tau^c - \tau$ larger than in the surroundings so that it can potentially exceed $\delta\tau$. Unfortunately, despite this established conceptual understanding, few detailed constraints are available to quantify the terms in equation (1), particularly over the short length scales known to govern earthquake nucleation (e.g., Dieterich, 1981; Li & Rice, 1983; Rice, 1983). Nevertheless, the observation that the largest earthquakes can only occur when equation (1) is satisfied over large areas suggests that a viable strategy for identifying those regions most prone to hosting such events might be based on an approximate analysis

<footnote n="1">¹Department of Earth Science, University of Oregon, Eugene, OR, USA, ²U.S. Geological Survey, Menlo Park, CA, USA</footnote>
of the mechanical conditions over long wavelengths. Fortunately, the available data on plate geometries, convergence directions, and stress field orientations make it possible to estimate subduction shear strengths $\tau^c$ and their variations over length scales relevant to mega-earthquake propagation (≥100 km). Accordingly, here we pursue the hypothesis that insights on key features of fault behavior—in particular on the maximum earthquake magnitude—may be inferred from the distribution of shear strength.

In an earlier paper we presented a simplified 2-D formulation, assuming hydrostatic pore pressure and neglecting cohesion, of the shear strength as a function of the friction coefficient $\mu$, bulk densities of the crust $\rho$ and water $\rho_w$, depth $h$, dip angle $\theta$, and the angle $\psi$ between the maximum principal component of stress and the horizontal (Bletery et al., 2016)

$$ \tau^c = \frac{gh\mu(\rho - \rho_w)(\sin 2\theta + \tan 2\psi \cos 2\theta)}{\sin 2\theta - \mu(1 - \cos 2\theta) + \tan 2\psi(\cos 2\theta - \mu \sin 2\theta)}. \quad (2) $$

Here we derive a similar expression for the shear strength on any point of a subduction fault in a more general 3-D case that takes into account the obliquity of plate convergence. Our expression for $\tau^c$ (equation (11)) depends on geometrical ($h, \theta$) and constitutive ($\mu$) parameters, in addition to other factors needed to describe the local state. After outlining our model formulation, we illustrate its implications by evaluating $\tau^c$ and its spatial variations along the Japanese subduction faults for which the controlling parameters are best constrained. We ultimately show that the 2011 $M_{w}9.0$ Tohoku-Oki earthquake ruptured an area of particularly low shear strength gradient, supporting the idea that shear strength is more likely to be exceeded simultaneously over large areas if its spatial distribution is homogeneous (Bletery et al., 2016).

## 2. Model

We consider the stress field that acts upon a fault increment on a subduction megathrust. In the $(x, y, z)$ coordinate system, where $x$ coincides with the convergence direction, $y$ is perpendicular to $x$ in the horizontal plane and $z$ points up (Figure 1), the stress tensor $\sigma$ can be written as

$$ \sigma = \begin{pmatrix} \sigma_{zz} + \Delta\sigma & \sigma_{xy} & \sigma_{xz} \\ \sigma_{xy} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_{zz} \end{pmatrix}, \quad (3) $$

where $\Delta\sigma = \sigma_{xx} - \sigma_{zz}$ is a stress difference. Let $\psi$ be the angle between the first principal component axis and the horizontal. Assuming, based on symmetry considerations, that $y$ is a principal component, the stress tensor can be rewritten as (Appendix A)

$$ \sigma = \begin{pmatrix} \sigma_{zz} + \Delta\sigma & 0 & \frac{\Delta\sigma}{2} \tan 2\psi \\ 0 & \sigma_{yy} & 0 \\ \frac{\Delta\sigma}{2} \tan 2\psi & 0 & \sigma_{zz} \end{pmatrix}. \quad (4) $$

Let $n$ be a unit vector normal to the megathrust and $s$ a unit vector in the direction of the plate convergence projected onto the megathrust (Figure 1). In the $(x, y, z)$ reference frame, $n = (-\sin \delta, 0, -\cos \delta)^T$ and $s = (\cos \delta, 0, -\sin \delta)^T$, with $\delta$ the apparent dip angle in the direction of the convergence (Figure 1). Based on simple geometrical considerations, $\delta$ may be expressed as a function of the strike angle $\phi$, the azimuth of the convergence direction $\lambda$, and the actual dip angle $\theta$ (Appendix B)

$$ \delta = \arctan(\sin(\phi - \lambda) \tan \theta). \quad (5) $$

The normal and shear stress $\sigma_n$ and $\tau$ are then

$$ \sigma_n = (\sigma \cdot n) \cdot n = \sigma_{zz} + \frac{\Delta\sigma}{2}(1 - \cos 2\delta + \tan 2\psi \sin 2\delta), \quad (6) $$

$$ \tau = \sqrt{(\sigma \cdot n) \cdot (\sigma \cdot n) - \sigma_n^2} = \frac{\Delta\sigma}{2}(\sin 2\delta + \tan 2\psi \cos 2\delta). \quad (7) $$

The Coulomb failure criterion gives the relationship between shear and normal stress at which the rupture initiates

$$ \tau = \mu(\sigma_n - p), \quad (8) $$
<figure><img src="image_22.png" /><figcaption>**Figure 1.** Sketch of the problem configuration. An incremental fault surface (of depth h, dip angle θ, and strike angle φ) is oriented with an apparent angle δ from the horizontal in the convergence direction of azimuth λ. The first principal component of the stress field σ₁ applied to this incremental fault surface forms an angle ψ with the horizontal. Inset: relationship between the apparent dip angle δ, the actual dip angle θ, the strike angle φ, and the horizontal convergence direction λ.</figcaption></figure>

where μ is the coefficient of friction and p is the pore fluid pressure, and we neglect cohesion for simplicity. Equations (6)–(8) enable us to express the critical stress difference required to initiate a rupture Δσ<sup>c</sup> as a function of μ, σ<sub>zz</sub>, p, δ, and ψ:

$$ \Delta\sigma^c = \frac{2\mu(\sigma_{zz} - p)}{\sin 2\delta - \mu(1 - \cos 2\delta) + \tan 2\psi(\cos 2\delta - \mu \sin 2\delta)}. \quad (9) $$

We can then derive from equation (7) (or equivalently equations (6) + (8)) an expression for the critical shear stress τ<sup>c</sup> at which a rupture initiates

$$ \tau^c = \frac{\mu(\sigma_{zz} - p) \sin(2\delta + 2\psi)}{\sin(2\delta + 2\psi) + \mu(\cos(2\delta + 2\psi) - \cos 2\psi)}. \quad (10) $$

### 3. Numerical Computation

Equation (10) is quite general. The only assumptions are that the fault strength is governed by the Mohr-Coulomb failure criterion (equation (8)) and that y is a principal component axis. In order to compute the distribution of shear strength, however, we need to make some stronger assumptions. A natural approximation for σ<sub>zz</sub> is to treat it as equal to the pressure due to the mass of the upper plate; that is, σ<sub>zz</sub> = ρgh, with ρ the bulk density of the crust, g the acceleration of gravity, and h the depth. This case is valid as long as the gravitational loading dominates the vertical state of stress, implying that any contribution of tectonic origin to the vertical stress is negligible compared to the lithostatic pressure. A hydrostatic fluid pressure distribution could produce p = ρ<sub>w</sub>gh, with ρ<sub>w</sub> the bulk density of water. Although commonly assumed and useful as a starting point for exploring the controls on strength variations, the distribution of fluid pressure is relatively poorly constrained and potential deviations of p from ρ<sub>w</sub>gh will be discussed the later sections. These two assumptions lead to

$$ \tau^c = \frac{\mu(\rho - \rho_w)gh \sin(2\delta + 2\psi)}{\sin(2\delta + 2\psi) + \mu(\cos(2\delta + 2\psi) - \cos 2\psi)}, \quad (11) $$

where δ = arctan(sin(φ - λ) tan θ) (equation (5)).
We compute the distribution of shear strength using equations (11) and (5) and the Generic Mapping Tools (GMT) software. We use the distributions of dip (θ) and strike (φ) angles along the Japanese megathrusts from the SLAB1.0 model (Hayes et al., 2012), the distribution of the convergence direction λ given by the NNR-MORVEL56 model (Argus et al., 2011), and the distribution of the stress field orientation ψ from Hardebeck (2015). To estimate h, we use the depth of the interface given by the SLAB1.0 model (Hayes et al., 2012) corrected from the 1 arc min bathymetry given by the ETOPO1 model. Other variables are treated as constant: g=9.8 m/s², ρ=2700 kg/m³, ρw=1000 kg/m³, and μ=0.6.

Though we recognize the contribution of the preconditioned dynamic stress perturbation (δτ) and the background shear stress (τ) highlighted in equation (1), we proposed that the spatial variations of shear strength might be a key parameter in controlling the global distribution of mega-earthquakes because shear strength is more likely to be exceeded simultaneously over large areas if its spatial distribution is homogeneous (Bletery et al., 2016). To support this hypothesis, we further calculate the gradient of shear strength. We use the GMT software to differentiate our map of shear strength along east (xE) and along north (xN) and project the shear strength gradient in the rake-parallel (dτc/ds) and rake-perpendicular (dτc/dt) directions (with t (= y) the normal to s in the fault plane (Figure 1)),

$$ \frac{d\tau^c}{ds} = \left( \frac{d\tau^c}{dx_E} \sin \lambda + \frac{d\tau^c}{dx_N} \cos \lambda \right) \cos \delta, \quad (12) $$

$$ \frac{d\tau^c}{dt} = \frac{d\tau^c}{dx_E} \cos \lambda - \frac{d\tau^c}{dx_N} \sin \lambda. \quad (13) $$

## 4. Results

Assuming σzz − p = (ρ − ρw)gh leads to a linear dependence of the shear strength on depth (equation (11)). This depth dependence is modified by variations in the stress field orientation with depth. Numerical computation of equation (11) using the stress field orientation of Hardebeck (2015) is shown in Figure 2. The pattern of variations in τc that we determine is qualitatively consistent with the dependence of stress drop on depth observed among small earthquakes in the Tohoku rupture area (Oth, 2013; Uchide et al., 2014): higher stress drops are favored at depth where the shear stress rupture threshold is larger. It is worth noting, however, that the observed stress drops are a factor of 15 smaller than the calculated shear strengths (inset of Figure 2). A significant part of this difference could be attributed to inaccuracies in assumed pore pressures and other parameter values. For instance, assuming a coefficient of friction μ of 0.1—such as reported by Fulton et al. (2013) in the shallow part of the Tohoku megathrust and in agreement with heat flow measurements (Gao & Wang, 2014)—instead of 0.6 reduces the difference to a factor of 2.5 (Figure S1 in the supporting information). Nevertheless, these offsets are consistent with the commonly accepted view that shear stress does not drop to zero following an earthquake; stress drops and shear strength clearly differ. In what follows, we use μ = 0.6, but we show all equivalent figures for μ = 0.1 in the supporting information.

The distributions of the rake-parallel dτc/ds and rake-perpendicular |dτc/dt| components of the shear strength gradient along the Japanese megathrusts are shown in Figures 3a and 3b. dτc/ds has a larger amplitude than dτc/dt and is positive as long as τc increases monotonically with depth (equation (11)). The norm of the shear strength gradient is dominated by the rake-parallel component, which has a higher amplitude (Figure 3c). The rake-parallel component dτc/ds indicates particularly low gradient offshore Nankai (red patch in the south of Japan in Figure 3a). On the other hand, the rake-perpendicular component |dτc/dt| indicates areas where the convergence direction λ relative to the fault orientation φ generates heterogeneities in the shear strength distribution. Such heterogeneities seem to be minimal in the Tohoku rupture area (red patch east of the main island) and south of Nankai (Figure 3b). Even though |dτc/dt| is of much smaller amplitude than dτc/ds, the distributions of the two components may be of similar importance as one (dτc/ds) influences along-dip rupture propagation, while the other (|dτc/dt|) influences the lateral propagation. The amplitude of the rake-parallel and rake-perpendicular components of the shear strength gradient, as well as the norm, is reduced by a factor of ~6 when computed with μ = 0.1 instead of μ = 0.6, but the relative distributions are rather insensitive to the choice of μ (Figure S2).
<figure><img src="image_23.png" /><figcaption>**Figure 2.** Shear strength $\tau^c$ along the Japanese megathrusts computed from equation (11). Inset: Median stress drops (colored squares) of small earthquakes in the Tohoku-Oki area (Uchide et al., 2014) superimposed on top of the shear strength distribution.</figcaption></figure>

The gradient maps are sensitive to the resolution of the different parameters in equation (11). The smaller wavelength features in Figure 3 reflect changes in the bathymetry, which is by far the best constrained parameter, with a resolution of 1 arc min. The slab1.0 model (constraining the dip and strike angles as well as the fault depth relative to sea level) has a nominal resolution of 0.5° but relies on uneven sampling in multiple data sets interpolated using Hermite splines of polynomial orders from 2 to 3, depending on depth (Hayes et al., 2012). The stress orientation is also interpolated from 0.5° wide pixels (Hardebeck, 2015). These uneven resolution limitations make our shear strength gradient maps overly sensitive to the spatial variations of the most well-resolved parameters. Though of comparable resolution, geometry dominates stress orientation effects because it varies more dramatically at the resolved scale (Figures 3a and 3c). Nevertheless, stress orientation and fault geometry may vary in different manners at much smaller scales and produce barriers at scales relevant to rupture nucleation that our long-wavelength imaging cannot capture. The rake-perpendicular component of the shear strength gradient is most sensitive to the obliquity term $\phi - \lambda$ because
<figure><img src="image_24.png" /><figcaption>**Figure 3.** (a) Rake-parallel $\frac{dr_c}{ds}$ and (b) absolute value of the rake-perpendicular $|\frac{dr_c}{dt}|$ components of the shear strength gradient along the Japanese megathrusts. White arrows show the directions of plate convergences relative to the Eurasian plate. (c) Norm of the shear strength gradient. The black contour shows the location of the largest historical rupture along the Japanese megathrusts: the 2011 M<sub>w</sub>9.0 Tohoku-Oki earthquake (Bletery et al., 2014). White contours show M ~ 8.5 historical earthquakes (Rhea et al., 2010; Simons et al., 2011; Furumura et al., 2011). Black dashed contours show recurrent M8.1–8.4 earthquakes that ruptured simultaneously during the 1707 M8.4–8.8 Hœi earthquake (Furumura et al., 2011; Stein & Okal, 2007). Note that slip contours of historical earthquakes are much less well constrained than for the Tohoku-Oki earthquake.</figcaption></figure>
<figure><img src="image_25.png" /></figure>

**Figure 4.** Average (black diamonds) and standard deviation (black vertical intervals) of the (a) rake parallel, (b) rake perpendicular, and (c) norm of shear strength gradient inside 17 slip contours corresponding to the (d) rupture area of the 2011 Tohoku earthquake ($M_w$ 9.0) distributed along the megathrusts. The red contour (n°7) corresponds to the actual slip contour of the Tohoku earthquake (Bletery et al., 2014). The other 16 black contours have the same shape and area than the Tohoku earthquake but are shifted along the megathrusts. The horizontal blue solid and dashed lines show the averages and standard deviations along the seismogenic zone (≤60 km) of the entire megathrusts.

the obliquity of the convergence makes the depth change along rake-perpendicular paths, inducing changes
in strength along such paths (Figure 3b).

We proposed in an earlier paper that mega-earthquakes are favored by low shear strength gradients (Bletery et al., 2016) because homogeneous shear strength distributions facilitate synchronized failure over large fault portions. Though limited by the different parameter resolutions, Figure 3 provides material to further investigate this hypothesis. We therefore consider the slip contour of the 2011 $M_w$ 9.0 Tohoku-Oki earthquake (Bletery et al., 2014), which is by far the largest rupture ever recorded in the region, and calculate the average shear strength gradient inside the rupture area (for both components and the norm). We then shift the Tohoku slip contour all along the Japanese megathrusts and compare the average shear strength gradient of Tohoku-sized areas. Our results indicate that the actual location of the Tohoku earthquake presents particularly homogeneous shear strength distributions (i.e., low shear strength gradient in all components) (Figure 4). As expected from Figure S2 relative to Figure 3, these results are insensitive to the choice of $\mu$ (Figure S3).

Though it seems difficult to find a Tohoku-sized area of more homogeneous shear strength distribution than
the actual area in which it took place, the smaller (still very large) 1963 $M_w$ 8.5 Kuril earthquake seems to
have ruptured a more heterogeneous fault portion (Figure 3c). Given the nature of our data sets, we are likely
insensitive to geometrical heterogeneities smaller than several tens of kilometers and frictional heterogeneities of any size, meaning that our shear strength gradient estimates essentially constitute lower bound constraints on the spatial variability of the static shear strength. In the case of the Kuril earthquake, this means that an average gradient in static shear strength of at least 4.4 MPa/km did not effectively halt a 245 km long rupture. In detail, this earthquake was actually composed of three distinct 40–60 km long subevents that ruptured successively (Beck & Ruff, 1987), consistent with a likely heterogeneous $\tau^c$ distribution at a scale we are not able to resolve and resulting in relatively high integrated gradients at the scale to which we are sensitive (Figure 3).

The other M ~ 8.5 historical earthquakes ruptured more favorable fault areas. The 1707 Hoei M8.4–8.8 earthquake ruptured the previously discussed red patch offshore Nankai in Figure 3a. Interestingly, this earthquake is thought to correspond to the simultaneous rupture of two patches that ruptured distinctively in 1944 and 1946 and before in 1854 (Furumura et al., 2011). The two patches seem to coincide with low shear strength gradient separated by a geometrical heterogeneity (Figure 3c).

## 5. Discussion and Conclusions

The rupture threshold $\tau^c$ depends upon geometry ($h$, $\theta$, $\phi$), rheology ($\mu$), hydrology ($p$), the orientation of the convergence ($\lambda$), and the orientation of the stress field ($\psi$), which is inherited from all the other parameters and past events (equations (10) and (5)). When the pore pressure and vertical normal stress both increase linearly with depth, the increase of $\tau^c$ with $h$ (equation (11) and Figure 2) suggests that larger values of $\delta\tau$ are to be expected when an earthquake initiates at depth. Remembering that a rupture propagates as long as $\delta\tau \ge \tau^c - \tau$ (equation (1)), it also suggests that it is much easier for a rupture to propagate in the updip direction than in the downdip direction. Indeed, $\tau^c$—and therefore likely ($\tau^c - \tau$)—decreases in the updip direction, meaning that the dynamic stress perturbation required to keep the propagation going is progressively reduced along an upwardly oriented rupture path. In special circumstances, this tendency may be counterbalanced by variations of the stress field orientation (equation (11)). However, the estimates of the stress field orientation (Hardebeck, 2015) included in our calculation of $\tau^c$ for Japanese subduction zones do not compensate for the expected linear increase of $\sigma_{zz} - p$ with depth (Figure 2).

The expected increase in $\tau^c$ with depth may also be counterbalanced by nonhydrostatic distributions of pore fluid pressure. It is likely that in some areas, the distribution of pore fluid pressures is close to lithostatic (e.g., Byerlee, 1990; Faulkner & Rutter, 2001; Rice, 1992; Streit & Cox, 2001), taking the form $p = \rho gh + C$ and breaking with the linear depth dependency of $\tau^c$ predicted by equation (11). Near-lithostatic pore fluid pressure appears, for instance, to be required to explain the characteristics of slow earthquakes (e.g., Audet & Bürgmann, 2014; Liu & Rice, 2007; Segall et al., 2010) and other slow-slip phenomena (e.g., Bletery et al., 2017; Shelly et al., 2006), in particular their sensitivity to extremely small tidal stress perturbations (Thomas et al., 2009). Such phenomena are well documented in Nankai both at depth and in the wedge (e.g., Nakajima & Hasegawa, 2016; Obara & Kato, 2016; Shelly et al., 2006), implying that the actual shear strength in the Nankai section (contours n°12 and n°13 in Figure 4) is probably locally much smaller than suggested in Figure 4. More generally, pore pressure is expected to follow complex variations that evolve with fluid flow along the megathrust and depend upon such poorly constrained variables as the permeability of the overriding plate (e.g., Nakajima & Hasegawa, 2016; Skarbek & Rempel, 2016). In the absence of firm constraints to the contrary, the assumption of hydrostatic pore pressures should be viewed as a strong limitation of our analysis.

The present numerical computation is exploring the influence of nonfrictional parameters—nominally seafloor bathymetry, fault geometry, stress orientation, and convergence direction—on the spatial distribution of the static shear strength along the Japanese megathrusts. Though we treat $\mu$ as constant, our approach does not require the assumption of constant friction and possible future constraints on the spatial distribution of $\mu$ could be included in the numerical computation. $\mu$ is also known to vary temporally as a function of the slip velocity and one or more state variables (e.g., review paper by Marone, 1998, and references therein). Although temporal variations in the friction coefficient due to rate-and-state effects are critical for our modern understanding of earthquake nucleation among other phenomena, the absolute magnitude of such changes is typically quite small (≤10%) and would be expected to have a correspondingly modest influence on our calculations of $\tau^c$. Dynamic weakening mechanisms, such as thermal pressurization (Lachenbruch, 1980), flash heating (Rice, 2006), or other material-dependent effects (e.g., Di Toro et al., 2011; Han et al., 2007), might
influence the evolution of shear strength more dramatically during the rupture process (e.g., Viesca & Garagash, 2015), possibly facilitating the propagation of rare large ruptures into domains that were strong enough (or rate strengthening) to withstand smaller seismic disturbances. Though beyond the scope of this study, our expression of $\tau^c$ could be coupled with time-dependent equations governing these different mechanisms in numerical simulations of the seismic cycle and/or earthquake rupture to introduce effects of geometry, stress orientation, and convergence obliquity in these approaches. Spatial variations of the upper plate density ($\rho$) could also, in future work, be included in the calculation of $\tau^c$ (equation (11)). Bassett et al. (2016) observed a correlation between contrasts in residual density anomaly and the rupture termination of the Tohoku-Oki earthquake. They proposed that the residual density anomaly reflects changes in frictional properties. As the dependency of shear strength on density is linear (equation (11)), Bassett et al.'s (2016) observation could be reinterpreted in terms of gradients in shear strength stopping rupture propagation.

In any case, the larger shear strength at depth predicted by equation (11) is counterbalanced by larger shear strength gradient (Figures 3a and 3c). This implies that earthquakes at depth are energetic (large $\delta\tau$) but spatially confined because their propagation is limited by the heterogeneity of the strength distribution. Such a mechanism would explain the high frequencies—thought to be radiated by small on-fault asperities (Madariaga, 1983)—detected at the downdip edge of the rupture of the Tohoku (Ide et al., 2011; Koper et al., 2011; Meng et al., 2011; Simons et al., 2011) and other large subduction earthquakes (Koper et al., 2012; Lay et al., 2012; Yao et al., 2013). Small-scale geometrical irregularities are themselves expected to influence rupture initiation and propagation (e.g., Dunham et al., 2011; Fang & Dunham, 2013) to the extent that such fault “bumps” are superimposed without a systematic dependence upon the large-scale geometry that we examine here. The variations in $\tau^c$ that determine which fault segments are most prone to large earthquakes are expected to be well represented by equation (10).

The shear strength formulation we present here expresses the rupture threshold as a function of parameters for which we have some estimates along actual faults. Present estimates of stress field orientation are still rather crude, but the development of seismic instrumentation should allow for more refined studies in the future. More generally, the small-scale variations on most parameters in equations (11) and (5) (e.g., $h, \theta, \phi$, and $\psi$) are not reflected in the broad-scale estimates we used and could generate small-scale heterogeneities relevant to the nucleation and propagation of small to moderate earthquakes. Though local variations of the background shear stress are also relevant to rupture propagation (equation (11)), we showed that insights on the maximum earthquake magnitude may be deduced from the shear strength distribution alone. Though this study focused on the Japanese megathrusts, the approach may be applied to any subduction zone and the idea may be transposed to other fault systems. In particular, in the case of strike-slip faults, several studies have shown that strong variations in azimuth ($\phi$) and lateral discontinuities coincide with brakes or stops in earthquake propagations (Biasi & Wesnousky, 2016; Lettis et al., 2002; Perrin et al., 2016; Wesnousky, 2006). Such variations in $\phi$ could be interpreted as the source of large gradients in shear strength producing barriers to ruptures. Finally, we note that improved constraints on the depth dependence of pore fluid pressure and variations in the coefficient of friction may, in future work, be introduced in equation (10) to explore their influences on the distribution of shear strength along large faults.

## Appendix A: Determination of the Nondiagonal Components of the Stress Tensor

If $y$ is aligned with one of the principal components of the stress field, then the stress tensor can be expressed in the principal component reference frame as

$$ \begin{pmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_3 \end{pmatrix} = \begin{pmatrix} \cos \psi & 0 & \sin \psi \\ 0 & 1 & 0 \\ -\sin \psi & 0 & \cos \psi \end{pmatrix} \cdot \begin{pmatrix} \sigma_{zz} + \Delta\sigma & \sigma_{xy} & \sigma_{xz} \\ \sigma_{xy} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_{zz} \end{pmatrix} \cdot \begin{pmatrix} \cos \psi & 0 & -\sin \psi \\ 0 & 1 & 0 \\ \sin \psi & 0 & \cos \psi \end{pmatrix}. \quad (A1) $$

The off-diagonal terms give a system of three equations that depend on $\sigma_{xy}$, $\sigma_{xz}$, and $\sigma_{yz}$

$$ 0 = \sigma_{xy} \cos \psi - \sigma_{yz} \sin \psi, \quad (A2) $$

$$ 0 = -(\sigma_{zz} + \Delta\sigma) \sin \psi \cos \psi + \sigma_{xz} \cos^2 \psi - \sigma_{xz} \sin^2 \psi + \sigma_{zz} \cos \psi \sin \psi, \quad (A3) $$

$$ 0 = \sigma_{xy} \sin \psi - \sigma_{yz} \cos \psi. \quad (A4) $$
Solving this system gives

$$\sigma_{xy} = 0, \qquad (A5)$$

$$\sigma_{xz} = \frac{\Delta\sigma}{2} \tan 2\psi, \qquad (A6)$$

$$\sigma_{yz} = 0. \qquad (A7)$$

The stress tensor can then be rewritten in the (x, y, z) reference frame as equation (4).

## Appendix B: Apparent Dip Angle

Considering the inset of Figure 1, we have

$$a = c \sin(\phi - \lambda), \qquad (B1)$$

$$b = a \tan \theta, \qquad (B2)$$

$$b = c \tan \delta. \qquad (B3)$$

Combining (B2) and (B3) gives

$$c \tan \delta = a \tan \theta. \qquad (B4)$$

Using (B1), we obtain

$$\tan \delta = \sin(\phi - \lambda) \tan \theta, \qquad (B5)$$

from which we can derive equation (5).

## Acknowledgments

This work was supported by the National Science Foundation (NSF) project EAR-1520238. We thank Takahiko Uchide for providing his stress drop estimates shown in Figures 2 and S1 and Joan Gomberg and Nicholas M. Beeler for thoughtful comments on the manuscript. The SLAB1.0 and ETOPO1 models are available online at the following addresses earthquake.usgs.gov/data/slab/ and ngdc.noaa.gov/mgg/global/, respectively. The convergence direction at the different points of the faults was calculated by the online motion calculator at this address: geoscience.wisc.edu/~chuck/MORVEL/. We used the Generic Mapping Tools to compute the distributions of shear strength and their gradients (gmt.soest.hawaii.edu/).

## References

Argus, D. F., Gordon, R. G., & DeMets, C. (2011). Geologically current motion of 56 plates relative to the no-net-rotation reference frame. *Geochemistry, Geophysics, Geosystems*, 12, Q11001. https://doi.org/10.1029/2011GC003751

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

Bassett, D., Sandwell, D. T., Fialko, Y., & Watts, A. B. (2016). Upper-plate controls on co-seismic slip in the 2011 magnitude 9.0 Tohoku-Oki earthquake. *Nature*, 531(7592), 92–96.

Beck, S. L., & Ruff, L. J. (1987). Rupture process of the Great 1963 Kurile Islands Earthquake Sequence: Asperity interaction and multiple event rupture. *Journal of Geophysical Research*, 92(B13), 14,123–14,138.

Biasi, G. P., & Wesnousky, S. G. (2016). Steps and gaps in ground ruptures: Empirical bounds on rupture propagation. *Bulletin of the Seismological Society of America*, 106, 1110–1124.

Bletery, Q., Sladen, A., Delouis, B., Vallée, M., Nocquet, J.-M., Rolland, L., & Jiang, J. (2014). A detailed source model for the $M_{w}9.0$ Tohoku-Oki earthquake reconciling geodesy, seismology, and tsunami records. *Journal of Geophysical Research: Solid Earth*, 119, 7636–7653. https://doi.org/10.1002/2014JB011261

Bletery, Q., Thomas, A. M., Rempel, A. W., Karlstrom, L., Sladen, A., & De Barros, L. (2016). Mega-earthquakes rupture flat megathrusts. *Science*, 354(6315), 1027–1031.

Bletery, Q., Thomas, A. M., Hawthorne, J. C., Skarbek, R. M., Rempel, A. W., & Krogstad, R. D. (2017). Characteristics of secondary slip fronts associated with slow earthquakes in Cascadia. *Earth and Planetary Science Letters*, 463, 212–220.

Byerlee, J. (1990). Friction, overpressure and fault normal compression. *Geophysical Research Letters*, 17(12), 2109–2112.

Di Toro, G., Han, R., Hirose, T., De Paola, N., Nielsen, S., Mizoguchi, K., ... Shimamoto, T. (2011). Fault lubrication during earthquakes. *Nature*, 471(7339), 494–8.

Dieterich, J. H. (1981). Constitutive properties of faults with simulated gouge. In N. L. Carter et al. (Eds.), *Mechanical behavior of crustal rocks: The Handin volume* (Vol. 24, pp. 103–120). Washington, DC: American Geophysical Union.

Dunham, E. M., Belanger, D., Cong, L., & Kozdon, J. E. (2011). Earthquake ruptures with strongly rate-weakening friction and off-fault plasticity. Part 2: Nonplanar faults. *Bulletin of the Seismological Society of America*, 101(5), 2308–2322.

Fang, Z., & Dunham, E. M. (2013). Additional shear resistance from fault roughness and stress levels on geometrically complex faults. *Journal of Geophysical Research: Solid Earth*, 118, 3642–3654. https://doi.org/10.1002/jgrb.50262

Faulkner, D. R., & Rutter, E. H. (2001). Can the maintenance of overpressured fluids in large strike-slip fault zones explain their apparent weakness? *Geology*, 29(6), 503–506.

Fulton, P., Brodsky, E. E., Kano, Y., Mori, J., Chester, F, Ishikawa, T., ... Toczko, S. (2013). Low coseismic friction on the Tohoku-Oki fault determined from temperature measurements. *Science*, 342(6163), 1214–1217.

Furumura, T., Imai, K., & Maeda, T. (2011). A revised tsunami source model for the 1707 Hoei earthquake and simulation of tsunami inundation of Ryujin Lake, Kyushu, Japan. *Journal of Geophysical Research*, 116, B02308. https://doi.org/10.1029/2010JB007918
<figure><img src="image_26.png" /></figure>

Gao, X., & Wang, K. (2014). Strength of stick-slip and creeping subduction megathrusts from heat flow observations. *Science*, **345**(6200), 1038–1041.

Han, R., Shimamoto, T., Hirose, T., Ree, J.-H., & Ando, J.-i. (2007). Ultralow friction of carbonate faults caused by thermal decomposition. *Science*, **316**(5826), 878–881.

Hardebeck, J. L. (2015). Stress orientations in subduction zones and the strength of subduction megathrust faults. *Science*, **349**(6253), 1213–1216.

Hayes, G. P., Wald, D. J., & Johnson, R. L. (2012). Slab1.0: A three-dimensional model of global subduction zone geometries. *Journal of Geophysical Research: Solid Earth*, **117**, B01302. https://doi.org/10.1029/2011JB008524

Ide, S., Baltay, A., & Beroza, G. C. (2011). Shallow dynamic overshoot and energetic deep rupture in the 2011 *M*<sub>w</sub> 9.0 Tohoku-Oki earthquake. *Science*, **332**(6036), 1426–1429.

Koper, K. D., Hutko, A. R., Lay, T., Ammon, C. J., & Kanamori, H. (2011). Frequency-dependent rupture process of the 2011 *M*<sub>w</sub> 9.0 Tohoku earthquake: Comparison of short-period p wave backprojection images and broadband seismic rupture models. *Earth, Planets and Space*, **63**(7), 16.

Koper, K. D., Hutko, A. R., Lay, T., & Sufri, O. (2012). Imaging short-period seismic radiation from the 27 February 2010 Chile (*M*<sub>w</sub> 8.8) earthquake by back-projection of *P*, *PP*, and *PKIKP* waves. *Journal of Geophysical Research: Solid Earth*, **117**, B02308. https://doi.org/10.1029/2011JB008576

Lachenbruch, A. H. (1980). Frictional heating, fluid pressure, and the resistance to fault motion. *Journal of Geophysical Research*, **85**(B11), 6097–6112.

Lay, T., Kanamori, H., Ammon, C. J., Koper, K. D., Hutko, A. R., Ye, L., ... Rushing, T. M. (2012). Depth-varying rupture properties of subduction zone megathrust faults. *Journal of Geophysical Research*, **117**, B04311. https://doi.org/10.1029/2011JB009133

Lettis, W., Bachhuber, J., Witter, R., Brankman, C., Randolph, C., Barka, A., ... Kaya, A. (2002). Influence of releasing step-overs on surface fault rupture and fault segmentation: Examples from the 17 August 1999 Izmit earthquake on the North Anatolian Fault, Turkey. *Bulletin of the Seismological Society of America*, **92**(1), 19–42.

Li, V. C., & Rice, J. R. (1983). Preseismic rupture progression and great earthquake instabilities at plate boundaries. *Journal of Geophysical Research*, **88**(B5), 4231–4246.

Liu, Y., & Rice, J. R. (2007). Spontaneous and triggered aseismic deformation transients in a subduction fault model. *Journal of Geophysical Research*, **112**, B09404. https://doi.org/10.1029/2007JB004930

Madariaga, R. (1983). High frequency radiation from dynamic earthquake. *Annales Geophysicae*, **1**, 17–23.

Marone, C. (1998). Laboratory-derived friction laws and their application to seismic faulting. *Annual Review of Earth and Planetary Sciences*, **26**(1), 643–696.

Meng, L., Inbal, A., & Ampuero, J.-P. (2011). A window into the complexity of the dynamic rupture of the 2011 *M*<sub>w</sub> 9 Tohoku-Oki earthquake. *Geophysical Research Letters*, **38**, L00G07. https://doi.org/10.1029/2011GL048118

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

Obara, K., & Kato, A. (2016). Connecting slow earthquakes to huge earthquakes. *Science*, **353**(6296), 253–257.

Oth, A. (2013). On the characteristics of earthquake stress release variations in Japan. *Earth and Planetary Science Letters*, **377**, 132–141.

Perrin, C., Manighetti, I., Ampuero, J.-P., Cappa, F., & Gaudemer, Y. (2016). Location of largest earthquake slip and fast rupture controlled by along-strike change in fault structural maturity due to fault growth. *Journal of Geophysical Research: Solid Earth*, **121**, 3666–3685. https://doi.org/10.1002/2015JB012671

Rhea, S., Tarr, A. C., Hayes, G. P., Villaseñor, A., Furlong, K. P., & Benz, H. (2010). Seismicity of the Earth 1900–2007, Kuril-Kamchatka arc and vicinity: US Geological Survey.

Rice, J. R. (1983). Constitutive relations for fault slip and earthquake instabilities. *Pure and Applied Geophysics*, **121**(3), 443–475.

Rice, J. R. (1992). Fault stress states, pore pressure distributions, and the weakness of the San Andreas Fault. *International Geophysics*, **51**, 475–503.

Rice, J. R. (2006). Heating and weakening of faults during earthquake slip. *Journal of Geophysical Research*, **111**, B05311. https://doi.org/10.1029/2005JB004006

Segall, P., Rubin, A. M., Bradley, A. M., & Rice, J. R. (2010). Dilatant strengthening as a mechanism for slow slip events. *Journal of Geophysical Research*, **115**, B12305. https://doi.org/10.1029/2010JB007449

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

Simons, M., Minson, S. E., Sladen, A., Ortega, F., Jiang, J., ... Webb, F. H. (2011). The 2011 magnitude 9.0 Tohoku-Oki earthquake: Mosaicking the megathrust from seconds to centuries. *Science*, **332**(6036), 1421–1425.

Skarbek, R. M., & Rempel, A. W. (2016). Dehydration-induced porosity waves and episodic tremor and slip. *Geochemistry, Geophysics, Geosystems*, **17**, 442–469. https://doi.org/10.1002/2015GC006155

Stein, S., & Okal, E. A. (2007). Ultralong period seismic study of the December 2004 Indian Ocean earthquake and implications for regional tectonics and the subduction process. *Bulletin of the Seismological Society of America*, **97**(1A), S279–S295.

Streit, J. E., & Cox, S. F. (2001). Fluid pressures at hypocenters of moderate to large earthquakes. *Journal of Geophysical Research*, **106**(B2), 2235–2243.

Thomas, A. M., Nadeau, R. M., & Bürgmann, R. (2009). Tremor-tide correlations and near-lithostatic pore pressure on the deep San Andreas Fault. *Nature*, **462**(7276), 1048–1051.

Uchide, T., Shearer, P. M., & Imanishi, K. (2014). Stress drop variations among small earthquakes before the 2011 Tohoku-Oki, Japan, earthquake and implications for the main shock. *Journal of Geophysical Research: Solid Earth*, **119**, 7164–7174. https://doi.org/10.1002/2014JB010943

Viesca, R. C., & Garagash, D. I. (2015). Ubiquitous weakening of faults due to thermal pressurization. *Nature Geoscience*, **8**(11), 875–879.

Wesnousky, S. G. (2006). Predicting the endpoints of earthquake ruptures. *Nature*, **444**(7117), 358–360.

Yao, H., Shearer, P. M., & Gerstoft, P. (2013). Compressive sensing of frequency-dependent seismic radiation from subduction zone megathrust ruptures. *Proceedings of the National Academy of Sciences*, **110**(12), 4512–4517.
# 4 Imager la distribution du couplage le long des grandes failles

Une approche plus classique pour identifier les zones de faille sur lesquelles de grands séismes sont susceptibles de se produire est d'estimer la distribution spatiale du couplage le long des grandes failles (voir Chapitre I.2.2). Pour cela, on inverse des données géodésiques long terme (GNSS ou InSAR) à la manière des inversions co-sismiques. Comme indiquée dans le Chapitre I.1.3, une faiblesse fondamentale de telles inversions est la non-unicité de leur solution. Dans l'étude qui suit, nous avons réalisé, à partir de données InSAR et GNSS, une inversion Bayésienne conjointe du couplage le long des failles Nord- et Est-Anatolienne ainsi que de la vitesse relative des plaques tectoniques Eurasienne, Anatolienne et Arabique. Cette approche nous permet d'estimer l'ensemble des solutions compatibles avec les observations. Nous avons ainsi pu estimer des incertitudes sur la distribution du couplage obtenue – et sur la vitesse long terme des plaques tectoniques – nous procurant une quantification probabiliste d'où s'accumulent les contraintes sur les failles.

Cette étude a été réalisée à Géoazur en collaboration avec Olivier Cavalié, Jean-Mathieu Nocquet et Théa Ragon. Elle a été publiée en 2020 dans le journal Geophysical Research Letters
[Bletery et al., 2020].
<figure><img src="image_27.png" /></figure>

# Geophysical Research Letters

## RESEARCH LETTER

10.1029/2020GL087775

**Key Points:**

* Distribution of interseismic coupling on the North and East Anatolian Faults

* Bayesian quantification of uncertainties on coupling and Euler poles

* The 2020 $M_w$ 6.8 Elazığ earthquake released 221.5 years (±26) of accumulated moment

**Supporting Information:**

* Supporting Information S1

**Correspondence to:**

Q. Bletery,
bletery@geoazur.unice.fr

**Citation:**

Bletery, Q., Cavalié, O., Nocquet, J. M., & Ragon, T. (2020). Distribution of interseismic coupling along the North and East Anatolian Faults inferred from InSAR and GPS data. *Geophysical Research Letters, 47*, e2020GL087775. https://doi.org/10.1029/2020GL087775

Received 4 MAR 2020
Accepted 5 AUG 2020
Accepted article online 17 AUG 2020

## Distribution of Interseismic Coupling Along the North and East Anatolian Faults Inferred From InSAR and GPS Data

Quentin Bletery¹ ID, Olivier Cavalié¹, Jean-Mathieu Nocquet¹,² ID, and Théa Ragon³ ID

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

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

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

¹Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, France, ²Institut de Physique du Globe de Paris, Université de Paris, CNRS, Paris, France, ³Seismological Laboratory, California Institute of Technology, Pasadena, CA, USA

**Abstract** The North Anatolian Fault (NAF) has produced numerous major earthquakes. After decades of quiescence, the $M_w$ 6.8 Elazığ earthquake (24 January 2020) has recently reminded us that the East Anatolian Fault (EAF) is also capable of producing significant earthquakes. To better estimate the seismic hazard associated with these two faults, we jointly invert interferometric synthetic aperture radar (InSAR) and GPS data to image the spatial distribution of interseismic coupling along the eastern part of both the NAF and EAF. We perform the inversion in a Bayesian framework, enabling to estimate uncertainties on both long-term relative plate motion and coupling. We find that coupling is high and deep (0–20 km) on the NAF and heterogeneous and superficial (0–5 km) on the EAF. Our model predicts that the Elazığ earthquake released between 200 and 250 years of accumulated moment, suggesting a bicentennial recurrence time.

**Plain Language Summary** Earthquakes are thought to occur on coupled fault portions, which are “locked” during the time separating two earthquakes, while tectonic plates are steadily moving. The spatial distribution of coupling has been imaged along numerous large faults in the world, but despite its considerable associated seismic hazard, not on the North Anatolian Fault (NAF). The recent $M_w$ 6.8 Elazığ earthquake (24 January 2020) has reminded us that the East Anatolian Fault (EAF) is also capable of producing large earthquakes. To better assess the seismic hazard associated with both the NAF and the EAF, we image the distribution of interseismic coupling along these faults. We find that the NAF is strongly coupled along most of the studied section. On the opposite, coupling is shallow and heterogeneous along the EAF. The initiation of the Elazığ earthquake coincides with a strongly locked but narrow (5 × 14 km) and superficial patch. The rest of the rupture extends over moderately coupled fault portions. We estimate that it took between 200 and 250 years to accumulate the moment released by the Elazığ event. Several fault segments along the EAF present similar coupling distributions, suggesting that, provided enough time, they could host earthquakes of similar magnitude.

## 1. Introduction

Earthquakes are thought to rupture fault portions that have previously accumulated a deficit of slip over tens to thousands of years (e.g., Avouac, 2015). Quantifying the spatial distribution of interseismic coupling—that is, the percentage of slip deficit with respect to the long-term drift of tectonic plates—along large faults is therefore crucial to anticipate earthquakes and better assess seismic hazard (e.g., Kaneko et al., 2010). The emergence of space geodetic techniques has allowed to infer interseismic coupling along a number of large faults during long quiescent periods of time separating one large earthquake to the next (e.g., Bürgmann et al., 2005; Jolivet et al., 2015; Loveless & Meade, 2011; Metois et al., 2016; Moreno et al., 2010; Nocquet et al., 2017; Protti et al., 2014). Though interseismic coupling models have been proposed to estimate the locking depth of the North Anatolian Fault (NAF) and East Anatolian Fault (EAF) (e.g., Aktug et al., 2013, 2016; Cavalié & Jónsson, 2014; Mahmoud et al., 2013; Tatar et al., 2012), none have quantified the lateral variations of coupling along these faults, which has limited the possibilities to study the spatial relationship between coupling and large earthquakes. The density of interferometric synthetic aperture radar (InSAR) observations (Cavalié & Jónsson, 2014) combined with sparser GPS measurements allows to infer these lateral variations of coupling on the eastern part of the NAF-EAF system (Figure 1).
<figure><img src="image_31.png" /><figcaption>**Figure 1.** (top) Tectonic setting of the NAF-EAF system (red lines) and study area (black rectangle). (bottom) Color maps show InSAR horizontal velocities (in a Eurasiafixed reference frame) in the satellite line-of-sight (LOS) direction (thick red arrows), ~103°N for Descending Tracks T264 and T493 (left) and ~77°N for Ascending Track T400 (right) (Cavalié & Jónsson, 2014). Black arrows show GPS measurements and their 95% ellipses of uncertainty (Nocquet, 2012; Ozener et al., 2010; Tatar et al., 2012). White diamonds indicate large (>100,000 people) cities.</figcaption></figure>
The eastern part of the NAF is known to produce large earthquakes (e.g., Ambraseys, 1971, 1989; Barka, 1996) and thought to be coupled from 0 to 15 km depth (Cavalié & Jónsson, 2014; Reilinger et al., 2006). On the other hand, simple back slip models showed that the EAF is weakly coupled and only in the first kilometers of the upper crust, from 0 to 5 km (Cavalié & Jónsson, 2014). This observation was in good agreement with the relative scarcity of large earthquakes recorded during the twentieth century (Burton et al., 1984; Jackson & McKenzie, 1988). For those reasons, the 24 January 2020 M<sub>w</sub> 6.8 Elazığ earthquake came as a surprise, on a segment that does not exhibit signs of past rupture (Duman & Emre, 2013) and in an area where the last earthquake of comparable magnitude (M<sub>s</sub> 6.8) occurred in 1905 (Nalbant et al., 2002). To understand this unexpected event, and more generally the seismicity in the region, we infer here the spatial distribution of interseismic coupling along the eastern part of the NAF-EAF system using InSAR (Cavalié & Jónsson, 2014) and GPS measurements (Nocquet, 2012; Ozener et al., 2010; Tatar et al., 2012) of the interseismic surface deformation.

Inferring spatially variable interseismic coupling along faults from geodetic observations—such as InSAR and GPS—of the Earth surface deformation requires solving an inverse problem that usually does not admit a unique solution (Nocquet, 2018; Tarantola & Valette, 1982). Most inversion techniques deal with this nonuniqueness by finding the solution that best fits the observations in a least squares sense, together with some roughness and/or damping penalty function. As a result, typical published coupling (or slip) models are the smoothest best fitting solutions among an infinity of possible models. We adopt here a Bayesian approach, which does not invert for a specific “ambiguously defined best solution” but explores the entire solution space, sampled with respect to the likelihood of each model. This approach—originally developed to invert for coseismic slip models (Minson et al., 2013)—also enables to reliably estimate uncertainties on coupling distributions (Jolivet et al., 2015, 2020).

## 2. Data

Our data set is composed of InSAR and GPS measurements in eastern Anatolia, all calculated in a stable Eurasia reference frame (Figure 1). Our InSAR data set is composed of two descending and one ascending tracks—all crossing both the NAF and EAF near their junction in eastern Turkey—processed by Cavalié and Jónsson (2014). Our GPS data set is composed of the horizontal components of 72 GPS stations located in the area (Nocquet, 2012; Ozener et al., 2010; Tatar et al., 2012).

InSAR data were derived from multiframe Envisat synthetic aperture radar (SAR) images provided by the European Space Agency. Each track includes between 16 and 19 SAR images acquired between 2003 and 2010. Interferograms were generated using the New Small BAseline Subset (NSBAS) processing chain (Doin et al., 2011). They were corrected for a ramp mostly due to a drift in the local oscillator on board the Envisat satellite (Marinkovic & Larsen, 2013). To avoid removing tectonic signals related to the motion of the Anatolian and Arabian plates, the ramps were estimated only on their Eurasian part that is considered as stable and orthogonal to the flight direction. All calculations were made considering stable Eurasia as a reference by setting the mean displacement of this area to 0, in the least squares sense. Surface displacement rates from the interferograms were derived using a small baseline time series approach, which maximizes coherence and the number of pixels to use in the analysis. A smoothing operator was applied to limit phase variations due to turbulent atmospheric delays. Finally, the linear component of the time series was extracted for each pixel in order to obtain the steady ground velocities. For a more detailed description of the InSAR processing, we refer the reader to the original study of Cavalié and Jónsson (2014).

Additionally, we compiled GPS data located between longitudes 38°E and 41°E and latitudes 35°N and 43°N from three independent studies. Velocity for 19 points were published by Tatar et al. (2012) derived from three surveys performed between 2006 and 2008. Another set of 19 points were published by Ozener et al. (2010) from three campaigns with 12-month interval. The remaining 34 points were originally published by Reilinger et al. (2006) and Reilinger and McClusky (2011) but recalculated in the continental-scale combination solution described in Nocquet (2012). The three data sets are expressed in a Eurasia-fixed reference frame. The lack of enough common sites shared among the three solutions prevents to properly combine them, but the few common sites and analysis of models residuals do not show any systematic pattern, suggesting that the three velocity fields are consistent within their uncertainties.
### 3. Bayesian Inversion of Rotation Poles and Interseismic Slip Deficit Rate Along Two Faults From InSAR and GPS Data

We invert the aforementioned InSAR and GPS measurements of the eastern Anatolia surface deformation to infer the distribution of interseismic slip deficit rate along the NAF-EAF system using a Bayesian sampling approach implemented in the AlTar1 package, originally developed by Minson et al. (2013) under the name of CATMIP. AlTar associates Markov chain Monte Carlo methods with a tempering process to explore the solution space, each step of the tempering being followed by a resampling to select only the most probable models. The probability density function (pdf) $p(\mathbf{m}|\mathbf{d})$ of a large number of likely models $\mathbf{m}$ given our data $\mathbf{d}$ is evaluated based on the ability of a model $\mathbf{m}$ to predict the data $\mathbf{d}$ (Minson et al., 2013):

$$p(\mathbf{m}|\mathbf{d}) \propto p(\mathbf{m}) \exp \left[ -\frac{1}{2}(\mathbf{d} - \mathbf{G}\mathbf{m})^T \mathbf{C}_{\chi}^{-1} (\mathbf{d} - \mathbf{G}\mathbf{m}) \right], \quad (1)$$

where $\mathbf{G}$ is the matrix of the Green's functions and $\mathbf{C}_{\chi}$ is the misfit covariance matrix. Vector $\mathbf{d}$ is composed of 144 GPS measurements (72 × 2 components) and a subset of InSAR pixels on the three tracks downsampled using the Quadtree algorithm (Jónsson et al., 2002).

Because the inferred distribution of coupling is presumably highly sensitive to the (usually) predetermined tectonic block motion, especially in a case involving three plates, we do not impose precalculated plate rotations but invert for them simultaneously with the interseismic slip deficit rate—similarly to the approach proposed by Meade and Loveless (2009) but adapted to a Bayesian framework. We discretize the eastern part of the NAF and EAF into 110 subfaults of depth-dependent sizes (Tables S1 and S2 in the supporting information) and invert for the model vector:

$$\mathbf{m} = \begin{pmatrix} \mathbf{w}^1 \\ \mathbf{w}^2 \\ \mathbf{S} \end{pmatrix}, \quad (2)$$

where $\mathbf{w}$ is the plate rotation vector expressed in Cartesian geocentric coordinates with unit of rad/year,¹ stands for Anatolia with respect to Eurasia,² for Arabia with respect to Eurasia, and $\mathbf{S}$ is the back slip on each subfault. Accordingly, we build $\mathbf{G}$ so that

$$\mathbf{G} = (\mathbf{A}, -\mathbf{G}_S), \quad (3)$$

where $\mathbf{A}$ is the matrix relating the plate rotation vectors to the horizontal velocities (see Appendix A) and $\mathbf{G}_S$ is the classical matrix of the Green's functions computed using the analytical solution of a shear finite fault embedded in an elastic half-space (Mansinha & Smylie, 1971; Okada, 1985).

$\mathbf{C}_{\chi}$ is the misfit covariance matrix, which translates data and epistemic uncertainties into uncertainties on the inverted model $\mathbf{m}$ (Bletery et al., 2016; Duputel et al., 2014; Ragon et al., 2018, 2019a, 2019b). Here, we only account for data uncertainties. For GPS records, we fill $\mathbf{C}_{\chi}$ with the (squared) standard deviations and covariances between the east and north components of a given station provided in the GPS solutions. For InSAR pixels, we first remove the tectonic signal from the unsampled interferograms using a preliminary model and calculate the covariance across the pixels of the residual interferograms as a function of their distances. We fit an exponential function (Figure S1) to the obtained cloud of points and express the covariance $C_{i,j}$ between two pixels as a function of their distance $D_{i,j}$:

$$C_{i,j} = a^2 \exp \left( \frac{-D_{i,j}}{b} \right), \quad (4)$$

by applying a regression to the parameters $a$ and $b$ independently on the three tracks (Jolivet et al., 2012, 2015, Sudhaus & Sigurjón, 2009). We then use Equation 4 to evaluate the covariance on the subsampled interferograms.

$p(\mathbf{m})$ is the pdf describing the prior information assumed on the different model parameters. We choose the less informative distributions for back slip parameters $\mathbf{S}$, that is, uniform distributions between 0 and the a priori long-term interplate velocities: 19.5 mm/year for the North Anatolian and 13 mm/year for the EAF (Cavalié & Jónsson, 2014). For the plate rotation vectors, we use the Euler poles and their associated
<table><caption>**Table 1**</caption><caption>A Priori (Le Pichon & Kreemer, 2010) and A Posteriori Euler Pole Coordinates and Angular Velocities With Respect to Eurasia</caption><thead><tr><th></th><th>Plate</th><th>Longitude (° E)</th><th>Latitude (° N)</th><th>Angular velocity (°/My)</th></tr></thead><tbody><tr><td rowspan="2">A priori</td><td>Anatolia</td><td>31.96 ± 0.10</td><td>32.02 ± 0.10</td><td>1.307 ± 0.083</td></tr><tr><td>Arabia</td><td>15.21 ± 0.10</td><td>28.31 ± 0.10</td><td>0.396 ± 0.010</td></tr><tr><td rowspan="2">A posteriori</td><td>Anatolia</td><td>34.22 ± 0.35</td><td>30.96 ± 0.60</td><td>1.087 ± 0.078</td></tr><tr><td>Arabia</td><td>16.13 ± 0.52</td><td>27.08 ± 0.37</td><td>0.386 ± 0.008</td></tr></tbody></table>

Note. A posteriori parameters are the mean and 2-σ standard deviation (95% confidence) of the posterior pdfs (Figure S6).

uncertainty from Le Pichon and Kreemer (2010) to derive a prior pdf. Plate rotation vectors (in Cartesian geocentric coordinates) **w**^p are related to Euler pole parameters through

$$ \mathbf{w}^p = \Omega^p \begin{pmatrix} \cos \phi^p \cos \lambda^p \\ \cos \phi^p \sin \lambda^p \\ \sin \phi^p \end{pmatrix}, \qquad (5) $$

where $\lambda^p$ and $\phi^p$ are the longitude and latitude of the Euler pole of a plate $p$ and $\Omega^p$ is its angular velocity (Bowring, 1985). Note that this change of coordinate system makes the problem linear (e.g., Maurer & Johnson, 2014; Meade & Loveless, 2009; Nocquet et al., 2001). We draw 100,000 sets of parameters ($\lambda^1, \phi^1, \Omega^1, \lambda^2, \phi^2, \Omega^2$) from normal distributions defined by means and standard deviations taken from previously published solutions (Le Pichon & Kreemer, 2010, summarized in Table 1). For each drawn set of parameters, we calculate the corresponding **w**^1 and **w**^2. We obtain Gaussian-like distributions for each component of **w**^1 (Figure S2) and **w**^2 (Figure S3). We extract the mean and standard deviation of these distributions and use them to define normal prior pdfs on $w_{x,y,z}^{1,2}$ in AlTar.

## 4. Results

We obtain a posterior marginal pdf for every inverted parameter in **m**, 110 fault slip parameters and 6 parameters describing the plate rotation vectors ($w_{x,y,z}^{1,2}$). The posterior pdfs on **w**^1 and **w**^2 parameters (Figure S4) appear uncorrelated (coefficients of correlation <0.013) with each other and—to a lesser extent—with fault slip parameters (coefficients of correlation <0.13) (Figure S5). Moderate anticorrelations are noticeable between fault slip parameters of patches located one beneath another (i.e., at the same location but different depth) (Figure S5a).

We convert the inverted pdfs on the rotation vectors ($\mathbf{w}^1, \mathbf{w}^2$) (Figure S4) into pdfs on the Euler pole coordinates and angular velocities (Figure S6). The means and 2-σ standard deviations of the inverted pdfs are summarized in Table 1. They are close to the previously published values we used as a prior (Le Pichon & Kreemer, 2010) but not equal (Figure S7). A possible explanation for this small discrepancy is that the plates are not strictly rigid (Aktug et al., 2013; England et al., 2016; Le Pichon & Kreemer, 2010; Nocquet, 2012), and thus, the rotations we invert from data in eastern Anatolia are slightly different from those obtained from data sampling on a larger area of the plate. Figure S8a shows the velocities corrected from plate motion using the Euler poles from Le Pichon and Kreemer (2010). It clearly shows a pattern of a residual rotation and unlikely large (5 mm/year) fault normal relative motion across both faults. Using our poles, residuals velocities appear to be consistent with the interseismic pattern (back slip) expected for strike-slip faults (Figure S8b). Our goal here is to infer the coupling distribution, and for that aim a refined estimate of the rotation parameters close to the fault is preferable to a plate-average solution, but one should be careful in using values in Table 1 for other purposes.

For each posterior Euler pole, we calculate the rotation predicted at the center of each patch and project the obtained vector along the fault strike direction to obtain posterior pdfs of the long-term slip rate along the faults (Figure S9). These pdfs are consistent with steady long-term slip rates of ~20 mm/year along the NAF and ~10 mm/year along the EAF (Figures 2 and S9 and Tables S1 and S2).

For each sampled model $\mathbf{m}_k = (\mathbf{w}_k^1, \mathbf{w}_k^2, \mathbf{S}_k)^T$, we divide the back slip parameters $\mathbf{S}_k$ by the long-term fault rate calculated at the center of each patch using the corresponding sampled Euler poles $\mathbf{w}_k^1$ and $\mathbf{w}_k^2$ to obtain
<figure><img src="image_32.png" /><figcaption>**Figure 2.** Interseismic coupling distribution inverted from InSAR and GPS data (mean of posterior pdfs in Figures S10 and S11). Black thick arrows indicate the long-term slip rate at depth derived from the inversion (mean and standard deviation of posterior pdfs in Figure S9). Focal mechanisms show *M* > 4.8 earthquakes (colors indicate event dates). Contours delineate the approximate rupture extent of the 1939 *M*<sub>S</sub> 8.0 Erzincan earthquake and of the 2020 *M*<sub>w</sub> 6.8 Elazığ earthquake (USGS finite fault solution). The light blue star indicates the epicenter of the Elazığ earthquake.</figcaption></figure>

the posterior marginal pdfs on the coupling coefficients (Figures S10 and S11). We show these pdfs in the form of their means (Figure 2) and standard deviations (Figure 3). Although restrictive, this representation gives an approximate view of the coupling spatial distribution and its associated uncertainties. Uncertainty is high (>25%) on the extreme west and—to a lesser extent—the extreme east parts of the fault system that are located outside of the InSAR tracks (Figure 1). The standard deviation on most parts of the faults is <20%, much lower on many subfaults (Figure 3). Note that standard deviation values are likely underestimated since we did not consider epistemic uncertainties here. The Earth structure is likely not homogeneous and the fault geometry not as simple as we modeled it, generating more uncertainties that we do not account for.

We calculate the GPS and InSAR measurements predicted for every posterior sampled model. We plot the predicted GPS means (red arrows) and 2-σ standard deviations (red ellipses) on Figure S12 and the residuals on Figure S13. For InSAR, we plot the mean predicted line-of-sight (LOS) displacements (Figures S14–S16)
<figure><img src="image_33.png" /><figcaption>**Figure 3.** Standard deviation of the coupling posterior pdfs. The extreme west and—to a lesser extent—the extreme east parts of the NAF-EAF system, presenting high (>25%) standard deviations, are located outside of the InSAR tracks.</figcaption></figure>

and standard deviations (Figure S17). The range of likely models that we found (Figures S10 and S11) is in very good agreement with both GPS and InSAR data. One way to quantify the relative amplitudes of residuals with respect to the observations is to calculate the ratio *r* of the mean of the absolute value of the residuals with the mean of the absolute value of the observations,

$$r = \frac{< |\mathbf{d} - \mathbf{d}_{\text{pred}}| >}{< |\mathbf{d}| >} \qquad (6)$$

This ratio is 15.9% for T264, 36.1% for T400, 24.3% for T493, and 21.6% for GPS. We attribute these reasonably small residuals—which do not exhibit coherent pattern (Figures S13 and S17)—to nontectonic sources. Furthermore, we find that every posterior sampled model predict very similar GPS and InSAR displacements; red ellipses are hardly visible on Figure S12, and the standard deviations of the predicted InSAR LOS displacements are very small (Figure S17). This highlights the limited resolution on the coupling model: If different models predict the same observations, discriminating between them is difficult.

## 5. Discussion

We show focal mechanisms of $M > 4.8$ earthquakes in the studied area from the Global Centroid Moment Tensor (GCMT) catalog (Dziewonski et al., 1981; Ekström et al., 2012) for events posterior to 1976 and from a compilation of historical earthquakes (Tan et al., 2008) for earlier events (1938–1976) (Figure 2). Focal mechanisms are represented at the location of their surface projections (i.e., at depth = 0). Colors indicate the dates of the events. The largest earthquake in the studied area is the 1939 $M_s 8.0$ Erzincan earthquake, which initiated near Erzincan and extended over the entire NAF segment west of Erzincan represented in
Figure 2 (Barka, 1996; Stein et al., 1997). We find that almost all of this section is strongly coupled, such as the rest of the studied NAF segment east of Erzincan. This easternmost segment of the NAF presents a moderate seismicity compared to the rest of the NAF. Our interseismic slip distribution suggests that it is as prone to generate large earthquakes as the rest of the NAF and as the Erzincan rupture segment in particular. In the middle of this overall strongly coupled (>75%) fault, we identify a few low-to-moderate coupling (10–50%) patches at depths between 5 and 10 km (Figure 2). These patches are associated with standard deviations between 5% and 25%, suggesting that these uncoupled patches are robust features. Interestingly, the most uncoupled patch coincides with the main step over of this section of the NAF. Step overs are thought to act as geometrical barriers that stop earthquake ruptures (e.g., Wesnousky, 2006). Although limited to one example, our results suggest that these geometrical features may also influence—or be influenced by—the interseismic behavior of the faults.

We find that locking on the EAF is much shallower with coupling values >50% limited to the first 5 km, consistently with previous studies (Cavalié & Jónsson, 2014). High coupling found at depth on the westernmost part of the fault is associated with standard deviations >20%, meaning that they are not reliable (Figure 3). Furthermore, we find that coupling also varies within the shallowest portion of the fault, alternating strongly coupled segments with weakly-to-moderately (0–60%) coupled ones (Figure 2). The most uncoupled shallow fault portion of the central EAF is located near Elazığ and coincides with the pull apart basin of Lake Hazar, as also observed on the Haiyuan fault (Jolivet et al., 2013). Different stress orientations around the basin could favor low coupling (Bertoluzza & Perotti, 1997; Wang et al., 2017; Van Wijk et al., 2017). This large reservoir of water may also provide the shallow part of the fault with fluids (although low resistivity associated to fluids is rather observed below 10 km depth, Türkoğlu et al., 2015) and locally weaken its mechanical friction, favoring asesimic slip. Such a behavior is observed both in laboratory and in situ (at the decametric scale) (Cappa et al., 2019). The mechanism invoked by the authors—consisting in an increase in nucleation length due to an increase in pore fluid pressure—may be at play at much larger scale here. On the other hand, the few earthquakes recorded on the EAF coincide with relatively high coupling. Before the recent Elazığ earthquake, the two largest events occurred near the localities of Bingol ($M_w$ 6.3, 2003) and Kovancilar ($M_w$ 6.1, 2010). The second one was followed by numerous aftershocks with magnitudes up to 5.6. All of these earthquakes occurred on >65% coupled fault portions, while fault segments with coupling <50% do not appear to have hosted $M > 4.8$ earthquakes.

According to the U.S. Geological Survey (USGS) finite-fault model (USGS, 2020), the Elazığ earthquake initiated between Elazığ and Malatya (light blue star in Figure 2) and propagated unilaterally westward (light blue contour in Figure 2). The early part coincides with a strongly locked (coupling coefficient: 100%) but narrow ($13.7 \times 5$ km) patch. The rupture seems to have then propagated throughout moderately coupled (coupling coefficient: 50–80%) fault segments. Although the USGS model is preliminary, its contours correlate fairly well with the coupling distribution, suggesting that the rupture stopped when reaching <25% coupled fault portions (Figure 2).

The last $M > 6.6$ earthquake in the approximate region dates back to 1905 ($M_s = 6.7$) (Nalbant et al., 2002). This event was located west of the recent Elazığ earthquake (38.6°E, 38.1°N) (Nalbant et al., 2002) but, given location uncertainties, could have ruptured the same fault portion. We calculate, for each sampled coupling model, the accumulated moment inside the rupture contour of the Elazığ earthquake since 1905. To simplify the problem, we assume that the earthquake ruptured the entire surface of the four main subfaults inside the rupture contour and not more, that is, the three shallowest subfaults plus the westernmost intermediate-depth one (Figure 3). We obtain a pdf of the seismic moment accumulated since 1905 (Figure 4a). The pdf mean is $7.3 \times 10^{18}$ N m, its standard deviation $0.8 \times 10^{18}$ N m. According to the USGS solution, the seismic moment released during the Elazığ earthquake is $13.87 \times 10^{18}$ N m—other solutions find even larger seismic moments (e.g., GCMT, Pousse-Beltran et al., 2020)—which is much larger than the $7.3 \pm 0.8 \times 10^{18}$ N m of moment deficit that we estimated since 1905. This seems to indicate that the recent Elazığ earthquake did not rupture the same fault portion than the 1905 earthquake. We further calculate the pdf of the time necessary to accumulate the seismic moment, which was released during the 2020 Elazığ earthquake (Figure 4b). The mean and standard deviation of the obtained pdf give a recurrence time for an Elazığ-type earthquake of $221.5 \pm 26$ years.
<figure><img src="image_34.png" /><figcaption>**Figure 4.** (a) Pdf of the accumulated seismic moment on the four patches inside the Elazığ rupture since 1905. The red vertical line indicates the seismic moment of the Elazığ earthquake according to the USGS solution (13.87 × 10¹⁸ N m). (b) Pdf of the time necessary to accumulate the seismic moment that was released during the Elazığ earthquake.</figcaption></figure>

## 6. Conclusion

We inverted InSAR and GPS observations to image interseismic coupling along the NAF and EAF in eastern Turkey. We adopted a Bayesian sampling approach in order to estimate posterior uncertainties on the coupling distributions and on the long-term fault rate. We did not impose a precalculated plate motion but inverted for the rotation of both the Anatolian and Arabian plates with respect to Eurasia, ensuring that the inferred coupling distribution is not biased in a systematic way by an inaccurate plate motion model. We found that the NAF is strongly coupled from 0 to 20 km depth, while the EAF is weakly coupled for the most part with high (>50%) coupling values limited to the shallowest part of the fault (0 to 5 km). Furthermore, we find that coupling is heterogeneous within this shallow portion, alternating seemingly creeping sections with strongly locked patches. Comparison between our interseismic coupling distribution and the preliminary finite-fault model of the USGS for the 2020 Mw 6.8 Elazığ earthquake reveals that this event likely initiated on one of this strongly locked (coupling coefficient: 100%) fault patches and then propagated into moderately coupled fault segments (coupling coefficient: 50–80%). Overall, we estimate that the Elazığ earthquake released 221.5 (±26) years of accumulated moment, suggesting a recurrence time ranging from 200 to 250 years.

## Appendix A: Rotation Matrix A

We build the rotation matrix **A** so that the motion due to the rotation of both the Anatolian and Arabian plates with respect to Eurasia equals **A** ⋅ **W**, where

$$ \mathbf{W} = \begin{pmatrix} \mathbf{w}^1 \\ \mathbf{w}^2 \end{pmatrix}. \qquad (A1) $$

Sorting all data points located on the Eurasian plate at the beginning of **d**, all data points located on the Anatolian plate in the middle and all data points located on the Arabian plate at the end, that is, writing **d** as

$$ \mathbf{d} = \begin{pmatrix} \mathbf{d}_0 \\ \mathbf{d}_1 \\ \mathbf{d}_2 \end{pmatrix}, \qquad (A2) $$

with **d**₀, **d**₁, and **d**₂ data points located on the Eurasian, Anatolian, and Arabian plates respectively, we can write **A** as a block matrix

$$ \mathbf{A} = \begin{pmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{A}' & \mathbf{0} \\ \mathbf{0} & \mathbf{A}' \end{pmatrix}, \qquad (A3) $$
so that $\mathbf{A} \cdot \mathbf{W}$ equals 0 for data points in Eurasia, $\mathbf{A}' \cdot \mathbf{w}^1$ in Anatolia, and $\mathbf{A}' \cdot \mathbf{w}^2$ in Arabia. $\mathbf{A}'$ is a transfer matrix relating the rotation vector in Cartesian geocentric coordinates $\mathbf{W}$ to the rotation block motion at each data point. It can be expressed at the location of an InSAR pixel or GPS station of longitude $\lambda$ and latitude $\phi$ as

$$ \mathbf{A}'_{\lambda, \phi} = \begin{pmatrix} -\sin \lambda & \cos \lambda & 0 \\ -\sin \phi \cos \lambda & -\sin \phi \sin \lambda & \cos \phi \\ \cos \phi \cos \lambda & \cos \phi \sin \lambda & \sin \phi \end{pmatrix} \begin{pmatrix} 0 & z & -y \\ -z & 0 & x \\ y & -x & 0 \end{pmatrix}, \quad (A4) $$

where

$$ \begin{pmatrix} x \\ y \\ z \end{pmatrix} = R_e (1 - \epsilon \sin^2 \phi)^{-1/2} \begin{pmatrix} \cos \phi \cos \lambda \\ \cos \phi \sin \lambda \\ (1 - \epsilon) \sin \phi \end{pmatrix}, \quad (A5) $$

with $R_e = 6,378.137$ km the Earth equatorial radius and $\epsilon = 0.00669438003$ the Earth eccentricity (Bowring, 1985).

Acknowledgments

The SAR data were provided by the European Space Agency through category1 project 6703, accessible at the following address (https://earth.esa.int/web/guest/data-access/). The inversion was performed using the AlTar software developed by Sarah Minson, Junle Jiang, Hailiang Zhang, Romain Jolivet, Zacharie Duputel, Michael Aivazis, James Beck, and Mark Simons at Caltech. We thank the Editor and Romain Jolivet and A. Ozgun Konca for thorough reviews.

References

Aktug, B., Dikmen, U., Dogru, A., & Ozener, H. (2013). Seismicity and strain accumulation around Karliova Triple Junction (Turkey). *Journal of Geodynamics*, **67**, 21–29.

Aktug, B., Ozener, H., Dogru, A., Sabuncu, A., Turgut, B., Halicioglu, K., et al. (2016). Slip rates and seismic potential on the East Anatolian Fault System using an improved GPS velocity field. *Journal of Geodynamics*, **94**, 1–12.

Ambraseys, N. (1971). Value of historical records of earthquakes. *Nature*, **232**(5310), 375–379.

Ambraseys, N. (1989). Temporary seismic quiescence: SE Turkey. *Geophysical Journal International*, **96**(2), 311–331.

Avouac, J.-P. (2015). From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle. *Annual Review of Earth and Planetary Sciences*, **43**, 233–271.

Barka, A. (1996). Slip distribution along the North Anatolian Fault associated with the large earthquakes of the period 1939 to 1967. *Bulletin of the Seismological Society of America*, **86**(5), 1238–1254.

Bertoluzza, L., & Perotti, C. R. (1997). A finite-element model of the stress field in strike-slip basins: Implications for the Permian tectonics of the Southern Alps (Italy). *Tectonophysics*, **280**(1-2), 185–197.

Bletery, Q., Sladen, A., Jiang, J., & Simons, M. (2016). A Bayesian source model for the 2004 great Sumatra-Andaman earthquake. *Journal of Geophysical Research: Solid Earth*, **121**, 5116–5135. https://doi.org/10.1002/2016JB012911

Bowring, B. R. (1985). The accuracy of geodetic latitude and height equations. *Survey Review*, **28**(218), 202–206.

Bürgmann, R., Kogan, M. G., Steblov, G. M., Hilley, G., Levin, V. E., & Apel, E. (2005). Interseismic coupling and asperity distribution along the Kamchatka subduction zone. *Journal of Geophysical Research*, **110**, B07405. https://doi.org/10.1029/2005JB003648

Burton, P. W., McGonigle, R., Makropoulos, K. C., & Üçer, S. B. (1984). Seismic risk in Turkey, the Aegean and the eastern Mediterranean: The occurrence of large magnitude earthquakes. *Geophysical Journal International*, **78**(2), 475–506.

Cappa, F., Scuderì, M. M., Colletti, C., Guglielmi, Y., & Avouac, J.-P. (2019). Stabilization of fault slip by fluid injection in the laboratory and in situ. *Science Advances*, **5**(3), eaau4065.

Cavalié, O., & Jónsson, S. (2014). Block-like plate movements in Eastern Anatolia observed by InSAR. *Geophysical Research Letters*, **41**, 26–31. https://doi.org/10.1002/2013GL058170

Doin, M.-P., Guillaso, S., Jolivet, R., Lasserre, C., Lodge, F., Ducret, G., & Grandin, R. (2011). Presentation of the small baseline NSBAS processing chain on a case example: The Etna deformation monitoring from 2003 to 2010 using Envisat data. In *Proceedings of the Fringe symposium*, pp. 3434–3437.

Duman, T. Y., & Emre, O. (2013). The East Anatolian Fault: Geometry, segmentation and jog characteristics. *Geological Society, London, Special Publications*, **372**(1), 495–529.

Duputel, Z., Agram, P. S., Simons, M., Minson, S. E., & Beck, J. L. (2014). Accounting for prediction uncertainty when inferring subsurface fault slip. *Geophysical Journal International*, **197**(1), 464–482.

Dziewonski, A. M., Chou, T.-A., & Woodhouse, J. H. (1981). Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. *Journal of Geophysical Research*, **86**(B4), 2825–2852.

Ekström, G., Nettles, M., & Dziewoński, A. M. (2012). The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. *Physics of the Earth and Planetary Interiors*, **200**, 1–9.

England, P., Houseman, G., & Nocquet, J.-M. (2016). Constraints from GPS measurements on the dynamics of deformation in Anatolia and the Aegean. *Journal of Geophysical Research: Solid Earth*, **121**, 8888–8916.

Jackson, J., & McKenzie, D. (1988). The relationship between plate motions and seismic moment tensors, and the rates of active deformation in the Mediterranean and Middle East. *Geophysical Journal International*, **93**(1), 45–73.

Jolivet, R., Lasserre, C., Doin, M.-P., Guillaso, S., Peltzer, G., Dailu, R., et al. (2012). Shallow creep on the Haiyuan Fault (Gansu, China) revealed by SAR interferometry. *Journal of Geophysical Research: Solid Earth*, **117**, B06401. https://doi.org/10.1029/2011JB008732

Jolivet, R., Lasserre, C., Doin, M.-P., Peltzer, G., Avouac, J.-P., Sun, J., & Dailu, R. (2013). Spatio-temporal evolution of aseismic slip along the Haiyuan Fault, China: Implications for fault frictional properties. *Earth and Planetary Science Letters*, **377**, 23–33.

Jolivet, R., Simons, M., Agram, P. S., Duputel, Z., & Shen, Z.-K. (2015). Aseismic slip and seismogenic coupling along the central San Andreas Fault. *Geophysical Research Letters*, **42**, 297–306. https://doi.org/10.1002/2014GL062222

Jolivet, R., Simons, M., Duputel, Z., Olive, J.-A., Bhat, H. S., & Bletery, Q. (2020). Interseismic loading of subduction megathrust drives long term uplift in northern Chile. *Geophysical Research Letters*, **47**, e2019GL085377. https://doi.org/10.1029/2019GL085377

Jónsson, S., Zebker, H., Segall, P., & Amelung, F. (2002). Fault slip distribution of the 1999 Mw 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements. *Bulletin of the Seismological Society of America*, **92**(4), 1377–1389.

Kaneko, Y., Avouac, J.-P., & Lapusta, N. (2010). Towards inferring earthquake patterns from geodetic observations of interseismic coupling. *Nature Geoscience*, **3**(5), 363–369.
Le Pichon, X., & Kreemer, C. (2010). The Miocene-to-Present kinematic evolution of the eastern Mediterranean and Middle East and its implications for dynamics. *Annual Review of Earth and Planetary Sciences*, **38**, 323–351.

Loveless, J. P., & Meade, B. J. (2011). Spatial correlation of interseismic coupling and coseismic rupture extent of the 2011 Mw = 9.0 Tohoku-oki earthquake. *Geophysical Research Letters*, **38**, L17306. https://doi.org/10.1029/2011GL048561

Mahmoud, Y., Masson, F., Meghraoui, M., Cakir, Z., Alchalbi, A., Yavasoglu, H., et al. (2013). Kinematic study at the junction of the East Anatolian fault and the Dead Sea Fault from GPS measurements. *Journal of Geodynamics*, **67**, 30–39.

Mansinha, L., & Smylie, D. E. (1971). The displacement fields of inclined faults. *Bulletin of the Seismological Society of America*, **61**(5), 1433–1440.

Marinkovic, P., & Larsen, Y. (2013). Consequences of long-term ASAR local oscillator frequency decay—An empirical study of 10 years of data. In *Living Planet Symposium*, European Space Agency Edinburgh, Frascati.

Maurer, J., & Johnson, K. (2014). Fault coupling and potential for earthquakes on the creeping section of the central San Andreas fault. *Journal of Geophysical Research: Solid Earth*, **119**, 4414–4428. https://doi.org/10.1002/2013JB010741

Meade, B. J., & Loveless, J. P. (2009). Block modeling with connected fault-network geometries and a linear elastic coupling estimator in spherical coordinates. *Bulletin of the Seismological Society of America*, **99**(6), 3124–3139.

Metois, M., Vigny, C., & Socquet, A. (2016). Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38–18 s). *Pure and Applied Geophysics*, **173**(5), 1431–1449.

Minson, S. E., Simons, M., & Beck, J. L. (2013). Bayesian inversion for finite fault earthquake source models I—Theory and algorithm. *Geophysical Journal International*, **194**(3), 1701–1726.

Moreno, M., Rosenau, M., & Oncken, O. (2010). 2010 Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone. *Nature*, **467**(7312), 198–202.

Nalbant, S. S., McCloskey, J., Steacy, S., & Barka, A. A. (2002). Stress accumulation and increased seismic risk in eastern Turkey. *Earth and Planetary Science Letters*, **195**(3-4), 291–298.

Nocquet, J.-M. (2012). Present-day kinematics of the mediterranean: A comprehensive overview of GPS results. *Tectonophysics*, **579**, 220–242.

Nocquet, J.-M. (2018). Stochastic static fault slip inversion from geodetic data with non-negativity and bound constraints. *Geophysical Journal International*, **214**(1), 366–385.

Nocquet, J.-M., Calais, E., Altamimi, Z., Sillard, P., & Boucher, C. (2001). Intraplate deformation in western Europe deduced from an analysis of the international terrestrial reference frame 1997 (ITRF97) velocity field. *Journal of Geophysical Research*, **106**(B6), 11,239–11,257.

Nocquet, J.-M., Jarrin, P., Vallée, M., Mothes, P. A., Grandin, R., Rolandone, F., et al. (2017). Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake. *Nature Geoscience*, **10**(2), 145–149.

Okada, Y. (1985). Surface deformation due to shear and tensile faults in a half-space. *Bulletin of the seismological society of America*, **75**(4), 1135–1154.

Ozener, H., Arpat, E., Ergintav, S., Dogru, A., Cakmak, R., Turgut, B., & Dogan, U. (2010). Kinematics of the eastern part of the North Anatolian Fault Zone. *Journal of Geodynamics*, **49**(3-4), 141–150.

Pousse-Beltran, L., Nissen, E., Bergman, E. A., Cambaz, M. D., Gaudreau, E., Karasözen, E., & Tan, F. (2020). The 2020 Mw 6.8 Elazığ (Turkey) earthquake reveals rupture behavior of the East Anatolian Fault. *Geophysical Research Letters*, **47**, e2020GL088136. https://doi.org/10.1029/2020GL088136

Protti, M., González, V., Newman, A. V., Dixon, T. H., Schwartz, S. Y., Marshall, J. S., et al. (2014). Nicoya earthquake rupture anticipated by geodetic measurement of the locked plate interface. *Nature Geoscience*, **7**(2), 117–121.

Ragon, T., Sladen, A., Bletery, Q., Vergnolle, M., Cavalié, O., Avallone, A., et al. (2019b). Joint inversion of coseismic and early postseismic slip to optimize the information content in geodetic data: Application to the 2009 Mw 6.3 L'Aquila earthquake, central Italy. *Journal of Geophysical Research: Solid Earth*, **124**, 10,522–10,543. https://doi.org/10.1029/2018JB017053

Ragon, T., Sladen, A., & Simons, M. (2018). Accounting for uncertain fault geometry in earthquake source inversions—I: Theory and simplified application. *Geophysical Journal International*, **214**(2), 1174–1190.

Ragon, T., Sladen, A., & Simons, M. (2019a). Accounting for uncertain fault geometry in earthquake source inversions—II: Application to the Mw 6.2 Amatrice earthquake, central Italy. *Geophysical Journal International*, **218**(1), 689–707.

Reilinger, R., & McClusky, S. (2011). Nubia-Arabia-Eurasia plate motions and the dynamics of Mediterranean and Middle East tectonics. *Geophysical Journal International*, **186**(3), 971–979.

Reilinger, R., McClusky, S., Vernant, P., Lawrence, S., Ergintav, S., Cakmak, R., et al. (2006). GPS constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions. *Journal of Geophysical Research*, **111**, B05411. https://doi.org/10.1029/2005JB004051

Stein, R. S., Barka, A. A., & Dieterich, J. H. (1997). Progressive failure on the north anatolian fault since 1939 by earthquake stress triggering. *Geophysical Journal International*, **128**(3), 594–604.

Sudhaus, H., & Sigurjón, J. (2009). Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: Application to the june 2000 Kleifarvatn earthquake, Iceland. *Geophysical Journal International*, **176**(2), 389–404.

Tan, O., Tapirdamaz, M. C., & Yörük, A. (2008). The earthquake catalogues for Turkey. *Turkish Journal of Earth Sciences*, **17**(2), 405–418.

Tarantola, A., & Valette, B. (1982). Inverse problems = Quest for information. *Journal of Geophysics*, **50**(1), 159–170.

Tatar, O., Poyraz, F., Gürsoy, H., Cakir, Z., Ergintav, S., Akpinar, Z., et al. (2012). Crustal deformation and kinematics of the eastern part of the north anatolian fault zone (Turkey) from GPS measurements. *Tectonophysics*, **518**, 55–62.

Türkoğlu, E., Unsworth, M., Bulut, F., & Çağlar, I. (2015). Crustal structure of the North Anatolian and East Anatolian Fault Systems from magnetotelluric data. *Physics of the Earth and Planetary Interiors*, **241**, 1–14.

USGS (2020). Finite fault model for the 01 24 2020 M 6.7 earthquake, Turkey. USGS finite-fault solution.

Van Wijk, J., Axen, G., & Abera, R. (2017). Initiation, evolution and extinction of pull-apart basins: Implications for opening of the Gulf of California. *Tectonophysics*, **719**, 37–50.

Wang, H., Liu, M., Ye, J., Cao, J., & Jing, Y. (2017). Strain partitioning and stress perturbation around stepovers and bends of strike-slip faults: Numerical results. *Tectonophysics*, **721**, 211–226.

Wesnosky, S. G. (2006). Predicting the endpoints of earthquake ruptures. *Nature*, **444**(7117), 358–360.
# 5 Perspectives

À la suite des études présentées dans les paragraphes II.2 et II.3, j’ai soumis un projet ERC Starting Grant (EarthS³Fault) à l’appel d’offre 2018 dans lequel je projetais d’utiliser les équations présentées dans ces deux études pour estimer, le long de l’ensemble des grandes de failles de la planète, la distribution de $\tau^c$, mais également de $\tau$. On peut en effet combiner les équations (7) et (A6) de l’article du paragraphe II.3 pour obtenir l’expression du taux de chargement en contraintes $\dot{\tau}$ suivante :

$$ \dot{\tau} = \frac{\sin 2\delta + \tan 2\psi \cos 2\delta}{\tan 2\psi} \dot{\sigma}_{xz}. \qquad (II-1) $$

La loi de Hooke donne

$$ \dot{\sigma}_{xz} = 2G\dot{\epsilon}_{xz}. \qquad (II-2) $$

$\dot{\epsilon}_{xz}$ étant une fonction du couplage, on peut exprimer $\dot{\tau}$ en fonction du couplage, du pendage apparent $\delta$ (lui-même fonction du pendage, de l’azimut et du rake, équation (5) de la partie II.3) et de l’orientation des contraintes $\psi$ [Hardebeck, 2015]. L’idée était ensuite de collecter un maximum d’informations sur la sismicité historique mondiale (temps depuis le dernier séisme $\Delta t$ et chute de contraintes $\Delta\sigma$ associée en chaque point de chaque faille) pour tenter d’estimer

$$ \tau \approx \tau^c - \Delta\sigma + \dot{\tau}\Delta t, \qquad (II-3) $$

puis d’introduire toutes ces contraintes dans des codes de simulations numériques pour tenter de prédire le comportement futur de chaque faille.

Le projet a été bien perçu et franchi la première phase. À l’issue de l’audition, le Jury a toutefois décidé de ne pas financer EarthS³Fault car il a jugé improbable que l’on puisse connaître assez bien les paramètres physiques en profondeur pour pouvoir utiliser des modèles à des fins prédictives. Après réception de cette décision, j’ai souscrit à ce point de vue. En particulier, la valeur de $\psi$ est d’une importance capitale dans les calculs et son estimation n’est pas résolue très finement spatialement [Hardebeck, 2015]. J’ai ainsi abandonné cette piste de recherche pour me focaliser sur d’autres, en particulier sur l’utilisation de l’IA pour l’alerte rapide et la détection de signaux précurseurs, pistes qui m’ont mené à soumettre un second projet ERC Starting Grant (EARLI) à l’appel d’offre 2020, second projet qui a lui été financé et dont les premiers résultats seront présentés dans le Chapitre IV (quelques idées du projet seront également esquissées dans la partie III.4).

Je continue, en revanche, à travailler sur des problèmes d’inversion – notamment à partir de
l'approche Bayésienne – pour imager le glissement lors des phases inter-sismique [Jolivet et al., 2020], mais également co-sismique [Ragon et al., 2021], post-sismique [Ragon et al., 2019b] et au cours d'épisodes transitoires [Bletery and Nocquet, 2020].
# Chapitre III

## Quand se produisent les grands séismes ?

### 1 Comment étudier la dynamique des grandes failles ?

Bien que nous ayons exploré des pistes intéressantes, nous avons vu qu'estimer où les grands séismes sont susceptibles de se produire est complexe. La question du « Où ? » est, toutefois, relativement accessible en comparaison de la question du « Quand ? ». Comment prédire quand vont se produire les prochains grands évènements ? Toutes les tentatives de réponse à cette question ont notoirement échoué.

La découverte des séismes lents – et de la micro-sismicité qui leur est associée – a ouvert la voie, à partir des années 2000 à un certain *monitoring* de l'activité des failles. L'étude de ces phénomènes permet en effet d'accéder à de l'information sur la dynamique des grandes failles. Lors de la phase inter-sismique, les failles ne sont pas « statiques », ou « figées » dans un état de couplage permanent. Elles présentent une dynamique complexe, observable par la géodésie (voir Chapitre I.2.3). Les propriétés des séismes lents renferment des informations sur les conditions physiques en profondeur. Les changements éventuels de ces propriétés au cours du temps indiquent une évolution des propriétés physiques en profondeur et pourraient potentiellement nous renseigner sur l'imminence d'un évènement. La micro-sismicité associée aux séismes lents apporte une information très riche, à haute résolution temporelle (les évènements micro-sismiques étant extrêmement nombreux). L'analyse de cette micro-sismicité peut ainsi servir à « mesurer » le glissement sur la faille en profondeur [Thomas et al., 2018]. Plusieurs études ont même proposé que des séismes lents ont précédé des grands séismes (voir Chapitre I.3.3).

Les séismes lents et la micro-sismicité associée constituent ainsi une grande opportunité d'étudier la dynamique des grandes failles. Depuis leur découverte, des séismes lents ont été identifiés dans de nombreuses régions du monde [Obara and Kato, 2016]. Qu'est-ce qu'un séisme lent ? Comment le glissement évolue lors de l'évènement qui peut durer plusieurs semaines ?