Léo, Joseph Alice Guillaume

Telly, Pierrick, Joseph Quentin Laurent

Keywords:

Il me semble que c'est lors de ma deuxième année de licence que la recherche, et de facto l'ambition de faire une thèse, m'a habité pour la première fois. Et ne m'a plus jamais quitté. Comment définir une thèse ? Une expérience scientifique, intellectuelle, humaine, une expérience de résilience surtout, et de découverte de soi. C'est une aventure solitaire -qui ne m'a pas déplu -et d'une richesse intellectuelle extraordinaire. Bien que solitaire, elle ne s'effectue pas seule ! Je pense évidemment tout d'abord à Jean-Mathieu. Nous avons commencé à travailler ensemble durant mon master puis en thèse. Merci de ta confiance durant ces cinq années. Ta rigueur scientifique et ton perfectionnisme dans le travail sont des qualités qui m'ont beaucoup inspirées et que j'admire. Tu m'as beaucoup appris, merci.

Merci également à Martin. Bien que nous n'ayons pas directement travaillé ensemble, tu t'es toujours montré disponible pour discuter de mon travail, et de mes opportunités de carrière.

Émilie et Christophe, merci de m'avoir permis de travailler sur le Chili. Merci également de m'avoir permis de découvrir le travail de terrain.

Je voudrais également remercier mes rapporteurs, Stéphane et Jeff, ainsi que Laëtitia, Marianne et Marie-Paule d'avoir examiné mon travail.

Ma vie de thésarde fut quelque peu nomade ! Entre l'IPGP, l'ENS et Géoazur. Le laboratoire de tectonique de l'IPGP m'accueille depuis le master 1. Yann, merci d'avoir fait attention à mon équilibre et de m'avoir renseigné et conseillé dans mes choix futurs. Ce fut précieux de comprendre comment construire sa carrière dans le monde académique. Je remercie également Raphaël, Arthur et Martine pour leur écoute et leur présence. J'ai beaucoup apprécié ces cinq années au sein de l'équipe, merci ! Je tiens également à remercier le laboratoire de géologie de l'ENS qui m'a accueilli durant la première moitié de ma thèse. Merci aux membres de l'équipe -dont Romain et Alexandre -et merci à mon cobureau Antoine, qui apportait toujours de la joie.

Je tiens à remercier Kristel, que j'ai retrouvé dans ces deux laboratoires. Les discussions avec toi ont été très enrichissantes. Ta gentillesse et générosité m'ont fait beaucoup de bien ! Merci.

Durant mes séjours à Géoazur, j'ai aussi été bien accueillie ! Un grand merci à Cédric. Tu as relu certains chapitres, tu m'as préparé à la soutenance. Au-delà du travail, ton amitié est précieuse. Merci également à François, Mohammed, Chao et Flo. Ce fut un plaisir d'être à Antibes en votre compagnie.

Nadaya, tu m'as donné l'opportunité de travailler avec toi en licence 3. J'ai découvert des questions passionnantes et le monde de la recherche académique. Je me suis beaucoup reposée sur les connaissances acquises avec toi ! Tu as été une encadrante incroyable et tu m'as conseillé d'aller à l'IPGP. Si j'ai réussi à faire cette thèse, c'est en partie grâce à toi. Merci également à Philippe, qui m'a suivi durant ce stage de licence 3. Tu m'as permis d'aller présenter mes travaux à Barcelone et tu m'as conseillé de travailler avec Jean-Mathieu. Merci de tes précieux conseils et de ta bienveillance. Je te dois beaucoup.

Merci à tous les enseignants que j'ai pu avoir depuis mes débuts. Je pense particulièrement à l'équipe pédagogique de la licence géologie à Paris 6, en particulier Laëtitia ! Tu es une enseignante exceptionnelle. Je pense aussi à Yves Noël et Philippe Huchon. En master à l'IPGP, j'ai aussi eu la chance d'avoir de bons enseignants. Je pense en particulier à Alexandre Fournier, Claude Jaupart et Alexandre Schubnel. Tant de professeurs inspirants durant ma scolarité m'ont permis d'avoir l'ambition et les moyens d'aller jusqu'à la thèse. Merci.

La thèse est une aventure -difficile -que l'on partage principalement avec d'autres thésards ! Spéciale dédicace à Sophie et Benjamin, mes cobureaux du 206*****, avec qui j'ai démarré l'aventure de la thèse et à qui j'ai permis de rendre les jeux de balles bien plus intéressants que d'accoutumé. Soso, merci pour ton écoute, ta sensibilité mais aussi ton humour noir ! Ben, ce fut un plaisir de découvrir quelqu'un d'aussi maniaque que moi ! Merci pour ta gentillesse et ton humour. La thèse aurait été bien plus triste sans vous ! Merci de votre soutien pendant cette aventure. Chloé, ta douceur et ta générosité sont des qualités qui m'ont beaucoup touché. Tu t'es toujours souciée de savoir si ça allait et tu ne te décourages pas pour essayer de me rendre moins dure ! Merci d'avoir été là toutes ces années. Merci aussi à Solène. Ces quelques mois dans le même bureau ont été très agréables ! Merci à Matthieu, avec qui j'ai eu de sacrés fous rires et qui me permettait de sortir un peu du travail quand on discutait. Javier, je suis très heureuse que tu aies décidé de venir faire une partie de ta thèse avec nous ! C'était un plaisir de partager ces quelques mois ensemble. Merci pour ta gentillesse et ton soutien. Je tiens aussi à remercier mon dernier cobureau, Nicolas, pour sa gentillesse dans les moments difficiles de fin de thèse ! Angèle, la rédaction du manuscrit de thèse je l'ai principalement (si ce n'est exclusivement !) vécu avec toi. Au rythme des jours, semaines et mois, matin, midi, soir et week-end, ta présence a été précieuse. Nous avons su nous épauler et rire de nos malheurs ! C'était aussi un plaisir d'avoir quelqu'un avec qui discuter de littérature, cinéma, théâtre. . . Je te remercie. Théo, tu t'es toujours soucié de savoir si ça allait. Tu n'as jamais failli et je t'en remercie. Ton soutien et ton calme (apparent) m'ont fait beaucoup de bien ! Tu as aussi toujours cru en moi, merci. Que ce soit à l'IPGP ou à l'ENS, je voudrais aussi remercier Alexandre, Manon, Thomas,

Vers une vision dynamique du glissement sur les failles

Les processus asismiques

Résumé

Pour mieux comprendre le comportement des failles, et donc peut-être un jour prévoir les grands séismes, il est nécessaire d'étudier tous les phénomènes qui s'y déroulent. Dans cette thèse, nous nous intéressons aux glissements asismiques, bien plus fréquents que les grands séismes, et qui contribuent à une partie importante du budget de glissement le long des failles. Le premier objectif de cette thèse était de mieux documenter la dynamique des glissements asismiques au cours du temps et d'identifier ses possibles interactions avec la sismicité. Le second objectif était de mesurer l'étendue et la fréquence des petits séismes lents, encore peu étudiés aujourd'hui. Avec une nouvelle méthode de modélisation spatio-temporelle des données GNSS, nous nous sommes intéressés au glissement post-sismique du grand séisme d'Illapel (M w 8.3) qui a eu lieu en 2015 le long de la subduction Chilienne, puis à une séquence mixant des séismes de taille intermédiaire et des glissements asismiques, qui a eu lieu en 2020 dans la région d'Atacama au Chili. Pour le séisme d'Illapel, nous montrons que le glissement post-sismique se développe à la périphérie de la zone de rupture. Cependant, une partie du glissement asismique a lieu dans la région rompue par le séisme, une observation qui remet en question les modèles classiques proposés pour décrire les propriétés frictionnelles des failles. Nous montrons également que le glissement post-sismique est très dynamique, avec plusieurs modulations dans le temps. En se développant le long de la faille, il est capable de conduire des aspérités à la rupture, à travers des répliques de magnitude 6.8-6.9. Pour la première fois, nous montrons aussi que le glissement post-sismique et un séisme lent peuvent avoir lieu simultanément et au même endroit. Ces résultats mettent en lumière une organisation spatiale et temporelle des glissements sismique et asismique dans les régions où le couplage inter-sismique est hétérogène. L'étude de la séquence d'Atacama illustre un processus où glissements sismique et asismique interagissent au cours du temps, expliquant la rupture de répliques de magnitudes proches de celle du choc principal pendant la première journée. Le glissement asismique qui s'ensuit est anormalement élevé. Il suggère aussi une réactivation du glissement dans une zone où un séisme lent a eu lieu en 2014. Enfin, je présente un outil de détection semi-automatique de petits séismes lents, que nous avons appliquons à la subduction Équatorienne. Des résultats préliminaires suggèrent que les petits séismes lents peuvent être bien plus fréquents que ce qui est documenté aujourd'hui, démontrant l'importance de les étudier pour corriger les déficits de glissements qui sont établis uniquement à partir des modèles inter-sismiques, afin de mieux évaluer l'aléa sismique.

Mots-clés

Zone de subduction, séismes, répliques, glissement post-sismique, séisme lent, don-nées GNSS, modélisation spatio-temporelle de glissement, micro-sismicité

Towards a dynamic view of slip along faults

The aseismic processes Abstract To improve our knowledge on the behavior of faults, and thus perhaps one day predict large earthquakes, it is necessary to study all the phenomena that take place on them. In this thesis, we focus on aseismic slip, which is much more frequent than great earthquakes, and contributes an important part of the slip budget along faults. The first objective of this thesis was to better document the dynamics of aseismic slip over time and to identify its possible interactions with seismicity. The second objective was to measure the extent and frequency of small slow slip events, which are still little studied today. With a new method of spatio-temporal modelling of GNSS data, we first studied the afterslip of the M w 8. [START_REF]The M w 8.3 2015 Illapel afterslip spatio-temporal modelling[END_REF] Illapel earthquake that occurred in 2015 along the Chilean subduction and then studied a sequence mixing intermediate-sized earthquakes and aseismic slip, that occurred in 2020 in the Atacama region of Chile. For the Illapel earthquake, we show that the afterslip develops at the periphery of the rupture zone. However, part of the aseismic slip occurs in the region ruptured by the earthquake, an observation that challenges the classical models proposed to describe the frictional properties of faults. We also show that afterslip is highly dynamic, with several modulations through time. As it develops along the fault, it is capable of driving asperities to rupture, through aftershocks of magnitude 6.8-6.9. For the first time, we also show that afterslip and a slow slip event can occur simultaneously and at the same location. These results highlight a spatial and temporal organization of seismic and aseismic slip in regions with heterogeneous interseismic coupling. The study of the Atacama sequence illustrates a process where seismic and aseismic slip interact over time, explaining the rupture of aftershocks of magnitudes close to that of the mainshock during the first day. The ensuing aseismic slip is anomalously high. It also suggests a reactivation of the slip in an area where a slow slip event occurred in 2014. Finally, I present a semi-automatic detection tool for small slow slip events, which we have applied to the Ecuadorian subduction. Preliminary results suggest that small slow slip events may be much more frequent than what is documented today, emphasizing the interest of accurately documenting them, to correct the slip deficit which are established only from interseismic models, in order to better assess the seismic hazard.

Keywords

Subduction zone, earthquakes, aftershocks, afterslip, slow slip event, GNSS data, slip time-dependent modelling, microseismicity

Résumé substantiel

Les failles de subduction produisent les plus grands séismes de notre planète. Ces derniers étant rares, notre compréhension du comportement des failles reste encore aujourd'hui incomplète. Si les séismes correspondent à des glissements rapides le long des failles, d'autres modes de glissement, dits "asismiques", ont lieu le long de ces failles et sont bien plus fréquents que les grands séismes. On en dénombre trois types différents à ce jour. Tout d'abord, les séismes lents, qui sont des évènements transitoires caractérisés par une faible vitesse de glissement, ponctuée d'une accélération puis d'une décélération. Le glissement post-sismique a lieu après un séisme et correspond à la réponse de la faille au changement de contraintes induit par le séisme. Ce glissement post-sismique est généralement observé en périphérie de la zone de rupture du séisme. Enfin, le glissement pre-sismique a lieu avant un grand séisme. Si plusieurs séismes lents ont été observés avant les séismes, il n'existe pas d'observation claire d'une accélération continue atteignant les vitesses de glissement sismique.

Ces différents types de glissements asismiques contribuent à une partie importante du budget de glissement le long des failles. Ainsi, dans cette thèse, je me suis intéressée aux séismes lents et aux glissement post-sismique avec deux objectifs. Le premier objectif était de mieux documenter la dynamique des glissements asismiques au cours du temps et d'identifier ses possibles interactions avec la sismicité. Le second objectif était de mesurer l'étendue et la fréquence des petits séismes lents, encore peu étudiés car proches du niveau de bruit des données GNSS. Pour cela, nous avons utilisé une nouvelle méthode de modélisation spatio-temporelle des données GNSS.

Le premier cas d'étude de cette thèse est la phase post-sismique du grand séisme de 2015 de M w 8.3 d'Illapel au Chili, présentée dans le chapitre 3. L'originalité de la phase post-sismique de ce séisme est que plus d'un mois et demi après le séisme, une grande réplique de M w 6.8 a lieu en profondeur de la rupture du choc principal. Elle est suivie quatre jours plus tard par deux grandes répliques de M w 6.9, qui se produisent à plus de 100 km au nord de la rupture principale. Une modélisation spatio-temporelle du glissement post-sismique journalier a permis une analyse fine de sa dynamique et de ses interactions avec la sismicité. Nous avons mis en évidence l'interaction entre les grandes répliques et le glissement post-sismique. D'une part, le développement du glissement post-sismique le long de la faille, ajouté à l'augmentation locale des contraintes due à la rupture du choc principal, a déclenché la rupture des deux répliques de M w 6.9 distante en temps et en espace du choc du principal. D'autre part, la vitesse de glissement postsismique accélère localement après la rupture des grandes répliques. Cette accélération correspond au glissement post-sismique des répliques qui s'ajoute à celui du choc principal. Nous avons également montré qu'un séisme lent se produit au sud de la rupture du choc principal, sans réplique majeure. L'approche spatio-temporelle a montré que ce séisme lent est corrélé avec une augmentation du taux de micro-sismicité. Cette étude montre que le glissement post-sismique est hautement dynamique, avec les effets des répliques majeures et des séismes lents qui se superposent à la décroissance générale des vitesses de glissement.

La deuxième cas d'étude est la phase post-sismique de la séquence de septembre 2020 dans la région d'Atacama au Chili, présentée dans le chapitre 4. Cette séquence a eu lieu quelques mois après l'initiation en mars 2020 d'un séisme lent dans une région de l'interface qui avait déjà hébergé un séisme lent en 2014 [START_REF] Klein | Return of the Atacama deep slow slip event: The 5-year recurrence confirmed by continuous GPS[END_REF]. La séquence a débuté par un séisme de M w 6.9, qui a été suivi le même jour par deux évènements de M w 6.3 et M w 6.4. Cette séquence est originale par la taille des principales répliques de magnitudes proches du séisme principal et en raison de sa proximité (environ 50 km) avec le séisme lent. Un rapide glissement post-sismique déduit de la migration rapide des répliques peut expliquer leur localisation et leur occurrence dans les heures suivant le choc principal. À partir du jour suivant, le glissement postsismique commence à une faible profondeur, chevauchant partiellement la rupture du choc principal. Simultanément, la zone du séisme lent subit un glissement asismique accéléré. L'incrément de contraintes induit par la séquence pourrait avoir accélérée ce séisme lent, et cette accélération a été capturée par la modélisation spatio-temporelle du glissement post-sismique de cette séquence.

La troisième étude de cette thèse porte sur la détection et la caractérisation des petits séismes lents. Depuis leur découverte en 1999 [START_REF] Hirose | A slow thrust slip event following the two 1996 Hyuganada earthquakes beneath the Bungo channel, southwest Japan[END_REF], de nombreux séismes lents de M w ≥ 6 ont été documentés le long des failles de subduction [START_REF] Lowry | Transient fault slip in Guerrero, southern Mexico[END_REF][START_REF] Dragert | A silent slip event on the deeper Cascadia subduction interface[END_REF]. Pourtant, si les séismes lents suivent une loi fréquence-magnitude, les petits séismes lents devraient être bien plus fréquents que les grands (Michel et al., 2019b). Afin de pouvoir vérifier cette hypothèse, j'ai développé une méthode de détection basée directement sur la modélisation des données GNSS, sans a priori concernant la région ou la période temporelle des séismes lents. Une première application de la méthode a été effectuée sur la subduction sud-équatorienne entre janvier 2010 et novembre 2022, où des séismes lents ont déjà été documentés. Les résultats préliminaires montrent une alternance entre des périodes "silencieuses" où le moment relâché est nul, durant lesquels il n'y a donc pas d'évènements, et des périodes où le moment relâché est non nul, synonyme d'un évènement asismique. Pour la première fois, l'approche utilisée détecte de petits évènements, de magnitudes équivalentes à M w ∼ 5.5. Sur toute la période de détection, l'outil trouve en moyenne deux évènements par an, avec une majorité de petits évènements (24 évènements dont le moment < 0.25 × 10 18 N.m, soit M w < 5.5), et 7 "grands" évènements (moment ≥ 0.25 × 10 18 N.m, soit M w ≥ 5.5). Ces résultats préliminaires corroborent l'hypothèse qu'il y a plus de petits séismes lents que de grands. Même si la technique développée détecte finement les évènements, elle ne permet pas encore de bien les caractériser, que ce soit pour leur durée ou leur localisation. Une exploration des paramètres de régularisation pour mieux modéliser les évènements et une détection appliquée à toute la subduction équatorienne sont nécessaires pour discuter de la fréquence et de la récurrence des séismes lents en Équateur.

Les différents projets réalisés pendant cette thèse ont permis de définir plus précisément la dynamique des glissements asismiques le long de la subduction sud-américaine. L'étude du glissement post-sismique du séisme d'Illapel (chapitre 3) a montré qu'il pouvait avoir lieu dans la région rompue par le choc principal. Cette observation a déjà été faite pour d'autres grands séismes (e.g. le séisme de 2010 M w 8.8 à Maule au Chili et le séisme de 2011 M w 9.0 à Tohoku-Oki au Japon), elle implique que certaines portions de la faille ne se rebloquent pas immédiatement et continuent de manière asismique.

L'étude jointe du glissement post-sismique du grand séisme d'Illapel et de ses grandes répliques ont montré qu'ils étaient étroitement liés. Le développement du glissement post-sismique le long de la faille, additionné au changement de contraintes induit par la rupture du choc principal, ont déclenché la rupture des répliques de M w 6.8 et 6.9. Ce phénomène a déjà été observé pour le grand séisme de M w 7.8 qui a eu lieu en 2017 près de Pedernales en Équateur. Deux grandes répliques de M w 6.7 et 6.9 ont eu lieu un mois après le séisme et corrélent spatiallement avec le développement du glissement post-sismique et la périphérie de la rupture du choc principal. L'interaction mis en lumière entre séisme, glissement post-sismique et grandes répliques montre que s'il est modélisé en temps réel, le glissement post-sismique pourrait déterminer où les grandes répliques vont rompre.

La fine modélisation spatio-temporelle du glissement post-sismique du séisme d'Illapel de 2015 au Chili a également permis de mettre en évidence un petit séisme lent qui a eu lieu après un mois et demi au sud de la rupture principale, en l'absence de réplique majeure. C'est la première fois qu'un séisme lent est observé dans une région de glissement post-sismique. Cette portion de faille a expérimenté plusieurs essaims sismiques avant le séisme d'Illapel et correspond à une région où il y a plusieurs failles intra-slab qui résultent de la subduction de la ride Juan Fernandez de la plaque Nazca (Poli et al., 2017). La géométrie complexe de cette région pourrait expliquer l'occurrence de ce séisme lent.

Les études réalisées dans cette thèse ont mise en évidence une organisation spatiale et temporelle particulière des glissements sismique et asismique. Les régions où le glissement post-sismique et les séismes lents superficiels se développent de manière privilégiée correspondent aux régions actives durant la phase inter-sismique, et ce à travers des essaims sismiques. Cette organisation semble être observée pour divers grands séismes de subduction parmi lesquels le séisme de 2011 de Tohoku-Oki au Japon (M w 9.0), le séisme de 2012 de Nicoya au Costa Rica (M w 7.6) ainsi que le séisme de 2017 de Pedernales en Équateur (M w 7.8). Ainsi, les régions expérimentant essaims sismiques et/ou séismes lents durant la phase inter-sismique, agissent comme des barrières à la propagation des grands séismes, et hébergent des séismes lents récurrents durant la phase post-sismique. Ainsi, bien modéliser les évènements asismiques et sismiques durant la phase inter-sismique devrait permettre d'identifier les régions à prédominance asismique, qui se situent à la périphérie des aspérités sismiques, et donc permettrait de les délimiter.

L'étude de la phase post-sismique de la séquence d'Atacama au Chili qui a eu lieu en 2020 a montré que cette séquence avait une forte composante asismique. En effet, la séquence a relâché un moment total équivalent à une M w 7.1, dont 60% est sismique et 40% est asismique. Les subductions japonaise et péruvienne ont également connu des séquences similaires [START_REF] Hirose | A slow thrust slip event following the two 1996 Hyuganada earthquakes beneath the Bungo channel, southwest Japan[END_REF]Villegas-Lanza, Nocquet, et al., 2016). L'origine du glissement asismique durant ces séquences n'est pas encore bien comprise, mais ces dernières relâchent une quantité non négligeable des contraintes le long des failles. Ces résultats soulignent l'intérêt de documenter précisément tous les glissements asismiques le long des failles, pour corriger les déficits de glissements qui sont établis uniquement à partir des modèles inter-sismiques, afin de mieux évaluer l'aléa sismique.

Chapter 1. Introduction

Earthquakes and plate tectonics

Earthquakes are not random, neither in space nor in time. They mostly occur along plate boundaries and are caused by the relative motion of tectonic plates. The plate tectonics concept arose from the continental drift, proposed over a century ago by Alfred [START_REF] Wegener | Die Entstehung der Kontinente und Ozeane: Braunschweig[END_REF] and provides a framework to explain numerous geological, paleontological, geophysical and geodetic observations. Later, the driver of plate tectonics was proposed with the first mantle convection model [START_REF] Holmes | Principles of physical geology[END_REF], eventually corroborated by the seafloor spreading model [START_REF] Hess | History of ocean basins[END_REF]. There are three different types of plate boundaries, each associated to different kinematics. A first type of plate boundary is the oceanic accretion zone, located at the oceanic ridges. They are regions of oceanic crust production by a basaltic volcanism. A second type of plate boundary is the subduction zone, where there is a large-scale sinking of an oceanic plate under another plate, often continental. It leads to the formation of a Wadati-Benioff zone -where the earthquake epicenters are located, an oceanic trench and a rich volcanism in the overriding plate. The third type of plate boundary is the transform fault, where there is no input or absorption of material, and plates are sliding parallel to the relative motion of the plates. The second type of plate boundary -subduction -is where the largest earthquakes occur. Large subduction earthquakes and processes occurring around them are the main focus of this thesis. Earthquakes have been documented for centuries, historically in writings and paintings [START_REF] Musset | Villes nomades du nouveau monde[END_REF][START_REF] Bernard | Pourquoi la Terre tremble[END_REF]. During the last century, a multitude of different types of observations has been used to study them. As a non-exhaustive list, one can cite the petrological (e.g. [START_REF] Sibson | Earthquakes and rock deformation in crustal fault zones[END_REF][START_REF] Fagereng | Mélange rheology and seismic style[END_REF][START_REF] Behr | What's down there? The structures, materials and environment of deep-seated slow slip and tremor[END_REF], paleoseismological (e.g. [START_REF] Sieh | Earthquake supercycles inferred from sea-level changes recorded in the corals of west Sumatra[END_REF][START_REF] Goldfinger | Superquakes and supercycles[END_REF]Weil-Accardo et 30 Chapter 1. Introduction al., 2020), seismological (e.g. [START_REF] Cotton | Frequency domain inversion of strong motions: Application to the 1992 Landers earthquake[END_REF][START_REF] Olsen | Three-dimensional dynamic simulation of the 1992 Landers earthquake[END_REF], geodetic observations (e.g. [START_REF] Massonnet | The displacement field of the Landers earthquake mapped by radar interferometry[END_REF]Freymueller et al., 1994;[START_REF] Avouac | The 2005, Mw 7.6 Kashmir earthquake: Sub-pixel correlation of ASTER images and seismic waveforms analysis[END_REF]. These different observations have made it possible to document earthquakes and their effects over different space and time scales, each data set having its own spatio-temporal resolution. These observations have led to the discovery of several phenomena, in particular before and after earthquakes, leading to the emergence of the concept of the seismic cycle [START_REF] Reid | The California earthquake of April 18, 1906: The mechanics of the earthquake[END_REF]. The interseismic strain grows linearly through time during decades to centuries, only disturbed by episodic slow slip events lasting days to a few years. The preseismic strain is millimetric to centimetric and can last seconds to months. The instantaneous coseismic strain induced by an earthquake is metric. The postseismic strain is millimetric to a few tens of centimeters, depending on the magnitude of the earthquake and can last weeks to years.

Time

The seismic cycle

The concept of seismic cycle is used to describe the behavior of faults through time, i.e. the temporal distribution of slip modes (seismic and aseismic), and their respective characteristics. During the last decades, the seismic cycle has been linked to friction laws, empirically determined from laboratory experiments on rocks. These laws allow us to infer the frictional characteristics driving each of the slip modes. The most commonly law used to study the seismic cycle is the rate-and-state friction law. These two concepts, i.e. the seismic cycle and the rate-and-state friction law, are detailed in this section and illustrated in Figures 1.1 and 1.2. For years to decades and even centuries, the fault and its surrounding medium accumulate stress due to the plate convergence and is progressively elastically deformed. This process corresponds to the loading period of the fault, referred to as the interseismic phase (see Figures 1.1 and 1.2a). The loading rate depends on the convergence velocity. Interseismic strain is often assumed to be constant in time. This vision, however, has been challenged by recent studies based on geodetic and seismological observations. These studies show that aseismic slip events, i.e. events that do not produce seismic waves, occur during the interseismic phase [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event[END_REF]Bartlow et al., 2014;[START_REF] Uchida | Periodic slow slip triggers megathrust zone earthquakes in northeastern Japan[END_REF][START_REF] Rousset | A geodetic matched filter search for slow slip with application to the Mexico subduction zone[END_REF]Michel et al., 2019a). A paleoseismological study on recent and fossil corals also showed that on several seismic cycles, there are decadal variations of the interseismic strain, suggesting that the loading rate of a fault is not constant in time [START_REF] Meltzner | Time-varying interseismic strain rates and similar seismic ruptures on the Nias-Simeulue patch of the Sunda megathrust[END_REF]. Thus, rather than accumulating stress in a continuous manner, the faults release part of the stress accumulated not only through seismic events but also through episodic aseismic slip.

When the stress applied to the fault reaches a certain threshold, all or most of the accumulated stress is suddenly released with a very high slip-rate, which leads to the generation of seismic waves. This is called an earthquake and corresponds to the coseismic phase (see Figures 1.1 and 1.2c).

Before an earthquake, the fault is believed to gradually deviate from its stationary state. This is the nucleation phase of an earthquake, otherwise called the preseismic phase (see Figures 1.1 and 1.2b). However, very few observations have been made for this phase, and it is therefore poorly understood. Two competing models are usually proposed to explain what happens during the preseismic phase: the preslip model and the cascade model [START_REF] Mclaskey | Earthquake initiation from laboratory observations and implications for foreshocks[END_REF]. The latter suggests that small to moderate earthquakes occur before large earthquakes, referred to as foreshocks. At some point, one will trigger the future large mainshock. The former suggests that large earthquakes nucleate from a progressively accelerating aseismic slip event. The 2011 M w 9.0 Tohoku-Oki earthquake in Japan experienced both a cascade of foreshocks and a long-term preslip [START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Uchida | Periodic slow slip triggers megathrust zone earthquakes in northeastern Japan[END_REF]. Two large earthquakes were triggered by a short-term preslip: the 2001 M w 8.4 Arequipa earthquake in Peru [START_REF] Ruegg | Co-, post-and pre (?)-seismic displacements associated with the Mw 8.4 Southern Peru earthquake of 23 June 2001 from continuous GPS measurements[END_REF][START_REF] Melbourne | Precursory transient slip during the 2001 Mw=8.4 Peru earthquake sequence from continuous GPS[END_REF] and the the 2014 M w 7.3 Papanoa earthquake in Mexico (Radiguet et al., 2016). The preseismic phase of the 1999 M w 7.6 Izmit earthquake in Turkey is still controversial. [START_REF] Ellsworth | Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks[END_REF] argue that there was a cascade of foreshocks while [START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF] suggest that there were repeating earthquakes, which they interpret as evidence of preslip. Another controversy exists on the 2014 M w 8.1 Iquique earthquake in Chile. [START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake[END_REF] propose that the mainshock was triggered by a cascade of foreshocks and preslip in the form of a slow slip event, while [START_REF] Schurr | Gradual unlocking of plate boundary controlled initiation of the 2014 Iquique earthquake[END_REF] argue that the area experienced a change of its frictional properties during the years prior to the mainshock, induced by several earthquake clusters and transient aseismic slip. [START_REF] Twardzik | Bayesian inference on the initiation phase of the 2014 Iquique, Chile, earthquake[END_REF] showed that just before the earthquake, the megathrust experienced both foreshocks and preslip but in distinct areas. [START_REF] Socquet | An 8 month slow slip event triggers progressive nucleation of the 2014 Chile megathrust[END_REF] showed that a long-term preslip took place during the months before the earthquake. However, for both models (i.e. the cascade model and the preslip), many other large earthquakes do not show any preseismic activity, or at least it has not been detected.

The interseismic and coseismic phases were introduced in a first concept of seismic cycle called the elastic rebound model, developed by Harry Fielding Reid and based on observations following the 1906 M w ∼ 7.8 San Francisco earthquake that ruptured a part of the San Andreas fault [START_REF] Reid | The California earthquake of April 18, 1906: The mechanics of the earthquake[END_REF]. But since the 1970s, another phase has been brought to light using the geodetic data from the 20 years of data following the 1906 San Francisco earthquake [START_REF] Thatcher | Strain release mechanism of the 1906 San Francisco earthquake[END_REF][START_REF] Thatcher | Strain accumulation and release mechanism of the 1906 San Francisco earthquake[END_REF].

After an earthquake, the surrounding medium continues to release stresses in part because of the stress increment induced by the mainshock. This is the postseismic phase, which lasts from a few hours to several years or decades, depending on the magnitude of the earthquake (see Figures 1.1 and 1.2d). During this phase, the stress continues to be released by several processes. A first process is the occurrence of aftershocks. While counting the number of aftershocks of the 1891 M w 7.5 Nōbi earthquake in Japan, the seismologist Fusakichi Omori noticed that they decrease over time in a certain way: the rate of aftershocks follows a power law as a function of time. This law predicts the number of aftershocks of an earthquake as a function of time (Omori, 1894), which was later amended by Tokuji Utsu (1957). The second process is the aseismic slip -called afterslip. [START_REF] Thatcher | Strain release mechanism of the 1906 San Francisco earthquake[END_REF][START_REF] Thatcher | Strain accumulation and release mechanism of the 1906 San Francisco earthquake[END_REF] put forward the concept of afterslip for the first time by studying geodetic data from 20 years after the 1906 M w ∼ 7.8 San Francisco earthquake. For years, it has been observed in the direct vicinity of the coseismic rupture, driving the aftershocks (Perfettini & Avouac, 2004), and increasing as a logarithmic function (Perfettini et al., 2010). The third process is the viscous flow in the lower crust and the upper mantle that can last several years [START_REF] Wang | Deformation cycles of subduction earthquakes in a viscoelastic Earth[END_REF][START_REF] Hu | Stress-driven relaxation of heterogeneous upper mantle and time-dependent afterslip following the 2011 Tohoku earthquake[END_REF][START_REF] Klein | Afterslip and viscoelastic relaxation model inferred from the large-scale post-seismic deformation following the 2010 Mw 8.8 Maule earthquake (Chile)[END_REF].

After the earthquake-induced stresses are released, the fault returns to its steady state and is again loaded by the plate convergence motion. This marks the beginning of a new interseismic phase and a new seismic cycle (see Figures 1.1 and 1.2a).

Whereas the vision of the seismic cycle seems simple at first glance, the increasing number of observations with increasingly better precision -together with new techniques to model them -have revealed a much more complex reality. In the description of the seismic cycle given above, it is necessary to consider that the stresses are accumulating heterogeneously along the fault, as a result of its spatially variable frictional properties. These frictional properties depend on the nature of the rocks, the presence of fluids, the depth (or equivalently the temperature and pressure), the topography of the subducted plate (seamounts, ridges...) or the faults being located either within the overriding plate or the subducted plate. The common friction law used to describe the framework of the faults is the rate-and-state friction law [START_REF] Dieterich | Time-dependent friction and the mechanics of stick-slip[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF][START_REF] Scholz | Earthquakes and friction laws[END_REF]. Two regimes emerge from this model, each associated with specific behaviors. The first one is called velocity-strengthening behavior: the friction increases with the slip rate and therefore the fault will resist even more to motion as the slip rate increases. According to [START_REF] Scholz | Earthquakes and friction laws[END_REF], this behavior allows for a steady continuous slip, corresponding to creep, i.e. aseismic slip taking place at the velocity of the convergence motion. The second one is velocity-weakening: the resistance to slip fault decreases as the slip rate increases, allowing for an acceleration of the slip rate. This ultimately leads to an instability: stress builds until sliding starts and gets released in an earthquake [START_REF] Scholz | Earthquakes and friction laws[END_REF]. A commonly accepted view is that the fault plane is predominantly velocity-strengthening but asperities of velocity-weakening patches are spatially distributed along the fault with various sizes (see Figure 1.2, top diagram) [START_REF] Lay | The surge of great earthquakes from 2004 to 2014[END_REF]Avouac, 2015;[START_REF] Bürgmann | The geophysics, geology and mechanics of slow fault slip[END_REF]. When these patches are small and far from each other, the slip remains aseismic (Avouac, 2015). But there also exist large velocity-weakening patches, which are mostly concentrated in what is called the seismogenic zone and which are usually embedded between two stable (velocity-strengthening) regions updip and downdip. The location of the transition between velocity-strengthening and velocity-weakening regions depends on the phase transitions of the minerals, thus on the temperature and pressure conditions at the fault [START_REF] Hyndman | The seismogenic zone of subduction thrust faults[END_REF]. Typically, for a continental strike-slip fault, the seismogenic zone is located between 3 and 20 km depth, whereas for a subduction fault, it is located between 15 and 45 km depth [START_REF] Scholz | Earthquakes and friction laws[END_REF][START_REF] Lay | The surge of great earthquakes from 2004 to 2014[END_REF].

To get an idea of the spatial distribution of the velocity-strengthening and velocityweakening behaviors along the fault, we use the notion of interseismic coupling. This concept was introduced by Hiroo [START_REF] Kanamori | Great earthquakes at island arcs and the lithosphere[END_REF] to explain the different sizes of earthquakes along subduction zones. The interseismic coupling coefficient Φ of the fault interface is defined as:

(

1 -Φ)V convergence = V interseismic
where V convergence is the plate convergence velocity and V interseismic is the interseismic slip rate [START_REF] Hyndman | The seismogenic zone of subduction thrust faults[END_REF][START_REF] Mccaffrey | Crustal block rotations and plate coupling. Plate Boundary Zones[END_REF]. When the coupling is equal to 1, the fault is fully locked and accumulates stress, while when the coupling is 0, the fault is unlocked and does not accumulate stress. Between these two values, the coupling is intermediate. Intermediate coupling could be explained by a spatially heterogeneous coupling at the fault. The presence of geometric complexities along the fault, such as seamounts, implies an abrupt change of the coupling at the fault, resulting in an intermediate coupling seen from the Earth's surface. Another explanation could be that the interseismic coupling is not stationary through time [START_REF] Bruhat | Deformation rates in northern Cascadia consistent with slow updip propagation of deep interseismic creep[END_REF]. Over several seismic cycles, [START_REF] Taylor | Rapid forearc uplift and subsidence caused by impinging bathymetric features: Examples from the New Hebrides and Solomon arcs[END_REF] and [START_REF] Thirumalai | Variable Holocene deformation above a shallow subduction zone extremely close to the trench[END_REF] have shown that subducted geometric complexities can change the interseismic strain in time. An alternative explanation could be that areas highly coupled, perhaps experience transient aseismic slip from time to time, resulting in an intermediate coupling when average over time. However, interseismic coupling models of subduction megathrusts enable to identify potentially seismogenic regions, i.e. fully locked areas. By estimating the locations of the earthquake epicenters along subduction zones, Tichelaar andRuff (1991, 1993) produced first estimations of the interseismic coupling. The continuous displacements recorded over several years by Global Navigation Satellite System (GNSS) sites were then used to estimate the interseismic coupling [START_REF] Mazzotti | Full interseismic locking of the Nankai and Japan-west Kurile subduction zones: An analysis of uniform elastic strain accumulation in Japan constrained by permanent GPS[END_REF][START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF]Metois et al., 2016).

Slow slip events, and more generally aseismic slips, have been documented during the interseismic phase (e.g. [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event[END_REF]Bartlow et al., 2014;Michel et al., 2019a), the proposed preseismic phase (e.g. [START_REF] Ruegg | Co-, post-and pre (?)-seismic displacements associated with the Mw 8.4 Southern Peru earthquake of 23 June 2001 from continuous GPS measurements[END_REF][START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF][START_REF] Socquet | An 8 month slow slip event triggers progressive nucleation of the 2014 Chile megathrust[END_REF] and the postseismic phase (e.g. Bedford et al., 2013;Hobbs et al., 2017) of the seismic cycle. Aseismic slips thus cover a wide spatio-temporal spectrum and release a significant part of the fault stress. Documenting them is therefore essential for a better assessment of the faults seismic hazard. This is also an opportunity to better document the evolution of fault slip and the relationship between the aseismic and seismic processes.

The different types of aseismic slip

To summarize, there are different types of aseismic slip taking place at faults: creep, slow slip event and afterslip.

Steady creep. The fault is unlocked and is slipping at a steady velocity equal to the plate convergence motion. Slow slip events. Slow slip events were discovered in 1999 along the Nankai subduction zone in Japan using GNSS data [START_REF] Hirose | A slow thrust slip event following the two 1996 Hyuganada earthquakes beneath the Bungo channel, southwest Japan[END_REF]. Shortly after, they were also spotted along the Guerrero subduction zone in Mexico [START_REF] Lowry | Transient fault slip in Guerrero, southern Mexico[END_REF] and along the Cascadia subduction zone [START_REF] Dragert | A silent slip event on the deeper Cascadia subduction interface[END_REF]. Compared to regular earthquakes, slow slip events are aseismic because they have a low slip velocity, preventing seismic waves to be generated. They are relatively long transient events (i.e. from hours to years) and include an acceleration phase and a deceleration phase, as for regular earthquakes. Slow slip events release a portion of the stress accumulated. Although the slip by itself is aseismic, different seismic signals often accompany slow slip events and can thus be used as proxies to identify them. Among these signals, one can mention an increase of the microseismicity, seismic swarms (i.e. increase of the seismicity rate, without a foreshockmainshock-aftershock sequence), repeating earthquakes (i.e. earthquakes with almost identical waveforms) and tremors (i.e. non-impulsive and emerging earthquakes, with a specific waveform) [START_REF] Bürgmann | The geophysics, geology and mechanics of slow fault slip[END_REF]. The tremors consist of numerous low-frequency earthquakes (short duration < 1 s and M w < 3) and sometimes, very-low-frequency earthquakes (long duration 10-200 s and M w of 3-4) [START_REF] Beroza | Slow earthquakes and nonvolcanic tremor[END_REF]. The scaling law that slow slip events follow is still being debated. Ide et al. (2007) suggest that they follow a different scaling law than earthquakes, where their moment proportional to their duration. Michel et al. (2019b) recently propose that they follow the same scaling law as earthquakes.

Afterslip. Aseismic slip can also occur after an earthquake, being this time referred to as afterslip. It corresponds to the response of the stress changes induced by the earthquake in its surroundings [START_REF] Hsu | Rapid afterslip following the 1999 Chi-Chi, Taiwan earthquake[END_REF]. Therefore, the afterslip distribution often complements the coseismic slip distribution and is assumed to occur preferentially in areas where the interseismic coupling is low. However, a small portion of the afterslip, and the aftershocks, have been observed within the coseismic rupture in a few cases (Bedford et al., 2013;Tsang et al., 2019). During the weeks that follow an earthquake, the afterslip rate decreases very quickly, as a function of the inverse of time (Perfettini et al., 2010). Afterslip and aftershocks have the same spatial and temporal evolution, and afterslip is thought to drive the occurrence of the aftershocks (Perfettini & Avouac, 2004;Hsu et al., 2006). Afterslip usually releases a moment equivalent to 25 % of the coseismic moment [START_REF] Pritchard | An aseismic slip pulse in northern Chile and alongstrike variations in seismogenic behavior[END_REF], but some specific cases showed that it can release a moment exceeding the coseismic moment (Villegas-Lanza, Nocquet, et al., 2016). The seismic moment released by the aftershocks is small compared to that of the afterslip [START_REF] Hsu | Rapid afterslip following the 1999 Chi-Chi, Taiwan earthquake[END_REF]Perfettini & Avouac, 2004). Nevertheless, a spatio-temporal study of very early afterslipi.e. the first 12 hours -showed that afterslip can be either seismic or aseismic (Twardzik, Vergnolle, Sladen, & Tsang, 2021).

Preslip. As we discussed earlier, preslip is aseismic slip during the preseismic phase. It is supposedly an aseismic slip that gradually accelerates over time. However, it seems to often take the form of a slow slip event (i.e., with an accelerating and decelerating phase). Slow slip events have been documented prior to several large megathrust earthquakes, and have been proposed to be part of a nucleation phase (e.g. [START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake[END_REF]Radiguet et al., 2016;[START_REF] Socquet | An 8 month slow slip event triggers progressive nucleation of the 2014 Chile megathrust[END_REF]. It has also been proposed that these slow slip events add periodic stress increments in their surrounding, which consequently modulate the occurrence of moderate-sized earthquakes [START_REF] Uchida | Periodic slow slip triggers megathrust zone earthquakes in northeastern Japan[END_REF].

In my thesis, I will focus on slow slip events and afterslip in subduction zones. In the next section, I introduce the general questions and explain why these questions are fundamental.

Motivations and objectives of this thesis

The purpose of my thesis is to provide some answers to two questions.

The first question concerns the relationship of the different types of aseismic slip with seismic processes. Afterslip occurs during the postseismic phase, while slow slip events occur during the interseismic phase (e.g. Michel et al., 2019a). The advantage of the afterslip is that we know where and when it takes place. It also has a large signal, allowing precise models to be derived. If afterslip general evolution is already rather well documented today, very few studies have investigated the possibility of short term modulations of slip superimposed to the overall evolution. A detailed study of the afterslip for earthquakes of various magnitudes would allow a better understanding of its dynamics and its relationship with the interseismic and coseismic phases. It will also allow us to assess if slow slip events can occur at the same time as the afterslip. Similarly, comparing the evolution of the afterslip -and potentially slow slip events -to the postseismic seismicity (i.e. large aftershocks and microseismicity) in space and time could inform us about possible interaction between seismic ruptures and aseismic slips. If so, the interactions would reflect in some way differences of the frictional properties along the fault. Systematically imaging the evolution of aseismic slip and seismic ruptures could thus provide some insights on the mapping of the velocity-weakening asperities, which host small to great earthquakes, and therefore could contribute to a better understanding of the faults' behavior.

The second question focuses on the importance of small slow slip events. Moderate to large slow slip events have been documented along several subduction zones (e.g. [START_REF] Hirose | A slow thrust slip event following the two 1996 Hyuganada earthquakes beneath the Bungo channel, southwest Japan[END_REF][START_REF] Dragert | A silent slip event on the deeper Cascadia subduction interface[END_REF][START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event[END_REF]Bartlow et al., 2011Bartlow et al., , 2014;;[START_REF] Rousset | A geodetic matched filter search for slow slip with application to the Mexico subduction zone[END_REF]Klein, Duputel, et al., 2018;Michel et al., 2019a). However, if slow slip events follow the [START_REF] Gutenberg | Seismicity of the world and associated phenomena[END_REF] law, which states that the frequency of earthquakes follows a power law as a function of their magnitude, small slow slip events should be more frequent than large ones. Many seismic sequences (i.e. seismic swarms or tremors), which generally go along with slow slip events, have been documented. However, these sequences are not always associated with any documented aseismic process. This lack comes from the difficulty to detect small slow slip events (as they generate surface displacement close to the GNSS data noise level) and to model them. Different techniques have been proposed to overcome this problem, using only GNSS data [START_REF] Rousset | A geodetic matched filter search for slow slip with application to the Mexico subduction zone[END_REF], or using specific seismic signals to target signals close to the noise level in the GNSS data. [START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes[END_REF]; [START_REF] Frank | Slow slip hidden in the noise: The intermittence of tectonic release[END_REF] used repeating earthquakes and very-low-frequency earthquakes respectively to target slow slip events in GNSS data. Yet, it is still difficult to build exhaustive catalogues of slow slip events, partly because there is no automatic approach with enough low numerical cost to allow a systematic search over large areas and large periods of time. Exhaustive catalogues would enable us to know if slow slip events follow the same scaling law as regular earthquakes (Michel et al., 2019b). They would also allow us to revisit our vision of interseismic coupling models taking into consideration the moment released by the slow slip events during the interseismic phase.

To address these two main questions, this thesis has two objectives. First, to image the spatio-temporal evolution of aseismic and seismic processes to study their relationships and dynamics. Second, to create a tool for semi-automatic detection of small slow slip events using geodetic data.

How to meet these challenges?

Seismological data provide precise imaging of seismic processes in time and space and enable to target aseismic processes through their seismic signature. Geodetic observations, on the other hand, enable us to image both aseismic and seismic (when large enough) processes. In particular, GNSS data highlighted the different aseismic slip processes. Although some other geodetic data, such as Interferometric Synthetic Aperture Radar (InSAR), can document aseismic processes [START_REF] Dalaison | The interplay between seismic and aseismic slip along the Chaman fault illuminated by InSAR[END_REF], the noise and resolution of SAR data do not allow the detection of small events when not close to the surface. It is worth noting that recently, using deep learning, millimetric surface displacements induced by small aseismic events have been highlighted with SAR data [START_REF] Rouet-Leduc | Autonomous extraction of millimeter-scale deformation in InSAR time series using deep learning[END_REF]. A dense GNSS network, however, provides sufficient resolution to detect small events and to document their temporal evolution. In this thesis, we therefore used GNSS data. But which modelling approach of GNSS data is the most adequate to model aseismic processes and study their dynamics and relationships with other processes? GNSS data modelling is often conducted today through static inversions. Static inversions have been used to study the coseismic slip model of an earthquake (e.g. Freymueller et al., 1994;[START_REF] Vigny | The 2010 Mw 8.8 Maule megathrust earthquake of central Chile, monitored by[END_REF][START_REF] Ozawa | Coseismic and postseismic slip of the 2011 magnitude-9 Tohoku-Oki earthquake[END_REF], its cumulative afterslip (e.g. [START_REF] Evans | Geodetic imaging of coseismic slip and postseismic afterslip: Sparsity promoting methods applied to the great Tohoku earthquake[END_REF][START_REF] Gualandi | Pre-and post-seismic deformation related to the 2015, Mw7.8 Gorkha earthquake, Nepal[END_REF] and the average interseismic locking of a subduction megathrust assumed to be steady in time (e.g. Metois et al., 2016;[START_REF] Graham | A global set of subduction zone earthquake scenarios and recurrence intervals inferred from geodetically constrained block models of interseismic coupling distributions[END_REF]. The advantage of this approach is that it requires little computing resources. To model a time evolution of slip, successive static inversions at selected dates are often used. The disadvantage of this approach is that it does not entirely describe the evolution of slip over time, preventing for instance a study of its dynamics and a detailed comparison with the seismicity rate through time. Yet, the comparison with the seismicity is particularly helpful to better characterize the small slow slip events that are at the limit of data resolution [START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes[END_REF][START_REF] Frank | Slow slip hidden in the noise: The intermittence of tectonic release[END_REF] and to investigate the relationship between aseismic slips and seismic ruptures. Moreover, the amount of available geodetic data covers large periods of time, which offers a great diversity of aseismic and seismic processes to study. Building a dynamic view of the slip is one of today's challenge in neotectonics. Kinematic inversions (also called time-dependent inversions) model the daily evolution of displacement, providing a slip model for each day. This approach is more robust than the static one since the inverted slip must fit the entire GNSS displacements history. The disadvantages are that it is often based on several assumptions and it requires more computing resources than the static inversions. However, the obtained kinematic model allows a direct comparison with the seismic processes (e.g. Bartlow et al., 2011;Rolandone et al., 2018), and a fine analysis of the slip dynamics. In this thesis, I used kinematic inversions of GNSS data to model aseismic slip along faults.

We are now focusing on the most suitable subduction zone. Beyond the numerous subduction zones in the world, the South American subduction zone appears to be among the best candidate to study aseismic and seismic processes. It experiences a great diversity of processes, has a heterogeneous coupling and it is monitored with several GNSS networks. Figure 1.3 gathers the large earthquake ruptures, the seismic swarms and the slow slip events that have taken place in Ecuador, Peru and Chile since 1995.

The Ecuadorian fault interface (Figure 1.3) is heterogeneously coupled [START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF][START_REF] Nocquet | Motion of continental slivers and creeping subduction in the northern Andes[END_REF] and is known to experience both slow slip events and seismic swarms. A deep slow slip event was detected at 60-80 km depth (Rolandone et al., 2018). South of the Pedernales region, the La Plata area experiences regular slow slip events and seismic swarms [START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF][START_REF] Segovia | Seismicity distribution near a subducting seamount in the Central Ecuadorian subduction zone, space-time relation to a slow-slip event[END_REF]. Moreover, slow slip areas were identified as barriers to the 2016 M w 7.8 Pedernales earthquake propagation (Nocquet et al., 2017;Vaca et al., 2018). Finally, a recent study of Chalumeau et al. (2021) also showed that afterslip stress perturbations trigger the occurrence of aftershocks.

Along the Peruvian fault interface (Figure 1.3), a special sequence of seismic swarm and aseismic slip occurred in 2009 (Villegas-Lanza, Nocquet, et al., 2016). In central Peru, at 13 • N-14 • N, the M w 8 Pisco earthquake ruptured in 2009 [START_REF] Sladen | Source model of the 2007 Mw 8.0 Pisco, Peru earthquake: Implications for seismogenic behavior of subduction megathrusts[END_REF]. Further south, at 16 • N, the 2001 M w 8.4 Arequipa earthquake ruptured at shallow depth [START_REF] Ruegg | Co-, post-and pre (?)-seismic displacements associated with the Mw 8.4 Southern Peru earthquake of 23 June 2001 from continuous GPS measurements[END_REF][START_REF] Perfettini | Geodetic displacements and aftershocks following the 2001 Mw= 8.4 Peru earthquake: Implications for the mechanics of the earthquake cycle along subduction zones[END_REF]. A transient event was documented and proposed as the preseismic phase of this earthquake [START_REF] Ruegg | Co-, post-and pre (?)-seismic displacements associated with the Mw 8.4 Southern Peru earthquake of 23 June 2001 from continuous GPS measurements[END_REF][START_REF] Melbourne | Precursory transient slip during the 2001 Mw=8.4 Peru earthquake sequence from continuous GPS[END_REF]. Along the Peruvian subduction interface, the interseismic coupling is also heterogeneous (Villegas-Lanza, Chlieh, et al., 2016).

Along the Chilean subduction zone (Figure 1.3), the 3500 km subduction interface is known to be a heterogeneous interface in terms of interseismic coupling (Metois et al., 2016). Seismic swarms and slow slip events prior to the 2014 M w 8.1 Iquique earthquake [START_REF] Duputel | The Iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty[END_REF] were documented, overlapping or located in the vicinity of the coseismic rupture [START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake[END_REF][START_REF] Socquet | An 8 month slow slip event triggers progressive nucleation of the 2014 Chile megathrust[END_REF]. At latitudes 23 • S-24 • S, a M w 7.8 ruptured at 30-40 km depth in 2007 in the Tocopilla area [START_REF] Peyrat | Kinematic rupture process of the 2007 Tocopilla earthquake and its main aftershocks from teleseismic and strong-motion data[END_REF]. South of this earthquake, in 1995, the M w 8 Antofagasta earthquake ruptured the megathrust at 30-55 km depth [START_REF] Pritchard | Co-seismic slip from the 1995 July 30 Mw=8.1 Antofagasta, Chile, earthquake as constrained by InSAR and GPS observations[END_REF]. A deep slow slip event and a seismic swarm were documented near Copiapó, along the Chilean subduction zone (Klein, Duputel, et al., 2018). A seismic and aseismic sequence occurred in the same area at shallow depth (i.e. 10-30 km) in 2020 (Klein et al., 2021). In 2015, at latitudes 30 • S-31.5 • S, the M w 8. [START_REF]The M w 8.3 2015 Illapel afterslip spatio-temporal modelling[END_REF] Illapel earthquake ruptured the megathrust from 40 km depth to the trench (Melgar et al., 2016;[START_REF] Grandin | Three-dimensional displacement field of the 2015 Mw8.3 Illapel earthquake (Chile) from across-and along-track Sentinel-1 TOPS interferometry[END_REF]Ruiz et al., 2016;Klein et al., 2017). Seismic swarms prior to this earthquake were documented downdip and south of the coseismic rupture (Poli et al., 2017). In 2017, a M w 6.7 earthquake occurred in Valparaíso, whose nucleation phase began with a slow slip event (Ruiz et al., 2017). The largest earthquake since 1960 was the 2010 M w 8.8 Maule earthquake. It ruptured a megathrust segment of approximately 400 km length [START_REF] Delouis | Slip distribution of the February 27, 2010 Mw=8.8 Maule earthquake, central Chile, from static and high-rate GPS, InSAR, and broadband teleseismic data[END_REF][START_REF] Moreno | Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone[END_REF]. This brief summary of the aseismic and seismic processes taking place along the South American subduction zone shows its rich diversity. This is the study area of this thesis.

This manuscript is divided into four Chapters. The first one introduces the data, the different kinematic inversion methods published to date together with different case studies and a detailed presentation of the method used in this thesis. The second Chapter presents a first case study, the postseismic phase of the 2015 M w 8. [START_REF]The M w 8.3 2015 Illapel afterslip spatio-temporal modelling[END_REF] Illapel earthquake (Chile) which was accepted in the Journal of Geophysical Research: Solid Earth in December 2022. The third Chapter presents a second case study, the Atacama sequence (Chile) in September 2020. The fourth Chapter introduces a semi-automatic detection tool of aseismic slip and its application on the Ecuadorian subduction zone. Finally, in the fifth and last Chapter, I summarize the elements brought by the different case studies of this thesis, before discussing some perspectives offered by this work. The fault has velocity-weakening asperities embedded in an overall velocity-strengthening area. a During the interseismic phase, the stress builds up on the velocity-weakening asperities. Slow slip events occur at the periphery of velocity-weakening asperities together with seismic swarms, or downdip of the seismogenic zone with tremors. b During the preseismic phase of an earthquake, the microseismicty increases and foreshocks occur. Slow slip events were also proposed during this phase. c During the coseismic phase, an earthquake release the stresses accumulated by the fault through seismic slip at velocity-weakening asperities. d During the postseismic phase, afterslip develops in the periphery of the earthquake rupture, and microseismicity mimics its evolution. Large aftershocks also occur.

Chapter 2. Data & method

Imaging the spatial and temporal evolution of aseismic and seismic processes along faults, and better detecting small aseismic signals are the main objectives of this thesis.

To do so, we use GNSS data and a kinematic inversion approach to model slip along faults. In the first section, I briefly introduce the Global Navigation Satellite Systems, the high-precision GNSS approach and an example of post-processing of time series. In the second section, I summarize existing methods for kinematic slip inversion and their application to different objects. Finally, in the last section of this Chapter, I develop the time-dependent inversion approach used in this thesis.

Data & processing

During the last decades, our comprehension of earthquakes and of the seismic cycle has been revolutionized by space geodesy. Monitoring of faults across the world with several space geodesy techniques has enabled recording of the deformation associated with the seismic cycle with millimeter precision. Among these space geodesy techniques, one can mention the Interferometric Synthetic Aperture Radar (InSAR) (e.g. [START_REF] Massonnet | The displacement field of the Landers earthquake mapped by radar interferometry[END_REF][START_REF] Grandin | Three-dimensional displacement field of the 2015 Mw8.3 Illapel earthquake (Chile) from across-and along-track Sentinel-1 TOPS interferometry[END_REF], optical imagery (e.g. [START_REF] Klinger | High-resolution satellite imagery mapping of the surface rupture and slip distribution of the Mw 7.8, 14 November 2001 Kokoxili earthquake, Kunlun fault, northern Tibet, China[END_REF][START_REF] Lauer | Fault geometry and slip distribution of the 2013 Mw 7.7 Balochistan earthquake from inversions of SAR and optical data[END_REF] and the various Global Navigation Satellite Systems (GNSS) (e.g. [START_REF] Bilham | GPS measurements of present-day convergence across the Nepal Himalaya[END_REF][START_REF] Mazzotti | Current tectonics of northern Cascadia from a decade of GPS measurements[END_REF]. InSAR is a method to measure Earth's surface changes, using radar signals. It consists of determining the variation in distance between an antenna carried by a satellite and a certain area of the Earth. By doing this before and after an event, seismic or aseismic, we can measure the deformation it has generated. Optical imagery, for its part, allows to measure the deformation of an event with satellite images.

On a pair of images, i.e. one acquired before the event (master scene) and one acquired after (slave scene), the pixels are cross-correlated. The distance between the position of a pixel in the master scene and the position of this same pixel in the slave scene represents the displacement of this pixel between the two satellite images, i.e. before and after the event. By doing this for all the pixels, we obtain the deformation related to the event. InSAR and optical imagery provide continuous spatial measurements of the Earth's surface. In contrast, GNSS provide discrete measurements, but with a better precision and time resolution. This latter technique is the one I used in my thesis, and is described in this first section.

Global Navigation Satellite System (GNSS): a brief overview

The Global Navigation Satellite System (GNSS) groups different satellite constellations that provide positioning for real-time navigation. The first GNSS was launched in 1973 by the United States of America and is called the Global Positioning System (GPS). It consists of a constellation of 24 operational satellites organized into groups of four, each group orbiting on 6 evenly distributed plans all with an inclination of 55 • . Later on, several countries or unions of countries launched their own satellite constellation, each of them providing either global or regional coverage. There are four global GNSSs including GPS. China is operating BeiDou (BDS, also called Compass) since 2000. It consists of a constellation of 35 satellites and is available worldwide since 2020. The European Union operates Galileo, a 24-satellite constellation completed in 2021. This constellation provides a better precision at high latitudes compared to other GNSSs. The Russian federation is operating since 1995 a constellation of 24 satellites called the Globalnaya Navigazionnaya Sputnikovaya Sistema (GLONASS). It was modernized between 2002 and 2011.

There are also GNSSs that provide regional coverage. The Government of India operates the Navigation Indian Constellation (NavIC, previously called Indian Regional Navigation Satellite System), a 8-satellite constellation completed in 2018. NavIC covers India and its surroundings. Japan is operating a 4-satellite constellation since 2018 called the Quasi-Zenith Satellite System (QZSS). It improves the coverage of East Asia and Oceania.

While monitoring geophysical objects, such as faults, landslides or glaciers, the GNSS antenna and receiver may detect some of all the GNSS satellites. Here, however, only GPS data has been used to produce the solutions presented in this work.

From GNSS raw data to time series

In this section, although I did not process the GNSS data myself, I provide a brief introduction to GPS data and a description of the data processing are given. In the first section, the concept of the Global Positioning is explained in details. In the second section, the steps required to obtain millimetric position time series from GNSS raw data is developed. The third section sets out the processing required to clean the time series and to extract the tectonic signal using the example of a deep slow slip event in Ecuador.

Satellite Global Positioning principle

Positioning by trilateration provides the position of the receiver with a meter precision using the travel time of electromagnetic waves -propagating through the atmosphere at the speed of light -between the satellite and the receiver. The unknown parameters here are the coordinates of the receiver: X r , Y r and Z r . The electromagnetic wave emitted by the satellite is coded with a starting time of emission t s and the satellite position. With this information and the reception time of the signal at the receiver t r , the receiver should be located on a sphere centered on the satellite with a radius equal to:

d rs = c(t r -t s ) (2.1)
with c the speed of light. Repeating this with a second satellite constrains the receiver position to be on the intersection with a second sphere. A third satellite is used to obtain three spheres that intersect at one point. We thus have three equations 2.1 to find the three unknowns (X r , Y r and Z r ). Figure 2.1 illustrates this approach.

Whereas the satellite is equipped with a very precise atomic clock, the receiver is not. This clock error is very large with respect to the propagation velocity of the wave, i.e. the speed of light. Thus, the position needs to be corrected from the clock error of the receiver called ϵ r . This is done using a fourth satellite.

The position must also be corrected from some delays of the wave propagation across the atmosphere, in particular across the layer called the ionosphere. The ionosphere is the ionized layer of the atmosphere and is a dispersive medium. Indeed, while crossing the ionosphere, the electromagnetic wave interacts with charged particles inducing a propagation delay proportional to the frequency. This delay can be corrected by combining two GPS frequencies. Let ϵ r be the clock error, ϵ i the ionosphere error and ϵ s the receiver state correction. Equation 2.1 is now: where p rs is the pseudo-range, i.e. the pseudo-distance (distance + clock error + delays) between the satellite and the receiver.

p rs = c(t r -t s ) + ϵ r + ϵ i (2.2) GNSS site X r , Y r , Z r X s1 , Y s1 , Z s1 X s2 , Y s2 , Z s2 X s3 , Y s3
This trilateration approach enables a real-time positioning with a precision of a few meters.

High precision GNSS

It is possible to obtain a higher precision than what can be achieved using trilateration -about a few millimeters. To achieve this, researchers have developed another approach. Rather than using the code of the signal, they use the phase of its carrier wave. If this alternative approach provides very precise positioning, it is at the expense of the loss of some information. Indeed, the carrier wave does not include information on the absolute distance satellite -receiver, but very accurately measures how that distance changes over time.

With this approach, at the first phase measurement, the number of full cycles between the satellite and the receiver is unknown. This is called the phase ambiguity. For high precision GNSS, and unlike the positioning by trilateration, the satellite clock error is important and must be corrected. To address these problems, many scientific softwares use an approach called the double-difference. The first simple difference consists in using two receivers which look at one satellite, to compute and remove the satellite clock error. The second simple-difference uses the same receivers but a different satellite to estimate and remove the receiver clock errors. The double-difference principle is shown on Figure 2.2. With the double-difference, we loose the absolute positioning and get a relative position of the receivers.

The carrier wave is also affected by several errors. After the double-difference, the positions are corrected for ionospheric and tropospheric delays. Because the measurement lasts approximately 24 hours, the positions should also be corrected from the Earth's tides and the Earth's response to oceanic tides. Finally, the apparent change of the antenna center position is also corrected (called the Phase Center Variation, PCV).

The resulting positions have a millimetric relative position with respect to other sites. It is called a free network solution.

The last step is to compute the positions in an appropriate reference frame by applying a Helmert transformation to the free network solution. An a priori external reference frame is required, i.e. the International Terrestrial Reference Frame (ITRF). Because the Earth's surface changes with time, an accurate reference frame must be defined on a regular basis. The ITRF is a cartesian geocentric reference frame particularly appropriate for the study of the Earth. It provides the position and velocity of specific points at the Earth surface. It is updated every 5 years by the Institut National de l'Information Géographique et forestière (IGN) by combining several spatial geodesy techniques: GNSS, Satellite Laser Ranging (SLR, operated by the NASA), Very Long Baseline Interferometry (VLBI, also operated by the NASA) and Doppler Interferometry and Radiopositioning Integrated by Satellite (DORIS, operated by the CNES).

The data obtained with the Helmert transformation is a position in the ITRF for each day. In our case, we obtain a GNSS displacement time series for each antenna (i.e. each GNSS site) and for all three components.: N in the north direction; E in the east direction and U in the vertical direction. When studying a specific region, these time series are converted in a local reference frame to provide the amount of displacement in the north, east and up components. The horizontal components (N and E) have a precision of approximately 2-3 mm. The vertical component (U ), has a larger error due to the absence of satellites below the antenna. The vertical precision is approximately 4-8 mm.

The first two projects conducted during this thesis, i.e. Chapters 3 and 4 were done using the solution of Klein, Vigny, Nocquet, and Boulze (2022). The last project, presented in Chapter 5, was done using the Ecuador data set described in (Nocquet et al., 2014) and regularly updated.

Post-processing of GNSS time series: common mode and Vondrák filters

When chasing for small modulations of the slip or small transient events (approximately a few millimeters), filtering the GNSS time series is key to obtain a clear signal. In this section, I present how I used the common mode and the Vondrák filters to clean the signal for the study of a deep slow slip event which occurred in 2015 in Ecuador. The first step is a challenging one and consists in analysing by eye the time series to identify the signal. I use the Vondrák filter to smooth the time series with different cut-off frequencies. The Vondrák filter [START_REF] Vondrák | A contribution to the problem of smoothing observational data[END_REF][START_REF] Vondrák | Problem of smoothing observational data II[END_REF] is a low pass filter, which consists in finding an equilibrium between a curve passing through all the observations and a smooth curve. The Vondrák filter is a general solution for time series with non-homogeneous uncertainties and is adequate for non-equidistant observations. This is particularly appropriate for our time series which often present several data gaps. The cut-off frequency can be adjusted depending on the target signal and its wavelength. A common mode error corresponds to small non-tectonic displacements in the time series common to all GNSS sites. The common mode filter is an empirical spatial technique developed by [START_REF] Wdowinski | Southern California permanent GPS geodetic array: Spatial filtering of daily positions for estimating coseismic and postseismic displacements induced by the 1992 Landers earthquake[END_REF]. A new automatic method was recently developed by [START_REF] Kreemer | Robust estimation of spatially varying common-mode components in GPS time-series[END_REF]. The common mode filter significantly reduces the highfrequency noise and the seasonal signal. This filter estimates the common displacement from several time series of a particular GNSS network and removes it. Ideally, the common mode should be computed on the entire time series, however, doing so brings a lot of complexity and the time series do not often span the same time period. The minimum requirement is that it should be computed on the period of the studied signal.

The method to compute the common mode is the following. First, having defined the GNSS network for the studied phenomenon, a subset of GNSS time series must be defined, from which the common mode will be estimated. This subset should be carefully selected regarding the location of the sites, the quality of the time series at the time of the event and must be uniformly distributed in space. For example, while studying a slow slip event with numerous sites located on the coast, a site located on a volcano is not relevant and should not be included in the subset. A low quality time series is not a good candidate for common mode either, but if during the event the quality is higher, it should be added to the sample. Second, outliers must be removed from the horizontal and vertical components of the time series. Third, the interseismic velocity is removed from the GNSS time series using the MIDAS method (Blewitt et al., 2016). The residual subset time series are stacked and averaged to obtain the common mode. This common mode is then removed from all the time series. Finally, in cases in which common mode solutions substantially vary as a function of the subset of sites used to compute the common mode, this variability is accounted for by increasing the uncertainties on the corresponding data in the kinematic inversion.

Let us take a practical example: a slow slip event occurred between 50 and 70 km depth at latitude 1 • S, in Ecuador between July and November 2015. It was first modelled by Rolandone et al. (2018). This event was studied at the beginning of my thesis and most of my work on this event consisted in clearing the time series. Figure 2.3a shows a selection of 23 raw GNSS sites. We can see about 8-10 mm westward displacement in August and September 2015 (red shaded area on Figure 2.3a). On the north component, the signal is unclear, although one should expect modulations of the displacements during the same period in association with the signal observed on the west component. To investigate this, after removing the interseismic velocity from the time series (following MIDAS method Blewitt et al. (2016)), we applied a Vondrák filter with a cut-off frequency f c = 10 shown in Figure 2.3b. The Vondrák filter helps to highlight the transient signal and to identify which sites recorded it. With this information, we can compute the common mode.

The reference subset of GNSS sites for the estimation of the common mode is composed of nine sites, surrounding homogeneously the signal area (shown with red triangles in Figure 2.4a). For this subset, I extracted the time series from 2014 to the day before the 2016 Pedernales earthquake that ruptured 100 km north of the transient signal. Then, I removed the interseismic velocity from the GNSS time series using the MIDAS method (Blewitt et al., 2016). Using a Vondrák filter, the outliers are also removed from the horizontal and vertical components. Figure 2.4b shows the resulting time series, from which we estimate the common mode. The time series show high-frequency noise, often correlated in time by several sites. They do not show any transient signal at the presumed occurrence of the transient event (red shaded area in Figure 2.4b). The estimated common mode is shown on Figure 2.5a. It is a pseudo-sinusoidal signal, close to 0. Finally, the common mode is removed from the time series shown in Figure 2.3a. The filtered time series are shown in Figure 2.5b. They are free from high frequency noise on the horizontal components and the transient signal is visible on the north component, while it was unclear on the detrended time series (Figure 2.3a). For a few time series, the transient signal reaches about 8 mm southward over more than a month. On the east component, the signal is also about 8 mm for the same period. The vertical component does not show any clear signal.

Considerable improvements on highlighting the transient signal have been achieved using the common mode filter. The displacements recorded by the GNSS, however, are very small. Modelling this deep transient event, with a displacement rate of about 0.2 mm.day -1 , was difficult to do.Therefore, I got myself involved in another project. This study has provided me with the opportunity to become familiar with the data and to learn the different approaches to clean the time series. This was a key step before performing time-dependent inversion of GNSS time series, described in section 2.3.

Existing methods for kinematic slip inversion from GNSS time series

Static inversions of GNSS data in neotectonics have significantly improved our knowledge and understanding of the seismic cycle. In particular, different slip modes, i.e. slow slip events, aseismic slip and afterslip following an earthquake have been well documented. With the improvement of GNSS precision and the increasing density of GNSS sites, the new challenge today is to focus on the smaller transient phenomena recorded by GNSS sites. If numerous slow slip events have been documented in several subduction zones, there must be more small slow slip events to be discovered. In this respect, time-dependent inversions of displacements in GNSS time series are valuable tools to model them, because they are difficult to see in the GNSS time series. These inversions also have the potential to model the spatial and temporal interactions between the different slip modes at the subduction megathrust during the seismic cycle. A few solutions for time-dependent inversion of GNSS displacements have been proposed, that I summarize in the following sections together with some applications to different transient events.

Principal Component Analysis Inversion Method (PCAIM)

and Variational Bayesian Independent Component Analysis Inversion Method (vbICAIM) approaches

The displacements seen in the GNSS time series find their origin in different physical processes. When studying phenomena that induce small displacements, it is crucial to identify and extract the right portion of displacement induced by the studied phenomenon (e.g. slow slip event, postseismic deformation...). The analysis of displacements recorded by GNSS sites can be conducted with a multivariate statistical technique, such as the Principal Component Analysis Inversion Method (PCAIM) developed by Kositsky and Avouac (2010) and the variational Bayesian Independent Component Analysis Inversion Method (vbICAIM), developed by [START_REF] Gualandi | Blind source separation problem in GPS time series[END_REF]. In a first section, I describe the general principle of signal decomposition through the PCA and vbICA. In a second section, I present an application of the vbICAIM to the Cascadia subduction zone (Michel et al., 2019a).

General principle

The PCAIM (Kositsky & Avouac, 2010) and the vbICAIM [START_REF] Gualandi | Blind source separation problem in GPS time series[END_REF] are both based on a data decomposition as a sum of several components, which can be inverted separately. Each component represents a particular displacement field and is combined with a time function. Each source is assumed to be stationary in space. These approaches aim to filter the GNSS time series and extract the relevant part of the data that explain the studied phenomenon.

On one hand, the PCA enables to decompose the signal into several components which corresponds to a combination of various sources signal. Because the data can be explained by a small part of the components, we can remove most of the components and obtain filtered time series. The PCA is thus a great method to filter GNSS time series. However, while seeking to model a transient signal, the PCA is not suitable. It does not focus on the sources of the components and assumes that the data are distributed following a normal law. However, as shown by [START_REF] Gualandi | Blind source separation problem in GPS time series[END_REF], some signal probability density functions are not Gaussian, or not even unimodal (e.g. seasonal signal, logarithmic decay signal...).

On the other hand, the vbICA focuses on the identification of the component sources. This approach is well-suited for all sources, which is not the case for the ICA approach (it is optimal when the sources have a unimodal distribution, that transient signals do not have). The vbICA approach assumes that the various processes that explain the data are mutually independent. This statement is strong while, for example, several studies showed that transient signals and middle-sized earthquakes are induced by the afterslip of a large earthquake [START_REF] Voss | Slow slip events in the early part of the earthquake cycle[END_REF]Rolandone et al., 2018). However, it is very efficient to separate postseismic signal from the volcanic signal. It is particularly well adapted for the South American subduction zone: in the Andes, there are many volcanoes that "disrupt" the signal from the subduction. Finally, in terms of resolution, the vbICA provides an a posteriori matrix which helps to assess the precision of the signal decomposition.

After the decomposition of the data into several components, the displacement field associated to each component is inverted for slip along the fault. The slip distribution is multiplied by the time function associated with the inverted component of the displacement field. We obtain the cumulative slip at each time step. The main disadvantages of the PCAIM and the vbICAIM are that (1) we do not directly compute for the slip rate at each time step, (2) they require observations for every date during the studied period and (3) they do not use any a priori constraint on the slip evolution.

An application to the Cascadia subduction zone

The Juan de Fuca plate converges beneath the North American plate at a velocity of about 38 mm.yr -1 . Despite a moderate convergence velocity, the Cascadia subduction megathrust has not experienced a large earthquake since 1700. This raises questions about the slip modes that occur along this megathrust. The Cascadia subduction zone is thus a good candidate to assess the potential of kinematic inversion and fully exploit its advantages. Michel et al. (2019a) used the vbICAIM [START_REF] Gualandi | Blind source separation problem in GPS time series[END_REF] to extract the tectonic loading along this subdiction and the transients signals. They produced a slow slip event catalogue with GNSS data between 2007 and 2017. They found 64 slow slip events with magnitudes 5.3 ≤ M w ≤ 6.8, exhibiting a variability in duration (14 days ≤ duration ≤ 106 days) and in propagation velocities (5.5-11 km.day -1 ). The kinematic models highlight that the propagation direction of the slow slip events can be either towards the north or the south. They also showed that several events propagated in the two directions. By comparing the cumulative slip evolution of the events to the tremor evolution, they were able to show a systematic spatio-temporal correlation between the slow slip events and the tremors (see Figure 2.2.1.2b).

Finally, by comparing the spatial distribution of slip released by the slow slip events to the interseismic coupling, they showed a clear segmentation of the slip mode along the megathrust. As shown on Figure 2.2.1.2a, the locked areas are at shallow depths whereas slow slip events always occur -between 2007 and 2017 -at 35-50 km depth. Therefore, they suggest that the area between the locked zone and the tremors/slow slip events must be creeping at a long-term slip rate.

The Network Inversion Filter (NIF)

The Network Inversion Filter (NIF) is another approach of full time-dependent inversion of GNSS time series. Rather than decomposing the signal, the time series are directly inverted at each time step. The first solution of time-dependent inversion is the Network Inversion Filter (NIF), developed by [START_REF] Segall | Time dependent inversion of geodetic data[END_REF] and then modified by [START_REF] Miyazaki | Spatial and temporal evolution of stress and slip rate during the 2000 Tokai slow earthquake[END_REF]. In a first section, we briefly describe the NIF approach and subsequently present in a second section two applications on the Cascadia and the Hikurangi subduction zones.

Methodology

The Network Inversion Filter aims to image the spatio-temporal evolution of slip rate along faults. [START_REF] Segall | Time dependent inversion of geodetic data[END_REF] developed a time domain filtering approach for GNSS data which analyses all the time series simultaneously and corrects the spatial and temporal correlated errors from the measurements. This approach avoids adding strong constraints on the slip rate evolution. The NIF solves for the cumulative slip and the slip rate at each time step using a Kalman filter approach. By working with the complete GNSS network, the NIF adjusts the prediction to fit a slip linear combination of all the previous time steps and the actual one. Hence, the time-dependent inversion must not only fit the data each day, but the entire GNSS time series. The NIF can handle a nonnegativity constraint on the slip rate as well as data gaps in the time series. [START_REF] Miyazaki | Spatial and temporal evolution of stress and slip rate during the 2000 Tokai slow earthquake[END_REF] added a covariance regularization with both spatial and temporal smoothing.

An application to a slow slip event in the Cascadia subduction zone

The Cascadia subduction megathrust is well-known for experiencing numerous slow slip events (cf. section 2.2.1.2). An episodic tremor and slip (ETS) event occurred in August 2009 and lasted a month and a half. In order to investigate the physical relationship between tremors and slip, Bartlow et al. (2011) inverted the GNSS time series for the spatio-temporal evolution of slip rate along the megathrust using the NIF approach.

The slip rate model reveals a strong correlation in space and time with the tremors (see Figure 2.7). The slow slip event started at latitude 46.5 • N on August 3 while the tremors appeared on August 6 in the same area. Both the tremors and the slip propagated bilaterally, towards north and south until August 16. From this day, the southern part of the slip propagated southward to latitude 44 • N, just like the tremors. Also noteworthy is that the tremors are located where the slip rate is maximum, suggesting that they are triggered by the slow slip event.

With the kinematic model of the slip rate, Bartlow et al. (2011) were able to suggest that tremors result from local heterogeneities in the slipping area of the ETS. Therefore, the tremors could be used as a proxy for active slip regions along the Cascadia subduction zone.

Slow slip event dynamics at the Hikurangi subduction zone

Between 2009 and 2012, a complex sequence of slow slip events and quiescence occurred along the Hikurangi subduction zone. This subduction zone being very well monitored with permanent seismic and geodetic networks, a thorough study of the rela-tionship between aseismic slip and seismicity is possible. Bartlow et al. (2014) used the NIF method to run a two-year time-dependent inversion of 113 GNSS sites. They found 12 slow slip events with very different characteristics in terms of magnitude (5.8 ≤ M w ≤ 6.9), duration (7 days ≤ duration ≤ 270 days), slip (1.2 cm ≤ total slip ≤ 27 cm) and dimension. The whole sequence moment is equivalent to a M w 7.4 event. They observed that when the stress drop is high, the slow slip event duration is longer. A global spatio-temporal comparison of the slip rate evolution and the seismicity did not show a clear correlation, suggesting a heterogeneous subduction megathrust in terms of frictional properties. However, at a more local scale, using the kinematic slip rate models, they highlighted three special features.

First, the 2010-2011 Manawatu slip event, located on the southwest deep part of the subduction zone does not show the classic acceleration and deceleration phases but rather shows several phases of acceleration of different time periods. These acceleration phases are separated by periods of quiescence with reduced slip (see Figure 2.8). Second, the slow slip events that ruptured the eastern part of the megathrust in 2010-2011 also present some specificities. They occurred simultaneously in different areas, however, in between these active areas, no significant event occurred. It shows that slow slip events, as well as afterslip following a large earthquake, occur on preferential areas, suggesting that the subduction megathrust have an heterogeneous distribution of frictional properties. Third, for a short time period, a seismic swarm occurred on the western part of the 2011 Cape Turnagain event during its migration. This is the only region at the Hikurangi subduction zone where the seismicity correlates with a slow slip event.

The Hikurangi subduction zone shows that slow slip events can be very different one from another, even in the same region. Unlike many other subduction zones, i.e. Cascadia and Japan (Bartlow et al., 2011), no clear correlation was found between slow slip event and tremors. Indeed, along the Hikurangi subduction zone, the tremors are rather shifted from the slow slip events areas.

Other approach: a semi-kinematic inversion inspired from seismology

In this section, I present another approach which is not strictly a kinematic inversion method but provide a semi-kinematic view of slip along faults. I then describe the study of [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event[END_REF] on the Guerrero subduction zone.

In 2006, a large slow slip event, which lasted one year, occurred along the Guerrero subduction zone. This event has been studied several times with a static approach. For a thorough study of its dynamics, [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event[END_REF] developed a spatio-temporal slip modelling approach, based on a frequency domain approach generally used in seismology [START_REF] Cotton | Frequency domain inversion of strong motions: Application to the 1992 Landers earthquake[END_REF]. Their inversion method is a two-step approach. First, a classic static inversion of the cumulative displacements extracted from the GNSS time series is performed to obtain the total slip distribution along the fault. Here, no constraint is imposed on the slip. Second, a quasi-static inversion of the GNSS time series in the frequency domain is performed using the Green function of the first step. Here, the final slip amplitude is constrained by the inversion from the first step. This inversion method does not provide the daily slip rate evolution, but the annual slip rate is calculated. They find that the slow slip event initiated in February 2006 at 42 km depth on the western part of the megathrust. At the initiation time of this slow slip event, two moderate-sized earthquakes ruptured in the same megathrust segment. It spreads on a 300 km × 150 km area with a maximum slip of 15 cm, and is equivalent to a M w 7.5. As shown on Figure 2.9, the semi-kinematic modelling provided 50-days cumulative slip snapshots to study the evolution of the slow slip event. It highlights that the event propagated eastward, 250 km away from its initiation location and ended in January 2007. They were also able to estimate the velocity of this event, which propagated at 0.8 km.day -1 , a low value compared to other subduction slow slip events.

The rise time, i.e. the duration needed to reach the total slip at a specific point, is key information to understand the rupture process of the slow slip event. The inversion approach used here provides it. They found a very large value, i.e. 50% of the total duration of the event, which is different from regular earthquakes (rise time usually is 10-15% of the total duration for a M w 7-7.5 event, [START_REF] Heaton | Evidence for and implications of self-healing pulses of slip in earthquake rupture[END_REF]).

This semi-kinematic inversion method provided key information about this one-year slow slip event regarding its dynamics (e.g. velocity, rise time, propagation). Compared to other subduction zones, the Guerrero megathrust geometry seems to control the propagation of slow slip events.

Various approaches to kinematic or semi-kinematic inversions have been briefly pre-2.3. PYEQ software: time-dependent slip inversion from GNSS time series 61 sented. We have seen that they are based on various assumptions and that they allow different slip outputs. The following section presents the approach used in my thesis.

PYEQ software: time-dependent slip inversion from GNSS time series

PYEQ is a software developed by Jean-Mathieu Nocquet (IRD, Géoazur & IPGP) for GNSS data analysis and modelling. PYEQ is built upon the PYACS software. PY-ACS implements the core library for producing time series from free solutions (see section 2.1.2.2). It also includes modules for handling, filtering, visualizing time series and producing velocity fields.

PYEQ estimates the slip evolution at faults through time from GNSS time series. PYEQ models the slip evolution as piecewise linear functions. It uses a linear formulation and least-squares inversion with non-negativity constraints on the incremental slip. The linear formulation allows fast computation of the solution and ensures the existence and uniqueness of the solution (Nocquet, 2018). Compared to the other approaches previously introduced, PYEQ does not assume that the signal sources are spatially stationary, and the linear formulation allows to handle gaps in the GNSS time series. It also allows to introduce both spatial and temporal smoothing of the slip. Additional tools are provided to discretize faults, compute the transfer functions in a homogeneous elastic half-space for rectangular and triangular dislocations. Finally, it implements a few approaches to assess resolution of the results.

In this section, I detail the approach implemented in the PYEQ software, starting with the forward problem. I then describe the regularization strategies. Finally, I provide some information on practical aspects of the software.

Forward problem

The forward problem consists in relating the 3D displacements observed by the GNSS sites at the Earth's surface to the slip taking place at the fault. Let d j (t) be the displacements vector at the time step t for site j for the three components (north N , east E and vertical U components) and S(t) the cumulative slip vector. For simplicity, let us consider the case with only two subfaults with indexes i 1 and i 2 , and only one component of slip. The observation equation is:

  g j,N i 1 g j,N i 2 g j,E i 1 g j,E i 2 g j,U i 1 g j,U i 2   S i 1 (t) S i 2 (t) =   d j,N (t) d j,E (t) d j,U (t)   (2.3) where g j,N i 1
is the transfer function relating the unit slip at subfault i 1 to the displacement component N (respectively E and U ) at site j. The PYEQ approach solves for the slip rate at user-defined discrete time steps ṡ(t). We first consider that for all time series at the model origin time t 0 , d(t 0 ) = 0. In other words, all time series are referenced with respect to the first observation time included in the inversion. At the first step t 1 , introducing slip rate into equation 2.3:

∆t 1   g j,N i 1 g j,N i 2 g j,E i 1 g j,E i 2 g j,U i 1 g j,U i 2   ṡi 1 (t 1 ) ṡi 2 (t 1 ) =   d j,N (t 1 ) d j,E (t 1 ) d j,U (t 1 )   (2.4) or ∆t 1 G ṡ(t 1 ) = d(t 1 ) (2.5)
where ∆t 1 = t 1t 0 and ṡi (t 1 ) is the slip rate at subfault i assumed to be constant over the ∆t 1 time period. For the second time step t 2 , the observation equation is:

∆t 1 G ṡ(t 1 ) + ∆t 2 G ṡ(t 2 ) = d(t 2 ) (2.6)
By extending the above equation for all time steps, the linear system for the forward problem is:

       ∆t 1 G 0 0 • • • • • • 0 ∆t 1 G ∆t 2 G 0 • • • • • • 0 ∆t 1 G ∆t 2 G ∆t 3 G 0 • • • 0 . . . . . . . . . • • • . . . 0 ∆t 1 G ∆t 2 G ∆t 3 G • • • • • • ∆t m G               ṡ(t 1 ) ṡ(t 2 ) ṡ(t 3 ) . . . ṡ(t m )        =        d(t 1 ) d(t 2 ) d(t 3 ) . . . d(t m )        (2.7) which corresponds to: G ṡ = d (2.8)
Formulating the problem in terms of slip rate (rather than in terms of cumulative slip) allows to explicitly force the slip rate to be positive.

Regularization

Regularization constraints are required for two reasons. First, the inverse problem is under-determined. Indeed, the model parameters vector ṡ has a dimension m × n (m time steps and n subfaults), larger than the displacements vector d with a dimension m × 3 × p (m time steps, 3 displacement components and p GNSS sites). Second, regularization is needed to avoid unrealistic and non-physical results, i.e. too large slip and large differences in the amount of slip between spatially and temporally adjacent subfaults. Therefore, two regularization constraints are added: (1) a condition of damping and ( 2) a spatial and a temporal smoothing. Their implementation in PYEQ are detailed in this section.

Damping

Regularization constraints are added to the linear system 2.8 as pseudo-observations. To avoid large slip in the model, we use a condition of damping, whose pseudo-observation states that slip is null within a given weight provided by C -1 D ṡ : ṡ(t) = 0, which can be rewritten as

     1 0 • • • 0 0 1 • • • 0 . . . . . . . . . . . . 0 0 0 1      Identity matrix I      ṡ(t 1 ) ṡ(t 2 ) . . . ṡ(t m )      =      0 0 . . . 0      (2.9) 2.3.2.

Spatial smoothing with a covariance regularization

For the spatial smoothing of a static inversion, PYEQ uses the formulation of [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF], also used in [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event[END_REF]. The covariance regularization matrix is:

C i,j s = σ 2 exp - d i,j d c (2.10)
where σ is the user-defined constraint on the slip rate, d i,j is the distance between subfaults with indexes i and j, and d c a user-defined distance of correlation.

Although this approach has been extended to both spatial and temporal domains, the resulting formulation is tedious and numerically expensive. For kinematic inversions over long time periods (i.e. more than ∼ 40-50 time steps), a Laplacian regularization is preferred. However, the covariance regularization has been used to model the 2020 Atacama sequence afterslip (see Chapter 4).

Spatial and temporal smoothing with a Laplacian regularization

For the sake of clarity, I describe here the principle assuming that subfault i has only two neighbours, the subfaults i -1 and i + 1. The Laplacian minimization is equivalent to an observation equation stating that slip rate at subfault i is the average of the slip rate of subfaults i -

1 and i + 1: ṡi = 1 2 ṡi-1 + ṡi+1 (2.11)
Equation 2.11 leads to the linear system:

      1 -1 • • • • • • • • • • • • • • • • • • • • • • • • • • • -1 2 -1 • • • • • • • • • • • • • • • • • • • • • • • • • • • -1 1       Discrete Laplacian Operator            ṡ1 . . . ṡi-1 ṡi ṡi+1 . . . ṡn            =            0 . . . 0 0 0 . . . 0            (2.12)
The temporal smoothing is added to the same discrete Laplacian operator. For the slip rate ṡi (t) at the time step t and at the subfault i, the temporal smoothing is the average of the slip rate at time steps t -1 and t + 1:

ṡi (t) = 1 2 ṡi (t -1) + ṡi (t + 1) (2.13)
Combining Equations 2.11 and 2.13, the slip rate at subfault i and time step t is the average of the slip rate of subfaults i -1 and i + 1 and at the time steps t -1 and t + 1:

4 ṡi (t) = ṡi-1 (t) + ṡi+1 (t) + ṡi (t -1) + ṡi (t + 1) (2.14)
The spatial and temporal smoothing are now in the same discrete Laplacian operator. However, PYEQ allows for different spatial and temporal smoothing constraints with two dimensionless constants, λ s and λ t . In that case, the observation equation states that the slip rate ṡi is the weighted average of spatial and temporal neighbours, with λ s and λ t , respectively the spatial and the temporal weights:

(2λ s + 2λ t ) ṡi (t) = λ s ṡi-1 (t) + ṡi+1 (t) + λ t ṡi (t -1) + ṡi (t + 1)
(2.15) Equation 2.15 leads to a discrete Laplacian operator D(λ s , λ t ) which allows to write the pseudo-observation equation:

D(λ s , λ t ) ṡ = 0 (2.16)

Towards the inverse problem

We obtain the regularized linear system by stacking vertically the Equations 2.8, 2.9 and 2.15. This gives:

  G I D(λ s , λ t )   ṡ =   d 0 0   (2.17)
Solving for this linear system is equivalent to find the minimum of the following cost function:

C( ṡ) = (G ṡ -d) T C -1 d (G ṡ -d) + λ D ṡ T ṡ + ṡ T D(λ s , λ t ) T D(λ s , λ t ) ṡ (2.18)

Resolution

From a practical point of view, PYEQ directly computes

A = G T C -1 d G and G T C -1
d d from the observations. It then adds λ D to A and finally D(λ s , λ t ) T D(λ s , λ t ), which is fast to compute because it is a sparse matrix. Without the non-negativity constraints, there is a classical direct solution to this problem which is:

ṡ = G T G + λ 2 D I + D(λ s , λ t ) T D(λ s , λ t ) -1 G T d (2.19)
The final linear system is solved by PYEQ using the optimized least-squares algorithm of Hansen (1992), to impose the non-negativity constraints on the slip rate.

Fault geometry and transfer functions

In order to calculate the transfer functions, we shall discretize the fault. In this thesis, we used a model of subduction interface geometry, the SLAB 2.0 model from Hayes et al. (2018) as input data. SLAB 2.0 is a curved surface obtained from a combination of tomography models, earthquakes catalogues and seismic data.

PYEQ discretizes the SLAB 2.0 fault model into quasi-equilateral triangular subfaults using successive division of an icosahedron, wich is a 20-face polyhedron. For this step, we need to define the area and depth ranges of the fault, as well as the size of the subfault, here provided as a number of subdivisions of an icosahedron face.

When the fault geometry is computed, we obtain n subfaults, from which we compute the transfer functions, using the formulation of Nikkhoo and Walter (2015) for triangular dislocations. For each center of the triangles, PYEQ calculates the strike and the dip. In this thesis, I used a fixed rake, determined from the relative motion of plates. The transfer function is therefore calculated in a unique direction.

Choice of model

With three varying parameters (i.e. damping, spatial and temporal smoothing), one can obtain very different models. In this section, we explore different ways to choose the best model.

L-curve

The L-curve approach allows to define the optimal trade-off between the overall model roughness and the misfit to the data. To do so, for each inversion, PYEQ calculates the misfit between the observed displacements vector d and the estimated slip rate vector ŝ, referred as a χ 2 :

χ 2 = (G ŝ -d) T C -1 d (G ŝ -d) (2.20)
PYEQ also calculates a reduced χ 2 r to evaluate if the misfit between the observed and estimated displacements is large (value > 1) or not (value close to 1). For m observations, the reduced χ 2 r is calculated as:

χ 2 r = χ 2 /m (2.21)
The L-curve is an empirical method to select the best model, i.e. the inversion that minimised the regularization constraints and has a good χ 2 (Hansen, 1992). In the case of a time-dependent inversion with spatial and temporal smoothing, the L-curve is a surface plot (Figure 2.11) of the misfit referred as the χ 2 as a function of the roughness and the temporal smoothing.

The best model is located at the inflexion point of the L-curve (the red star on Figure 2.11). An other approach to select the best model is the cross-validation, which is implemented in PYEQ and described in the next section.

Cross-validation

A more quantitative approach to select the optimal regularization parameters is the cross-validation. This is a re-sampling method widely used in machine learning problems and also used for static slip inversions of geodetic data (Matthews & Segall, 1993;Freymueller et al., 1994).

The K-Fold cross validation consists in separating in time and space the data into K subsets of the same size. For a given model, the subset K is removed from the observed displacements vector d. The ability of the remaining data to correctly predict the data of subset K is evaluated by the Mean Square Error (MSE), defined as:

M SE K = (G K ŝK -d K ) T C -1 d K (G K ŝK -d K ) n d K (2.22)
where d K is the observed displacements vector for subset K with a length n d K . ŝK is the modeled displacements estimated without including observations d K . G K is the sub-matrix of G whose lines correspond to the displacements contained in d K .

To avoid too much dependence of the M SE on the chosen subset, an inversion is computed for all subsets K and its M SE K . We then compute a final M SE which is the average of all M SE K . In the case of GNSS data, K = 9 subsets are used (see Figure 2.12). The GNSS sites are split in three sub-networks, and the GNSS time series are split into three time windows.

Too smooth models give a bad fit to all the data and a large M SE K value. Models that are too rough, even if they give a good fit to the data used in the inversion, will give a poor prediction to the d K subset data removed from the inversion.

We have introduced two ways to evaluate which regularization parameters to choose in order to have the best model. We will now see how to evaluate the resolution of kinematic models.

Resolution tests

The full time-dependent inversion approach of PYEQ does not provide the a posteriori covariance matrix because it is very tedious considering today's computational power. However, we performed two different resolution tests to assess the model resolution. I will describe here static and kinematic synthetic tests.

Static synthetic tests

The static synthetic tests aim to assess the spatial resolution of the model. They are usually used for static inversions. These tests consist in building a synthetic model, to generate the corresponding displacements and to invert them. The differences between the model built and the inverted result are proxies to the resolution of the GNSS network and the inversion approach.

Building the synthetic model can be done with different approaches. The first one consists in building a checkerboard model with different sizes of slip patches, regarding the model to assess. The second one is to build a model similar to the one we want to evaluate. For example, if studying the postseismic deformation of a megathrust earthquake, we can build a slip patch surrounding the mainshock rupture as a circle arc.

After building the synthetic model, we calculate the product of the transfer functions (without added noise) and the synthetic model. We obtain the synthetic displacements of the synthetic model. We noise the data with a normal distribution. Finally, we compute the inversion of these synthetic displacements. Here, we shall apply the L-curve approach (section 2.3.6.1) or the cross-validation (see section 2.3.6.2) to select the optimal regularization parameters. When testing different size of slip patches, we can assess what kind of pattern the GNSS network and the inversion approach can resolve. It also provides the spatial resolution, but not the temporal resolution. In the next section, I describe kinematic tests that provide spatial and temporal resolution.

Kinematic synthetic tests

The kinematic synthetic tests are more appropriate for kinematic inversions because they enable to assess how the inversion can retrieve the synthetic model in space and time. The approach to be followed is identical to that of the static tests, except that a synthetic kinematic model is built. As for the static tests, the kinematic synthetics are multiplied by the Green functions which provide the synthetic time series. We invert these synthetic time series and perform a cross-validation to select the best regularization parameters.

Compared to the static tests, besides the added value of the temporal resolution, the kinematic synthetic tests enable to compare the cumulative slip evolution, and also the slip rate evolution for the synthetic and the inverted results.

Considering the previous example of the postseismic deformation, we can assess different aspects of the model. First, a homogeneous slip surrounding the coseismic rupture as a circle arc, following a logarithmic increase through time. To go into more details, we can investigate if small slip modulations through time and space can be retrieved by the GNSS network and the kinematic inversion approach.

Conclusions

Throughout this Chapter, we described how the positioning measurements were obtained, and how these measurements were processed to obtain the time series that we use to make kinematic models. I then briefly described several kinematic inversion approaches and application examples, before presenting the approach used in this thesis, implemented in the PYEQ software. Finally, I detailed the different techniques we use to choose the right model and evaluate its resolution. I will now present the case studies performed during this thesis.

GNSS network

divided in three equivalent subset GNSS time series divided into three equivalent sub-time series

K = 1 K = 2 K = 3 K = 5 K = 6 K = 4 K = 7 K = 8 K = 9
. . .

Inversion without K = 1 Inversion without K = 2
Calculation of the Mean Square Error (MSE) between the remained observations of subset K = 1 and the predicted data for K = 1 computed with K = 2 to 9

Calculation of the Mean Square Error (MSE) between the remained observations of subset K = 2 and the predicted data for K = 2 computed with K = 1 and K = 3 to 9 . . . . . .

. . .

The cross-validation is computed for several regularization parameters. The best prediction of the remained observation K subset gives the optimal regularization parameters.

K = 1 K = 2 K = 3 K = 5 K = 6 K = 4 K = 7 K = 8 K = 9 K = 1 K = 2 K = 3 K = 5 K = 6 K = 4 K = 7 K = 8 K = 9 K = 1 K = 2 K = 3 K = 5 K = 6 K = 4 K = 7 K = 8 K = 9 prediction prediction Figure 2
.12: Illustration of the K-Fold cross-validation using K=9. The GNSS network is divided into three randomly chosen subsets and the time series are divided into three periods of same duration. Each K subset is successively removed for the inversion. The Mean Square Error (MSE) is calculated between the predicted data (red arrows) and the removed K subset and summed for all K. The cross-validation is computed for several regularization parameters. The minimum MSE gives the optimal regularization parameters.

Figure 2.13: Illustration of a static resolution test using a checkerboard model. The checkerboard is build with a certain size of slip patch. The product of the checkerboard and the transfer functions gives the synthetic displacements (blue arrows), noised with a normal distribution. We invert these synthetic displacements for a slip distribution with different regularization parameters (inverted displacements are the green arrows). The best model (blue star) is chosen using the L-curve.

Chapter 3. Afterslip of the M w 8.3 2015 Illapel earthquake imaged through a time-dependent inversion of continuous and survey

GNSS data

After studying time-dependent slip inversion approaches and becoming familiar with PYEQ, I started work on the postseismic sequence of the great M w 8.3 2015 Illapel earthquake in Chile. Survey data acquired after the earthquake had not yet been exploited and the time series showed complexities that had not yet been discussed. In a first section, I present the seismo-tectonic context of the Illapel region. I then briefly present the Illapel earthquake and the interest to study its postseismic phase, and the questions associated with it. Finally, the paper was published in Journal of Geophysical Research: Solid Earth is in section 3.4 and the Supporting Information are in section 3.5.

Seismo-tectonic context of the Illapel region

Description of the regional tectonics

The Chilean subduction zone is approximately 3500 km long, making it one of the largest in the world. It is also one of the fastest: the Nazca plate subducts beneath the South American plate at a rate of 65.5 ± 0.8 mm.yr -1 [START_REF] Jarrin | Current motion and deformation of the Nazca plate: new constraints from GPS measurements[END_REF]. The Juan Fernandez ridge, located on the Nazca plate at latitude 33 • S (see Figure 3.1), subducts beneath the South American plate and coincides with the southern edge of the 2015 M w 8. [START_REF]The M w 8.3 2015 Illapel afterslip spatio-temporal modelling[END_REF] Illapel earthquake rupture, where the subduction megathrust interface is highly coupled (Metois et al., 2016). One specificity of this region is that the Nazca plate flattens at a depth of approximately 80 km, implying that the subduction interface is shallow [START_REF] Barazangi | Spatial distribution of earthquakes and subduction of the Nazca plate beneath South America[END_REF].

Intraslab earthquakes, i.e. thrust events, were documented in this region. In 1997 and 1998, a sequence of earthquakes with magnitudes higher than 6 has brought to light the existence of the Punitaqui thrust fault [START_REF] Lemoine | Evidence for earthquake interaction in Central Chile: the July 1997-September 1998 sequence[END_REF]). An other sequence of intraslab earthquakes occurred in 2019 northwest of Punitaqui, at Coquimbo, along another intraslab normal fault [START_REF] Ruiz | The January 2019 (Mw 6.7) Coquimbo earthquake: insights from a seismic sequence within the Nazca plate[END_REF].

Historical and recent earthquakes

Chile experienced the largest earthquake ever recorded, the 1960 Valdivia M w 9.5 earthquake. As far as we can go back in history, the Chilean subduction zone has frequently hosted M w ∼ 8 events that either ruptured at intermediate depths (35 km in the seismogenic zone) or at deeper depths (> 60 km). Very large events (M w > 8.4) ruptured the whole megathrust, inducing large tsunamis, and have occurred approximately twice per century (Ruiz & Madariaga, 2018).

The Illapel segment (i.e. latitudes 30 • S-32 • S) has ruptured numerous times through large megathrust earthquakes (orange lines in Figure 3.1). In 1730, a M ∼ 9 event ruptured approximately 800 km of the subduction interface, which includes the Illapel segment. It ruptured through a M ∼ 8 event in 1880, and 60 years later, in 1943, with a M w ∼ 7.9 earthquake. In 1971, part of the Illapel segment hosted a M w 7.8 event. Recently, in 2015, the entire Illapel segment experienced a M w 8.3 earthquake.

North of the Illapel segment (at latitudes 28 • S-30 • S, Figure 3.1), two earthquakes partially ruptured the subduction interface. Again, the giant M ∼ 9 earthquake ruptured the southern half of this segment. And in 1922, almost the entire segment ruptured with a M w 8.4-8.7 earthquake.

South of the Illapel segment, at latitudes 32 • S-33.5 • S (Figure 3.1), four large earthquakes ruptured the megathrust. The first time was through the giant 1730 M ∼ 9 earthquake. Almost a century later, it was partially ruptured with the 1822 M ∼ 8 -8.5 event.

The same portion of the megathrust hosted the 1906 M w ∼ 8.2 earthquake and the 1985 M w 8 event.

The M w 8.3 2015 Illapel earthquake

A M w 8.3 earthquake occurred September 16 th , 2015 and ruptured a 200 km long segment of the Chilean subduction zone near Illapel, in northern Chile. This earthquake has been largely studied, using different data sets and methodology. The coseismic models find between 8 and 10 meters of slip, with 2-3 meters of slip that reached the trench, inducing a large tsunami (e.g. Ruiz et al., 2016;Melgar et al., 2016;[START_REF] Grandin | Three-dimensional displacement field of the 2015 Mw8.3 Illapel earthquake (Chile) from across-and along-track Sentinel-1 TOPS interferometry[END_REF][START_REF] Fuentes | A study of the 2015 Mw 8.3 Illapel earthquake and tsunami: Numerical and analytical approaches[END_REF]Klein et al., 2017;An & Meng, 2017). The rupture extended downdip to 55 km through a narrow portion of the megathrust.

The interseismic phase of the Illapel region was also well studied, both regarding the seismicity (Poli et al., 2017) and the interseismic coupling (Métois et al., 2014;Metois et al., 2016). Poli et al. (2017) studied the seismicity during the 20 years preceding the Illapel earthquake and identified regular seismic swarms. These swarms occurred in three distinct areas close to the megathrust area that will surround the coseismic rupture. According to the interseismic coupling model of Metois et al. (2016), the Illapel earthquake ruptured a highly coupled area of the subduction interface.

The postseismic deformation associated with the Illapel earthquake

The postseismic phase of the M w 8.3 earthquake has been well documented. Two studies focused on the very early afterslip (Twardzik et al., 2021;Liu et al., 2022) and several other studies focused on the mid-term (approximately one month) afterslip (Barnhart et al., 2016;Shrivastava et al., 2016;Feng et al., 2017;Huang et al., 2017;[START_REF] Xiang | Spatial-temporal properties of afterslip associated with the 2015 Mw 8.3 Illapel earthquake, Chile[END_REF]. These studies used geodetic and/or seismological data together with different methodologies. The aftershock sequence was also studied and precise aftershock catalogues were produced (Frank et al., 2017;Huang et al., 2017;Poli et al., 2017).

Although the postseismic phase of the Illapel earthquake has already been thoroughly studied, several elements show that additional work is required. First, the continuous GNSS time series of the Illapel area show interesting complexities in the postseismic displacements during the months following the earthquake. These complexities were not modelled nor discussed in previous studies. Second, a large GNSS survey was conducted a few days after the earthquake, by Christophe Vigny and Émilie Klein (CNRS, ENS) to monitor the postseismic deformation of this large earthquake. These additional data, so far not used for afterslip modelling, gave us the opportunity to develop a more detailed model of the postseismic deformation in this region. Third, a full time-dependent inversion of the postseismic displacements using the PCAIM approach (Shrivastava et al. (2016), see section 2.2.1 for the PCAIM description) was published for the Illapel earthquake and represents an opportunity to compare results obtained with PYEQ (see section 2.3). As a matter of fact, only one other study used a time-dependent inversion to image the postseismic deformation. It was for the 2016 M w 7.8 Pedernales earthquake in Ecuador (Rolandone et al., 2018). Finally, a kinematic model of the afterslip would allow a thorough comparison with the available aftershock catalogues (Frank et al., 2017;Huang et al., 2017;Poli et al., 2017).

The sequence of aftershocks of the Illapel earthquake presents some particularities. On November 7, 51 days after the mainshock, a large M w 6.8 aftershock ruptured downdip of the coseismic rupture. Four days later, on November 11, two large M w 6.9 aftershocks ruptured approximately 100 km north of the mainshock area. Why did these aftershocks occur so late in the postseismic sequence and so far from the mainshock rupture? Does the afterslip experience pulses of slip after their occurrence? Are there any other dynamics at play, such as a slow slip event? To study the spatio-temporal evolution of the different types of slip, i.e. the seismic rupture of the M w 8. [START_REF]The M w 8.3 2015 Illapel afterslip spatio-temporal modelling[END_REF] Illapel earthquake, its afterslip and the aftershocks that followed (the large ones and the microseismicity) would provide a first answer to the questions outlined in this thesis. Namely, to study the interaction of aseismic slip processes with seismicity to better understand the behavior of faults, possible with a time-dependent inversion of the postseismic deformation. The wide range of studies conducted on the entire seismic cycle of the Illapel segment is also a great opportunity to discuss the relationship between the different interseismic, coseismic and postseismic phases.

Abstract

We use continuous and survey GNSS data to image the spatial and temporal evolution of afterslip during the two-months following the M w 8.3 2015 Illapel earthquake. Our approach solves for the incremental daily slip at the subduction interface using non-negative least-squares with spatial and temporal Laplacian regularization constraints. We find that afterslip developed at three specific areas at the megathrust, surrounding the coseismic rupture. In addition, well resolved afterslip also occurs within the coseismic rupture area that experienced ∼4 meters of seismic slip. Our afterslip model shows striking correlations with the spatial distribution of aftershocks and repeating earthquakes.

We capture the local afterslip triggered by a M w 6.8 and two 6.9 aftershocks that ruptured downdip and north of the coseismic rupture respectively. The latter ones were possibly triggered by the afterslip that developed north of the rupture. We also find a pulse of enhanced aseismic slip lasting a few days south of the rupture that spatially and temporally correlates with a seismicity burst. We finally find that areas of enhanced afterslip spatially correlates with areas having experienced regular seismic swarms observed during the years prior to the Illapel earthquake. This correlation supports the view of localized fluid high pore pressure areas behaving aseismically and surrounding a highly locked asperity, preventing the seismic rupture to propagate into them.

Plain Language Summary

After a large earthquake, aseismic slip occurs in the vicinity of the seismic rupture.

In this study, we focus on the aseismic slip that occurred after a large subduction earthquake of magnitude 8.3 in 2015 in central Chile. Using GPS time series, we obtain a spatiotemporal view of the aseismic slip during the two months following the earthquake. We find that aseismic slip occurred at three preferential areas located at the periphery of the rupture, as observed for many earthquakes. However, we also identify a patch of aseismic slip within the rupture area. We could image the daily afterslip rate and compare it to the daily seismicity and found striking correlations with the spatial distribution of aftershocks. In particular, we find that the largest aftershocks with magnitude 6.8-6.9 are tightly related to the evolution of afterslip. This result suggests that precisely monitoring the evolution of afterslip after a large earthquake might help to define preferential areas of large aftershocks. We also find a potential slow slip event south of the mainshock rupture which correlates with a local increase of seismicity. Finally, we find that the areas of enhanced afterslip had experienced seismic swarms decades before the earthquake, possibly reflecting areas of the fault with high fluid pore pressure. These results highlight new aspects about the dynamics of afterslip.

Introduction

During the weeks following a large earthquake, transient aseismic slip, referred as afterslip, develops at the fault and dominates the deformation observed at the Earth' surface. The development of continuous GNSS networks at subduction zones for almost three decades has allowed to capture the main pattern of the afterslip that followed some major megathrust events: afterslip initiates immediately after the earthquake, with a slip rate decreasing through time as 1 t leading to a logarithmic growth of the cumulative slip (e.g. Perfettini et al., 2010). Afterslip takes place at areas of the fault surrounding the coseismic rupture (Perfettini et al., 2010), with little if any slip inside the coseismic rupture. The evolution of the cumulative number of aftershocks mimics the afterslip evolution both in time and space, and the moment released through aftershocks is only a small fraction of the afterslip equivalent moment (Perfettini & Avouac, 2004;Hsu et al., 2006).

Afterslip has been interpreted as the frictional response of the fault to the stress increment induced by the earthquake at the area surrounding the coseismic slip. In a rate-and-state friction framework, areas experiencing afterslip are supposed to follow a ratestrengthening regime where friction increases with the sliding velocity [START_REF] Lay | An asperity model of large earthquake sequences[END_REF]Avouac, 2015). Specifically, rate-strengthening law successfully predicts the logarithmic evolution of the cumulated afterslip. First order comparison between the location of areas that are locked during the interseismic phase, the location of the seismic rupture and afterslip support a view where the fault is composed of locked patches accumulating slip deficit released through seismic ruptures, embedded in an overall slipping and aseismic fault plane [START_REF] Lay | An asperity model of large earthquake sequences[END_REF]Avouac, 2015).

However, despite the overall success of the above "rate-state asperity model" in explaining the main pattern of afterslip, there are several observations that suggest a more complex view. For a few cases, afterslip and aftershocks have been found to be located within the coseismic rupture area (Johnson et al., 2012;Agurto et al., 2012;Bedford et al., 2013;Tsang et al., 2019). Some afterslip has also been found at areas locked during the years prior the earthquake (Rolandone et al., 2018). Finally, afterslip develops at specific locations rather than as a rim surrounding the coseismic rupture as would be expected from the coseismic stress increment (Perfettini et al., 2010;Bedford et al., 2013;Rolandone et al., 2018).

Finally, some studies argued that Slow Slip Events (SSEs) during the interseismic period and areas of afterslip take place at the same locations (Hobbs et al., 2017;Rolandone et al., 2018). Those results question the friction law and the processes that could simultaneously explain the spontaneous occurrence of SSE during the interseismic phase and afterslip. A related open question is to know whether we can observe some modulation of the afterslip rate, that would be superimposed to the overall long-term afterslip decay. Such a behavior has been proposed by Bedford et al. (2013) for the 2010 M w 8.8 Maule earthquake, where the inversion of GNSS data shows several pulses of aseismic slip during 2.5-9 months after the mainshock. However, such pulses have not been documented for others large megathrust events. (i.e. 1880, 1943 and 2015), and is thought to be the northern segment of the M∼9 1730 great earthquake with a proposed rupture length exceeding 500 km (Ruiz & Madariaga, 2018). Coseismic slip models consistently find a maximum slip of 8-10 m, with several meters of slip reaching the trench, attested by a significant tsunami (Melgar et al., 2016;Ruiz et al., 2016;Barnhart et al., 2016;Shrivastava et al., 2016;Klein et al., 2017). All coseismic slip models also find that large slip (> 5m) remained confined within the first shallowest ∼35 km. They however also highlight that several meters of seismic slip extended downdip to ∼55 km depth along a 50 km wide portion of the megathrust at latitude 31 • S (Figure 1).

Previous studies used different data sets, time periods and methodology to estimate the afterslip. The early stage of afterslip (hours to 3 days) was studied by Twardzik et al. (2021) and Liu et al. (2022), using kinematic GNSS data and successive static inversions through time. Both studies found two patches of afterslip, surrounding the northern and the southern deep edges of the mainshock rupture, but highlighted some peculiar behaviors of the early afterslip. Liu et al. (2022) found that the rapid displacement evolution observed at GNSS site CNBA (Figure 1) evolves as a power law, rather than the usual logarithmic model. However, Twardzik et al. (2021) found that the coseismic slip associated with large aftershocks (up to M w 7.1) for the southern patch (close to CNBA) can explain all the geodetically determined slip during the first 12 hours for that area.

This result highlights that locally, seismic slip through large aftershocks can significantly contribute to the total slip and certainly impact the afterslip evolution.

After the early stage, several studies based on daily GNSS solutions and InSAR results concurrently find that afterslip develops from two discrete areas located along the northern and southern edges of the coseismic rupture. The two areas of afterslip are separated by the deepest portion of the coseismic rupture (Barnhart et al., 2016;Shrivastava et al., 2016;Huang et al., 2017;Feng et al., 2017). Detailed studies of the aftershock sequence show that aftershocks are mainly located at or close to areas experiencing large afterslip (Huang et al., 2017). Frank et al. (2017) and Liu et al. (2022) further find that the aftershock and afterslip areas are spreading along strike through time, an observation which is consistent with the prediction for rate-strengthening areas response to a coseismic Coulomb stress increment.

Although these observations agree with the 'rate-state asperity' model, there are also several characteristics that challenge that view. A peculiar feature is the absence of significant afterslip downdip the deepest part of the rupture (at latitude 31.2 • S, 50 km depth). There is also no afterslip at shallow depth (<10 km) from latitude 30.8 • S along the southern edge of the coseismic rupture. In addition, for most models, part of the afterslip overlaps with areas having slipped seismically by 3-6 m during the mainshock (Barnhart et al., 2016;Shrivastava et al., 2016).

In this study, we aim at precisely characterizing the spatial evolution of afterslip following the Illapel earthquake. For that purpose, we perform a full time-dependent inversion of the daily slip from GNSS time series, including both continuous and survey measurements. While previous studies have described the cumulative afterslip evolution, our approach focuses on the dynamic of afterslip over periods of a few days. In particular, we investigate whether afterslip smoothly increases through time as a logarithmic function or whether slip pulses or modulations can be identified. We specifically compare our model with the seismicity rate evolution and investigate the relationship between afterslip and the occurrence of the largest aftershocks.

3 Data

Data

Our input data set is composed of 28 continuous GNSS (cGNSS) sites and 17 survey (sGNSS) sites (Figures 1 andS1). cGNSS sites in the Illapel area were installed by the joint international laboratory (LIA) Montessus de Ballore [START_REF] Vigny | Upper plate deformation measured by GPS in the Coquimbo Gap, Chile[END_REF] 2018)). We use the cGNSS time series from the processing SOAM GNSS solENS solution described in Klein et al. (2022). Postseismic time series are obtained by removing the interseismic velocity before the Illapel earthquake estimated from four years of data with the MIDAS method of Blewitt et al. (2016).

The selected cGNSS network used in this study (see Figure 1) spreads ∼150 km along-strike north and south of the coseismic rupture and ∼ 400 km inland from the trench.

22 cGNSS sites recorded data during the whole period of our study from September 17

to November 31 (74 days), while 6 of them provided approximately a month and a half of data. Two large aftershocks (M w 7.1 and M w 6.8) ruptured within the 6 hours that followed the mainshock south of the epicenter (Figure 1) causing a few centimeters displacements (Twardzik et al., 2021). Our first observation corresponds to the 24 hours average position during day September 17, mitigating the potential bias induced by the aftershocks on the first position. The cGNSS data are complemented by 17 sGNSS sites (Ruiz & Madariaga, 2018). The dashed lines are the Slab2.0 iso-depth contours every 20 km (Hayes et al., 2018).

re-occupied after the Illapel earthquake and provide 4 to 9 days of measurements, between 2015/09/27 and 2015/11/30. sGNSS sites bring additional information because postseismic daily displacements were still of the order of several millimeters per day during the survey period (Figure 2b), larger than the daily position uncertainty. Among the sGNSS data used, 5 sites are located in the immediate vicinity of the mainshock rupture and 7 others are located less than 50 km from the rupture edge, providing dense measurements above the downdip limit of the coseismic rupture. Compared to previous studies, these sGNSS sites provide additional constraints on the afterslip evolution during their observation period, contributing to make the model more robust.

Time series patterns

After the Illapel earthquake, cGNSS times series show logarithm-like, rapid, postseismic displacements on the east component during the first two weeks following the mainshock (Figure 2). The westward displacement reaches 10 cm at PFRJ and 5 cm at CNBA close to the rupture (Figure 1). Located one hundred kilometers north of the rupture area, cGNSS site TOLO shows almost as much displacement as CNBA, suggesting a large amount of afterslip north of the rupture.

As an example, Figure 2 shows the time series at survey site POBR which was installed on 2015/09/30, two weeks after the earthquake north of the rupture. POBR shows almost 2 cm of westward displacement during a week period, demonstrating the interest of including sGNSS data in afterslip modelling.

The GNSS network also recorded the coseismic displacements of the largest aftershocks of the sequence. The first one was a M w 6.8 earthquake and ruptured on November 7 (51 days after the mainshock) at ∼40 km depth. The others were two M w 6.9 earthquakes that occurred on November 11 (55 days after the mainshock) at shallow depth close to the trench, ∼100 km north of the rupture. In order to investigate the relationship between the afterslip and the occurrence of those earthquakes, for both dates, we extracted the coseismic offsets and invert for the coseismic slip spatial distribution (see Supporting Information, Figure S6).

4 Time-dependent modeling of the postseismic deformation

Method

The afterslip of the Illapel earthquake has been studied using static inversions of GNSS displacement observed at successive dates (Barnhart et al., 2016;Feng et al., 2017;[START_REF] Xiang | Spatial-temporal properties of afterslip associated with the 2015 Mw 8.3 Illapel earthquake, Chile[END_REF]. Although static inversions are useful to determine the location of afterslip, it has several drawbacks: it does not warrant that the modeled slip from two successive models is meaningful, and could for instance lead to non-physical negative slip.

This would be critical for successive dates where the observed displacement is small, close to the data noise level. Successive static inversions cannot account for possible differential displacement among sites during the period, that carry information of the path followed by the slip. Also, new GNSS sites installed after the earthquake cannot be included in the inversion because their displacement since the earthquake is unknown. In that case, a different time window can be selected, but with the difficulty of inserting the results in a consistent way in the inversion of the whole time period.

Our approach overcomes these limitations by simultaneously solving for the spatially variable daily slip rate using all available time series. It therefore includes the information of the cumulated displacement at a given date as static inversions would do, but also uses the information of the trend of the observed displacement. It further includes the constraint that incremental daily slip must be positive. Shrivastava et al. (2016) apply the PCAIM method of Kositsky and Avouac (2010) with a time function. Compared to PCAIM, our approach does not include the assumption that slip is a combination of spatially stationary sources. Here, we summarize the method together with the relevant information for the kinematic inversion of the Illapel afterslip.

The plate interface is discretized into 1014 quasi-equilateral triangular subfaults with 14 km long edges (see Figure S2), following the curved surface from the Slab2.0 model (Hayes et al., 2018). The fault area extends from latitude 32 • S to 28 • S and from the trench down to 90 km depth. The cumulative slip s i (t) at each subfault i at time t is parametrized as a monotonic increasing piecewise linear function of time:

s i (t) = t k ≤t t k =0 ṡi (t k )∆ k (1)
where ṡi (t k ) is the constant slip rate during day k with duration ∆ k = 1 day. The displacement d jl (t) observed for component j at site l is:

n i=1 g ijl s i (t) = d jl (t) ( 2
)
where n is the number of subfaults, g ijl is the Green function, here calculated in an homogeneous elastic half space for triangular dislocation (Nikkhoo & Walter, 2015). The rake is kept fixed at a value derived from the direction of the Nazca/South America plate convergence from [START_REF] Vigny | Upper plate deformation measured by GPS in the Coquimbo Gap, Chile[END_REF]. Combining Equations ( 1) and ( 2) gives the ob-servation equation:

t k ≤t t k =0 n i=1 g ijl ṡi (t k ) = d jl (t) (3)
Stacking ṡi (t k ) for all dates in a single vector ṡ, and all time series d jl (t) in vector d leads to the linear system:

G ṡ = d (4)
The inverse problem associated with Equation ( 4) is highly under-determined. We add regularization constraints in the form of minimum Laplacian in space and time. For that, we build a discrete Laplacian operator by setting that every ṡi (t k ) equals (in the least squares sense) the weighted average of its spatial neighbor at t k and its time neighbors at time t k-1 and t k+1 . Although applied simultaneously in the time and space domains, this approach allows for different spatial and temporal smoothing constraints, controlled by the dimensionless constants λ space and λ time respectively. As input data, we use both horizontal and vertical displacements at GNSS sites together with their associated uncertainty re-scaled to account for the observed daily scatter. Because of larger uncertainty and sensitivity to the Earth' elastic structure, the vertical component is further down-weighted by a factor of 3 in the inversions. This value has been empirically chosen so that the sum of the squared normalized residuals (residual divided by the standard deviation) is similar for the three components east, north and up.

Our method offers the advantage of explicitly imposing a non-negative constraint on the inverted daily slip. The method also accepts that some data might be missing. This is particularly useful for postseismic modelling where additional sites are often deployed several days or weeks after the mainshock. Here, for instance, it allows us to include 4 to 16 days of survey data in our analysis. Finally, it allows a simple handling of aftershocks without the need to correct the time series for coseismic offsets. Indeed, Equation ( 2) can also be written as:

n i=1 g ijl s i (t -t a ) = d jl (t -t a ) (5)
where t a is the origin time of an aftershock. So handling aftershocks only requires to consider a time series as two separate time series before and after an aftershock. This approach was applied in our inversion for the M w 6.8 aftershock of November 7 and the M w 6.9 doublet aftershocks of November 11, considered as a single event.

We apply our method to solve for the daily incremental slip from September 17, 2015 to November 30, 2015. We selected this period because it ensures that afterslip will be the dominant process with only a small visco-elastic contribution (e.g. [START_REF] Tian | Spatio-temporal variations of afterslip and viscoelastic relaxation following the mw 7.8 gorkha (nepal) earthquake[END_REF]).

This period also includes the occurrence of several large aftershocks, allowing to investigate their potential relationship with afterslip.

Regularization and choice of model

We used two approaches to select optimal regularization parameters. The so-called L-curve approach defines the optimal trade-off between the overall model roughness and misfit to the data (Hansen, 1992). Plotting the misfit as a function of our regularization parameters defines a surface where misfit decreases as weaker regularization parameters are used. An inflexion point is searched for which smaller regularization value only result in a marginal reduction of misfit. Here, we explored a range of combination of spatial and temporal smoothness parameters and compute the misfit defined as χ 2 = (Gŝ- -d) for each inverted model ŝ (Figure S3). The obtained surface shows sharper variation with λ space leaving the choice of λ time more subjective. We selected a model with λ space = 0.8 and λ time = 5 that we refer as the smooth model. As an alternative method to select an optimal model, we use a cross-validation approach. The cross-validation has been used to select the optimal spatial smoothing parameters for static slip inversion from geodetic data (Matthews & Segall, 1993;Freymueller et al., 1994) using the

d) T C -1 d (Gŝ
Leave-Out-One approach. We generalize this to the time-and-space case using K-Fold cross-validation, nowadays widely used in Machine Learning problems. In that approach, the data is splitted into K subsets of equivalent size. For a given model, the K th subset is removed from the observations included in the inversion and the ability of the remaining data to correctly predict data of the K th subset is evaluated by computing the Mean Square Error (M SE) defined by: Compared with the result of the L-curve, the cross-validation approach suggests that optimal models require lower temporal smoothing, allowing rougher models in time.

M SE k = (G k ŝk -d k ) T C -1 d k (G k ŝk -d k )/n d k (6
We select a model with λ space = 1 and λ time = 0.1 that we refer as the rough model (Fig- ure 4b).

Resolution tests

Resolution of the kinematic inversion is difficult to assess since it would require to evaluate the model resolution both in space and time. Semi-analytical solutions for the posterior covariance and marginal probability density functions have been proposed for the case of a linear inverse problem with covariance regularization and non-negative constraint (Nocquet, 2018). However, such approach remains computationally out of reach for the present problem which includes 75 000 parameters (1014 subfaults × 74 days).

In the Supporting Information (see Figure S4), we provide static inversion checkerboard tests to evaluate the spatial resolution of our models. These tests show that patches with size of 50 × 50 km 2 can be resolved close to the coastline, while the resolution decreases close to the trench and to a lesser extent at depth. A kinematic inversion is expected to provide a better resolution than a static inversion at a final date of slip because the kinematic approach requires the whole history of displacement to be fitted. Hence, the tests presented in the Supporting Information are probably pessimistic indicators of the true spatial resolution of the kinematic inversion results.

To get a sense of both the temporal and spatial resolution, we perform dynamic resolution tests by generating time evolving synthetic models and evaluate how well they can be retrieved in space and time. Details are provided in the Supporting Information We find that the kinematic inversion is able to retrieve with a good accuracy the spatial distribution of afterslip distribution, with only reduced slip close to the trench at 31.5 • S (Figure S5a). The synthetic test reveals only mild method artifacts which would erroneously concentrate the slip [START_REF] Page | Constraining earthquake source inversions with GPS data: 1. resolution-based removal of artifacts[END_REF]. The overall time evolution is very well retrieved, as demonstrated by the similarity of the moment rate through time, with only slight underestimation of moment (Figure S5a).

Given the good results obtained, we generate a new synthetic model by adding a slip pulse to the previous synthetic afterslip model. The synthetic slip pulse has a gaussian spatial distribution and a gaussian slip rate evolution through time. This model aims at evaluating whether some modulation of slip rate, small compared to the afterslip, could be retrieved by the inversion. Figure S5b shows the resulting inversion. Aside the overall decreasing rate, the inversion images a patch of enhanced slip, slightly spread compared to the synthetic model, but occurring at the time of the pulse.

Cumulative afterslip evolution

As a result of our time-dependent inversion, we simultaneously obtain the spatial distribution of cumulative afterslip through time and its time derivative, the afterslip rate.

We first describe the evolution of the cumulative afterslip and discuss its correlation with the seismicity. Then in section 6, we focus on the high frequency dynamics of afterslip and its relation with the large aftershocks.

Cumulative afterslip

Both the smooth and the rough models find that afterslip initiates immediately after the mainshock at three distinct areas (Figure 5a At a first order, the three areas of slip take place along the edges of the coseismic rupture, but do not perfectly surround it. More specifically, our results, with denser data, confirm that no afterslip developed downdip of the coseismic rupture at depth greater than 60 km. Overall, the pattern of afterslip shows significant north-south asymmetry in the amount of afterslip. Indeed, as also found in previous studies, the northern patch (Box A in Figure 5g) appears to be the largest, both in size and amount of slip (Barnhart et al., 2016;Huang et al., 2017;[START_REF] Xiang | Spatial-temporal properties of afterslip associated with the 2015 Mw 8.3 Illapel earthquake, Chile[END_REF]. Its size is 80 km along dip and 70 km along strike, about twice the size of the deeper patches (Boxes B & C in Figure 5g).

The cumulative slip at Box A reaches 50 cm after a month and exceeds 80 cm after 74 days, while it is 50 cm and 40 cm for the Boxes B and C respectively at the end of the studied period. Hence, slip in Box A contributes to 70% of the moment released through aseismic slip during the 74 days. Our inversion finds slip occurring at shallow depth, reaching the trench while previous studies found it to be restricted deeper than 10 km. Given the large distance from the GNSS sites, the shallowest afterslip cannot be well resolved.

The observed difference among afterslip models, aside the different data sets, most probably reflects the impact of different regularization constraints. It however shows that a model involving shallow afterslip reaching the trench in Box A is allowed by the data.

As also found in previous studies (Barnhart et al., 2016;Shrivastava et al., 2016;Feng et al., 2017;[START_REF] Xiang | Spatial-temporal properties of afterslip associated with the 2015 Mw 8.3 Illapel earthquake, Chile[END_REF]Twardzik et al., 2021;Liu et al., 2022), the shallowest part of Box B appears to overlap with an area of significant (2-4 m) coseismic slip.

This latter area is close to the coast which hosts several GNSS sites, making this result robust. Box C shows similar size and slip as Box B, with a propagation of slip southward away from the rupture occurring during the first 10 days (see movie S1). Box C is located where the two early aftershocks (M w 7.1 & 6.8) occurred within the first 6 hours after the mainshock. The early afterslip imaged by the inversion at this area must also include the afterslip of these early aftershocks.

Overall, after 74 days, taking a shear modulus of 30 GPa, the moment released by the afterslip is 3.53×10 20 N.m, equivalent to a M w 7.6 (Figure 5h), which represents approximately 13% of the seismic moment released by the M w 8. [START_REF]The M w 8.3 2015 Illapel afterslip spatio-temporal modelling[END_REF] Illapel earthquake. The cumulative moment and cumulative number of aftershocks show a similar time evolution (Figure 5h), as commonly found for other earthquakes (e.g. Hsu et al., 2006;Perfettini et al., 2010).

Comparison of cumulative afterslip with seismicity

The aftershocks sequence following the Illapel earthquake has been studied in great details, offering the opportunity for further tests of the dynamics of slip at specific locations along the plate interface. For comparison with our model, we used two different studies of seismicity. Huang et al. (2017) used a template matching approach to produce an improved earthquake catalogue during a one month period following the mainshock (2015/09/16 to 2015/10/16). From this catalogue, they extracted 291 sequences of repeating earthquakes. They used the empirical relation between repeaters seismic moment and aseismic slip from [START_REF] Nadeau | Seismological studies at Parkfield VI: Moment release rates and estimates of source parameters for small repeating earthquakes[END_REF] to derive estimates of incremental aseismic slip occurring along the fault. These increments of slip are then summed to provide a map of the cumulative aseismic slip, hereafter referred as repeater afterslip model. terslip derived from geodetic data (Barnhart et al., 2016;Shrivastava et al., 2016;Huang et al., 2017).

Our geodetic afterslip models and the repeater afterslip model consistently find that both the slip amount and the moment release are significantly larger in Box 1 north of the coseismic rupture than in Boxes 2 & 3 (Figure 6a,b). The density of aftershocks from Frank et al. (2017) catalogue does not highlight a significantly higher aftershock activity at Box 1 (Figure 6c). On the contrary, it finds more aftershocks taking place at Box 3, perhaps because of a large number of aftershocks in this area during the first hours following the mainshock (Twardzik et al., 2021). Some spatial correlations between our models and the repeater afterslip model for some details are striking (Figure 6a There are also some noticeable differences. In terms of relative slip between the three areas of afterslip, Huang et al. (2017) model find that the Box 3 experiences more slip than the Box 2 while they are similar in our model (Figure 6a,b). This is perhaps also due to the contribution of aftershocks during the first hours following the mainshock (Twardzik et al., 2021), a pattern that our inversions cannot capture because they use 24 hours averaged data, starting the day following the mainshock.

The overall good spatial correlation found between the geodetic afterslip models with aftershocks confirms the view that afterslip drives the seismicity evolution and the repeated failure of small asperities (Figure 6). However, the comparison of seismicity with our afterslip models suggests that areas of higher aftershock activity are located at the edges of high afterslip areas rather than centered at the maximum slip areas (Figure 6c).

This pattern can clearly be seen for the southern patch (Box 3 in Figure 6c) where patches with higher density of aftershocks are located up-dip of high afterslip patches. To a lesser extent, a similar pattern is observed for the northern and deeper patches (Figure 6c).

The comparison with the repeater afterslip model also suggests that repeaters were preferentially activated at the periphery of high afterslip patches (Figure 6b). We acknowledge that the observed spatial correlation between slip gradient and aftershock/repeater activity cannot be demonstrated given ( 1

Slip rate evolution

Most studies of afterslip so far have focused on the cumulative slip evolution and little attention has been paid on the dynamics over short periods. Our approach which solves for the slip rate offers the opportunity to test whether slip rate decreases smoothly in time or experiences pulses of slip at some areas of the fault. Figure 7a-f shows the spatial distribution of the afterslip rate at selected dates, from our smooth model. Figure 7g shows the afterslip moment rate together with the daily number of aftershocks for the smooth model. Because the cross-validation test suggests that rougher time evolution is resolved in the data, Figure 7h shows the afterslip moment rate from our rough model.

Velocity-strengthening response to stress increment predicts that high slip rate should decrease rapidly as 1 t . The overall moment rate function shown in Figure 7g,h shows such a rapid decay, with a moment rate decreasing by a factor of 3 after ten days. Fitting a simple ṡ0 /(1+t/τ ) function to the moment rate function of Figure 7g, we find a characteristic decay time of τ ∼ 5 days. Afterslip rate appears to be the largest at the northern patch (∼50 mm/day versus 25 mm/day for the deep and southern patches). The slip rate at the three patches decreases through time down to a few mm/day after 2.5 months.

However, superimposed to the overall 1 t decay, an increase of the afterslip moment rate is observed around the period of occurrence of the largest aftershocks (Figure 7g,h), between 50 and 60 days after the mainshock. This period is also marked by an increase of the microseismicity daily rate (Figure 7g,h). In the following, we investigate in details the potential relationship between the afterslip and those large aftershocks.

6.1 Afterslip behavior and the M w 6.8, November 7 deep aftershock

The M w 6.8 aftershock occurred on November 7, 51 days after the mainshock. The inversion of the static offset finds that it ruptured a 20 × 20 km 2 patch at 40 km depth, north of the deepest extent of the coseismic rupture (see Supporting Information and Figure S6). Our cumulative afterslip model finds that this aftershock occurred in one of the three areas experiencing enhanced afterslip (Figure 5e). For that area, the afterslip rate model indicates that the slip rate decreases from ∼35 mm/day on the day after the mainshock down to less than 5 mm/day on October 29, that is 42 days after the mainshock (Figure 7a-c). This value of slip rate is similar at the three areas of afterslip.

However, from day 48 onward, while the slip rate continues to decrease at Boxes A & C (respectively north and south of the coseimic rupture), slip rate increases at the deepest part of Box B (the 31 • S patch, Figure 7d). A maximum slip rate of ∼10 mm/day is found the day following the M w 6.8 aftershock (Figure 7e), before decreasing during the next following 10 days.

The aftershock coseismic offset does not impact our inversion because times series of nearby GNSS sites have been split at the time of the aftershock to form independent observations and hence does not contribute to the modeled afterslip. However, the temporal smoothing constraints required to stabilize our inversion would tend to spread any daily sharp slip rate increment over a period of a few days. In order to evaluate this potential 'time spreading' effect in our inversion and to better characterize the slip rate evolution around November 7, we perform an inversion starting on September 17 but ending on November 7, the day before the aftershock. This inversion does not show any evidence of slip acceleration before the aftershock (Figure 8a). We then perform a second inversion starting on November 7 and ending on November 10. This latter inversion highlights localized rapid slip (∼10 mm/day), encompassing the aftershock rupture area and rapidly decreasing to 4 mm/day the next day (Figure 8b,c).

We conclude that ( 1) the M w 6.8 November 7 aftershock occurred in an area of enhanced afterslip compared to the neighbour area at the fault, ( 2) this aftershock generated its own local afterslip during a few days that superimposed to the overall evolution of the mainshock afterslip, (3) our inversion approach could capture the local slip rate increase from the small signal in the time series, that had not been seen in previous studies. However, the local slip acceleration is blurred in time because of the temporal constraints we imposed, whose optimal values were adjusted over the whole period of two months and all GNSS sites. Specifically 'tuned' inversions can however overcome this limitation and can successfully image the local afterslip triggered by the aftershock.

6.2 Afterslip behavior and the two M w 6.9, November 11 shallow aftershocks

The doublet of moderate size M w 6.9 earthquakes occurred 55 days after the mainshock on November 11, at a distance of ∼100 km north of the mainshock rupture area (Figure S6). Both the remote location of these aftershocks with respect to the mainshock rupture area and the 55 days delay question whether they can be considered as aftershocks or not. Our kinematic inversion first reveals large afterslip developing north of the rupture and possibly at shallow depth. Our smooth geodetic afterslip model further finds that afterslip would have reached the area ruptured by the November 11 events (Fig-

ures 5e and 7e

). This result is in agreement with Frank et al. (2017) who showed that the aftershocks front propagated fast and on a large distance at shallow depth, north of the rupture. Using the same procedure as for the deep November 7 aftershock, we tested whether the slight slip rate increase noted in Figure 7e actually occurred during the days preceding the two November 11 events. The inversion starting on November 7 (after the M w 6.8 event) and ending on November 10 (before the M w 6.9 events) does not show ev- (Dziewonski et al., 1981;Ekström et al., 2012). The mainshock slip is contoured in blue every 2 m (Klein et al., 2017). The blue star is the Mw 8.3 2015 Illapel epicenter. The dashed lines are the Slab2.0 iso-depths every 20 km (Hayes et al., 2018).

idence for significant slip acceleration before November 11 (Figure 9a). We therefore believe that the small slip increase in Figure 7e rather reflects the local afterslip induced by the November 7 earthquake. A small amount of afterslip is also noticed north of the November 11 earthquakes (Figure 7f) likely triggered by these earthquakes.

In order to confirm this view, we performed an inversion starting on November 11 (Figure 9b-d), a few hours after the M w 6.9 aftershocks. Despite some noise, this inversion clearly captures the shallow afterslip triggered by the two M w 6.9 November 11 aftershocks. The aftershock-triggered afterslip area appears to be relatively large given the moderate magnitude of the two aftershocks (equivalent to a M w 7.1). It is certainly not well resolved considering the large distance to the nearest GNSS sites, but the location of the found afterslip correlates pretty well with the location of the aftershocks triggered by the November 11 earthquakes (Figure 9b,c). The obtained model further clearly shows a slip rate decrease, for instance from 8 mm/day during the day of the aftershock down to 2 mm/day close to the aftershocks, a decrease rate consistent with the overall decay time estimated for the mainshock.

Interestingly, we also note that additional slip is suggested by our model for the deep patch. It also correlates with some seismicity activity possibly reflecting ongoing local afterslip following the November 7 aftershock (Figure 9a-d).

In summary, we find that the largest aftershocks during the 74 days following the Illapel earthquake are linked with the afterslip evolution. The deepest one on November 7 probably corresponds to a locked asperity brought to failure by the surrounding afterslip. The 100 km remote distant and 2 months delayed aftershocks were likely triggered by the large amount of afterslip north of the coseismic rupture.

A Slow Slip Event at the southern patch?

The rougher afterslip model shows additional modulations of the moment rate, which are not related to the occurrence of any significant aftershock, but correlates with an increase of the seismicity rate (Figure 7h). By inspection of movie S4, we do not find a clear spatial correlation between the seismicity rate and the local slip rate increase neither around day 19 nor day 25. However, from day 67 onward, the slip rate increases at the southern patch of afterslip at ∼25 km depth. The slip rate culminates with more than 5 mm/day of slip observed on day 69. After this day, the slip rate slowly decreases until it vanishes on day 72 (Figure 10). We find that this potential 6-day slow slip correlates with a local seismicity increase at the tip of the transient slip area. The synthetic test shown in movie S6 suggests that the data can resolve very small transient signal over a few days.

Here, the equivalent moment released through the suspected transient is M w 6.1-6.2. This is indeed very small, but the correlation with the increase of seismicity simultaneously in space and time is a strong support for the hypothesis of a slow slip event superimposed to the overall afterslip evolution. The comparison with the interseismic regime described in Section 7.3 indicates that this area of the subduction interface is prone to regular slow slip events and seismic swarms.

Discussion

Overall pattern of afterslip

Overall, our results confirm previous findings of immediate relocking of the coseismically ruptured area and afterslip developing at the edge of the coseismic rupture. Our kinematic inversion finds the logarithm evolution of the cumulative slip and the 1 t decrease of the slip rate with a characteristic decay time of ∼5 days. The cumulative moment of afterslip is 2.37×10 20 N.m after a month, equivalent to a M w 7.5, a value similar to Barnhart et al. (2016) and Huang et al. (2017). After 74 days, the cumulative mo- 2 M w 6.9 2 M w 6.9 2 M w 6.9 2 M w 6.9 M w 6.8 M w 6.8 M w 6.8 M w 6.8 M w 6.8

2 M w 6.9

Slip rate (mm/day) The blue star is the Mw 8.3 2015 Illapel epicenter. The dashed lines are the Slab2.0 iso-depths every 20 km (Hayes et al., 2018).

ment of afterslip is 3.53×10 20 N.m equivalent to a M w 7.6, corresponding to 10% of the coseismic moment. Compared to previous studies, our model finds that, rather than taking place homogeneously around the rupture, afterslip develops from relatively localized areas. No afterslip developed downdip the deep edge of the coseismic rupture and the afterslip located north of the rupture is the largest. Aside this overall behaviour, there are specific patterns that we discuss in the next paragraphs.

Afterslip within the coseismic rupture area

An average of ∼3-4 meters of coseismic slip is imaged by all coseismic models (Melgar et al., 2016;Shrivastava et al., 2016;Barnhart et al., 2016;Ruiz et al., 2016;Klein et al., 2017) at a 30-45 km deep patch. Our model infers significant afterslip at latitude 31.2 • S (∼30 km depth) reaching 40 cm, well inside the coseismic rupture (western half of Box B in Figure 5g and southwest of Box 2 in Figure 6a). Given its dimension (∼ 40 × 20 km 2 ), the vicinity of the coastline and the direct observation from dense GNSS data for both coseismic and afterslip models, this overlap is certainly resolved. In addition, that area hosts a significant density of aftershocks (Figure 6c, Frank et al. ( 2017)) and among them significant repeating earthquakes (Figure 6b, Huang et al. ( 2017)), providing additional evidences for aseismic slip at this location. Overlap of afterslip with coseismic slip is also found at shallow depth at the northern edge of the coseismic rupture, (southwestern part of Box A in Figure 5g and west of Box 1 in Figure 6a). This is an area where a relatively large distance from the coast prevents the precise location and amount of slip to be retrieved from the inversion, making this overlap speculative. We therefore conclude that although most of the afterslip occurred at the edge of the coseismic rupture area and most of the seismically rupture appears to have fully relocked quickly, a ∼ 40 × 20 km 2 region (western half of Box B in Figure 5g and southwest of Box 2 in Figure 6a) slip seismically and continued to slip aseismically during the weeks following the mainshock. The area of afterslip within the coseismic rupture is adjacent to the area found to have slipped seismically through large aftershocks during the first 6 hours following the mainshock by Twardzik et al. (2021). It might therefore also include the afterslip generated by these aftershocks. However, the time evolution of afterslip at that area (Box 2 in Figure 6e) does not show slower or faster decay with respect to the other areas of afterslip (Boxes 1 & 3, Figure 6d,f).

The observation of significant afterslip taking within the coseismic area of the mainshock is somehow at odd with the asperity model described in the introduction of the paper, where seismic rupture occurs at velocity-weakening areas of the fault, while afterslip witnesses velocity-strengthening friction. Both earthquake rupture dynamic simulations [START_REF] Noda | Stable creeping fault segments can become destructive as a result of dynamic weakening[END_REF][START_REF] Kaneko | Towards inferring earthquake patterns from geodetic observations of interseismic coupling[END_REF] and observations [START_REF] Harris | Large earthquakes and creeping faults[END_REF] show that seismic ruptures can propagate into velocity-strengthening areas of faults, through dynamic weakening induced by high shear stress as the seismic rupture front arrives at the edge of a velocity-strengthening area. This mechanism possibly happened during the Illapel earthquake. [START_REF] Aochi | Early stage and main ruptures of the[END_REF] proposed that the Illapel earthquake started with a foreshock and that the rupture immediately propagated northward, similarly as observed in the back-projection and the kinematic slip inversions (e.g. Melgar et al. (2016)).

The epicentral area experienced 3-4 m of slip during the first 30s of the rupture, before reaching the future afterslip area which broke between 30-40s, before propagating further northward and updip. We speculate that this area is velocity-strengthening that would explain the observed afterslip and that it resisted the seismic rupture propagation during the Illapel earthquake, as in the models proposed by [START_REF] Kaneko | Towards inferring earthquake patterns from geodetic observations of interseismic coupling[END_REF]. al., 2014), much smoother than the afterslip distribution imaged from our inversions (Fig-

ure 11). The deepest areas of afterslip correlate with the downdip limit of high interseismic coupling. The northern patch of afterslip also correlates with weaker interseismic coupling for our smoother model. However, the afterslip area that occurred within the seismically ruptured area (white rectangle on Figure 11) shows high coupling, perhaps because it lies in the stress shadow of both updip and downdip locked area. The southern patch of afterslip is located in a large area of high coupling. This area is located close to the coastline so it is expected to be well resolved for both our afterslip models and the interseismic coupling model. This observation is at odd with the asperity model. A similar observation has been made for the Ecuador M w 7.8 2016 earthquake for one patch of afterslip (Rolandone et al., 2018). In Rolandone et al. (2018), the patch that was locked during the years before the earthquake was found to host regular Slow Slip Events (SSE)

and Seismic Swarms (SS).

Along the Illapel segment of the Chilean subduction zone, no SSE has been geodetically evidenced yet, possibly due to a lack of a dense permanent monitoring network.

However, Poli et al. (2017) studied the seismicity during the 20 years preceding the Illapel earthquake and identified regular seismic swarms. Seismic swarms occur at three distinct areas at the megathrust showing a striking correlation with both our afterslip (Poli et al., 2017). The green dots are the repeater earthquakes locations of Poli et models and the aftershocks locations (Figure 12). The first area of seismic swarm at latitudes 30.6 • S-31 • S matches the patch of enhanced afterslip where the November 7 aftershock occurred. The second area at latitudes 31.3 • S-31.9 • S is at the edges of the afterslip in our models. Finally, the last one located at latitudes 31.9 • S-32.6 • S is at the periphery in our afterslip model. Poli et al. (2017) further show that repeating earthquakes occur during seismic swarms, witnessing associated SSEs. The location of the repeaters found Poli et al. (2017) is found to occur at the periphery of the afterslip areas of our rough geodetic afterslip model. Such correlation between the location of enhanced afterslip and areas experiencing SSEs and SS is very similar to the one observed for the Pedernales 2016 earthquake in Ecuador (Rolandone et al., 2018). It confirms a key-role of SSEs and SS areas in stopping the seismic rupture for large earthquakes (Vaca et al., 2018) and in promoting afterslip. Afterslip and SSEs occurring at the same location have also been observed at the Nicoya peninsula, Costa Rica (Hobbs et al., 2017).

A highly dynamic afterslip

While the cumulative afterslip evolution shows a smooth increase, the afterslip rate evolution shows higher dynamics in space and time. Our modeling suggests that the data not only carry information about the overall location of afterslip, but also includes information about the high frequency dynamics of afterslip. This high frequency dynamics is also suggested by the days-long modulation in the microseismicity rate (movies S4 and S6).

Large aftershocks appear to occur at area or at the edge of enhanced afterslip. These large aftershocks in turn trigger their own, local afterslip that superimposes to the overall mainshock-induced afterslip. Our approach made possible to identify such modulations, but with the limitation of the smoothing in time of the actual slip rate. Nonetheless, specific inversions at selected time windows were able to capture the local afterslip of the November 7 M w 6.8 and the two November 11 M w 6.9 earthquakes and their time evolution.

Aside aftershocks, our rough model highlights a small afterslip rate modulation south of the coseismic rupture. We suggest that this afterslip modulation is a small slow slip event (SSE). Indeed, the observed slip shows a progressive slip acceleration and decrease as SSEs do. The proposed magnitude of M w 6.1-6.2 for a duration of 5-6 days agrees with scaling laws proposed for SSE (Ide et al., 2007). Finally, it occurred at a location of the interface where regular seismic swarms, most probably associated with SSEs have been documented prior the 2015 Illapel earthquake (Poli et al., 2017).

These modulations remain small compared to the afterslip, and do not challenge the general view that the mainshock coseismic stress increment drives the overall pattern of afterslip. They however suggest that some other processes are superimposed to it. 10 km scale locked asperities locally resist to the aseismic slip around or at the edge of them, until the shear stress from both the mainshock coseismic slip and afterslip leads them to break. Then these aftershocks locally trigger their own afterslip, possibly helping the overall afterslip to propagate farther. Interpreting SSE is more challenging since they are often associated with velocity-weakening behaviour (e.g. Rubin (2008)). Perhaps, following Poli et al. (2017) suggestion made for the interseismic period, relatively -25-sudden pore pressure change due to fluid injection might provide an explanation to be explored in future.

Conclusions

A full time-dependent inversion of both cGNSS and sGNSS provides a kinematic view of the 74-days afterslip that developed during two months after the Illapel M w 8.3 megathrust earthquake. It highlights that afterslip developed at preferential areas at the subduction interface, which had experienced regular seismic swarms and likely slow slip events during the decades preceding the earthquake. Those areas possibly correspond to fractures within the slab that promote high pore-pressure and are favorable to aseismic slip. Such areas prevented the rupture to extend along strike at their locations. Afterslip and SSEs are usually thought to witness different friction laws at the megathrust.

Here, they appear to take place at the same location, certainly obeying the same friction law. This observation is not accounted for in the 'rate-state asperity' model, but it might be key in anticipating the extent of future larges ruptures.

Our detailed modelling further suggests a close link between afterslip and large aftershocks. 10 km scale locked asperities might be brought to failure by the afterslip developing in their vicinity, in addition to the stress increment induced by the mainshock.

This process might explain remote and delayed moderate size aftershocks as our results

suggest for the two M w 6.9 events on November 11, 2015. Relationship between a large aftershock and afterslip propagation was recently observed for a M w 6.9 earthquake in central Chile (Klein et al., 2021). Hence, survey of afterslip evolution after a large earthquakes might therefore help to define preferential location for large aftershocks to occur.

For the first time, our study images the afterslip rate evolution and highlights that the afterslip is highly dynamic, like the seismicity. GNSS sites are able to capture the afterslip of the large aftershocks and a slow slip event, south of the mainshock rupture.

This slow slip event correlates with an increase of microseismicity rate.

Those observations are similar to the one made after the M w 7.8 2016 Ecuador earthquake (Nocquet et al., 2017). Two large M w 6.7-6.9 aftershocks occurred a month after the mainshock in an area of enhanced afterslip having experienced regular seismic swarms and slow slip events prior the earthquake (Rolandone et al., 2018;Vaca et al., 2018). A similar pattern was also observed after the M w 7.6 2012 Nicoya earthquake at least updip the rupture (Dixon et al., 2014;Hobbs et al., 2017;[START_REF] Voss | Slow slip events in the early part of the earthquake cycle[END_REF]. These similarities suggest a common behaviour, perhaps corresponding to the seismic behaviour of the highly heterogeneous subduction interface.

Appendix: Supporting Information

Static resolution tests

In order to assess the robustness of the afterslip model, we perform several resolution tests. We started with static tests to estimate the spatial resolution of our model, using continuous and survey GNSS sites.

We build checkerboard models with different size of slip patches regarding our model: The static tests show that the GNSS network resolve all the 100 km × 100 km patches, until the trench. The resolution of smaller patches, i.e. 75 km × 75 km patches, is slightly lower than for larger ones. Onshore, it is still very good but offshore, we lose resolution between the trench and 15 km depth on almost the whole segment. However, resolution remains high between 32 • S and 33 • S. The data can solve smaller patches, i.e. 50 km × 50 km size, but resolution is significantly lower than the other tests. From 15 km to 60 km depth, the resolution is good except at La Serena bay (between 30 • S and 30.5 • S), where the coast is ≈ 100 km away from the trench. At shallower (< 15 km) and greater (> 60 km) depths, the GPS network is unable to solve the small asperities.

September 25, 2022, 10:48pm

Dynamic resolution tests

First, we assess the accuracy of the kinematic inversion to model a classical postseismic deformation. We made synthetic displacements corresponding to a homogeneous afterslip rate following a logarithmic decay over two weeks, surrounding as a circle arc the coseismic rupture. We noised the data with a normal distribution, with 1.4 mm on the horizontal components and 5 mm for the vertical component. We performed a cross-validation to select the optimal model, that is with regularization parameters λ space = 0.1 and λ time = 0.01 (Figure S5a). The kinematic inversion is able to retrieve with a good accuracy the synthetic cumulated afterslip distribution in space, although the shallowest parts of the afterslip -poorly constrained -show weak or no slip. The inversion found the maximum slip amount at the same areas as for the synthetic model, except at 31.5 • S. The maximum cumulative slip, highly dependent of the regularization parameters, is lower after the inversion (974 mm vs. 782 mm). The cumulated source time functions show that the kinematic inversion find the same geodetic moment (M w 7.5) for the synthetic model and its inverted result.

Second, we carried out a specific tests regarding our Illapel afterslip model. We assess the robustness of the afterslip pulse located on the southern patch on the last days of the kinematic inversion. To that end, we made synthetic displacements as a sum of 1) a homogeneous afterslip rate following a logarithmic decay over two weeks (like the first step) and 2) a slow slip event which maximum slip rate of 1500 mm/day occurred on day 9 and located on the southern part of the afterslip. We noised the data with a normal distribution, with 1.5 mm on the horizontal components and 5 mm for the vertical component. We performed a cross-validation to select the optimal model, that is September 25, 2022, 10:48pm

with regularization parameters λ space = 0.1 and λ time = 0.1 (Figure S5b). The kinematic inversion is able to model the spatial distribution of the synthetic cumulated slip and especially the southern extension of slip induced by the slow slip event. The maximum cumulated slip (991 mm for the synthetic model and 792 mm for the inverted result), located on the whole length of the circle arc is only found at three distinct areas by the kinematic inversion. However, the synthetic slip rate is modelled with great accuracy over the 14 days period, and above all the slow slip event slip rate is modelled with high precision in space and time by the inversion, as shown by the source time function. The geodetic moment found by the inversion (M w 7.5) is similar to the synthetic's.

The kinematic inversion can solve a classic afterslip rate both in time and space, with slight vagueness regarding the location of the maximum slip. While solving an afterslip with a modulation, the kinematic inversion is able to solve with great precision the modulation but with a lowest quality of the maximum slip location. We conclude that the modulation of the afterslip seen when the largest aftershocks occurred (Figure 5) is real and well resolved.

September 25, 2022, 10:48pm

Large aftershocks inversions

We model the coseismic slip induced by the largest aftershocks, using the same approach as described in the main text restricting the time period to two successive days around the aftershock date. In that case, no temporal smoothing applies. For the M w 6.8 aftershock that ruptured on November 7 th , we use the time series from 9 cGNSS sites (BTON, CMBA, CNBA, EMAT, LVIL, OVLL, EDR, PFRJ, SLMC) and discretize the subduction interface into 256, 5 km long edges, equilateral triangles following the Slab2 model (Hayes et al., 2018). We find a optimal spatial smoothing parameter of λ space = 0.1 from a L-curve.

For the two M w 6.9 aftershocks that ruptured on November 11 th , we use 7 cGNSS sites (BTON, CRZL, LSCH, OVLL, PEDR, PFRJ, TOLO) and discretize the megathrust into 1168 subfaults, 5 km long edges, equilateral. We find an optimal spatial smoothing parameter λ space = 0.4.

September 25, 2022, 10:48pm

Movie S1. Cumulative afterslip slip evolution contoured in orange from 20 cm and then every 10 cm. The black dots are the seismicity per day and the gray dots are the cumulative past seismicity. Slip inversion results for the largest aftershocks on November 7 and November 11 are contoured in brown every 5 cm. Their focal mechanisms is from Global CMT (Dziewonski et al., 1981;Ekström et al., 2012). The dashed lines are the Slab2.0 iso-depth contours every 20 km (Hayes et al., 2018). The mainshock slip is contoured in blue every 2 m (Klein et al., 2017). The blue star is the M w 8.3 2015 Illapel epicenter.

Movie S2. Daily slip rate at selected dates. The slip rate is contoured in green every 5 mm/day (the first 3 mm of slip are hidden). Black dots are the seismicity per day.

Slip inversion results for the largest aftershocks on November 7 and November 11 are contoured in brown every 5 cm. Focal mechanisms are from GlobalCMT (Dziewonski et al., 1981;Ekström et al., 2012). The dashed lines are the Slab2.0 iso-depths every 20 km (Hayes et al., 2018). The mainshock slip is contoured in blue every 2 m (Klein et al., 2017). The blue star is the M w 8.3 2015 Illapel epicenter.

Science, 362 (6410), 58-61.

Klein, E., Vigny, C., Fleitout, L., Grandin, R., Jolivet, R., Rivera, E., & Métois, M. (2017). A comprehensive analysis of the Illapel 2015 Mw8. The dashed lines are the Slab2.0 iso-depth contours every 20 km (Hayes et al., 2018).

September 25, 2022, 10:48pm (Hayes et al., 2018). The fault area extends from latitude 32 • S to 28 • S and from the trench down to 90 km depth. The dashed lines are the Slab2.0 iso-depth contours every 20 km (Hayes et al., 2018). On September 9, 2020, a M w 6.9 earthquake ruptured near the town of Totoral (at latitude 28 • S), along the Chilean subduction zone (Figure 4.1a). It was followed by a M w 6.3 earthquake twenty minutes later, itself followed by a M w 6.4 event later that day. This sequence is intriguing because of the relatively large and rapid occurrence of the aftershocks, compared to the magnitude of the mainshock. Such a scenario has less than 5% probability to occur (Klein et al., 2021). Indeed, according to the Omori law (Omori, 1894;Utsu, 1957) and the Reasenberg and Jones (1989) approach, the occurrence probability of aftershocks of M w ⩾ 6.3 within 24 hours after a mainshock of M w 6.9 is between 0.3 and 4%. Furthermore, [START_REF] Båth | Lateral inhomogeneities of the upper mantle[END_REF]'s law predicts that the difference in magnitude between the mainshock and its aftershock must average 1.2. However, in this sequence, the observed magnitude difference is only 0.5. The question then arose as to whether this sequence was classical (i.e. mainshock, aftershocks) or if a slow slip event induced the early and large aftershocks. In a classical sequence, the afterslip corresponds to about 25% of the coseismic moment [START_REF] Pritchard | An aseismic slip pulse in northern Chile and alongstrike variations in seismogenic behavior[END_REF]. But for some sequences departing from mainshock-aftershocks scheme, such as the 2009 northern Peru sequence, the aseismic slip represents more than 1000% of the coseismic moment. This process rather showed an interlaced seismic and aseismic sequence that lasted 7 months, with a predominant aseismic mode (Villegas-Lanza, Nocquet, et al., 2016). Therefore, 128 Chapter 4. Afterslip kinematic modelling of the 2020 Atacama sequence (Klein et al., 2021) documenting afterslip helps to identify if the sequence is classical, or if other processes are involved. Documenting sequences that are both aseismic and seismic is the objective of the Agence Nationale de la Recherche (ANR) project "Synchronous Slow Slip & Seismic Swarm" (S5) coordinated by Jean-Mathieu Nocquet (IRD, Géoazur & IPGP). An S5 process is defined as a slow slip event simultaneously accompanied by an intense increase in microseismicity, i.e. a seismic swarm. Several regions are targeted for this project, such as island La Plata in Ecuador, northern Peru, and the regions of Copiapó and La Serena in Chile. Documenting precisely these processes provides insights about the periodic stress redistribution induced along the subduction megathrust. Furthermore, S5 have been spotted before large earthquakes and could also be a precursor.

The occurrence of this sequence during my thesis was an opportunity to work on its postseismic phase. I interrupted my work on the Illapel postseismic phase and carried out this study, presented in this Chapter. This work was included in a collaborative study led by Émilie Klein (CNRS, ENS). In the first section, I introduce the seismo-tectonic context of the region before presenting the S5 already documented in Chile. I then present the 2020 Atacama sequence, and detail the work I have done on the associated postseismic phase. Finally, in the last section, I discuss the results and conclude this Chapter.

Seismo-tectonic context of the Atacama region

The Atacama region (i.e. latitudes 26 • S-29 • S) megathrust segment ruptured in 1819, through a M ∼ 8.3 earthquake. It experienced its largest major earthquake a century ago, in 1922 with a M w 8.4-8.7 event (Ruiz & Madariaga, 2018). The Copiapó ridge subducts beneath the South America plate at latitude 27.3 • S and has a direction identical to the plate convergence motion (Figure 4.1a).

The interseismic coupling model of the region is heterogeneous (Figure 4.1b, Metois et al. (2016). It is characterized by two highly coupled segments, the Chañaral and the Atacama segments, separated by a Low Coupling Zone (LCZ) called the Baranquilla LCZ. The Chañaral segment (latitudes 27 • S-27.5 • S) shows an intermediate coupling from the trench down to 30 km depth. At 30-45 km depth, Figure 4.1b shows a highly coupled area, and the 60-90 km depth range shows a low to intermediate coupled area. The Baranquilla LCZ (latitudes 27.5 • S-28.1 • S) is a low coupled area except for depths between 45 and 65 km. The Atacama segment (latitudes 28.1 • S-29.1 • S) shows a large highly coupled area between 15 and 45 km depth. Downdip, the segment is uncoupled for 15 km.

The Atacama region, and more precisely the Copiapó area (latitude 27 • S) experienced several seismic and aseismic events before the 2020 sequence. These alternating (Dziewonski et al., 1981). The black triangles represent the GNSS sites used in this study. The red ellipse show where the 1973,1976,2006 and 2015 seismic swarms occurred (Comte et al., 2002[START_REF] Comte | Copiapo Chile seismic swarm analysis: mapping the interplate contact[END_REF]Holtkamp et al., 2011). The blue curves indicate the Copiapó 2014 slow slip every 10 cm (Klein, Duputel, et al., 2018). The orange line indicates the proposed lateral extent of the historical 1819 and 1922 megathrust earthquakes (Lomnitz, 2004;Ruiz & Madariaga, 2018). b) Interseismic coupling model from Metois et al. (2016). The grey lines delineate three segments. The dashed lines are the Slab2.0 iso-depth contours every 20 km (Hayes et al., 2018). seismic and aseismic processes (S5) are also thought to regularly occur in the La Serena area (latitude 30 • S). Both these areas show heterogeneous interseismic coupling (Metois et al., 2016), which raises the question of whether a link between heterogeneous interseismic coupling zone and S5 can be established. The Copiapó and La Serena S5 are described in section 4.2.

Seismic and aseismic processes (S5) in Chile

The Chilean subduction zone is approximately 3500 km long, making a good GNSS coverage challenging. This could explain the poor documentation of aseismic events. One of the objectives of the S5 project is to densify the seismological and GNSS networks in areas where S5 might occur to refine the analysis of S5 processes. The Atacama region is well monitored thanks to the efforts of Christophe Vigny and Émilie Klein (CNRS, ENS) for the GNSS network (CNRS-INSU Tellus program) and Zacharie Duputel (IPGP, Volcanological observatory of Piton de la Fournaise) for the broadband seismological network (European Research Council starting grant PRESEISMIC). The Centro Sismológico Nacional (CSN, Universidad de Chile, Baez et al. (2018)) also made it possible to densify the network and maintain it. In sections 4.2.1 and 4.2.2, I present the S5 that took place in the Copiapó and La Serena areas.

The Copiapó area

The Copiapó region corresponds to the megathrust segment located at latitudes 27-28 • S. This segment experienced four seismic swarms during the last 50 years and, as far as we know, one slow slip event in 2014. These processes occurred in the heterogeneous coupled Chañaral segment (Figure 4.1b).

A slow slip event was spotted in Copiapó (at latitudes 26.8 • S-27.5 • S) between 40 and 60 km depth, downdip of the seismogenic zone (blue contours in Figure 4.1a, Klein, Duputel, et al. (2018)), at the transition between the coupled and the uncoupled part of the Chañaral segment. It lasted almost two years, between September 2014 and June 2016 and was equivalent to a M w 6.9 earthquake. Transient signals were also spotted in 2005 and 2009 on a single time series, lasting approximately one year. Klein, Duputel, et al. (2018) therefore proposed a recurrence time of approximately 4-5 years. Assuming that the 2009 slow slip event ruptured the same area as the 2014 event, released all the stress, and considering the interseismic coupling of the region Metois et al. ( 2016); [START_REF] Klein | Bridging the gap between North and Central Chile: insight from new GPS data on coupling complexities and the andean sliver motion[END_REF], the subduction megathrust could not have accumulated the stresses necessary to produce the 2014 slow slip event.

The four seismic swarms of the Copiapó area occurred 80 km west of Copiapó, between 10 and 30 km depth along the subduction megathrust (red shaded area in Figure 4.1a, Comte et al. (2002, 2006); Holtkamp et al. (2011)). The first one was discovered in 1973, with at least 72 earthquakes. It occurred where the 1922 M w 8.4-8.7 earthquake ruptured. Six years later, in 1979, a smaller seismic swarm occurred, with fewer earthquakes of lower magnitudes. Another swarm took place in 2006 and is the largest documented in the area (180 earthquakes). The last seismic swarm occurred in 2015, but so far no study has been published on this event.

The Copiapó area, with its heterogeneous interseismic coupling, therefore experienced recurrent seismic swarms as well as one and probably more slow slip events. The La Serena area shares some of these characteristics.

The La Serena area

The La Serena area is located south of Copiapó, at latitudes 27.5 • S-30.5 • S. The interseismic coupling at the shallow part of the megathrust in this area is not well resolved due to the large distance of the GNSS sites from the trench (i.e. approximately 80-100 km).

The La Serena area has experienced significant seismicity three months before the 1997 M w 7.3 Punitaqui intrasblab earthquake, a seismicity that has not decreased since [START_REF] Lemoine | Evidence for earthquake interaction in Central Chile: the July 1997-September 1998 sequence[END_REF][START_REF] Pardo | The October 15, 1997 Punitaqui earthquake (Mw=7.1): a destructive event within the subducting Nazca plate in central Chile[END_REF]Pardo, Comte, & Monfret, 2002). [START_REF] Gardi | Modeling of stress transfer in the Coquimbo region of central Chile[END_REF] suggest that aseismic slip could be at the origin of this increased seismicity. The displacements recorded by GNSS sites for 3 years showed that the megathrust in the La Serena area is partially coupled while most of the stress is released through aseismic slip processes [START_REF] Vigny | Upper plate deformation measured by GPS in the Coquimbo Gap, Chile[END_REF]. More recent interseismic coupling models showed that La Serena is a Low Coupling Zone (LCZ) where the 2015 M w 8.3 Illapel earthquake stopped. Therefore, the La Serena LCZ is likely to experience shallow slow slip events, hardly detectable due to the low data resolution in the area (Metois et al., 2016). North of the La Serena LCZ, the megathrust is highly coupled at shallow depths (5-40 km) and uncoupled at greater depths (below 40 km).

With its heterogeneous interseismic coupling and the likely occurrence of aseismic slip, the La Serena area is conducive to seismic swarms and probably to processes that are both seismic and aseismic. The GNSS network should be densified in the future to search for slow slip events. We now return to the Copiapó area to discuss the 2020 Atacama sequence.

The 2020 Atacama sequence

The September 2020 Atacama sequence is unusual because of the large size and timing of the aftershocks, in comparison with the mainshock. In this section, I summarize the study conducted on this sequence. It should be noted that, given the 20 minutes between the M w 6.9 mainshock and the first large aftershock (M w 6.3), it was not possible to differentiate them in the GNSS data. Therefore, the study first focused on the mainshock and the second large aftershock (M w 6.4) that occurred sixteen hours later. Seismological and geodetic data were used to investigate this sequence. A catalogue of earthquakes was built and magnitudes were estimated. From these results, the probability of the M w 6.3 aftershock occurrence within the first 24 hours after the mainshock was estimated to range from 0.3% to 4%. The mainshock dynamic rupture was modelled using the seismic waveforms. Figure 4.2c shows the dynamic slip model of the mainshock. With a maximum slip of 1.1 meters, the rupture extended on a surface of 25 × 25 km 2 and the model is equivalent to a M w 6.7.

High-rate GPS was calculated to extract the coseismic displacements of the mainshock and the M w 6.4 aftershock. It was subsequently used to model the coseismic slip of these events using static inversions. Figure 4.2b,c shows the coseismic slip distributions of the mainshock and the M w 6.4 aftershock. The mainshock ruptured a 80 × 40 km 2 area between 20 and 40 km depth, a rather large area for a M w 6.9 earthquake. The maximum slip is reached at latitude 28 • S and is equal to 1.3 meters. The M w 6.4 aftershock ruptured updip of the mainshock, at ±20 km depth and latitude 28 • S. The coseismic model shows a maximum slip of almost 1 meter and a rupture area of 30 × 30 km 2 . The updipward propagation, between the mainshock rupture area and the second aftershock rupture area, suggests that another phenomenon plays a role in the sequence. To investigate this, the postseismic slip of the mainshock, until the occurrence of the second aftershock, was modelled.

To obtain the postseismic displacements of the mainshock, the high-rate GPS and the daily GNSS solutions were used. The coseismic displacements estimated from the high-rate GPS of the mainshock and the second aftershock were first summed, and subsequently removed from the total displacements estimated from the daily GNSS time series between before the mainshock and after the second aftershock. The obtained displacements correspond to the coseismic displacements of the first aftershock (M w 6.3) and to the postseismic displacements of the mainshock. They have been modelled using a static inversion and are shown in Figure 4.2d. The slip is embedded between the two coseismic ruptures (mainshock and second aftershock) and reached more than 80 centimeters. This result confirms that large afterslip occurred updip the mainshock rupture, leading rapidly to the occurrence of the second large aftershock.

The postseismic displacements following the second aftershock were modeled using a time-dependent inversion of the daily GNSS time series, as for the Illapel postseismic deformation. The next section detailed this work that I have done.

4.4

The time-dependent inversion of the postseismic deformation of the 2020 Atacama sequence

Time series and methodology

In this study, the data set is composed of 12 continuous GNSS time series, shown in Figure 4.1a. I used the time series of the processing solution from Klein, Vigny, Nocquet, and Boulze (2022). I removed the interseismic velocity from this solution, estimated using the MIDAS approach of Blewitt et al. (2016) over a period of three years preceding the sequence.

The observed time series after the sequence show logarithmic-like increasing trenchward postseismic displacements on the east component over almost 25 days (Figure 4.3). The site BAR2, located approximately 80 km north of the sequence (located on Figures 4.1a and 4.4), shows 13 mm of westward displacements after three weeks. The north By inspecting the time series one by one, I decided to perform an inversion including 22 days, starting on the day after the second aftershock (September 2, 2020) until September 24, 2020. I followed the methodology described in section 2.3. A difference here is that we used a covariance matrix for the regularization of the inverse problem, instead of the Laplacian regularization. The covariance regularization is described in section 2.3.2.2 and Equation 2.10. Several models were calculated with different smoothing parameters and correlation distances. The model achieving the best fit to the time series was obtained with a smoothing parameter of σ = 20 mm. √ day (average slip) and a correlation distance between subfaults of D corr = 35 km. In the following section, I describe the results obtained.

Results

Figure 4.5a-f shows the cumulative slip and the seismicity evolution at selected dates. The afterslip initiates where the second aftershock ruptured, at approximately 20 km depth (Figure 4.5a). It propagates northward, where the mainshock ruptured. On September 5, after three days, the afterslip reaches more than 15 mm (Figure 4.5b). On September 7, five days after the second aftershock, a deeper afterslip patch initiates north of the first one and at greater depth (Figure 4.5c). Interestingly, this patch overlaps the area of the 2014 Copiapó slow slip event (Klein, Duputel, et al., 2018). Slip is required at this location since several GNSS sites record a displacement of about 5 mm. Besides, the agreement between observed and modelled displacement is good at this location (see Figure 4.6). From this day onward, the shallow and deep patches of afterslip develop independently. After 22 days, the shallow patch reaches almost 60 mm of slip, while the deep one accumulates 15 mm of slip (Figure 4.5f). The deep patch grew at a significantly slower rate than the shallow one. The seismicity is located on the southern edge of the shallow afterslip patch, between 10 and 30 km depth, and curiously did not migrate during the 22 days (black and grey dots in Figure 4.5). The 22-day cumulative afterslip model released a moment of 1.23 × 10 19 N.m, equivalent to a M w 6.8 event (by taking a shear modulus of 30 GPa). Figure 4.7a-f shows the daily afterslip rate and the seismicity rate (black dots) at selected dates. On the first day (Figure 4.5a), the shallow patch shows more than 6 mm.day -1 of slip rate over a 60 × 30 km 2 patch. It is the largest in terms of size and slip rate. Indeed, the deeper patch shows 1 mm.day -1 of slip rate. The afterslip rate and the seismicity rate rapidly decrease over the first few days, especially at the shallow patch, from September 8 onwards (Figure 4.5c), the shallow patch shows slip rate of 3 mm.day -1 whereas the deep one still shows 1 mm.day -1 of slip rate, but on a smaller area. As a consequence of the decreasing slip rate in the shallow patch, the seismicity decreases. The deeper patch of slip disappears on September 15 (Figure 4.5e), while the shallow one disappears after September 17, (Figure 4.5f). Considering that they did not have the same slip rate (6 mm.day -1 versus 1 mm.day -1 ), the afterslip rate decay is significantly slower for the deep patch. The seismicity rate decreases for the most part over the first six days, and slightly decreases with time afterwards. Surprisingly, it does not follow the spatial evolution of the afterslip and remains in the same small region (i.e. between 10 and 30 km depth and at latitude 20 • S). Figure 4.8 shows the source time function of the afterslip rate model. We can observe quasi-steady moment rate decrease through time rather than the 1/t typical decay usually observed for afterslip. After this description of the time-dependent models of the afterslip, I discuss the 2020 Atacama sequence in the context of the South American subduction zone, then I discuss its relationship with the 2014 deep slow slip event. The study of the Atacama sequence contributes to documenting the diversity of processes by which stresses are released along the South America subduction zone. Large earthquakes release stress accumulated at large velocity-weakening asperities and generate afterslip surrounding their rupture, together with aftershocks that rupture small to intermediate size velocity-weakening asperities, as we described for the Illapel earthquake in Chapter 3. This process also occurs more frequently for moderate-sized earthquakes (i.e. M w 5-7 events). Because afterslip generated by those earthquakes is small, they appear as mainshock-aftershocks sequences in seismological data. The South American subduction zone also experiences slow slip events. At shallow depth (less than 45 km depth), slow slip events are synchronous with intense microseismicity, as for example the 2010 and 2013 La Plata slow slip events in Ecuador [START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF][START_REF] Segovia | Seismicity distribution near a subducting seamount in the Central Ecuadorian subduction zone, space-time relation to a slow-slip event[END_REF]. We define this process as S5 for Synchronous Seismic Swarm and Slow Slip. During S5, seismicity accounts for at most a few percent of the total moment released and as a consequence slow slip appears to drive the evolution of the process. Deep slow slips (> 40 km depth) were found at Copiapó in Chile (Klein, Duputel, et al., 2018) and in central Ecuador (Rolandone et al., 2018). Neither microseismicity nor tremors were found, also witnessing that aseismic slip appears to dominate the process.

Discussion

The 2020 Atacama sequence does not belong to any of these two classes of stress release process, but shares some characteristics of the mainshock-aftershocks sequence and some of the S5 process. The Atacama sequence mixes ruptures of medium-sized velocity-weakening asperities and aseismic slips. Specifically, seismic ruptures and aseismic slips interact and trigger each other. The contribution of aseismic slip to the overall moment released during the sequence is significant. The spatial and temporal organization of seismic and aseismic slips (see Figure 4.9) shows that the mainshock triggered the afterslip, which migrated updip, where it triggered the second aftershock. In turn, this second aftershock triggered its own afterslip, added to that of the mainshock and the first aftershock. This 22-day afterslip released a moment equivalent to a M w 6.8, which represents 80% of the M w 6.9 mainshock. This is abnormally large since usually we expect about 25% of the coseismic moment after one month [START_REF] Pritchard | An aseismic slip pulse in northern Chile and alongstrike variations in seismogenic behavior[END_REF]. The entire sequence released a moment equivalent to a M w 7.1, of which 60% is seismic and 40% is aseismic. Regarding the spatial extent of slip, the aseismic part appears significantly larger than that of seismic. The aseismic slip is therefore almost as important as the seismic slip in this sequence.

The 2020 Atacama sequence fuels a growing set of similar intermediate-sized mixed seismic-aseismic sequences along the South American subduction zone. In Peru, a 7month sequence of medium-sized seismic ruptures and aseismic slip occurred in 2009 (Villegas-Lanza, Nocquet, et al., 2016). The aseismic slip following the second earthquake of the sequence released a moment 10 times that of the earthquake moment and can therefore not be considered as afterslip anymore. The whole sequence was predominantly aseismic, where the seismic moment was 25% of the overall moment. In Chile, a similar sequence occurred in 2006 near Copiapó, but at shallow depth [START_REF] Ojeda | Seismic and aseismic slip during the 2006 copiapó swarm in northcentral chile[END_REF]. This 1.5-month sequence hosted aseismic slip, seismic swarm and two medium-sized seismic ruptures. The aseismic slip was also predominant and released 60% of the whole sequence moment. As for the 2020 Atacama afterslip and the 2009 Peru afterslip, the 2006 Copiapó afterslip was very large compared to the coseismic ruptures. It is thus a new type of sequence, which presents a certain pattern, with the mutual interaction between large aseismic slips and ruptures of some medium-sized earthquakes. Finally, the afterslip of these sequences is abnormally large compared to the seismic ruptures.

Possible interaction between shallow and deep aseismic processes

An interesting feature of the 22-day afterslip I derived for the Atacama sequence is that, despite the sequence having occurred at shallow depth, a patch of aseismic slip also developed at depth, about 70 km from the area of the seismic slip. This area overlaps with that of the deep slow slip event documented by Klein, Duputel, et al. (2018). After the Atacama 2020 paper was published, [START_REF] Klein | Return of the Atacama deep slow slip event: The 5-year recurrence confirmed by continuous GPS[END_REF] found that a slow slip event started in March 2020 on the southern part of the 2014 event slip area, where the deep patch of afterslip occurred (see Figure 4.9a,b). The slow slip developed from March 2020 until the September 2020 Atacama sequence. The occurrence of the Atacama sequence probably disrupted the slow slip event evolution. This new result raises the question to know if the deep patch of afterslip, with its particular dynamics (see section 4.4.2, second paragraph), could actually be the continuation of this slow slip event. Possibly, due to the increment of stress induced by the sequence, the slow slip might have accelerated and this acceleration is captured by the time-dependent modelling of the afterslip. It would be interesting to investigate whether the slow slip event played a role in the occurrence of the sequence, since the slow slip event and the Atacama sequence occurred within about 75 km of each other. Performing a time-dependent inversion between March 2020 to the end of the sequence and study the seismicity over the same period would allow to know if the slow slip event triggered the sequence. Finally, a time-dependent modelling of the slip between 2014 and 2020 would enable a further study of the 2014 and the 2020 slow slip events and to know if it is not in fact a single one.

This study demonstrates the very dynamic behavior of a heterogeneous subduction megathrust and the particularly decisive role of aseismic slip in the occurrence of earthquakes.

The detailed study of the 2020 Atacama sequence was published in the July 2021 issue of the Earth and Planetary Science Letters journal. The paper is in Appendix 4.6. An earthquake sequence occurred in the Atacama region of Chile throughout September 2020. The sequence initiated by a mainshock of magnitude M w = 6.9, followed 17 hours later by a M w = 6.4 aftershock. The sequence lasted several weeks, during which more than a thousand events larger than M l = 1 occurred, including several larger earthquakes of magnitudes between 5.5 and 6.4. Using a dense network that includes broad-band, strong motion and GPS sites, we study in details the seismic sources of the mainshock and its largest aftershock, the afterslip they generate and their aftershock, shedding light of the spatial temporal evolution of seismic and aseismic slip during the sequence. Dynamic inversions show that the two largest earthquakes are located on the subduction interface and have a stress drop and rupture times which are characteristics of subduction earthquakes. The mainshock and the aftershocks, localized in a 3D velocity model, occur in a narrow region of interseismic coupling (ranging 40%-80%), i.e. between two large highly coupled areas, North and South of the sequence, both ruptured by the great M w ∼8.5 1922 megathrust earthquake. High rate GPS data (1 Hz) allow to determine instantaneous coseismic displacements and to infer coseismic slip models, not contaminated by early afterslip. We find that the total slip over 24 hours inferred from precise daily solutions is larger than the sum of the two instantaneous coseismic slip models. Differencing the two models indicates that rapid aseismic slip developed up-dip the mainshock rupture area and down-dip of the largest aftershock. During the 17 hours separating the two earthquakes, micro-seismicity migrated from the mainshock rupture area up-dip towards the epicenter of the M w 6.4 aftershocks and continued to propagate upwards at ∼ 0.7 km/day.

The bulk of the afterslip is located up-dip the mainshock and down-dip the largest aftershock, and is accompanied with the migration of seismicity, from the mainshock rupture to the aftershock area, suggesting that this aseismic slip triggered the M w = 6.4 aftershock. Unusually large post-seismic slip, equivalent to M w = 6.8 developed during three weeks to the North, in low coupling areas located both up-dip and downdip the narrow strip of higher coupling, and possibly connecting to the area of the deep Slow Sleep Event detected in the Copiapo area in 2014. The sequence highlights how seismic and aseismic slip interacted and witness short scale lateral variations of friction properties at the megathrust.

Introduction

The Atacama region (26 • S-30 • S) is one of the long lasting seismic gaps of Chile (Lomnitz, 2004;[START_REF] Métois | Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 • S)[END_REF]Ruiz and Madariaga, 2018). In this region, the last megathrust earthquake occurred in 1922, a M w 8.6 event that stroke North-Central Chile and triggered a transpacific tsunami [START_REF] Willis | Studies in Comparative Seismology: Earthquake Conditions in Chile[END_REF][START_REF] Beck | Source characteristics of historic earthquakes along the central Chile subduction zone[END_REF][START_REF] Soloviev | Catalog of Tsunamis on the Western Coast of the Pacific Ocean[END_REF]Ruiz and Madariaga, 2018;[START_REF] Kanamori | New constraints on the 1922 Atacama, Chile, earthquake from historical seismograms[END_REF]. After 1922, the largest earthquake that occurred in the area was in 1983 with a magnitude 7.7 [START_REF] Pacheco | Seismic moment catalog of large shallow earthquakes, 1900 to 1989[END_REF][START_REF] Comte | Crustal seismicity and subduction morphology around Antofagasta, Chile: preliminary results from a microearthquake survey[END_REF]. More recently, in 2013, an event of magnitude 6.8 located around 50 km depth occurred, probably at the bottom-end of the seismogenic zone along the plate interface. A decade of survey GPS measurements conducted in this region revealed two large highly coupled zones, the Atacama and the Chañaral segments, separated by a relatively large intersegment of intermediate to low coupling, named the Baranquilla low coupling zone (LCZ) [START_REF] Métois | Revisiting the North Chile seismic gap segmentation using GPS-derived interseismic coupling[END_REF][START_REF] Métois | Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 • S)[END_REF]Klein et al., 2018a). Additionally, a 1.5 year-long, M w ∼7, slow slip event (SSE) was also detected in the region in 2014, but occurred deeper (40-60 km) than usual seismogenic depths (10-40 km) (Klein et al., 2018b). A detailed analysis of the only continuous GPS site in the region at this time also revealed two episodes of transient deformation, prior to the 2014 event, in 2005 and 2009, suggesting a possible recurrence of about 5-years for deep slow slip events in the region (Klein et al., 2018b).

Here, we study a large seismic sequence that occurred in the Atacama region throughout September 2020 (Fig. 1), South of an area where seismic swarms have occurred several times in the past, i.e. in 1973, 1976 and 2016, offshore the town of Caldera (27 • S, Fig. 1, Comte et al., 2002;Holtkamp et al., 2011). The 2020 sequence initiated on September, 1st, at 04:09 UTC, with an earthquake of magnitude 6.9. It was followed 20 minutes later by an event of magnitude 6.3, close to the mainshock epicenter and 17 hours later, at 21:09 UTC, by another event of magnitude 6.4, the largest aftershock of the entire sequence, 20 km updip the mainshock. Overall, the sequence lasted several weeks with more than a thousand events and includes several large earthquakes of magnitude larger than 5. We use a complete set of seismological sites deployed in the area prior to the sequence that includes broad-band, strong motion and GPS to monitor the spatio-temporal evolution of this sequence (Fig. 1). Thanks to this dense network we greatly improve the threshold detection down to magnitude 1 (with a magnitude of completeness of 2.5) and the precision of the localization through a 3-D refined velocity model. Focusing on the first day, we compare the high rate and the daily GPS solutions to quantify the amount of seismic and aseismic deformation that took place after the mainshock. Finally, we discuss how this sequence takes place in the earthquakes history of this area and how it may alter the potential seismic hazard of the nearby highly coupled zones.

Seismic analysis

The Atacama region is poorly covered by the national seismic network (CSN, Centro Sismológico Nacional, University of Chile, Santiago) with only 2 broad-band stations at less than 100 km from the sequence. Since 2013, less than 2000 earthquakes have been located in Chile between latitudes 30 • S and 26 • S.

Building the sequence catalog

We built a catalog using data from 14 broad-band stations of the CSN in a 300 km radius around the sequence, completed by data from three semi-permanent stations in the Copiapo region (30-150 km North), 10 temporary stations between Vallenar and Ovalle (100-300 km South), and 30 stations from the national strong-motion network of the CSN [START_REF] Barrientos | The seismic network of Chile[END_REF][START_REF] Leyton | Geophysical characterization of the Chilean seismological stations: first results[END_REF] providing data only for the 16 largest events (Fig. 1). Event detection was performed by STA/LTA method using the six closest stations, with two constraints: firstly, one of the three closest stations had to be first in triggering a detection and secondly, each event had to trigger detections on at least 5 stations to be considered. These criteria geographically restricted the area of study and filtered out the smallest local events and the hundreds of earthquakes happening everyday in Chile. Between the 25th of August and the 25th of September included, 1354 events have been detected, out of which 50% happened within the first four days of the sequence. No significant raise in seismic activity was detected prior to the main event: 1 to 9 events/day occurred between the 25th and the 31st of August (Fig. 1-B). Manual P-and S-wave arrival-times readings were performed, leading to 916 earthquake locations out of which 843 events belong to the dense core of the sequence and 74 correspond to surrounding activity that may or may not be related to the sequence. Specifically, half of these (35 events) occur up North in a 80 x 80 km 2 area, 11 are located further inland, 15 are poorly located beyond the trench or very deep below the contact and the last 13 correspond to quarry blasts.

Earthquake locations were determined by a double-difference approach in a regional 3D velocity model obtained by regional tomography [START_REF] Potin | Tomography of Chile[END_REF]. Fig. 2 represents a trenchperpendicular vertical cross-section across the sequence, with Pwave velocities and P-over S-wave velocity ratios (based on earthquakes arrival times, [START_REF] Potin | Tomography of Chile[END_REF]. The seismicity associated with the sequence is located at the interface, mainly between 15 km and 40 km deep, with some events scattered within the first 15 km of the upper plate. The background seismicity visible on Fig. 2, located within a 50 km range on both sides of the cross- section, appears to extend within the plunging oceanic plate and can be interpreted as the double seismic zone observed in several places along the Chilean coast [START_REF] Bloch | High-resolution image of the North Chilean subduction zone: seismicity, reflectivity and fluids[END_REF][START_REF] Bloch | Watching dehydration: seismic indication for transient fluid pathways in the oceanic mantle of the subducting Nazca slab[END_REF][START_REF] Comte | An inverted double seismic zone in Chile: evidence of phase transformation in the subducted slab[END_REF][START_REF] Sippl | Seismicity structure of the northern Chile forearc from >100,000 double-difference relocated hypocenters[END_REF], although these events are poorly located due to the lack of local observations. P-wave velocities and P-over S-wave velocity ratios for the upper plate, the interface and the upper oceanic mantle are consistent with others local tomographic models obtained in northern Chile [START_REF] Pastén-Araya | Fluids along the plate interface influencing the frictional regime of the Chilean subduction zone, northern Chile[END_REF], 2021).

Fig. 3 shows the spatio-temporal evolution of the seismicity over the first 72 h following the mainshock. Immediately after the mainshock, seismicity spread over a 20 x 20 km 2 region, a size roughly consistent with the rupture area (Fig. 3-C and 3-D). This initial spatial extension shows the area of influence of the stress increase due to the mainshock. Throughout these first 72 hours, both the North-South and Eastern (downdip) boundaries of seismicity remain stable. On the contrary, seismicity slowly spreads updip (westward), with an average velocity of approximately ∼ 0.7 km/hour (considering a dip of 20 • ; red dashed line on Fig. 3-C), resulting in almost doubling the initial area of aftershock.

Moment magnitudes

To constrain the magnitude of the largest events of the sequence, we perform regional W-phase source inversions [START_REF] Duputel | W phase source inversion for moderate to large earthquakes (1990-2010[END_REF] combined with a bootstrap analysis [START_REF] Efron | An Introduction to the Bootstrap[END_REF]. We use broad-band velocimetric data from the Federation of Digital Seismic Networks (FDSN) (C, C1 (doi .org /10 .7914 /SN/ C1), CX (doi :10 .14470 /PK615318), G (doi :10 .18715 /GEOSCOPE .G), GT (doi .org /10 .7914 /SN /GT) and II (doi .org /10 .7914 /SN /II) networks) within 26 degrees of epicentral distance. To improve the homogeneity of the data coverage, we select one station per cell in a 100 km × 100 km grid in the vicinity of the source. The used time window starts at the P-wave arrival time. Its duration is 300 s for epicentral distances smaller than 12 • and grows with distance (15 × s/•) for farther stations. Waveforms are filtered using a frequency band-pass that varies with the Global CMT magnitude. Here we filter between 50-80 s and 120-250 s. The average M w and ±2σ uncertainties are: 6.87 ± 0.07, 6.29 ± 0.04, 6.42 ± 0.07 for the events that occurred on 2020/09/01 at 04:09 UTC, 04:30 UTC and 21:09 UTC. The bootstrap histograms are shown on Fig. 4 and estimated parameters are gathered in the supporting information.

Characterization of the sequence: Mainshock-Aftershock sequence or Seismic swarm?

To evaluate the difference of the 2020 Atacama seismicity compared to a standard mainshock-aftershock sequence, we analyze earthquake sizes and temporal distribution in the area. Considering seismic events in the epicentral area since 2017 in the CSN catalog, we estimate a b-value of b = 0.8 ± 0.2 using the [START_REF] Aki | Maximum likelihood estimate of b in the formula log N = a-bM and its confidence limits[END_REF] approach (consistently, we estimate b = 0.7 ± 0.1 for the 2020 Atacama sequence using the catalog presented in section 2.1). The time of aftershocks relative to the M w = 6.9 mainshock is consistent with the Omori-Utsu law r(t) = K (t + c) -p with p = 1.0, c = 0.1 days and K = 16.3 (see Fig. S1 of the Supporting Information; Omori, 1894;Utsu, 1957). Looking independently into the magnitude and temporal distribution of the earthquakes, the sequence does not seem different from a classical mainshock-aftershock sequence. However, what seems anomalous is the occurrence of two M w > 6 aftershocks within 24 hours after the mainshock. Using a simple approach similar to Reasenberg and Jones (1989), we forecast the number of aftershocks of magnitude M w ≥ 6.3 within 24 hours after the mainshock using b = 0.8 and Omori-Utsu parameters mentioned above. Results shown in Fig. S1 indicate that there is only a probability of 0.3% of having at least two aftershocks of magnitude M w ≥ 6.3 shortly after the mainshock. However, this estimate depends on the assumed b-value. If we consider b = 0.7 as for the Atacama sequence, the aforementioned probability increases to about 4%.

GPS data analysis

Early 2019, in order to densify the CSN network [START_REF] Báez | The Chilean GNSS network: current status and progress toward early warning applications[END_REF], we installed 5 continuous GPS (cGPS) stations in the Atacama region. Three of them were collocated with broad-band seismometers (see section 2.1). Overall, we benefit from 12 cGPS stations located in the area of the sequence complemented by 5-6 stations further away for the reference (Fig. 1). In this study, we use both the stations positions obtained from 24-hours daily solutions throughout the whole duration of the sequence and the high rate (1 Hz,hereafter HRGPS) data that allow to decipher the successive displacements during the first day.

24-hours daily solutions

In addition to the data from the national Chilean network [START_REF] Báez | The Chilean GNSS network: current status and progress toward early warning applications[END_REF] and from the 5 additional stations, we include data of the Argentinian RAMSAC network (Piñón et al., 2018) and of the Brazilian RBMC network. We also include all the IGS stations available on the South American continent. This dataset is processed using the GAMIT/GLOBK software following the classical MIT methodology (Herring et al., 2010a,b). In a second step, we produce daily time series by constraining continental stations to their well-known coordinates in the ITRF2014 [START_REF] Altamimi | Itrf2014 plate motion model[END_REF] with the PYACS toolbox.

A specific difficulty needs to be addressed when several large earthquakes occur during the same day. If a single coordinate is calculated for the entire day, it will end up being anywhere between the pre-and post-earthquakes coordinate, depending on different parameters: when exactly the earthquakes occur during the day, which data segment (before, between and after the earthquakes) is the longest, and how the filter will handle data that does not fit the obtained average position of the day. In order to eliminate the pre-seismic observation (before 4:30 UTC) and to separate the two events in the data (see Fig. S2), we consider at which time the two main events occurred (M w = 6.9 at 4:09 UTC and M w = 6.4 at 21:09 UTC) and the day of the earthquake was processed using only the observations acquired between 4:30 UTC and 21:00 UTC. Therefore, this day's position corresponds to an averaged position of the station after the first event (M w = 6.9) and before the second event (M w = 6.4). Note that the selected time window also allows us to exclude the M w = 6.3 aftershock. Because only 25 min separates this aftershock from the mainshock, the potential deformation generated by this aftershock is most likely impossible to differentiate from the mainshock, using daily GPS solutions.

Time series reveal significant displacements on at least 7 stations (Fig. 5). Steps between days 244 and 245 (resp. 245 and 246) correspond to the coseismic displacements generated by the first (resp. the second) event, both occurring during day 245 (September 1st). The typical curvature of the time series of the stations nearest the events during the days (possibly weeks) following the mainshock also reveals postseismic deformation. This deformation seems unusually large (∼100% in only a couple of days) at the nearest station (TTRL). The estimation of the coseismic displacement of the first event of M w = 6.9 at 04:09 (the mainshock) is obtained by differentiating between the position at midday 245 (between 4:30 and 21:00 UTC) and the position of the day before (244) (Fig. 6-A, vectors in light red). It includes part, but not all, of the post-seismic deformation occurring during the 15 hours time span between the mainshock and the large aftershock at 21:09, which is potentially a combination of rapid after-slip and a-seismic deformation, but also potential deformation due to the M w = 6.3 aftershock of 04:30. The estimation of the co-seismic displacement of the second event (M w = 6.4) is obtained by difference between the position of the following day (246 -02/09/2020) and the position of the day of the 2 earthquakes previously described (midday 245, between 4h30 and 21:00 UTC). In a similar way, it also includes a combination of rapid after-slip and potential a-seismic deformation that might have occurred after both events (Fig. 6-B, vectors in light blue).

High rate GPS observations

High rate data are processed with Track software from MIT (T. Herring) which is a double-difference software, meaning that we compute the motion of a "rover" station relative to a "fixed" station. In this processing we use 5 "fixed" stations surrounding the area of interest (represented by black diamonds on Fig. 1): LSCH (La Serena) and LHOR (LosHornos) to the South; PAZU (Pan de Azucar) to the North; ALUM in Argentina and MRCG (Maricunga) to the East and North-East. We use the tropospheric zenithal delays (ZTD) generated by the 24 h static solution (one delay estimated every 2 hours at every site) to constrain the tropospheric delay in the kinematic processing to the static value. For the three largest events, we generate motograms (high rate evolution of position with time, from the Latin word "moto" for motion) of one hour spanning the events (see Fig. S5 for the mainshock at 4:09 UTC, Fig. S6 for the largest aftershock at 21h09 UTC and Fig. S7 for the smaller aftershock at 4:30 UTC). For all motograms, we built a sidereal filter by simply stacking the 1-hour data segments, of 3 (or 6) days before the earthquake with a 4 m 7 s time delay everyday following [START_REF] Choi | Modified sidereal filtering: implications for high-rate gps positioning[END_REF]. We then filter the co-seismic motogram, by simply subtracting this common mode to the original data. Then, the co-seismic jump is simply estimated as the offset between the 3-minutes data segment before and after the time of the earthquake (Fig. 6). Uncertainties are estimated visually from the motograms and range between 1 and 5 mm for the horizontal components and 5 and 10 mm for the vertical component. We are able to identify clear co-seismic jumps at most stations for the mainshock, small but discernable jumps at several stations for the largest aftershock, but nothing for the smaller aftershock of 4:30 UTC. This is an indication of the threshold detection of our current cGPS network: between magnitude 6.3 and 6.4.

Comparing coseismic offsets extracted from both the daily solution and from the HRGPS solution offers some confidence. Although the HRGPS is associated with larger uncertainties (5 mm) than daily solutions (1-2 mm), both solutions appear very con- sistent and show very similar offset. Specifically, stations located more than 50 km from the epicenter compare very well (BING, MMOR, UDAT, TAMR, TOT5, TRST). However, for both events, nearfield stations (TTRL, BAR2, and LLCH) exhibit a smaller HRGPS coseismic offset (smaller by 50%) than the daily solution one. This is very significant and indicates additional deformation is present immediately after the earthquake occurrence.

Analysis of major earthquakes

Coseismic slip static inversions

We built a fault geometry with triangular patches based on Slab2.0 (Hayes et al., 2018) between 26.5 • S and 29 • S and down to 60 km depth. We evaluated the slip distributions generated by the two largest earthquakes by inverting the coseismic displacements estimated from the HRGPS. We compute constrained least squares inversions using the CSI toolbox [START_REF] Gombert | Impulsive source of the 2017 mw= 7.3 Ezgeleh, Iran, earthquake[END_REF]. For both models, we apply as little smoothing as possible and we forbid back slip in the thrust direction. We assume only one slip component which direction is fixed parallel to the plate convergence (convergence vector from Klein et al., 2018a). Green functions are calculated at each node of the fault plane, assuming a homogeneous elastic half-space [START_REF] Meade | Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space[END_REF].

10 to 14 stations spanning the area were used in the inversion (Fig. 6). Resolution tests are fully described in the Supplemental material. They show that (1) a good recovery for ∼40 x40 km patterns is found between 15-55 km depth even with conservative noise budget for co-seismic offsets; (2) a very good (1-2 km) ability to locate the area of maximum slip; (3) peak-slip amounts are recovered within 10-30% and magnitude by 0.1; (4) extent of slip might be smeared by a few km.

For the mainshock, we find a slip distribution spreading over a rather large surface of 80 x 40 km 2 , between 27.5 • S and 28.5 • S. This surface seems too large for a M w 6.9 earthquake, but the bulk of the slip is concentrated in a much smaller area of only about 25 x 20 km 2 (Fig. 7-A). There is a trade-off between the quantity of slip and the size of the rupture zone. We test several models in which we concentrate larger slip amount in a narrower zone (for ex. within the region currently yielding more than 60 mm, or more than 80 mm of slip, see Fig. S8). Southward offsets can be reproduced by a larger amount of slip in the north (see Fig. S8-B).

But reducing the rupture zone to ∼30 x 30 km 2 leads to significantly larger residuals at closest stations (BAR2 and TTRL, Fig. S8-C). Therefore, the extension of the rupture zone to the north is required by the observations at more than 50 km, yielding significant westward coseismic offsets which are not converging toward a pin point. The best fit model includes a narrow strip of slip, elongated below the coastline south of the high slip area. This feature depicts only several cm of slip and is requested only by millimetric variations at a few stations. It may be beyond the resolution of our data and modeling. The deep extension of slip, reaching 40 km down, observed at 28 • S seems required both by the large coseismic displacement measured at station TOT5 located some 75 km away from the epicenter, and by the coseismic uplift measured at BAR2 and LLCH. Although vertical data do not appear essential since an inversion considering only the horizontal coseismic displacements produces similar slip pattern. We tested models with pure dip slip direction perpendicular to the trench, and models with two slip directions, but neither provides satisfying results (see supporting information for more details). We estimate a seismic potency of 4.14 • 10 8 m.m 2 , which corresponds to a moment of 2.01 • 10 19 N.m (M w = 6.8) using a shear modulus of 4.9 • 10 10 Pa (which is the value used for the W-phase). The geodetic moment appears slightly smaller than the seismic moment re-estimated at long-period using the W-phase but still lies within the error bar. Considering the size of this event, we made the approximation of a homogeneous half-space for all our inversions, which could account for part of the difference.

For the largest aftershock, because it generates smaller displacements than the mainshock at many stations, we dispose of less well determined co-seismic vectors. In particular we do not use any vertical displacement in the inversion of the aftershock slip distribution. Also, considering that we have very few observa- tions, we decreased the uncertainties of non-zero vectors to 1 mm, in order to strongly encourage the model to fit these. We find a circular slip distribution, significantly smaller with about 30 x 30 km overall (only 10 x 10 km for the bulk of the slip), with a peak slip at 95 mm (Fig. 7-B). For this event as well, the geodetic moment also appears slightly larger than the seismologic one, with 6.04 • 10 18 N.m (M w = 6.5, corresponding to a seismic potency of 1.24 • 10 8 m.m 2 ). Finally, the epicenter of the mainshock is located on the updip-western edge of the rupture zone, suggesting a downdip-bilateral propagation. The aftershock slip distribution is located updip the mainshock rupture zone and shows a striking complementary (Fig. 7-C). The aftershock lies in the hole left by the bean-shaped mainshock. Together they homogenize the slip over a larger and rounder area.

Dynamic inversions

We used seismic waveforms from strong motion stations deployed in the area (Fig. 1) to estimate the dynamic properties of the coseismic rupture. The low-frequency source properties (e.g., average slip and stress drop) of the largest event were estimated using an elliptical patch approach (e.g., [START_REF] Ruiz | Determination of the friction law parameters of the mw 6.7 Michilla earthquake in northern Chile by dynamic inversion[END_REF][START_REF] Herrera | Dynamic inversion of the 2015 Jujuy earthquake and similarity with other intraslab events[END_REF]). In this model, the rupture nucleates within a circular area and then propagates through a larger elliptical area. This rupture process is controlled by the friction law proposed by [START_REF] Ida | Cohesive force across the tip of a longitudinal-shear crack and Griffith's specific surface energy[END_REF]. Hence, in addition to the geometric parameters defining the circle and the ellipse, this dynamic model also includes the stress drop (T e ), the yield stress (T μ ), and the slip-weakening distance (D c ). We used strong-motion records integrated to displacement and band-pass filtered in low frequency (0.02-0.2 Hz for the mainshock). The AXITRA code [START_REF] Bouchon | A simple method to calculate Green's functions for elastic layered media[END_REF][START_REF] Coutant | Programme de simulation numerique AXITRA[END_REF] was used to simulate the source-to-receiver wave propagation via an appropriate 1-D velocity model for the area, which was extracted from [START_REF] Potin | Tomography of Chile[END_REF]. The inversion of the dynamic model was performed using the Neighborhood Algorithm (Sambridge, 1999), which finds the model that best fits the observed waveform data. The misfit between observed and modeled waveforms was calculated using an L2 norm.

The best solution for the mainshock converged toward an elliptical rupture of 24.4 km by 26 km (Fig. 8), with a minimum misfit of 0.24 (Figs. 8 andS11), a maximum slip of 1.1 m and a M w =6.7, which is similar to the solution obtained from HRGPS (Fig. 7). Also, the associated dynamic parameters are T e = 5.3 MPa, T μ = 5.59 MPa, and D c = 0.72 m. These dynamic parameters are similar to those obtained for inter plate events along Chilean subduction (Ruiz et al., 2017;[START_REF] Otarola | Dynamic rupture of subduction earthquakes located near the trench[END_REF] and the stress drop parameters are in the average of thrust earthquakes occurring on a subduction interface [START_REF] Kanamori | Theoretical basis of some empirical relations in seismology[END_REF]. Isodepth from Slab2.0 (Hayes et al., 2018).

Interplay between seismic and aseismic slip

Static inversion of afterslip on the day of the mainshock

On one hand, we compute the total co-seismic motion due to both events, both quantified by daily cGPS between the 30/08/2020 and the 02/09/2020 (shown in light red on Fig. 9-A). This calculation includes both events and the total amount of aseismic slip that occurred over the two days. On the other hand, we compute the total displacements measured by HRGPS (shown by dark red arrows on Fig. 9-A). Considering that the HRGPS allows to extract the pure co-seismic motion over a couple of minutes around the earthquakes, the difference between the total daily cGPS co-seismic estimates and the total HRGPS estimates (Fig. 9-B) should highlight the amount of early afterslip during the day of the earthquakes. Indeed, this difference shows a significant westward motion at TTRL and BAR2, similar to the post-seismic motion observed over the following days (Fig. 9-C).

Using the same methodology and parameters as previously (section 4.1), we compute static inversions of the different displacements fields. Unsurprisingly, the slip distribution inferred from the total daily cGPS displacements (Fig. 9-Ai, noted in following Ai) shows significantly more slip than the slip distribution inferred from the total HRGPS displacements (Fig. 9-Aii, noted in following Aii). In particular the peak slip of (Ai) reaches 17 cm compared to only 10 cm for (Aii). But both distributions show very similar patterns over a somehow circular area extending from 27.4 • S to 28.5 • S. The distribution of early postseismic shows slip occurring on a significantly smaller, narrow peanut-shape area elongated along a roughly NS direction (Fig. 9-B). Part of this slip could be coseismic slip due the M w = 6.3 aftershock which occurred at 4:30 UTC and that we were not able to extract from HRGPS. Small amount of slip observed at greater depth is most likely unresolved.

Time-dependent inversion of the postseismic deformation

Significant displacements are observed on the cGPS time series during a 22-days period, between the 2nd and the 24th of September. In order to quantify the slip evolution after the second large aftershock, we perform a kinematic inversion of the cGPS times series. We invert for slip on the subduction interface (following Rolandone et al. (2018) and [START_REF] Bletery | Slip bursts during coalescence of slow slip events in Cascadia[END_REF]). We find that the best fit to the time series is obtained with a smoothing parameter σ = 20 mm. day and a correlation distance between subfaults of D corr = 35 km. The total slip after 22 days is equivalent to M w = 6.8. Overall, it spreads over roughly the same area as the area ruptured by the mainshock and its largest aftershock, between 27.5 • S and 28 • S (Fig. 9-C). A static inversion of the cumulative postseismic displacement (following the same methodology as the coseismic static inversions) over the same period yields a very similar same pattern (Fig. S12). Regarding its spatio-temporal evolution, the post-seismic slip begins offshore and starts developing onshore and deeper after 6 days (Fig. 9-C1 & 9-C2). At a later stage, on the ninth day, a dissociated smaller patch begins more to the North, between 27 • S and 27.5 • S (Fig. 9-C3). It is deeperat a depth of approximately 35 to 55 km -and localized in the updip vicinity of the 2014 slow slip event (Klein et al., 2018b). The northward migration of post-seismic slip is associated with a northward rotation of post-seismic vectors wrt. co-seismic vectors at several stations near the epicenter area (BAR2, LLCH, TTRL) and the development of Westward vectors North of the epicenter area (BING, MMOR, UDAT). The source time function associated to this inversion shows a quasi-steady decrease in the slip's intensity. Then, negligible slip is found to occur after approximately 18 days. A movie of the postseismic slip evolution is provided in the supporting information.

Discussion and conclusions

General agreement and small discrepancies

Concerning the mainshock, the different types of modeling presented here are in good agreement, with some discrepancies regarding the magnitudes, the size of slip distributions and the peak slip. The dynamic model yields a smaller magnitude (M w = 6.7) than the one inferred from the W-phase (M w = 6.9). This is common and due to the simple elliptical geometry used for the dynamic inversion, which can therefore not fully capture the correct slip distribution and concentrate the solution. GPS constrained slip models yield a magnitude of 6.8, slightly smaller than the W-phase magnitude, but the difference is within the error bar (cf section 2.2 and table S1, same observation for the difference in magnitude of the largest aftershock between the static inversion and the Wphase analysis).

Slip models inferred from GPS show a quite larger rupture zone, which could have several origins. First, it could be an artefact imposed by wrongly detected small displacements at stations located farther away from the epicenter, although this should mostly be taken into account by the uncertainties. Second, the geodetic models might be contaminated by inaccuracies in the Earth model as we assume an homogeneous half-space and neglect topography (e.g., [START_REF] Duputel | Accounting for prediction uncertainty when inferring subsurface fault slip[END_REF][START_REF] Langer | Impact of topography on earthquake static slip estimates[END_REF]. Finally, the model resolution is limited by the number of observations, resulting in a trade-off between the amount of slip and the size of the rupture. Eventually, from both analyzes, we are confident that the greatest slip is well concentrated in an area of 30 x 30 km 2 , associated with the rupture of a single asperity. The HRGPS inversion shows a more extended rupture area, the lesser slip regions probably being at the resolution limit of our data.

Relation with coupling on the interface

We compare the two slip distributions with the coupling distribution proposed in the region by Klein et al. (2018a) (Fig. 10. The whole September sequence takes place in between the highly coupled Atacama segment (South of 28 • S), and the Chañaral segment (North of 27 • S). There, in the so-called Baranquilla inter-segment, we observe a narrow strip highly coupled connecting the 2 segments with significantly lower coupling on both the shallower and deeper part of the interface. We find that most of the slip due to the 01/09/2020 mainshock (dark blue contour) occurred downdip of its epicenter (dark blue dot), mostly overlapping the narrow strip of higher coupling. The largest aftershock at 21:09 UTC (light blue contour) shows a striking complementarity with the mainshock, occurring updip and extending in the low coupled region (Fig. 10). Early afterslip that occurs during the 17 hours between the mainshock and the largest aftershock (Fig. 9-B), is located mostly between the rupture zones of the two earthquakes, in a peanut-like shape (Fig. 10). Part of the obtained slip could be co-seismic due to the 4:30 UTC M w = 6.3 aftershock, and part indeed due to aseismic slip.

Interplay of seismic and aseismic slip in an area of heterogeneous coupling

We showed that the probability of having at least two aftershocks of magnitude M w > 6 within 24 hours is quite low. This leads us to question whether it is a simple mainshock-aftershock sequence or a seismic swarm, which is commonly defined as an increase of seismicity rate without a clear mainshock earthquake (Holtkamp et al., 2011). It could also mean that there is room for other processes that could have triggered these earthquakes so shortly after the mainshock.

During the first 17 hours, seismicity spread updip the mainshock epicenter and outside its rupture zone, into what later became the largest aftershock rupture zone, (Fig. 7-C). The asymmetry observed between the updip and downdip propagations of seismicity over the first 72 h (Fig. 3-C) is most likely driven by a specific source. Incidentally, the quantity of aseismic slip occurring directly after the mainshock, and over the following 20 days, is abnormally high. There is some overlap between the afterslip distribution, and the co-seismic slip distributions of the mainshock and the largest aftershock (Fig. 10). But the bulk of the distributions are disconnected and the overlap lies within the regions of lesser slip. We suggest that slow slip could be responsible for increasing shear stress at the front of the slip zone, propagating updip at approximately 0.7 km/hour ((which is within the range of slow slip propagation speeds observed elsewhere, [START_REF] Gao | Scaling relationships of source parameters for slow slip events[END_REF] -Fig. 3-C), until surrounding a locked asperity which eventually triggered the M w = 6.4 aftershock, 17 hours after the mainshock. Such a relation between seismicity at the front of the slip has been proposed in various context, during the interseismic phase but also during SSE, associated or not with non-volcanic tremors (Bartlow et al., 2011(Bartlow et al., , 2014;;Vaca et al., 2018;[START_REF] Bletery | Slip bursts during coalescence of slow slip events in Cascadia[END_REF], and is consistent with numerical models of seismicity driven by slow slip (e.g. [START_REF] Ariyoshi | Migration process of very low-frequency events based on a chain-reaction model and its application to the detection of preseismic slip for megathrust earthquakes[END_REF][START_REF] Yingdi | Preprint: Tremor Migration Patterns and the Collective Behavior of Deep Asperities Mediated by Creep[END_REF][START_REF] Wynants-Morel | Stress perturbation from aseismic slip drives the seismic front during fluid injection in a permeable fault[END_REF] The equivalent moment released over a period of 22 days, following the mainshock, reaches more than 80% of the coseismic moment, spreading in a much broader region than the coseismic rupture zone where the coupling is lower, as well as in a broader region than the aftershocks area. Usually, postseismic deformation reaches around 25% of the co-seismic deformation after a month. However, several cases have been documented where moderate size earthquakes are followed by abnormally large afterslip in Japan [START_REF] Yagi | Co-seismic slip, post-seismic slip, and aftershocks associated with two large earthquakes in 1996 in Hyuga-Nada, Japan[END_REF][START_REF] Suito | Interplate fault slip along the Japan trench before the occurrence of the 2011 off the Pacific coast of Tohoku earthquake as inferred from gps data[END_REF] and northern Peru where moment released through aseismic slip during a sequence was several time larger (3 to 14) than the moment released through earthquakes [START_REF] Villegas Lanza | A slow slip and seismic swarm sequence in a weakly coupled subduction zone in northern Peru[END_REF]. For the latter, it has been suggested that additional processes -i.e. not only an earthquake but also, for example, one or several slow slip eventswere involved to explain such a large amount of afterslip. A similar hypothesis was proposed to explain the abnormally rapid and large early afterslip following the 2016 M w 7.8 Ecuador earthquake (Rolandone et al., 2018). Complex sequence with large afterslip occurring very close a recurrent SSE patch was also observed in Mexico (Radiguet et al., 2016). The Atacama region seems propitious to slow slip events, while such an event was observed in the region in 2014 (Klein et al., 2018b). Here, overall, we estimated from the geodetic models that the sequence released a total moment of 4.94 • 10 19 N.m (M w = 7.1), with close to 60% through earthquakes and 40% through aseismic slip. Slip occurred spread over an area of ∼100 x 100 km 2 , much larger than expected for M<7 earthquakes, also highlighting the role of aseismic slip during the sequence. Postseismic slip migrates to greater depth 6 days after the mainshock, reaching eventually the 2014 slow slip area. Therefore, the Baranquilla LCZ seems prone for aseismic processes, potentially recurrent at depth as observed in the past, and favors large postseismic slip.

Considerations on seismic hazard in the area

Considering the historical seismicity in the region, i.e. the M w 8.5 mega-earthquakes of 1819 and 1922, and the high coupling imaged in the Atacama and Chañaral segments, we previously suggested that a joint rupture of these two segments was highly plausible in the future (Klein et al., 2018b). Both segments have indeed accumulated enough deformation since 1922 to gen-erate a M w 8 earthquake (Klein et al., 2018a). What is the impact of this sequence regarding scenarios for future megathrust ruptures in the region? Different scenarios seem plausible. On one hand, the whole September 2020 sequence is likely to have increased the stress at the edges of the highly locked Atacama and Chañaral segments, promoting future rupture(s) there. In particular the whole sequence occurred very near the northern edge of the Atacama segment. Could this initiate the destabilization of this highly locked patch and trigger a rupture of this segment already? And would a rupture of the Atacama segment trigger in turn the rupture of the Chanaral segment, initiating a 1819 or 1922 like megathrust earthquake? On the other hand, this same sequence may have released a significant amount of stress in the Baranquilla LCZ, which could in turn decrease the potential for a joint rupture of the Atacama and Chañaral segments by reinforcing its ability to act as a barrier for megathrust rupture propagation. In this scenario, Atacama and Chañaral segments could rupture independently, at different times and with smaller earthquakes than in 1819 and 1922. It is difficult to decipher between these scenarios, but the occurrence of a seismic sequence between two highly locked patches identified to be responsible for devastating earthquakes 100 and 200 years ago is a clear sign that this region should be monitored closely in the next future. Chapter 5. Semi-automatic detection tool for small slow slip events: application to Ecuador

CRediT

In this chapter, I focus on the systematic detection and characterization of small slow slip events (SSEs). For that, I provide the rationale for searching for small SSEs. Then, I describe a novel approach based on kinematic slip inversion of sliding time windows. Subsequently, I present preliminary results for the Ecuadorian subduction zone, showing the potential of the proposed approach. Finally, I propose several ideas to improve the detection.

Motivations

The development of an automatic SSE detection tool is grounded in the need to build a comprehensive catalogue of small SSEs. Since their discovery, numerous intermediate to large SSEs have been documented along subduction zones. The first slow slip event documented was evidenced by GNSS data along the Bungo Channel (Japan) and had a moment estimated to be M w 6.6 [START_REF] Hirose | A slow thrust slip event following the two 1996 Hyuganada earthquakes beneath the Bungo channel, southwest Japan[END_REF]. Along the Cascadia subduction zone, today known for experiencing frequent and large Episodic Tremor and Slip (ETS), the first slow event documented released a moment equivalent to a M w 6.7. Nowadays, every subduction zone monitored with a relatively dense GNSS networki.e. Alaska [START_REF] Fu | Repeated large slow slip events at the southcentral Alaska subduction zone[END_REF], Cascadia (Bartlow et al., 2011), Japan [START_REF] Obara | Episodic slow slip events accompanied by non-volcanic tremors in southwest Japan subduction zone[END_REF], Costa Rica (Dixon et al., 2014), Guerrero [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event[END_REF], Ecuador (Vaca et al., 2018), Peru (Villegas-Lanza, Nocquet, et al., 2016), Chile (Klein, Duputel, et al., 2018), Hikurangi (Bartlow et al., 2014) -experiences SSEs. Continental strike-slip faults, for instance the San Andreas fault [START_REF] Michel | Searching for transient slow slips along the San Andreas fault near Parkfield using independent component analysis[END_REF] and the North Anatolian fault [START_REF] Rousset | Slow slip events in the roots of the San Andreas fault[END_REF], are also experiencing SSEs. Even if the resolution of GNSS data is better than 20 years ago, the events documented today are still equivalent to M w ≥ 6. Michel et al. (2019b) showed that SSEs that occurred between 2007 and 2017 along the Cascadian subduction megathrust follow a frequency-magnitude distribution similar to that of earthquakes (see Figure 5.1a), the Gutenberg-Richter law. Moreover, Vaca et al. (2018) showed that along the Ecuadorian megathrust, near Esmeraldas (see location on Figure 5.3), the little and frequent modulations observed by the GNSS sites are correlated to the modulations observed in the seismicity rate (Figure 5.1b). These two observations suggest that smaller SSEs must be much more frequent than large ones.

Building comprehensive catalogues of SSEs would allow us to know if they occur only at the same place of already documented intermediate-sized SSEs or at other areas along the subduction interface. If so, we would identify new regions with a predominant aseismic behavior. We could also see which scaling law follow SSEs over a broader range of magnitude. Ide et al. (2007) suggest that SSEs follow a different scaling law than earthquakes, where their moment is proportional to their duration. In contrast, Michel et al. (2019b) propose that SSEs follow the same scaling law as earthquakes. A comprehensive catalogue of SSEs of different sizes would allow the re-examination of this question. In fine, it would also be an opportunity to investigate the origin of intermediate coupling regions. Several origins could explain intermediate coupling: the presence of geometrical complexities along the fault [START_REF] Taylor | Rapid forearc uplift and subsidence caused by impinging bathymetric features: Examples from the New Hebrides and Solomon arcs[END_REF][START_REF] Thirumalai | Variable Holocene deformation above a shallow subduction zone extremely close to the trench[END_REF], a non-stationary interseismic coupling [START_REF] Bruhat | Deformation rates in northern Cascadia consistent with slow updip propagation of deep interseismic creep[END_REF] or repeated short transient aseismic events of otherwise highly coupled areas.

The apparent scarcity of small SSEs can be partly explained by the fact that detecting and modelling them is challenging, since signals in data are close to the noise level. Another explanation can be found in the fact that the amount of GNSS data available today is very large, which is an exciting opportunity but also a real challenge because we simply do not know where to look. Two different approaches have been developed to detect small SSEs, either guessing their occurrence from the seismic signals that often accompany SSEs, or using GNSS time series. [START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes[END_REF]; Frank (2016) used the repeating low-frequency earthquakes and the tremors activity to target transient signals in the GNSS time series that are hidden because of the noise level. They applied this approach to the Guerrero and the Cascadia subduction zones and found many small slow slip events (i.e. corresponding to displacements < 1 mm in the GNSS time series). They suggested that the interseismic loading phase is punctuated by large and small slow slip events which release a part of the accumulated stress. This method allows the detection of submillimetric signals of slow slip events in the GNSS times series, but the disadvantage is that it requires a preliminary work on the seismological data. [START_REF] Rousset | A geodetic matched filter search for slow slip with application to the Mexico subduction zone[END_REF] developed a matched filter for GNSS time series. They built synthetic time series including surface displacements induced by SSEs that they correlated to the Guerrero subduction zone time series, known for experiencing SSEs. They detected numerous events that had not yet been documented, suggesting a highly dynamic subduction interface. This approach enables retrieving the location, duration and magnitude of each event, but it does not allow a good characterization of these parameters for events with M w < 5.6. Both approaches briefly described here are based on GNSS time series. In this thesis, I have developed a technique based on kinematic inversions of GNSS time series, because it requires a spatial coherence of displacements as well as a direction compatible with a slip on the subduction interface. This approach is described in the next section.

A semi-automatic detection tool based on kinematic inversions of GNSS time series

While the previously presented approaches do detect small events with seismological or GNSS data [START_REF] Frank | Uncovering the geodetic signature of silent slip through repeating earthquakes[END_REF][START_REF] Frank | Slow slip hidden in the noise: The intermittence of tectonic release[END_REF][START_REF] Rousset | A geodetic matched filter search for slow slip with application to the Mexico subduction zone[END_REF], the kinematic inversion approach applied to the 74-day afterslip of the 2015 M w 8. [START_REF]The M w 8.3 2015 Illapel afterslip spatio-temporal modelling[END_REF] Illapel earthquake in Chile, showed that it was able to identify a SSE that was invisible to the eye in the GNSS time series. The SSE signal was found within an already slipping patch of afterslip, and the SSE slip rate was estimated to be ≤ 5 mm.day -1 (see Chapter 3). We therefore propose here a method based on the direct modelling of GNSS time series, without preconceived notions of a particular region nor time period where SSEs occur. With this method, the SSE detection is robust because it requires that two conditions are met. First, a spatial coherence between small displacements not visible to the eye in the GNSS time series. Second, this set of coherent displacements must be consistent with slip along the subduction interface multiplied by the transfer functions. The interseismic velocity is not (Nocquet et al., 2017). The blue star shows its epicenter from Global CMT (Dziewonski et al., 1981). The dashed lines are the Slab2.0 iso-depth contours every 20 km (Hayes et al., 2018). removed from the time series. They are computed in the reference frame of the overriding plate, leading to a mean velocity in the opposite direction to that produced by reverse slip along the subduction interface. Thus, with the slip non-negativity condition implemented in PYEQ (see 2.3), this implies that when there is no trenchward displacement within GNSS sites in the overriding plate, there is no slip found in the inversion. This approach offers other advantages. The small computing resources required by PYEQ for performing time-dependent inversions enables to apply it to a wide spectrum of space and time and it directly provides a first estimation of occurrence and magnitude of small SSEs.

Figure 5.2 shows the general principle of the approach. For a given data set of GNSS time series, we perform successive kinematic inversions. The successive time windows are overlapping by the half of their duration (windows 1 and 2 on Figure 5). For instance, here we chose a 2-month long window to cover the longest possible period and to reduce the possibility of an event straddling two inversions, and we select a 1-month step for sliding the time windows. In other words, every month, a 2-month kinematic inversion is performed. If the region is too large, in order to keep the inversion efficient, we do not perform the inversion on the full subduction interface. Instead, we divide the subduction interface into several equivalent subdivisions and we perform an independent set of inversions for each subdivision (subdivisions A, B and C in Figure 5). Therefore, we obtain m x n slip distributions, where m are the number of subdivisions of the subduction interface and n the number of time windows. It is possible to define a sliding window for time and not for space, and vice versa. Since the interseismic velocity is not removed from the time series, the time window of the kinematic inversions must be relatively small so that the displacements induced by the plate convergence are small compared to potential SSEs induced displacements. It is to be judged according to the convergence rate, but in the case of Ecuador where the Nazca plate subducts beneath the North Andean Sliver at 46.4 mm.yr -1 (Nocquet et al., 2017;[START_REF] Jarrin | Current motion and deformation of the Nazca plate: new constraints from GPS measurements[END_REF], displacements over two months are reasonable since the interseismic displacement should be smaller than than 2 mm.

Slow slip events along the Ecuadorian subduction zone

The Ecuadorian subduction zone is well suited for applying the proposed methodology for several reasons. First, it experiences small and large SSEs. There are other indications that it would experience small ones very frequently. Second, the continuous GNSS network is dense since 2011 and provides data since 2009. I present the context of the Ecuador megathrust and then I discuss the choice made for the GNSS time series data set. Finally, I show some preliminary results. (Nocquet et al., 2017;[START_REF] Jarrin | Current motion and deformation of the Nazca plate: new constraints from GPS measurements[END_REF]. The Ecuadorian fault interface is heterogeneously coupled [START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF]Nocquet et al., 2017). The Colombian-Ecuadorian part of the subduction interface (i.e. at latitudes 1 • N-2.5 • N) is highly coupled from the trench down to 20 km depth. From 20 to 60 km depth, the interseismic coupling gradually diminishes. Further south, at latitudes 0.5 • N-0.5 • S, the fault is coupled in a southwest-northeast orientation, which is the direction of the Carnegie ridge subducted at the same area. At the latitude of La Plata island (1 • S), the fault is also highly coupled from the trench down to 25 km depth. Finally, at Salinas (latitude 2 • S), a small part of the megathrust is intermediately coupled and correlates with a subducted geometric complexity of the Nazca plate.

Context of the Ecuadorian subduction zone

Along the Ecuador megathrust, SSEs and seismic swarms occur concurrently (blue rectangles at shallow depth in Figure 5.3). In 2010, a one-week SSE with a moment release equivalent to a M w 6.0-6.3 took place below the La Plata island, where the fault is highly coupled [START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF]. A large seismic activity was reported at the same area and over the same time period, demonstrating a clear interaction with the SSE. In 2013, another one-week SSE accompanied by a seismic swarm was reported in the same area [START_REF] Segovia | Seismicity distribution near a subducting seamount in the Central Ecuadorian subduction zone, space-time relation to a slow-slip event[END_REF]. Synchronously to the SSE, seismicity occurred along intraslab faults, which are witnesses of subducted topographic features. In 2015, a 2month SSE took place at depth, between 50 and 70 km (Rolandone et al., 2018). At latitude 2.25 • S, below Salinas, recurrent small SSEs are suspected to occur but no study has been published to date. Nonetheless, Figure 5.4 shows the time series of SALN, located at Salinas (see location on Figure 5.5), from which we can observe a lot of complexity on the north and the east components. Regular millimetric or centrimetric southwestward displacements are observed since 2010. For instance, the east component shows several transient displacements including large ones (≥ 1 cm) and smaller ones (< 1 cm) almost every three months. This leads us to think that at this place the subduction interface is highly dynamic and regularly releases small fractions of stress. Further north, at latitudes 0.5 • N-0.5 • S, in 2016, a M w 7.8 earthquake ruptured a 110 km long segment of the megathrust near Pedernales (Nocquet et al., 2017). The large event ruptured two large areas that are highly coupled between 15 and 30 km depth, and stopped at the Punta Galera Mompiche region. In this region, regular SSEs concomitant to seismic swarms occur (Vaca et al., 2018(Vaca et al., ), in 2007(Vaca et al., -2008(Vaca et al., , 2011(Vaca et al., -2012(Vaca et al., and then in 2013(Vaca et al., -2014. These events released periodically a portion of the stress accumulated by the fault in this region, which might explain why the 2016 M w 7.8 Pedernales earthquake stopped and did not propagate into that part of the megathrust. The one-month afterslip triggered by this large earthquake took place in the regions where the SSEs took place prior the earthquake (Rolandone et al., 2018). The stress perturbation added by the mainshock also triggered a SSE at La Plata, located 100 km south of the mainshock rupture. The Ecuadorian subduction megathrust is therefore highly dynamic, has experienced several SSEs and is suspected of experiencing many more of smaller size SSEs. I now discuss the data set before moving on to detection.

Data set

Ecuador has many volcanoes. There are 17 active volcanoes in Ecuador and 16 potentially active (visit https://igepn.edu.ec/red-de-observatorios-vulcanologicos -rovig, last accessed: January 2023). Their activity shows signals that are different depending on the region and the time period, which adds a lot of complexity to the time series. First tests taking the whole Ecuador GNSS network showed that the tool detected events that were not SSEs. And when the events detected were SSEs, their model was not well constrained as a result of the signal recorded at sites located near the volcanoes. Therefore, it was decided to perform the inversions by removing all the sites nearby the volcanoes. In addition to its many volcanoes, Ecuador also contains part of the Amazon basin, east of the Andes. Therefore, GNSS time series record strong seasonal and hydrological signals. Filtering the time series with a model based on GRACE data (i.e., measurements of the gravity field anomalies of the Earth) would be a possible solution to remove these signals. But, for the preliminary results presented in this Chapter, we have simply excluded the GNSS sites located in the Amazon basin. The GNSS data set used here is shown in Figure 5.5. As for a classical slip inversion, we build the geometry of the fault where we want to look for small SSEs and we compute the transfer functions of the chosen data set (see section 2.3.5).

Preliminary results

Initially, our goal was to study the dynamic of SSEs at the scale of the whole Ecuadorian segment of the subduction. This study could not be achieved during the time of my thesis. Even if the tool was applied to the whole Ecuador subduction zone, we focus here on the detection results of shallow (depth < 20 km) SSEs in the Salinas area (see Figure 5.3 for location), where the GNSS time series presented in Figure 5.4 show many modulations that might be induced by small SSEs. In addition, this area is still poorly known compared to other areas of SSEs in Ecuador.

We show the result of the proposed methodology restricted to a 70 km long segment between latitudes 1.8 • N and 2.6 • S and from the trench down to 25 km. Figure 5.6 shows the moment released per day between January 2010 and November 2022 from successive inversions, with a time window of two months sliding every month. In order to assess how the found moment rate relates to small modulation in the observed displacements, the east component of time series of two coastal sites is superimposed to the moment rate (SALN as blue dots and SEEC as cyan dots). We have chosen these sites because they are the closest to the area where the SSEs are expected.

We first notice a large moment release -2.12 × 10 18 N.m/day -in March 2019. On March 31, 2019, a M w 6.2 earthquake occurred at latitude 1.99 • S and longitude 81.05 • W, 10 to 20 km from SALN site. The coseismic offset has been corrected in the time series so the large moment rate increase observed here is generated by the afterslip following this earthquake. The successive inversions capture the 1 t decay of the afterslip over a week following the mainshock.

Outside this period, the method highlights marked alternations between weeks to months long quiescent periods with null moment release, and well-defined periods lasting from a few days to three weeks, where the moment rate of the inversions is enhanced. The smallest signals for this period are of the order of 3 × 10 16 N.m/day, equivalent to M w 4.9. The methodology thus seems to be successful in identifying small SSEs. The moment peaks between 2010 and November 2022 correspond to small but clear modulations in the time series. For instance, the methodology highlights a moment release of 0.6 × 10 18 N.m/day (equivalent to a M w 5.8) in July and August 2012 (Zoom 1 in Figure 5.6). This moment peak corresponds to a transient event in SALN time series that induced a 5 mm westward displacement. In April and September 2015, three moment peaks are observed (Zoom 2 in Figure 5.6). In April, the moment release is about 1 × 10 18 N.m/day (equivalent to a M w 6) and corresponds to a 1-month westward displacement of 7 mm observed on SALN and SEEC time series. Similarly, in August and September, two distinct moment release peaks of 1.7 × 10 18 N.m/day and 0.5 × 10 18 N.m/day are observed (respectively equivalent to M w 6.1 and 5.7). These peaks correspond to a transient event on the time series. More precisely, between these two peaks, the time series show that the event seems to stop before starting again, demonstrating the fine sensitivity of the approach used. In August 2019, another moment release of 1.8 × 10 18 N.m/day is observed and correlates with a transient event both captured by SALN and SEEC sites through a 5 mm displacement (Zoom 3 in Figure 5.6). There are also subtle modulations in the time series that are not considered as SSEs by the approach used, probably because they are not consistent with modulations seen at other GNSS sites. Indeed, the methodology in that sense acts as a filter finding signal only when subtle transient displacements are consistent with slip multiplied by the transfer functions (see section 5.2). For instance, the SALN time series shows a transient event in August 2010 which does not correlate with an enhanced moment release, because it does not have a spatial consistency with other sites nearby.

Finally, between January 2010 and November 2022, we find an average of about two events per year. More specifically, we identify 7 large events (moment rate ≥ 0.25 × 10 18 N.m/day) and 24 smaller events (moment < 0.25 × 10 18 N.m). It seems that there are more small SSEs than large ones, and that the SSEs follow, at least in this preliminary study, a frequency-magnitude law as proposed by Michel et al. (2019b). It would be useful to better characterize the detected events to investigate this question. A smaller number of events is detected between 2010 and 2014, and at that time, there was a scarcer GNSS network.

Perspectives

The preliminary result shown in this Chapter was performed without temporal smoothing and with a weak spatial regularization. This allows a fine detection of SSEs, but does not allow a good characterization of their duration nor location on the subduction interface. To better characterize them, we must explore in the future the optimal regularization parameters. It would also be useful to compute the integral of the moment rates to obtain the total moment of each event. We could then discuss the significance of the SSEs.

Over almost 13 years, the SSEs detected with our methodology have released a moment equivalent to a M w 6.8. This is significant, especially since there is no record of a major historical earthquake in this southern Ecuador. Assuming that there has been no great earthquake since 1900 in this region, and that the subduction interface accumulates a slip deficit constant in time, it would have accumulate a M w 7.4 since 1900. This shows the importance of looking at small SSEs along the faults, at least along the southern Ecuadorian subduction.

The frequency of SSEs detected along the Ecuadorian subduction shows that it is active throughout the interseismic phase. The numerous modulations observed in the coastal time series of Figure 5.6 calls into question the way the interseismic velocity is estimated. Indeed, the estimation of a constant slope is not adequate since we observe an average of two SSEs per year. The work presented in this thesis highlights several aspects of the dynamics of slip along the South America subduction zone. As a conclusion, I summarize the most salient results of my work and discuss their importance in the global context of subduction zones. Finally, I conclude this manuscript by discussing the opportunities opened by the work done during this thesis for future research.

Conclusions

Afterslip within the coseismic rupture

The study of the Illapel earthquake afterslip showed that it mostly developed at the periphery of the rupture (Chapter 3). However, a non-negligible part, i.e. an area of about 40 × 20 km 2 , is located in the deep part of the coseismic rupture. According to the rate-and-state friction laws and its mechanical interpretation [START_REF] Dieterich | Time-dependent friction and the mechanics of stick-slip[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF][START_REF] Scholz | Earthquakes and friction laws[END_REF], seismic slip occur at velocity-weakening asperities whereas afterslip slip occur in velocity-strengthening regions. An earthquake breaks one or more velocityweakening asperities, which relock immediately after the rupture since they no longer have any seismogenic potential. The stress builds-up again on these velocity-weakening asperities with the plate convergence motion. Therefore, it is rather surprising that a part seismically broken during the Illapel earthquake has not relocked and has continued to slowly slip aseismically after having slipped rapidly. One can obviously question the resolution of the different models, and the existence of this overlap. However, most of the coseismic models find the deep propagation of the mainshock rupture (Melgar et al., 2016;Ruiz et al., 2016;[START_REF] Tilmann | The 2015 Illapel earthquake, central Chile: A type case for a characteristic earthquake[END_REF][START_REF] Fuentes | A study of the 2015 Mw 8.3 Illapel earthquake and tsunami: Numerical and analytical approaches[END_REF]Klein et al., 2017), just as postseismic models find the overlap of afterslip in the deepest part of the mainshock rupture (Barnhart et al., 2016;Shrivastava et al., 2016;Huang et al., 2017;Twardzik et al., 2021;Liu et al., 2022). Moreover, the resolution tests performed in Chapter 3 show that the afterslip is well resolved at this location. Along the Chilean subduction, it has already been observed that seismic and aseismic slip can take place at the same area. The 2010 M w 8.8 Maule earthquake ruptured 500 km of the megathrust, at latitudes 33.5 • S-38.5 • S. Three studies that used different data and methodologies found that afterslip occurred within the coseismic rupture at latitudes 36 • S-37 • S [START_REF] Vigny | The 2010 Mw 8.8 Maule megathrust earthquake of central Chile, monitored by[END_REF]Agurto et al., 2012;Bedford et al., 2013). In a global context of subduction zones, coseismic slip and afterslip were also proposed to have taken place in the same area, along the Japanese subduction megathrust. The 2011 M w 9.0 Tohoku-Oki earthquake afterslip took place in areas where moderate to large earthquakes have occurred prior the Tohoku-Oki event.

A long-term slow slip event was also observed prior to the 2011 M w 9.0 Tohoku-Oki earthquake and within its future coseismic rupture [START_REF] Mavrommatis | A decadal-scale deformation transient prior to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Mavrommatis | Long-term acceleration of aseismic slip preceding the Mw 9 Tohoku-Oki earthquake: Constraints from repeating earthquakes[END_REF]. What do these observations imply for the frictional anatomy of faults? A possibility is that velocity-weakening friction allows hybrid behavior with seismic and aseismic slip. In terms of seismic hazard, it implies that regions where occasional aseismic slip occurs does not mean that they do not have seismogenic potential.

The afterslip triggers large aftershocks

The Illapel afterslip model also highlighted a close relationship between the afterslip development and the largest aftershocks (Chapter 3). Two relationships between afterslip and large aftershocks have been identified. The first one is that the M w 6.8 aftershock ruptured a small 25 × 25 km 2 velocity-weakening asperity north of the deep part of the coseismic rupture as a result of the stress increment induced by the mainshock and the afterslip that developed in its vicinity. The second phenomenon is that the two M w 6.9 aftershocks, located 100 km north of the mainshock, were triggered prematurely by the afterslip developing north of the mainshock which probably increased the shear stress around very shallow velocity-weakening asperities. Both phenomena show that afterslip is able to put the clock forward for rupture of velocity-weakening asperities, already loaded by the stress increment induced by the mainshock rupture. During the 2020 Atacama sequence, a similar phenomenon was observed. While the first aftershock occurred in the immediate vicinity of the mainshock rupture, rapid afterslip then propagated updip toward the location of the second large aftershock. Together with the stress increment of the mainshock-first aftershock rupture, afterslip likely contributed to trigger the early and abnormally large aftershock later that day (Chapter 4). A similar process has also been observed for the postseismic phase of the 2016 M w 7.8 Pedernales earthquake along the Ecuadorian subduction. The one-month afterslip has developed in three different regions, including one north of the mainshock where two aftershocks of M w 6.7 and M w 6.9 ruptured occurred a month after the mainshock. These aftershocks were also likely triggered by the stress increment of the mainshock and of the afterslip (Rolandone et al., 2018). If confirmed for many earthquakes, this cascade of interactions between mainshock, afterslip and large aftershocks shows that modelling the afterslip in near-real time would make it possible to anticipate the regions where large aftershocks will occur.

A slow slip event within the afterslip

For the first time along a subduction megathrust, we were able to highlight a slow slip event during the afterslip following the Illapel earthquake (Chapter 3). The slow slip event occurred south of the mainshock rupture where large afterslip occurred for 47 days and where no significant aftershock took place. It was a small-sized event, with a maximum slip rate of 5 mm.day -1 but it correlates in space and time with an increase of the microseismicity rate. The event lasted 6 days and had a moment release equivalent to a M w 6.1-6.2. How could this slow slip event have taken place? The mainshock induced a stress increment in its periphery, which was released by afterslip. The slow slip event marks a re-acceleration of aseismic slip which was decreasing in rate. Poli et al. (2017) identified regular seismic swarms during the interseismic phase in this region, where there is also intraslab faults and where the Juan Fernandez ridge is subducting. Perhaps this geometric complexity induces a complex stress distribution along the fault, possibly favorable to slow slip event occurrence. The mainshock must have loaded the intraslab faults in addition to the subduction megathrust. The afterslip must have released some of the stress increment from the mainshock along the subduction fault, but it also added stress increment around the subducted objects and on the intraslab faults. Since slow slip events are assumed to occur on small velocity-weakening asperities (Rubin, 2008) and/or with fault geometrical complexities [START_REF] Romanet | Fast and slow slip events emerge due to fault geometrical complexity[END_REF], it must have occurred in the immediate vicinity of a subducted object. If this scenario is a hypothesis, the slow slip event found in the afterslip of the Illapel earthquake proves that the afterslip is not always steady but can experience short term slip modulation, even two months after the mainshock.

A specific spatial organization of seismic and aseismic slip

In Chapter 3, we identified three distinct afterslip patches and one slow slip event around the mainshock rupture, which spatially correlate with seismic swarms documented during the interseismic phase prior the Illapel earthquake. This particular spatial and temporal pattern has already been identified along several subduction zones worldwide. Along the Japanese subduction zone, a long-term slow slip event was spotted during the decade prior and updip to the 2011 M w 9.0 Tohoku-Oki megathrust earthquake [START_REF] Mavrommatis | A decadal-scale deformation transient prior to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Mavrommatis | Long-term acceleration of aseismic slip preceding the Mw 9 Tohoku-Oki earthquake: Constraints from repeating earthquakes[END_REF]. The afterslip developed updip of the mainshock rupture, in the same region of this decennial slow slip event (Johnson et al., 2012;[START_REF] Periollat | Transient brittle creep mechanism explains early postseismic phase of the 2011 Tohoku-Oki megathrust earthquake: Observations by high-rate GPS solutions[END_REF]. It was also observed for the subduction zone of Costa Rica. Regular slow slip events were identified between 2007 and 2012 at three preferential areas (Dixon et al., 2014). In 2012, the M w 7.6 Nicoya earthquake ruptured the megathrust in between the areas that experienced slow slip events during the interseismic phase. The 2.5 years afterslip following this large event also occurred both updip and downdip of the coseismic rupture (Hobbs et al., 2017), where the slow slip events were documented prior to the Nicoya earthquake. Then in 2014 and 2015, slow slip events occurred again, in the identical areas as for the afterslip and the interseismic slow slip events. Along the Ecuadorian subduction fault, the same pattern was identified. Several slow slip events and seismic swarms have been documented prior to the 2016 M w 7.8 Pedernales earthquake (Holtkamp et al., 2011;[START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF][START_REF] Segovia | Seismicity distribution near a subducting seamount in the Central Ecuadorian subduction zone, space-time relation to a slow-slip event[END_REF]Vaca et al., 2018) and acted as barriers for the propagation of its rupture. The afterslip developed both updip and downdip of the coseismic rupture, and the stress increment of the mainshock triggered a slow slip event 100 km of its rupture (Rolandone et al., 2018). These aseismic slip regions during the postseismic phase of the Pedernales earthquake are identical to the ones that hosted slow slip events and seismic swarms during the interseismic phase. A similar spatial and temporal organization of the slow slip events, the large earthquake rupture and the afterslip is thus found for four different regions. Regions that experience regular slow slip events during the interseismic phase are barriers to the propagation of seismic rupture of large earthquakes. These same aseismic regions will be the preferential zones for afterslip and slow slip events during the postseismic phase. Detecting and modelling slow slip events during the interseismic phase is a good way to identify active aseismic slip regions, which could potentially delineate large velocity-weakening asperities, that rupture through large earthquakes. 6.1.5 Interactions between aseismic and seismic slip: an indication to map the frictional properties along a fault

It is known that the rupture of a large earthquake controls the afterslip and aftershocks (Perfettini et al., 2010) and that aftershocks follow the evolution of the afterslip (Perfettini & Avouac, 2004). In this thesis, a global interaction between seismic slip and aseismic slip has been demonstrated for sequences of different sizes. Figure 6.1 shows the interaction model proposed after the Illapel afterslip study and its similarities with the Pedernales earthquake and the Nicoya earthquake as discussed in section 6.1.4. During the interseismic phase, regular slow slip events and seismic swarms occur at specific areas, releasing the stress on a regular basis, which implies that the fault is weakly to moderately coupled (Figure 6.1a). At some point, a large earthquake ruptures, whose rupture propagation is stopped by the areas that are aseismically active during the interseismic phase (i.e. areas experiencing slow slip events, Figure 6.1b). This coseismic rupture induces stress at its periphery, where afterslip initiates. It develops at preferential areas around the coseismic rupture, adding itself stress increment in its vicinity (Figure 6.1c). The mainshock and the afterslip stress increment triggers large aftershocks and slow slip events, where the afterslip developed (Figure 6.1d). A time-dependent model of the afterslip made it possible to highlight these interactions. Aside from the nature of these interactions, their location and their timing provide information on the size and distribution of the velocity-weakening asperities along the fault. Let us take the example of the Illapel earthquake. The M w 8. [START_REF]The M w 8.3 2015 Illapel afterslip spatio-temporal modelling[END_REF] Illapel earthquake, because of its size, broke a large velocity-weakening asperity. Around this large asperity, very small asperities were broken by the afterslip through microseismicity. The timing, the size of the large aftershocks (one M w 6.8 and two M w 6.9) and their location demonstrate that the fault is composed of sparse and rare medium-sized asperities located further away from the large one. The 2020 Atacama sequence showed different interactions, imaging a different map of the frictional properties along the fault in the region. The M w 6.9 mainshock ruptured a medium-sized asperity. At its immediate vicinity, smaller asperities ruptured through a M w 6.3 and a M w 6.4 events, explaining their timing and location. The Atacama sequence therefore broke a narrow network of medium-sized asperities. In both cases, we were able to demonstrate that seismic and aseismic slip interact with each other, and that the delay and distance between each interaction reflect the size and arrangement of the velocity-weakening asperities along the subduction megathrust.

Predominant aseismic sequences

Large earthquakes, slow slip events and seismic swarms, the latter two being often associated, are the different phenomena observed along faults. More original sequences with predominant seismic slip but still significant aseismic slip, such as the 2020 Atacama sequence (Chapter 4), seems to be typical sequences of the South American subduction zone. The Atacama sequence was a succession of three moderate events, followed by afterslip that released a moment equivalent to a M w 6.8. The entire sequence released a moment equivalent to a M w 7.1, of which 60% is seismic and 40% is aseismic. Aside from the 2020 Atacama sequence, a predominant aseismic sequence occurred in 2009 along the Peruvian subduction zone. It lasted 7 months and was a succession of mediumsized earthquakes (M w 5. 3-6.0) between which abnormally large aseismic slip occurred (Villegas-Lanza, Nocquet, et al., 2016). This kind of sequence was also observed along the Japanese subduction zone. The 1996 Hyuganada sequence started with two M w 6.7 earthquakes, which triggered a long and large afterslip [START_REF] Hirose | A slow thrust slip event following the two 1996 Hyuganada earthquakes beneath the Bungo channel, southwest Japan[END_REF]. The seismic slip released a moment equivalent to a M w 6.9, while the aseismic slip released a moment equivalent to a M w 6.6. The origin of such a large aseismic slip in this kind of sequence is still not understood, but they release a significant portion of the stresses along the fault. It seems therefore necessary to be able to detect and characterize any aseismic event along the faults, to better evaluate the slip deficit and thus the seismic hazard.

Perspectives

Aseismic processes along faults are closely related to seismic processes in space and time. Even if we manage to describe the major events (i.e. large earthquakes, large slow slip events) that take place along faults, these events are rare and do not allow us to study the diversity of processes and their interactions. Yet, there are several indications that along the South American subduction zone, smaller events occur more frequently and show new interactions between seismic and aseismic processes. Indeed, improved data resolution and time-dependent modelling approaches, will allow us to think towards smaller sequences to better understand the mechanism of larger ones. Smaller sequences should be thoroughly studied, as well as aseismic events. But for that, we have to be able to detect them.

A semi-automatic detection tool of small slow slip events has been developed based on time-dependent inversions of GNSS time series (Chapter 5). We applied it to the Ecuadorian subduction zone and the obtained results, though preliminary, are promising. This tool can be used to detect and target small events in order to model them. Since small slow slip events appear to be more frequent than large ones, detecting them would help to better assess seismic hazard in two ways. They will first help to better decipher the frictional properties variations along the subduction interfaces. They will also help to move toward the next generation of interseismic model and could include the fraction of aseismic slip accounted by small slow slip events. Rupture of an earthquake whose propagation is stopped at predominantly aseismic areas, where the stresses were regularly released. c

Afterslip occurrs in the immediate vicinity of the earthquake, with a possible overlap, and is triggered by the stress increment of the mainshock rupture. d Large aftershocks rupture, at the periphery of the afterslip, triggered by the stress increment induced by the mainshock and the afterslip. Repeating earthquakes are located at the periphery of the afterslip, triggered by its stress increment. Slow slip events can occur as a response of the stress increment induced by the mainshock, the afterslip and the aftershocks. Then, it is the beginning of a new interseismic phase with recurrent slow slip events, seismic swarms and repeating earthquakes.

Figure 1 . 1 :

 11 Figure 1.1: Strain at the Earth's surface throughout the seismic cycle. The interseismic strain grows linearly through time during decades to centuries, only disturbed by episodic slow slip events lasting days to a few years. The preseismic strain is millimetric to centimetric and can last seconds to months. The instantaneous coseismic strain induced by an earthquake is metric. The postseismic strain is millimetric to a few tens of centimeters, depending on the magnitude of the earthquake and can last weeks to years.

Figure 1 . 2 :

 12 Figure1.2: Summary of the processes at a subduction interface during the seismic cycle. The fault has velocity-weakening asperities embedded in an overall velocity-strengthening area. a During the interseismic phase, the stress builds up on the velocity-weakening asperities. Slow slip events occur at the periphery of velocity-weakening asperities together with seismic swarms, or downdip of the seismogenic zone with tremors. b During the preseismic phase of an earthquake, the microseismicty increases and foreshocks occur. Slow slip events were also proposed during this phase. c During the coseismic phase, an earthquake release the stresses accumulated by the fault through seismic slip at velocity-weakening asperities. d During the postseismic phase, afterslip develops in the periphery of the earthquake rupture, and microseismicity mimics its evolution. Large aftershocks also occur.

Figure 1 . 3 :

 13 Figure 1.3: Compilation of the large earthquakes approximate ruptures, seismic swarms and slow slip events documented since 1995. A Vaca et al. (2018), B Nocquet et al. (2017), C Rolandone et al. (2018), D Vallée et al. (2013); Segovia et al. (2018), E undergoing work, F Villegas-Lanza, Nocquet, et al. (2016), G Sladen et al. (2010), H Ruegg et al. (2001); Perfettini et al. (2005), I Ruegg et al. (2001); Melbourne and Webb (2002), J Duputel et al. (2015) K Ruiz et al. (2014), L Socquet et al. (2017), M Peyrat et al. (2010), N Pritchard et al. (2002), O Comte et al. (2002, 2006); Holtkamp et al. (2011), P Klein, Duputel, et al. (2018), Q Melgar et al. (2016); Grandin et al. (2016); Ruiz et al. (2016); Klein et al. (2017), R Poli et al. (2017), S Ruiz et al. (2017), T Delouis et al. (2010); Moreno et al. (2010).

Figure 2 . 1 :

 21 Figure 2.1: Principle of real-time positioning by trilateration. Three satellites are required to locate the receiver and a fourth satellite is used to correct the clock error of the receiver.

Figure 2 . 2 :

 22 Figure 2.2: Principle of double-difference. A first simple-difference enables to estimate and remove the satellite clock error using two receptors. A second simple-difference allows to estimate and remove the receiver clock errors using the same two receptors but a different satellite.

Figure 2 . 3 :

 23 Figure 2.3: Selection of GNSS time series that recorded the 2015 deep slow slip event in Ecuador. a) Detrended time series. b) Time series filtered with a Vondrák (cut-off frequency of 10 cycles per year). The red shaded area is the presumed occurrence of the deep slow slip event.

Figure 2 . 4 :

 24 Figure 2.4: a) GNSS network used to study the 2015 deep slow slip event. GNSS sites are shown as black triangles. GNSS sites used to estimate the common mode are shown as red triangles. The dashed lines are the Slab2.0 iso-depth contours every 40 km (Hayes et al., 2018). b) Time series used to estimate the common mode. The red shaded area is the presumed occurrence of the deep slow slip event.

Figure 2

 2 Figure 2.5: a) Common mode estimated from the time series shown in Figure 2.4b. b) Selection of GNSS time series filtered with the common mode. It is the same selection as in Figure 2.3. The red shaded area is the presumed occurrence of the deep slow slip event.

Figure 2 . 6 :

 26 Figure 2.6: a) Interseismic coupling model and tremors represented as grey dots. The blue contours are the regions where slow slip events occurred. b) Slip rate of the events as a function of time. The black dots are the tremors. a) and b) are at the same scale. These Figures are from Michel et al. (2019a).

Figure 2 .

 2 Figure 2.7: Slip rate evolution of the 2009 slow slip event along the Cascadia megathrust averaged over two-day intervals. The black dots show the tremors. This Figure is from Bartlow et al. (2011).

Figure 2 .

 2 Figure 2.8: Slip rate evolution of the Manawatu slow slip event along the Hikurangi megathrust, averaged over sixteen day intervals. This Figure is from Bartlow et al. (2014).

Figure 2 .

 2 Figure 2.9: Slip evolution through 50-day snapshots of the 2006 slow slip event along the Guerrero megathrust. This Figure is from Radiguet et al. (2011).

Figure 2 .

 2 Figure 2.10: Illustration of the general principle for GNSS data modelling for the case of a subduction megathrust.

  a l r o u g h n e s s S p a t ia l ro u g h n e s s

Figure 2 .

 2 Figure 2.11: Surface plot of the L-curve. χ 2 as a function of the model defined as the spatial roughness and the temporal roughness. The red star corresponds to the best model.

Figure 3 . 1 :

 31 Figure 3.1: Seismo-tectonic context of the Illapel region. The orange lines indicate the proposed lateral extent of historical and recent megathrust earthquakes (Ruiz & Madariaga, 2018). The green stars represent the large aftershocks epicenters of the Illapel earthquake (CSN catalogue). The red curves indicate the Illapel 2015 M w 8.3 coseismic slip distribution contoured every 2 m(Klein et al., 2017). The dashed lines show the Slab2.0 iso-depth contours every 20 km(Hayes et al., 2018).

78

  Chapter 3. The M w 8.3 2015 Illapel afterslip spatio-temporal modelling(Tissandier et al., 2023)

3. 4 .

 4 Afterslip of the M w 8.3 2015 Illapel earthquake imaged through a time-dependent inversion of continuous and survey GNSS data 79

2

  Previous works on the 2015 Illapel earthquake With a good GNSS data coverage and a coastline located ∼70 km away from the trench, the 2015/09/16 M w 8.3 Illapel earthquake offers a good opportunity for a detailed study of afterslip time evolution following a large megathrust earthquake. It ruptured a 200 km long segment of the megathrust along the Nazca/South America subduction in Chile. The Illapel segment experiences regular M8+ earthquakes every 60-80 years

Figure 1 .

 1 Figure 1. GNSS network used in this study and past megathrust earthquakes. Black and white triangles are permanent and survey GNSS sites respectively. The green stars are the large aftershocks epicenters (Mw ≥ 6.8) of the Illapel earthquake (CSN catalogue). The red curves indicate the Illapel 2015 Mw 8.3 coseismic slip distribution contoured every 2 m (Klein et al., 2017). The red star is the Illapel earthquake epicenter. The orange lines indicate the proposed lateral extent of historical and recent megathrust earthquakes(Ruiz & Madariaga, 2018). The

Figure 2 .Figure 3 .

 23 Figure 2. Observed and modeled time series at selected sites located on Figure 1. a) cGNSS sites PFRJ, CNBA, TOLO and b) sGNSS site POBR. Left column is the North component, middle column the East component and right column the vertical component. The gray dots indicate the daily displacement since the first date of observation together with the 1-σ error bar. The red curve shows the modeled time series according to our preferred afterslip model.

  ) where d k is the vector of the observation for subset k with length n d k and ŝk is the estimated model estimated without including observations d k and G k the sub-matrix of G whose lines correspond to the observations in d k . Too smooth models provide bad fit to all data and a large value of M SE k . Too rough models, although providing good fit to the data used in the inversion, will provide bad prediction at their spatial and temporal neighbors that were not included in the inversion. In order to reduce the dependence of M SE on the chosen tested data subset, the inversions are repeated for all k subsets, and the final M SE score is simply the mean of all M SE k . Most studies use 5 to 10 subsets, ensuring that results do not depend on the chosen division of the data. Here we use K=9, splitting the analyzed period into 3 consecutive time windows and dividing the GNSS site into 3 separate subsets. Thus for each K, one third of the sites have one third of their observations removed, which are used to compute M SE k . Results showing the MSE score as a function of λ space and λ time are shown in Figure 4a.

(Figure 4

 4 Figure 4. a) Cross-validation and b) L-curve. The black star is the smooth afterslip model (λspace = 0.8 and λtime = 5) and the black square is the rough afterslip model (λspace = 1 and λtime = 0.1).

  and movie S1). A first patch is located north of the coseismic rupture at relatively shallow depth (∼0-20 km, Box A, Figure 5a-g). The two other patches of slip (Boxes B & C) develop at greater depth (∼40 km) on the interface, on each side north and south of the deepest part of the coseismic rupture (Figure 5a-g).

Figure 5

 5 Figure 5. a) to g) Cumulative afterslip evolution (smooth model) contoured in red from 20 cm and then every 10 cm. Slip inversion results for the largest aftershocks on November 7 and November 11 are contoured in brown. Their focal mechanisms is from Global CMT (Dziewonski et al., 1981; Ekström et al., 2012). The mainshock slip is contoured in blue every 2 m (Klein et al., 2017). The blue star is the Mw 8.3 2015 Illapel epicenter. The dashed lines are the Slab2.0 iso-depth contours every 20 km (Hayes et al., 2018). h) Cumulative source time function in green and cumulative number of events (CSN catalogue) as blue bars. The vertical dashed lines in red show the date of the largest aftershocks.

Frank

  et al. (2017) used the waveform of earthquakes from the CSN catalogue(Baez et al., 2018) and applied a matched-filter search for additional earthquakes. They obtained a catalogue including more than 16 000 events spanning 19 months, from 2015/01/01 to 2016/06/27, covering the period analyzed in our study.At a first order, our geodetic afterslip model, the repeater afterslip model and the density map of aftershocks all show very similar spatial patterns, outlining a crescentshaped area surrounding the coseismic rupture with a strip of slip and aftershocks crossing the deep extent of the coseismic rupture (Figure6a-c). The observed spatial correlations between the geodetic afterslip models, the repeater afterslip model and the aftershocks spatial distribution are much better than any previous published models of af-

Figure 6 .

 6 Figure 6. Comparison of a) the 30 days cumulative afterslip smooth model contoured in red every 20 cm with b) Huang et al. (2017) repeater afterslip model and c) the first 30 days of Frank et al. (2017) earthquakes catalogue (following the mainshock). The Boxes are identical to Huang et al. (2017) and differ from Figure 5g. The green dots are the repeater earthquakes locations of Huang et al. (2017). The mainshock slip is contoured in blue every 2 m (Klein et al., 2017). The dashed lines are the Slab2.0 depths every 20 km (Hayes et al., 2018). Comparison of the slip rate evolution of our afterslip models with Huang et al. (2017) repeater afterslip model within d) Box 1, e) Box 2 and f ) Box 3.

  ,b), like the southernmost extension of afterslip in Box 3 or the overall southwest-northeast orientation of the slip area of Box 3.

  Figure 7. a) to f ) Daily slip rate (smooth model) contoured in green every 5 mm/day (the first 3 mm/day of slip are hidden). The black dots are the seismicity occurring at the selected dates. Slip inversion results for the largest aftershocks on November 7 and November 11 are contoured in brown. The focal mechanisms are from GlobalCMT (Dziewonski et al., 1981; Ekström et al., 2012). The mainshock slip is contoured in blue every 2 m (Klein et al., 2017). The blue star is the Mw 8.3 2015 Illapel epicenter. The dashed lines are the Slab2.0 iso-depths every 20 km (Hayes et al., 2018). Source time function of g) the smooth model and h) the rough model in green and events per day (CSN catalogue) as blue bars. The vertical dashed lines in red are the aftershocks occurrence.

Figure 8 3 2015

 83 Figure 8. a) Daily slip rate of an inversion ending on November 6. b) and c) Daily slip rate of inversion starting on November 7 and ending on November 10. The slip rate is contoured in green every 2 mm/day. The black dots are the seismicity occurring at the selected dates (CSN catalogue). Slip inversion results for the Mw 6.8 on November 7 is contoured in brown. The mainshock slip is contoured in blue every 2 m (Klein et al., 2017). The blue star is the Mw 8.3 2015 Illapel epicenter. The focal mechanisms are from GlobalCMT(Dziewonski et al., 1981; Ekström et al., 2012). The dashed lines are the Slab2.0 iso-depths every 20 km(Hayes et al., 2018).

Figure 9

 9 Figure 9. a) Daily slip rate of an inversion starting on November 7 and ending on November 10. b) to d) Daily slip rate of an inversion starting on November 11. The slip rate is contoured in green every 2 mm/day. The black dots are the seismicity occurring at the selected dates (CSN catalogue). Slip inversion results for the Mw 6.8 on November 7 and November 11 are contoured in brown. The focal mechanisms are from GlobalCMT(Dziewonski et al., 1981; Ekström et al.,

Figure 10

 10 Figure 10. a) to e) Daily slip rate (rough model) of the afterslip modulation south of the mainshock at selected dates. The slip rate is contoured in green every 5 mm/day (the first 3 mm/day of slip are hidden). The black dots are the seismicity occurring at the selected dates (CSN catalogue). Slip inversion results for the largest aftershocks on November 7 and November 11 are contoured in brown. The focal mechanisms are from GlobalCMT(Dziewonski et al., 1981; Ekström et al., 2012). The mainshock slip is contoured in blue every 2 m(Klein et al., 2017).

7. 3 Figure 11 .

 311 Figure 11. Comparison of the 74 days cumulative afterslip smooth model contoured in gray from 20 cm and then every 10 cm with the interseismic coupling of Métois et al. (2014). The white rectangle is the overlapping area with the coseismic rupture. The mainshock slip is contoured in blue every 2 m (Klein et al., 2017). Slip inversion results for the largest aftershocks on November 7 and November 11 are contoured in brown. The dashed lines are the Slab2.0 depths every 20 km(Hayes et al., 2018).

Figure 12 .

 12 Figure 12. Comparison of the 74 days cumulative afterslip smooth model with the interseismic seismicity(Poli et al., 2017). The green dots are the repeater earthquakes locations ofPoli et

  Poli et al. (2017) noticed that the spatial clustering of the seismic activity before the Illapel earthquake shares similar orientation with the main fracture zones observed on the outer rise of the subducting Nazca plate. Therefore, they propose that areas of the megathrust experiencing seismic swarms correspond to fractures within the subducted Nazca slab, characterized by high pore-fluid pressure and frequent fluid releases. High pore pressure reduces the normal stress, favours stress release through transient aseismic slip during the interseismic phase and through afterslip after a nearby large rupture.

  50 km × 50 km, 75 km × 75 km and 100 km × 100 km. The slip of the checkerboard is 898 mm corresponding to M 0 = 3.16e+20 N.m, which is the moment determined after 74 days. We noised the data with a normal distribution, with 2 mm on the horizontal components and 5 mm for the vertical component. For each checkerboard, we perform several inversions using different smoothing parameters and the χ 2 as a function of model roughness. The chosen model corresponds to the best trade-off between model roughness and misfit indicated by the L-curve.

TISSANDIER

  Figure S1. GNSS network used in this study. Black and white triangles are permanent

Figure S2 .

 S2 Figure S2. Fault geometry used in the kinematic inversion. The fault is discretized into

Figure

  Figure S3. χ 2 as a function of the roughness of the model defined as 1/λ space (1/spatial

Figure S4 .Figure S5 .Figure S6 .

 S4S5S6 Figure S4. Spatial resolution tests. a) is 100 km × 100 km patches, b) is 75 km × 75 km and c) is 50 km × 50 km. The left plot is the checkerboard and the right one is the inverted result. Blue and green arrows are respectively the synthetic and the modeled

Figure 4 . 1 :

 41 Figure 4.1: a) Atacama seismo-tectonic context. The focal mechanisms are the 2020 Atacama mainshock and aftershocks from Global CMT(Dziewonski et al., 1981). The black triangles represent the GNSS sites used in this study. The red ellipse show where the1973, 1976, 2006 and 2015 seismic swarms occurred(Comte et al., 2002[START_REF] Comte | Copiapo Chile seismic swarm analysis: mapping the interplate contact[END_REF] Holtkamp et al., 2011). The blue curves indicate the Copiapó 2014 slow slip every 10 cm(Klein, Duputel, et al., 2018). The orange line indicates the proposed lateral extent of the historical 1819 and 1922 megathrust earthquakes(Lomnitz, 2004; Ruiz & Madariaga, 2018). b) Interseismic coupling model fromMetois et al. (2016). The grey lines delineate three segments. The dashed lines are the Slab2.0 iso-depth contours every 20 km(Hayes et al., 2018).

Figure 4 . 2 :

 42 Figure 4.2: a) Dynamic slip model of the mainshock. b) Coseismic slip model of the mainshock. c) Coseismic slip model of the second large aftershock. The slip contours start at 40 mm and are represented every 20 mm for b)and c). The grey lines are the Slab2.0 iso-depth contours every 20 km(Hayes et al., 2018). d) First aftershock and afterslip slip model between the mainshock and the second aftershock. The slip contours are represented every 25 mm. The Figures were modified afterKlein et al. (2021).

Figure 4 . 3 :

 43 Figure 4.3: BAR2 and TTRL continuous GNSS time series from the day after the mainshock (2020/09/02) to the 2020/09/29 with respect to the South American plate. The blue dots indicate the daily displacements and the grey vertical bars are the 1-σ error bar. Their location can be found on Figures 4.1aand 4.4.

4. 4 .

 4 The time-dependent inversion of the postseismic deformation of the 2020 Atacama sequence 135

Figure 4 . 4 :

 44 Figure 4.4: Observed cumulative horizontal displacements field on 2020/09/24, 22 days after the mainshock, with respect to the South American plate. The error ellipses are 1-σ confidence. The black triangles are the location of the GNSS sites time series in Figure 4.3. The dashed lines are the Slab2.0 iso-depth contours every 20 km(Hayes et al., 2018).

  and 4.4) shows 12 mm of westward displacements after one month (Figure4.3), while the north component does not display any clear displacement. The vertical component of the two sites is noisy and does not show a clear signal. Thus, I did not use it.

Figure 4 .

 4 4 shows the observed horizontal displacements field, with an overall westward motion. The greatest displacements are recorded at sites BAR2 and TTRL. Almost 5 mm of southwestward displacement is also visible on the sites located at latitudes 27.25 • W-28.25 • W and longitudes 70.25 • W-71 • W.

  4.5.1 The 2020 Atacama sequence in the context of stress release process along the South America subduction zone

Figure 4 .

 4 Figure 4.5: Cumulative afterslip evolution contoured in red from 10 mm and then every 5 mm. The black and grey dots are respectively the seismicity occurring at the selected dates and the cumulative past seismicity (CSN catalogue). The focal mechanisms correspond to the 2020 Atacama mainshock and aftershocks from Global CMT(Dziewonski et al., 1981). The blue curves indicate the Copiapó 2014 slow slip every 10 cm(Klein, Duputel, et al., 2018). The dashed lines are the Slab2.0 iso-depth contours every 20 km(Hayes et al., 2018).

Figure 4 . 6 :

 46 Figure 4.6: The grey and red arrows are respectively the observed and modelled cumulative horizontal displacements. The error ellipses are 1-σ confidence. The cumulative afterslip model on the 2020/09/24 is contoured in red in red from 10 mm and then every 5 mm.The dashed lines are the Slab2.0 iso-depth contours every 20 km(Hayes et al., 2018).

Figure 4 .

 4 Figure 4.7: Daily slip rate evolution contoured in green every mm. The black dots show the seismicity occurring at the selected dates (CSN catalogue). The focal mechanisms correspond to the 2020 Atacama mainshock and aftershocks from Global CMT (Dziewonski et al., 1981). The blue curves indicate the Copiapó 2014 slow slip every 10 cm(Klein, Duputel, et al., 2018). The dashed lines are the Slab2.0 iso-depth contours every 20 km(Hayes et al., 2018).

Figure 4 . 8 :

 48 Figure 4.8: Source time function of the afterlip rate model in green and events per day (CSN catalogue) as blue bars.

Fig. 1 .

 1 Fig. 1. Overall context of the sequence of September 2020 in the Atacama region of Chile. The relocated earthquakes catalog is plotted as a function of time since the mainshock (in days since the mainshock). Events represented with white contours were relocated outside of the core sequence. Mechanisms and M w of the 3 largest events are the re-estimated one. The different observation networks used in this study are represented. A. Cross section of the relocated catalog of the core sequence, as function of depth, with the same color scale function of time. B. Local magnitudes M l of the relocated catalog of the core sequence as function of time. Violet stars show swarms locations(Holtkamp et al., 2011). Slab isodepth fromHayes et al. (2018). The dashed red lines illustrate the approximate length 1819 and 1922 earthquakes rupture zones. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 2 .

 2 Fig. 2. Vertical cross-section of tomography model for both P-wave velocity and V P /V S velocity ratio. Cross-section oriented perpendicular to the trench, across the sequence. Blue dots represent background seismicity in a 50 km range from the cross-section (CSN catalog, 2013 -07/2020). Red dots represent the seismic sequence between August 25th and September 25th included. All earthquakes were relocated in the local 3D tomography model presented. Seismicity of the sequence spreads along subduction contact between the trench and about 40 km at depth.

Fig. 3 .

 3 Fig. 3. Analysis of the relocated catalog over the first 72 h after the mainshock; A. map view, distances are in km, coordinate (0,0) correspond to the mainshock and colors represent time; B. number of events/h, time is relative to the mainshock origin time. Bins are centered on the hour; C. distance to mainshock in the West-East direction vs. time, in km; D. distance to mainshock in the South-North direction vs. time, in km; E. horizontal distance to main shock in km. The 3 main events are highlighted by the red diamonds on subplots C, D & E, the mainshock is black contoured. The red dashed lines on subplots C and D depict the seismicity boundaries.

Fig. 4 .

 4 Fig. 4. Histograms of the bootstrap analysis for the moment magnitude of the three largest events with 10 5 inversions. Average M w and ±2σ uncertainties are given in the legend.

Fig. 5 .

 5 Fig. 5. Time series of GPS daily positions from stations in the region of the sequence on the 3 components. The vertical black lines flag the exact time of the 2 events of September 1st, at the beginning of the sequence.

Fig. 6 .

 6 Fig. 6. Comparison of 24 hours CGPS and HRGPS static co-seismic offset estimations; A. mainshock of 04:09 (reddish vectors) and B. aftershock of 21:09 (blueish vectors). Horizontal top row, vertical bottom row. Earthquakes' locations from the relocated catalog and mechanisms from the W-phase analysis.

Fig. 7 .

 7 Fig. 7. Slip distributions of A. the mainshock at 4:09 UTC inverted from the HRGPS (Fig. 6-A); B. the aftershock at 21:09 UTC inverted from the HRGPS (Fig. 6-B); Distributions are represented as the blue color scale (in mm), blue isolines are represented every 20 mm; Horizontal coseismic displacements are depicted by arrows: Observations (red) vs predictions (pink); Vertical coseismic displacements are depicted by colored dots: Observations (big circles) vs predictions (small circles) with amplitude represented with the polar color scale; C. Zoom in to compare both slip distributions and the relocated catalog of aftershocks occurring between the 2 events represented with the color scale. Isodepth from Slab2.0 (Hayes et al., 2018).

Fig. 8 .

 8 Fig. 8. Coseismic model of the mainshock obtained from the dynamic inversion. A. Geographic context of the mainshock rupture and stations used for modeling. The moment tensor was obtained from GCMT. B. Dynamic slip model on the fault plane and waveform misfit convergence colored with the stress drop. The bottom plot shows the E-W observed (blue) and modeled (red) waveforms of the best dynamic model.

Fig. 9 .

 9 Fig. 9. Slip history over the sequence: A. Total coseismic: vectors show the total coseismic displacement on September, 1st (including both events and the aseismic slip that occurred during that period) measured by CGPS (light red) and the corresponding Slip distribution (i), compared with the total coseismic displacement due to the 2 events measured by HRGPS (dark red) and the corresponding slip distribution (ii); B. Early afterslip estimated from the difference between CGPS and HRGPS estimates and the corresponding slip distribution; C. Slip-time dependent inversion of the postseismic deformation 22 days with 3 snapshots of the cumulative slip distribution. Yellow and red arrows are respectively model-predicted and observed displacements for CGPS sites recorded since the mainshock. Postseismic slip contours are every 10 mm. Gray lines are Slab2.0 isodepth from Hayes et al. (2018).

Fig. 10 .

 10 Fig. 10. Slip distributions of the M w = 6.9 mainshock (01/09/2020-4:09 UTC, dark blue contours every 20 mm starting at 60 mm), the M w = 6.4 aftershock (01/09/2020-21:09 UTC, light blue contours every 20 mm starting at 60 mm); the rapid afterslip between the 2 events (red contours every 20 mm starting at 40 mm), and 1 month of postseismic slip (yellow contours every 15 mm starting at 15 mm). The epicenter of the M w = 6.2 aftershock (4:30 UTC) is depicted by the orange dot. Comparison with the coupling distribution in the region (Klein et al., 2018a) and the 2014 SSE distribution (Klein et al., 2018b, represented by the dark green contours every 50 mm starting at 200 mm)). Background seismicity from relocated catalog depicted by white dots. Slab2.0 isodepth from Hayes et al. (2018) every 10 km.

  authorship contribution statement E. Klein: Conceptualization, Data curation, Formal analysis, Investigation, Validation, Writing -original draft, Writing -review & editing. B. Potin: Data curation, Formal analysis, Writing -review & editing. F. Pasten-Araya: Data curation, Formal analysis, Writing -review & editing. R. Tissandier: Formal analysis, Investigation, Writing -review & editing. K. Azua: Data curation, Formal analysis, Writing -review & editing. Z. Duputel: Formal analysis, Funding acquisition, Investigation, Writing -review & editing. C. Herrera: Formal analysis, Investigation, Writing -review & editing. L. Rivera: Formal analysis, Investigation, Writing -review & editing. J.M. Nocquet: Formal analysis, Funding acquisition, Writingreview & editing. J.C. Baez: Resources, Writing -review & editing. D. Zigone: Writing -review & editing. R. Madariaga: Validation, Writing -review & editing. J.P. Ampuero: Resources, Writing -review & editing. S. Ruiz: Data curation, Formal analysis, Funding acquisition, Investigation, Validation, Writing -review & editing. C. Vigny: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Validation, Writing -original draft, Writing -review & editing.

  Figure 5.1: a) Frequency-magnitude distribution of SSEs along the Cascadia subduction zone from Michel et al. (2019b). b) Number of events per month (red bars) and GNSS time series ESMR (blue dots) for the Punta Galera in Ecuador (see location on Figure 5.3) from Vaca et al. (2018).

Figure 5 . 2 :

 52 Figure5.2: Illustration of the principle for the semi-automatic detection of SSEs using PYEQ. The time series are divided into successive time windows which overlap each other by the half of their duration. We can divide the subduction interface into several equivalent subdivisions and perform an independent set of inversions for each of them. We obtain m x n kinematic inversions, where m are the number of subdivisions and n the number of time windows.

Figure 5 . 3 :

 53 Figure 5.3: Seismo-tectonic context of Ecuador. Interseismic coupling model from Nocquet et al. (2017). The blue rectangles show the regions where slow slip events and seismic swarms were documented. The blue curves indicate the M w 7.8 Pedernales earthquake slip contoured every 2 m(Nocquet et al., 2017). The blue star shows its epicenter from Global CMT(Dziewonski et al., 1981). The dashed lines are the Slab2.0 iso-depth contours every 20 km(Hayes et al., 2018).

Figure 5 .

 5 Figure 5.3 shows the seismo-tectonic context of Ecuador. The Nazca plate subducts beneath the North Andean Sliver at 46.4 mm.yr -1(Nocquet et al., 2017;[START_REF] Jarrin | Current motion and deformation of the Nazca plate: new constraints from GPS measurements[END_REF]. The Ecuadorian fault interface is heterogeneously coupled[START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF] Nocquet et al., 2017). The Colombian-Ecuadorian part of the subduction interface (i.e. at latitudes 1 • N-2.5 • N) is highly coupled from the trench down to 20 km depth. From 20 to 60 km depth, the interseismic coupling gradually diminishes. Further south, at latitudes 0.5 • N-0.5 • S, the fault is coupled in a southwest-northeast orientation, which is the direction of the Carnegie ridge subducted at the same area. At the latitude of La Plata island (1 • S), the fault is also highly coupled from the trench down to 25 km depth. Finally, at Salinas (latitude 2 • S), a small part of the megathrust is intermediately coupled and correlates with a subducted geometric complexity of the Nazca plate.

Figure 5 . 4 :

 54 Figure 5.4: SALN observed time series located at Salinas. The time series is in the North Andean Sliver referenceframe. The interseismic velocity was removed following the MIDAS approach(Blewitt et al., 2016). The gray dots indicate the daily displacement together with the 1-σ error bar.

  Figure 5.5: GNSS network used for the Ecuador detection shown in Figure 5.6. Black triangles are the GNSS sites.

Figure 5 . 6 :

 56 Figure 5.6: Geodetic moment rate of 2-month long kinematic inversions every month from January 2010 to November 2022 as red bars. The blue and cyan dots are respectively the SALN and SEEC east component time series with respect to the overriding plate. The vertical blue scale is the eastern displacement of the time series.

Figure 6 . 1 :

 61 Figure 6.1: Illustration of the behavior suggested by several earthquakes. a Recurrent slow slip events, seismic swarms and repeating earthquakes occurring in several areas, predominantly aseismic. b

  [START_REF] Vigny | Upper plate deformation measured by GPS in the Coquimbo Gap, Chile[END_REF]. Upper plate deformation measured by GPS in the Coquimbo Gap, Chile. Physics of the Earth and Planetary Interiors, 175 (1), 86-95. Voss, N. K., Malservisi, R., Dixon, T. H., & Protti, M. (2017). Slow slip events in the early part of the earthquake cycle. Journal of Geophysical Research: Solid Chapter 3. The M w 8.3 2015 Illapel afterslip spatio-temporal modelling (Tissandier et al., 2023)

	Earth, 122 (8), 6773-6786.
	Wessel, P., Luis, J., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H., & Tian, D.
	(2019).	The generic mapping tools version 6.	Geochemistry, Geophysics,
	Geosystems, 20 (11), 5556-5564.
	Xiang, Y., Yue, J., Jiang, Z., & Xing, Y. (2021). Spatial-temporal properties of af-
	terslip associated with the 2015 Mw 8.3 Illapel earthquake, Chile. Earth, Plan-
	ets and Space, 73 (1), 1-15.

-30-110

Remerciements

. The dashed lines are the Slab2.0 depths every 20 km (Hayes et al., 2018). -24-manuscript submitted to JGR: Solid Earth Acknowledgments RT PhD was funded by the Ministère de l'Enseignement supérieur, de la Recherche et de l'Innovation. This research was supported by the French Agence National de la Recherche (ANR) through the project S5 (Grant Number ANR-19-CE31-0003-01). JO acknowledges support from the ANID scholarship ANID-PFCHA/Doctorado Nacional/2020-21200903. SR thanks ANID/FONDECYT; project no. 1200779. Numerical computations were partly performed on the S-CAPAD platform, IPGP, France. Figures were made using Global Mapping Tools (GMT) version 6 (Wessel et al., 2019). We thank two anonymous reviewers for their suggestions that helped to improve the manuscript.

Acknowledgements

We would like to warmly thank all CSN staff, for their precious help for the fieldwork. We also are thankful to the CSN, IRIS and IPOC for providing the data used in this study, as well as to the Institut National des Sciences de l'Univers (INSU-CNRS) and the Réseau Sismologique & Géodésique Français (RESIF, as part of the "Investissements d'Avenir" program, ANR-11-EQPX-0040, and the French Ministry of Ecology, Sustainable Development and Energy) for providing the geodetic instruments for campaigns and 3 broadband seismometers. Part of the computations presented in this paper were performed using the GRICAD infrastructure (https:// gricad .univ-grenoble -alpes .fr), which is supported by Grenoble research communities. Finally, we would like to thank our reviewer for the constructive comments which helped us improving our manuscript.

Fundings

This work was supported by the Agence Nationale de la Recherche (projects ANR-19-CE31-0003), the CNRS-INSU Tellus program and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement N • : 805256). SR and JCB were supported by Fondecyt project (N • 1200779, ANID, Chile). JCB was also supported by ANID PIA (ACT192169).

Figures have been made using Generic Mapping Tools GMT (Wessel et al., 2013) using topography from ETOPO5. Python toolboxes used: PYACS+ PYEQ (https://github .com /JMNocquet /pyacs36, CSI (http://www.geologie .ens .fr /~jolivet /csi) and Pyrocko (https:// pyrocko .org/).

Availability statement

The GNSS data used in this study are in the Supplemental Material. The continuous GNSS data solution is regularly updated, for more information please refer to Klein et al. (2022). The kinematic slip inversion code [START_REF] Nocquet | PYEQ: a software for fault slip kinematic inversion from GNSS time series[END_REF] is available on Zenodo at https://zenodo.org/record/7467256#.Y6Ljb7KZNb8.

Data availability

Seismic data collected are available through the Incorporated Research Institutions for Seismology (IRIS) Data Management Center. The results of the W-phase analysis are available in the supporting information. The relocated catalog is available in the supporting information.

The coseismic offset tables extracted from daily and HRGPS presented in the study are in the supporting information. Position time series of the sequence can be made available upon request.

Supporting Information for "Afterslip of the M w 8.3 2015 Illapel earthquake imaged through a time-dependent inversion of continuous and survey GNSS data" R. Tissandier 1, * , J.-M. Nocquet 1,2 , E. Klein 3 , C. Vigny 3 , J. Ojeda 1,4 , S. Ruiz The focal mechanisms correspond to the 2020 Atacama mainshock and aftershocks from Global CMT (Dziewonski et al., 1981). The blue curves indicate the Copiapó 2014 slow slip every 10 cm (Klein, Duputel, et al., 2018). The dashed lines are the Slab2.0 iso-depth contours every 20 km (Hayes et al., 2018). b) Comparison of the slip distributions of the mainshock in blue, the early afterslip in red, the second aftershock in cyan and the 22-day afterslip in yellow, the 2014 slow slip in green with the interseismic coupling model of [START_REF] Klein | Bridging the gap between North and Central Chile: insight from new GPS data on coupling complexities and the andean sliver motion[END_REF] and the 2020 slow slip [START_REF] Klein | Return of the Atacama deep slow slip event: The 5-year recurrence confirmed by continuous GPS[END_REF]. The white dots show the seismicity from a relocated catalogue. The black lines are the Slab2.0 iso-depth contours every 10 km (Hayes et al., 2018).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary material

Supplementary material related to this article can be found online at https://doi .org /10 .1016 /j .epsl .2021.117081.