Gaelle

David

Malcon, Josipa, Marquito, Thibaut Julia Hugo

René

Andrés Nour

Andrés Hernandez

Estelle Barajas

Hester

Ana Lara

Ginno, Dilrubaer Arianna Giuseppe

Anne Socquet

Jean Paul Ampuero

María Leoncio Cabrera
email: leoncio.cabrera@univ-grenoble-alpes.fr

Sergio Ruiz

Piero Poli
email: pieropoli85@gmail.com

Eduardo Contreras-Reyes

Axel Osses

Renzo Mancini

Constanza Flores

Juan Carlos Báez

Jean B Ammirati

Joaquín Vásquez

María Constanza Flores

Josquín Vasquez

William B Frank

Chapter 3 -Depth Earthquakes Controlled by

Keywords: Intraslab Intermediate Depth Earthquakes, Aftershock, Chile, Subduction, Slab Earthquake dynamics. Earthquake source observations. Seismicity and tectonics. Seismic tomography

my apologies if I have left anyone out, but you know you are part of this). Thank you for all the conversations, hugs, words of encouragement and love. Thank you for every scientific and non-scientific discussion as well. Many thanks also to my friends in Chile, you have been part of my family for several years. Thank you for your support in spite of the distance, different time zones, and so many other difficulties.

Many thanks to my

iii

Abstract

Seismic sequences are usually classified into three types: Mainshock-Aftershocks, Swarm and Foreshocks-Mainshock-Aftershocks. However, which are the physical processes that control them is still not well understood (e.g., static/dynamic stress transfer, fluids, aseismic slip, or a combination of processes). By studying three types of sequences in different seismotectonic settings using high-resolution seismic catalogues in combination with statistical seismology, modelling, and geodetic observations among others, we aim to better understand the physical process driving the seismicity and their role during fault slip. In the first part, we analyze the variations of the source properties and aftershock activity of six ~Mw 6.3 intermediate depth earthquakes in the subduction zone of northern Chile to characterize the mainshock-aftershock process. We show that while the mainshocks exhibit similar rupture geometry and stress drop, the aftershock productivity systematically decreases for the deeper events within the slab, especially below the 400 450°C isotherm depth. We propose that this isotherm separates highand low-hydrated zones, thus controlling the aftershock productivity. Subsequently, we study a seismic swarm that occurred in the Antarctica. We create a seismic catalog of ~36,000 events (Aug/2020-Jun/2021). In addition, we observe a prominent geodetic deformation signal at a nearby GNSS station. Based on the dynamics of the seismicity and the geodetic deformation, we infer a volcanic origin for this swarm, which occurred close to a ridge axis and the Orca Volcano. In the third part we start the study of crustal seismic sequences in Central Italy by a seismic catalogue of ~5,000 events starting ~3.5 months before the mainshock. We observe that the precursory phase experiences multiple accelerations of the seismicity rate that we divide into two main sequences with different features. While the first part is characterized by weak earthquake interactions, smooth moment release, slow spatial migration patterns, and a lower effective stress drop, the second sequence exhibits the opposite behavior. We interpret these differences as distinct physical processes (aseismic slip in the first part and stress transfer in the second one) that are controlled by different physical properties of the fault system. In the fourth part we focus on the rupture initiation of the L'Aquila earthquake. We characterize a ~0.6-s signal preceding the large dynamic rupture. Our results indicate that the geometry of the rupture initiation is in agreement with that of the main event and the fault system. Furthermore, we estimate that the rupture velocity of this rupture initiation ( km/s) is slower than for conventional earthquakes, with a seismic efficiency of 0.24. This latter indicates that ~80% of the energy budget corresponds to fracture energy. Based on this, we interpret that the rupture struggled to start, likely due to the complexities of the medium. Finally, we explore the use of unsupervised learning (clustering) in the study of a dense seismic catalogue. We start by creating the first and longest seismic catalog (~280,000 events) that connects both spatially and -Visso-Norcia sequence (Mw 6.0-5.9-6.5). By using different physical parameters of the seismicity and clustering algorithms, we can classify the seismicity in an automated and data-driven way. We obtain different seismicity groups (clusters) along the fault system and others scattered in the region with different physical properties (e.g., interactions, durations). Thus, our methodology presents a powerful tool in the study, classification and analysis of massive seismic catalogs. v

Résumé

Il y a trois types de séquences sismiques: des séries évenements séisme principal-réplique, des essaims et des séries précurseurs-séisme principal-répliques. Malheureusement, le mécanisme physique qui les contrôle n'est encore compris (e.g., glissement asismique). En étudiant trois types de séquences dans différents contextes géologiques, en suivant une approche sismologique, et en incorporant des observations géodésiques et la modélisation nous cherchons à mieux comprendre le processus physique. Premièrement, nous analysons les variations des propriétés de six séquences séisme principal-répliques de profondeur intermédiaire de magnitudes ~Mw 6.3, localisées dans la partie nord de la zone de subduction chilienne. Nous montrons que lorsque les séismes principaux présentent une géométrie de rupture et une chute de contrainte similaires, le taux des répliques diminue systématiquement pour les événements les plus profonds, notamment en dessous de la profondeur de l'isotherme 400-450 °C. Nous proposons que cet isotherme sépare les zones fortement et faiblement hydratées, contrôlant ainsi la productivité des répliques. Ensuite, nous étudions un essaim qui s'est produit dans le bassin de Bransfield, en Antarctique. Nous avons construit un catalogue avec ~36000 événements (août/2020-juin/2021). Nous observons aussi un signal géodesique à une station GNSS proche. Nous expliquons cet essaim par une source volcanique, localisée près d'un ride et du volcan Orca. Dans la troisième partie, nous conduisons une étude des séquences sismiques en Italie centrale en analysant la phase précurseur du séisme de (Mw6.1), en 2009. Pour ce propos, nous avons créé un catalogue sismique de ~5000 événements, qui commence ~3,5 mois avant le séisme principal. Nous observons que la phase précurseur subit de multiples accélérations dans le taux de sismicité, que nous pouvons diviser en deux séquences principales. Alors que la première partie est caractérisée par de faibles interactions entre les séismes, une libération continue du moment, des schémas de migration spatiale lents, et une baisse de la contrainte effective plutôt faible, la deuxième séquence présente le comportement opposé. Nous interprétons ces différences comme des mécanismes physiques distincts (glissement asismique dans la première partie et transfert de contraintes dans la seconde), contrôlés par des propriétés physiques différentes du système de failles. Dans la quatrième partie, nous nous concentrons sur l'initiation de la rupture du séisme de L'Aquila. Nous caractérisons un signal de ~0.6-s précédant la grande rupture dynamique enregistrée dans 15 stations de mouvement fort. Nos résultats indiquent que la géométrie de l'initiation de la rupture est en accord avec celle du séisme principal et du système de failles. Nous estimons que la vitesse de rupture de cette initiation de rupture (0,9±0,2 km/s) est plus lente que pour les séismes conventionnels, avec une efficacité sismique de 0,24. Cette dernière indique que ~80% du bilan énergétique correspond à l'énergie de rupture. Sur cette base, nous interprétons que la rupture a démarré lentement, très probablement à cause de l des propriétés de la faille. Enfin, nous explorons l'utilisation de l'apprentissage non supervisé (clustering) dans l'étude d'un catalogue sismique dense. Nous commençons par créer le premier et le plus long catalogue sismique (~280000 séismes) qui relie à la fois spatialement et temporellement le séisme de 2009 à L'Aquila et la séquence de 2016 à Amatrice-Visso-Norcia. Nous montrons qu'en utilisant différents paramètres physiques de la sismicité et des algorithmes de clustering, nous pouvons classer la sismicité de manière automatisée. Nous obtenons différents groupes de sismicité (clusters) le long du système de failles et en dehors avec des propriétés physiques différentes (comme interactions et durées). Nous pouvons alors suggérer, que notre méthodologie permet de séparer automatiquement la sismicité correspondant à différents processus. vi vii 3.1 Resume of kinematic parameters and number of aftershocks obtained for each mainshock Table 6.1 Parameters retrieved for the rupture initiation considering a line source model xii

List of Figures and Tables

Chapter 1 Introduction

Since its origins, one of the main tasks of seismology has been to characterize earthquakes. One of the fundamental elements for this has been the creation of seismic catalogs, which basically correspond to a compilation of earthquake parameters such as origin time, location and magnitude, among others. The analysis of these seismic catalogs has multiple applications in seismology, ranging from scientific purposes such as the study of the precursory phase of earthquakes (e.g., [START_REF] Kato | The generation of large earthquakes[END_REF] to risk assessment (e.g., Gerstenberger et al., 2016). The analysis of seismic catalogs has also shown a fundamental property of earthquakes: they cluster in time and space (e.g., [START_REF] Beaucé | Systematic detection of clustered seismicity beneath the Southwestern Alps[END_REF][START_REF] Gardner | Is the sequence of earthquakes in southern California, with aftershocks removed, poissonian?[END_REF][START_REF] Knopoff | The statistics of earthquakes in southern California[END_REF][START_REF] Zaliapin | Earthquake clusters in southern California II: Classification and relation to physical properties of the crust[END_REF][START_REF] Utsu | A statistical study on the occurrence of aftershocks[END_REF]. In addition, different types of earthquake sequences (which are described below) have been identified based on their spatiotemporal characteristics. Thus, seismic catalogs and the observations about the clustering of earthquakes are a key tool to provide us valuable insights into the physical processes driving seismicity (e.g., [START_REF] Beaucé | Systematic detection of clustered seismicity beneath the Southwestern Alps[END_REF][START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas[END_REF].

With the recent increment of continuous seismological data (e.g., USArray [START_REF] Kerr | Geophysical exploration linking deep Earth and backyard geology[END_REF]; AlpArray in the Alps [START_REF] Hetényi | The AlpArray seismic network: A large-scale European experiment to image the Alpine orogen[END_REF]), development of advanced data analysis techniques (e.g., [START_REF] Mousavi | Machine Learning in Earthquake Seismology[END_REF] in combination with novel computational ways (e.g., GPUs), we now enter the era of high-resolution catalogs [START_REF] Arrowsmith | Big Data Seismology[END_REF]). These new catalogs can improve our observations by providing better spatial and temporal resolution, and through its analysis we can study the dynamics of seismicity and better understand the physical process (e.g., aseismic slip, fluids, stress transfer) that drives earthquake sequences.

Patterns of earthquake sequences

Seismic sequences are classified from the combination of two main factors: the occurrence or not of an event whose magnitude is prominently larger than the majority of the rest of the events with respect to the mainshock (e.g., occurrence of seismicity after the mainshock, before the mainshock, or both). If we plot the time evolution of the different earthquakes occurring in a sequence, three types of earthquake sequences can be distinguished: Mainshock(MS)-Aftershocks, Foreshocks-Mainshock-Aftershocks and Swarms, which illustrated in Fig. 1.1 and described below [START_REF] Mogi | Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (third paper)[END_REF][START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF].

Mainshock-Aftershocks

Aftershocks is the label for earthquakes occurring after a mainshock and their number and magnitude decrease gradually over time [START_REF] Mogi | Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (third paper)[END_REF]. Aftershocks had traditionally been considered a manifestation of only stress redistribution once mainshock took place given the common spatial overlap between the location of aftershocks and the regions of increased stress after mainshock (e.g., [START_REF] Stein | The role of stress transfer in earthquake occurrence[END_REF]. Figure 1.2 shows an example of the spatial overlap between regions of increased Coulomb stress (red colors) after the occurrence of the 1994 Mw=6.7

Northridge earthquake and the reported aftershocks (black dots). However, the emergence of better-quality aftershock catalogs shed light on the occurrence of other physical processes such as afterslip (e.g., Frank et al., 2017;[START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF][START_REF] Perfettini | A model of aftershock migration driven by afterslip[END_REF]. Mainshock-Aftershocks-type sequences are probably the most common and wellstudied seismic sequences [START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF]. In this type of sequences, a general approach is to characterize the source of the mainshock, and the spatial and temporal distribution of the aftershocks with respect to that source (e.g., [START_REF] Freed | Earthquake triggering by static, dynamic, and postseismic stress transfer[END_REF] in order to understand the physical properties of the rupture (e.g., stress drop, amount and distribution of the slip, among others), and the physical mechanisms that control the subsequent aftershocks. An interesting observation that has been reported regarding the interaction between mainshocks and aftershocks is the difference in the aftershock productivity (e.g., [START_REF] Marsan | How variable is the number of triggered aftershocks?[END_REF]. This observation becomes even more evident when shallow earthquakes are compared with deeper earthquakes, the latter being the ones with the lowest aftershock productivity [START_REF] Frohlich | Aftershocks and temporal clustering of deep earthquakes[END_REF][START_REF] Wiens | A deep earthquake aftershock sequence and implications for the rupture mechanism of deep earthquakes[END_REF][START_REF] Zhan | Supershear rupture in a Mw 6.7 aftershock of the 2013 Sea of Okhotsk earthquake[END_REF][START_REF] Houston | Deep earthquakes[END_REF].

The above-mentioned differences might reflect differences in the physical properties of the medium where the aftershocks take place (e.g., temperature, stresses, presence of fluids), which has direct implications for hazard assessment after a mainshock has occurred (see for example [START_REF] Dascher Cousineau | What controls variations in aftershock productivity?[END_REF]. Therefore, understanding which properties of the medium control these differences is an essential task.

Foreshocks-Mainshock-Aftershocks

Foreshocks are small earthquakes that sometimes precede a large mainshock [START_REF] Omori | On the Fore-shocks of Earthquakes[END_REF][START_REF] Bouchon | The long precursory phase of most large interplate earthquakes[END_REF][START_REF] Dodge | Detailed observations of California foreshock sequences: Implications for the earthquake initiation process[END_REF]. Thus, Foreshocks-Mainshock-Aftershocks corresponds to earthquake sequences where the seismicity in the epicentral region of the incoming mainshock increases prior to the earthquake, and it is also followed by aftershocks once the mainshock takes place (Mogi, 1963, Fig. 1.1b). They have been observed in different seismotectonic settings such as strike-slip faults (e.g., [START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF][START_REF] Chen | California foreshock sequences suggest aseismic triggering process[END_REF][START_REF] Dodge | Detailed observations of California foreshock sequences: Implications for the earthquake initiation process[END_REF][START_REF] Dodge | Foreshock sequence of the 1992 Landers, California, earthquake and its implications for earthquake nucleation[END_REF][START_REF] Durand | A Two Scale Preparation Phase Preceded an M w 5.8 Earthquake in the Sea of Marmara Offshore Istanbul, Turkey[END_REF][START_REF] Ellsworth | Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks[END_REF][START_REF] Shelly | A high resolution seismic catalog for the initial 2019 Ridgecrest earthquake sequence: Foreshocks, aftershocks, and faulting complexity[END_REF][START_REF] Tape | Earthquake nucleation and fault slip complexity in the lower crust of central Alaska[END_REF][START_REF] Yoon | Foreshocks and mainshock nucleation of the 1999 M w 7.1 Hector Mine, California, Earthquake[END_REF], subduction zones [START_REF] Bouchon | The long precursory phase of most large interplate earthquakes[END_REF]Gardonio et al., 2020;[START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake[END_REF][START_REF] Ruiz | Nucleation phase and dynamic inversion of the Mw 6.9 Valparaíso 2017 earthquake in Central Chile[END_REF], and extensional regimes [START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF]. An example of the foreshocks-mainshock-aftershocks sequence for the 2017 Mw 6.9 Valparaiso earthquake (Chile) is presented in Fig. 1.3. We can see how the seismicity preceding the mainshock started in the northern part of the epicenter, and after the mainshock occurred a larger segment (-33.3° to -32.9°) is covered by aftershocks. In addition, foreshocks have also been observed in laboratory experiments (e.g., [START_REF] Acosta | Can precursory moment release scale with earthquake magnitude? A view from the laboratory[END_REF][START_REF] Mclaskey | Foreshocks during the nucleation of stick slip instability[END_REF][START_REF] Mclaskey | Preslip and cascade processes initiating laboratory stick slip[END_REF][START_REF] Mclaskey | Earthquake initiation from laboratory observations and implications for foreshocks[END_REF]. Thus, the existence of foreshocks has led to the idea that the occurrence of foreshocks could be related to process occurring before a large earthquake.

Different physical models have been proposed to explain the processes that lead to large seismic events. One of them is the cascade model (Ellsworth & Beroza, 1995), in which an earthquake is the consequence to stress perturbations triggered by previous earthquakes. In this model, the earthquakes themselves provide the triggering mechanism that initiates failure of the next event [START_REF] Gomberg | Unsettled earthquake nucleation[END_REF]. On the other hand, the pre-slip model proposes that the occurrence of a large earthquake starts with an expanding slow slip event, which eventually reaches a critical size after which the mainshock occurs. In this case, foreshocks are a consequence of the slow slip event (Ellsworth & Beroza, 1995). Other factors such as localization of deformation have been also proposed [START_REF] Kato | The generation of large earthquakes[END_REF]. However, which of these mechanisms (if any of them) best represents the physics of the precursory phase of earthquakes is still under debate, with examples in which the same earthquake is studied using different approaches, and opposite conclusions are obtained (e.g., [START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF]Ellsworth &[START_REF] Ellsworth | Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks[END_REF].

Seismic Swarms

Earthquake swarms corresponds to sequences of earthquakes for which the number of events often increases with time, and then gradually decreases after a certain period. In addition, no single earthquake dominates in size [START_REF] Mogi | Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (third paper)[END_REF]Scholz, 2019, Fig. 1.1c). An example of a seismic Swarm in the Corinth Gulf, Greece is presented in Figure 1.4. We can see how the number of events (Fig. 1.4b) and the magnitudes (Fig. 1.4c) gradually increase and decrease over time.

Earthquake swarms have been observed in different tectonic settings such as volcanic and geothermal regions (e.g., [START_REF] Cattania | A nonplanar slow rupture episode during the 2000 Miyakejima dike intrusion[END_REF][START_REF] Klein | The Reykjanes Peninsula, Iceland, earthquake swarm of September 1972 and its tectonic significance[END_REF], transform faults (e.g., [START_REF] Roland | Earthquake swarms on transform faults[END_REF][START_REF] Vidale | A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers[END_REF], subduction zones (e.g, [START_REF] Valenzuela-Malebrán | Seismicity clusters in Central Chile: investigating the role of repeating earthquakes and swarms in a subduction region[END_REF] and extensional regimes (e.g., [START_REF] Essing | Spatiotemporal evolution of the Seismicity in the Alto Tiberina Fault System revealed by a Template Matching Catalog[END_REF]. Unlike aftershocks which usually occur immediately all over the fault plane and along the edges of mainshock rupture, earthquake swarms show spatial migrations starting from a small region [START_REF] Fischer | The Growth of Earthquake Clusters[END_REF]. Based on this migratory pattern, earthquake swarms are thought to be primarily driven by external processes such as transient increases in fluid pressure and/or aseismic slip (e.g., Chen et al., 2012;[START_REF] De Barros | Imbricated aseismic slip and fluid diffusion drive a seismic swarm in the Corinth Gulf, Greece[END_REF][START_REF] Vidale | A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers[END_REF]. However, the primary processes driving seismic swarms are still under debate [START_REF] Valenzuela-Malebrán | Seismicity clusters in Central Chile: investigating the role of repeating earthquakes and swarms in a subduction region[END_REF]. Therefore, understanding the driving physical processes of swarms is fundamental for the characterization of the seismic cycle and earthquake hazard mitigation [START_REF] De Barros | Imbricated aseismic slip and fluid diffusion drive a seismic swarm in the Corinth Gulf, Greece[END_REF].

Multiple Physical Processes Controlling Seismicity

The features of each type of earthquake sequence described above show how each sequence is controlled by different physical processes or combinations of them. Moreover, sequences of the same type, but occurring in different tectonic settings show differences (e.g., aftershock productivity), indicating that, in addition to the process, the environment in which the sequence occurs plays a fundamental role. Thus, characterizing these earthquake sequences allows us to better understand the different physical processes and factors at work.

For example, foreshocks are one of the most powerful observational tools to study the physical processes before an earthquake (e.g., [START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF]. Therefore, if we understand what drives foreshocks, we can better understand the process before an earthquake and thus, infer whether or not the processes are explained by the current conceptual models (e.g., Ellsworth & Beroza, 1995). As well as for seismic swarms, we can determine if the process is controlled by the interaction of earthquakes, or if for example aseismic slip also takes place. This allows us to understand the interplay between seismic and aseismic processes during fault slip.

For all of the above-mentioned reasons, we need more and better observations as well as a robust analysis of these observations to obtain quantitative results and provide new insights about the physical mechanisms and the properties of the medium that controls the earthquake sequences previously described. In the following section we briefly describe how we improve our observations by creating high-resolution seismic catalogs, which are then analyzed and complemented with other observations and tools.

Towards High-Resolution Seismic Catalogs

Seismology has a strong observational component. Without going any further, the classification of earthquake sequences such as those described in the previous section emerged from the pure observation of the temporal evolution of the magnitudes of earthquakes, within a specific region and registered in a seismic catalog (e.g., [START_REF] Mogi | Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (third paper)[END_REF]. This highlights the relevance of seismic catalogs, which when analyzed in detail can reveal valuable information about the physical processes taking place and the characteristics of the medium. For example, [START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF] improved the catalog of foreshocks before the 2011 Mw 9.1 Tohoku earthquake using a waveform correlation technique (template matching, see more detail in Chapter 2). Thus, they identified two distinct sequences of foreshocks migrating at rates of 2 to 10 km/day along the trench axis toward the epicenter of the mainshock (Fig. 1.5), suggesting slow-slip transients during the precursory phase of the mainshock. Using the same technique, Frank et al. (2017) analyzed the aftershocks of the 2015 Mw 8.3 Illapel earthquake to map the rheology of the Central Chile subduction zone. Their results suggested that afterslip acted as the main driver of aftershock activity, and they were also able to estimate the rate and state as a function of depth. As a last example, De Barros et al. ( 2020) tracked the spatiotemporal evolution of a seismic swarm in the Corinth Gulf, Greece (Fig. 1.4). Based on a detailed analysis of the swarm migration, they conclude that the swarm was the result of a combination of multiple driving processes such as fluid diffusion, aseismic slip and stress transfer. Earthquake migration toward the rupture initiation point of the mainshock. Space-time diagram of all detected events between 13 February and the mainshock origin time, with earthquake origin locations indicated in terms of the distance along the trench axis (blue circles scaled to magnitude). Red dashed lines, approximate locations of the fronts of earthquake migration; red stars, repeating earthquakes in the JMA catalog; green stars, newly detected events that were found to resemble those repeating events; black star, Mw 9.0 mainshock; yellow star, Mw 7.3 largest foreshock. (Inset) Time variations in seismicity rates inside and outside the EMZ after the Mw 7.3 largest foreshock. The blue dashed curve denotes the least-squares fitting of the modified Ohmori law. Modified from [START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF].

The above-mentioned examples show that knowing in detail the location, magnitude and time of occurrence of as many earthquakes as possible leads to a more detailed and comprehensive understanding of the different seismic sequences. Therefore, in this thesis we use high-resolution seismic catalogs as the main input for the further analysis. To create such catalogs, we use already existing seismic catalogs and a method called template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF]. The basic idea of template matching is to correlate template waveforms (of earthquakes reported in the initial catalog) with continuous data in a sliding window that preserves the seismic moveout, in order to detect new events that were not reported in the initial catalog (e.g., De Barros et al., 2019;Frank et al., 2017;[START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF]Shelly, 2017). Since template matching plays a fundamental role in this work, it is described in more detail in the next chapter. Once the new and denser catalog is obtained, we proceed to analyze it with different tools (e.g., statistics, modelling) and complementary observations (e.g., geodetic signals). We thus obtain new and independent results for the different seismic sequences studied here.

Goals and Structure of the Thesis

The scope of this work is to provide new observations and analysis methods to better understand the physical processes associated with different types of earthquake sequences. In each chapter described below, particular scientific questions are addressed. Each chapter has the format of a scientific article (except Chapter 2), since Chapters 3 to 5 (and Appendix 1) correspond to articles already published. On the other hand, Chapter 6 is already submitted, and Chapter 7 is soon to be submitted.

Chapter 2 since this is the main method used throughout the thesis for the creation of new high-resolution seismic catalogs. Likewise, each subsequent chapter summarizes the method and the particular specifications used (e.g., number of stations, frequency band, templates, etc.). In addition, the template matching results are complemented in each chapter with other tools and observations, which are pertinently described in each chapter.

In [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF], in which we take advantage of three similar earthquakes occurring at the same latitude of the Chilean subduction zone, but at different locations along-dip. Since the magnitudes of these events are higher (~Mw 6.8), we also performed source dynamic inversions.

In Chapter 4 -Lived Swarm in the Central Bransfield Basin, smic swarm that started in August 2020 and was located in a very complex back-arc basin at the transition from rifting to ocean spreading, and close to the submarine Orca volcano. In this work, we address the question of whether the swarm is a product of tectonic extension or a volcanic process. In addition, given the limited availability of data in such a remote region, our work shows how to use single-station approaches to obtain conclusive results from the complement between seismic and geodetic data.

In Chapter 5 -Temporal Evolution of Foreshocks Preceding the seismic sequences in Central Italy. In particular, this chapter deals with the question of how we took advantage of the permanent network of broadband stations in Central Italy, which has good spatial and temporal coverage. Thus, we created a dense seismic catalog of foreshocks, which we analyzed using different statistical tools. Based on this we can observe differences between the behavior of the seismicity occurring on two different faults before the mainshock, which we attribute to different physical mechanisms and different fault rheologies. Chapter 6 deals with the question of how was the early rupture initiation of the L'Aquila earthquake. To do so, we take advantage of a dense and well azimuthally distributed network of strong motion stations. These stations are located in the near field, and recorded with a high sampling rate of 200 Hz a small emergent signal (~0.6s duration) that we attribute to the onset of the rupture. Through modeling and analysis of the source of this signal, we propose that the onset of the rupture was slow and complex, which is represented by low seismic efficiency. Based on the results of this chapter and Chapter 5, as well as previous work by other authors, we propose a conceptual model to explain the rupture initiation.

Chapter 7 Clustering): 8 Years of Seismicity in Central Italy (2009to analyze a dense seismic catalog in an automated and data-driven way. We explore the use of unsupervised learning (clustering) to study the first and longest seismic catalog (~280,000 Amatrice-Visso-Norcia sequence (Mw 6.0-5.9-6.5) created by us using template matching. By using different physical parameters of the seismicity and clustering algorithms, we obtain different seismicity groups (clusters) along the fault system and others scattered in the region. Based on this, we interpret that our methodology allows to automatically separate seismicity corresponding to different processes (e.g., background seismicity, mainshock-aftershocks).

Finally, a summary of the main findings and conclusions are presented in the last Chapter (8) "Conclusions and Outlook", as well as future perspectives.

Chapter 2 Template Matching

As mentioned in the introduction, the general strategy we follow in this thesis to study earthquake sequences is based on the generation of a new high-resolution seismic catalog, which is then analyzed in detail with different tools (e.g., statistics, modelling). To achieve higher spatio-temporal resolution we need to detect earthquakes that have not been reported in the local seismic catalogs. This brings us to the problem of earthquake detection, which consists of determining when an instrument (e.g., seismometer) or a network of instruments registers an earthquake-related signal, and not just ambient noise. In addition, once a new earthquake is detected, its origin time, location and magnitude are required for the compilation of the catalog.

Over the years, different methods have been used in seismology to address the problem of earthquake detection (for a recent review see [START_REF] Arrowsmith | Big Data Seismology[END_REF]. The most classic approach would be the visual detection of events in continuous records by an analyst. This corresponds to the work routinely performed in monitoring centers and seismological observatories. In general, this approach provides a good spatio-temporal coverage of earthquakes over a certain magnitude of completeness that will depend among other things on the instrumental region coverage (e.g., [START_REF] Barrientos | & National Seismological Center (CSN) Team[END_REF]. However, many earthquakes of smaller magnitude are missed due to the difficulty of identifying them within the noise.

The next level for earthquake detection could be using power-based detection algorithms, such as the short-term average/long-term average (STA/LTA) detector (Allen, 1982). However, such methods have difficulties to discriminate earthquakes from transient noise sources, and the detection gets even more complicated when the amplitudes of target seismic events are on the same order of magnitude as the ambient seismic noise, which is measured by the signal-to-noise ratio (SNR).

One way to overcome the SNR issue is to analyze several seismic stations at the same time, taking advantage of the coherent information in the network. One method that does exactly this is template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF], which has demonstrated amazing performance in detecting new events and densifying seismic catalogs (e.g., [START_REF] Chamberlain | Illuminating the and Phases of the 2016 M7. 8 Earthquake With 10 Years of Seismicity[END_REF]Shelly, 2017). This method has two elements as necessary input: an initial catalog of events and continuous records from several stations simultaneously. Although single-event initialization and then recursive expansion of the detections have also been explored, see [START_REF] Wimez | A recursive matched-filter to systematically explore volcanic long-period earthquake swarms[END_REF]. Since our study areas already have existing catalogs (see detail in each chapter) and continuous records at several stations, we used template matching for the densification of the initial catalogs.

The General Approach of Template Matching

When two collocated (or very close) earthquakes occur with the same focal mechanism, their waveforms will exhibit similarity and therefore high correlation coefficient (CC). We can thus detect new events in the vicinity of known ones by using template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF], an algorithm that computes the average CC between the waveforms of a reference event (called template) considering their respective moveouts, and the continuous data at multiple stations (Fig. 2.1). The similarity between templates and continuous data recorded by a seismic array can be expressed as: (2.1) In which and are the indexes for the station, the channel, and the time sample, respectively (Fig. 2.1 shows an example where nine stations with one channel each are used);

is the user-designated weight associated with station and component . If the sum of the weights is normalized to one, then attributing equal weights to every station and component results in computing the arithmetic mean. is the length of the template in samples (Fig. 2.1b,c show examples of templates). and are the template and the continuous seismic data, respectively (Fig. 2.1b); and is the moveout on station and channel . The moveout describes the delay with which an event is observed on every station relative to the beginning of the earliest 2.1b,c). It is important to note that in equation ( 2.1) it is assumed that the waveforms are centered around zero, which allows for a large speed up of the algorithm but makes it inappropriate for template matching on raw data [START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF]. Therefore, a pre-processing of the raw data (which is described in the next section) is needed. Once the pre-processing is done, the general process can be summarized as in Figure 2.1. is not sufficient to declare a detection, while the case of (f) is (see details on how to define the threshold in the text).

Using equation ( 2.1) with continuous data and templates results in a detection trace as illustrated in Figure 2.2 for each template. This trace corresponds to a time series of the average correlation coefficient ( ) between the template and continuous data. Subsequently, a threshold value of the must be defined to declare the new detections done by the template (Fig. 2.2b). When using several templates from a catalog, the same event can be detected by more than one template. Therefore, a declustering process must be carried out to eliminate multiple detections. For this purpose, a time window is usually defined and the event with the highest average CC within the window is retained. We describe more details performed after the detections are declared in the next section (post-processing). Now that we have seen how template matching works, we will mention some steps of the different stages prior to the calculations (Sect. 2.2), during the calculation stage (Sect. 2.3) and after the calculations (Sect. 2.4). These steps must be considered as part of the process, and although they are covered in detail in each further chapter, we summarized them here so that the reader can become more familiar with the whole process.

Pre-processing

The most laborious stage during the creation of seismic catalogs using template matching is the preparation of continuous data. Once the raw data is collected, the following processing is applied:

Check the days with data available for each station. Days with no data can be ignored in the processing by assigning a weight of 0 in equation (2.1). Deconvolution of the instrumental response to account for instrument changes.

Remove mean, trend and filter in the desired frequency range (In general, filtering aims to remove the noisiest frequency bands in the region). Check and fill data gaps. Downsampling. Although optional, in cases where it is desired to scan using thousands of templates, several stations and years of continuous data, it can greatly decrease the computational load. Save processed data. e.g., 24-hour files.

Create the templates waveforms from the continuous data. To do this, it is necessary to know the times at which each waveform must be cut. Some examples of how to obtain such information are: manual picks from a monitoring center, estimation of theoretical arrival times for a local velocity model (e.g., ObsPy TauP Toolkit [START_REF] Beyreuther | ObsPy: A Python toolbox for seismology[END_REF][START_REF] Crotwell | The TauP Toolkit: Flexible seismic traveltime and ray-path utilities[END_REF]), or using some automatic picking algorithm (e.g., Phasenet, Zhu & Beroza, 2018). From these same times the moveout information is created. Assess the quality of the templates. In general, only templates that meet some criteria are retained e.g., being registered on a minimum number of stations/channels, having an SNR greater than a chosen threshold. It is also important to check that the templates correctly sample the study area, since the main limitation of template matching is that it will improve the catalog spatially represented by the templates. The following chapters detail how each initial catalog was spatially selected.

Once the continuous data, templates, moveouts and weights have been constructed, it is possible to proceed to use the FMF algorithm.

Running Template Matching

After data processing, we can proceed to the stage that generally consumes the most computational time: the calculation of the correlations. For this task, we used the Fast Matched Filter (FMF) algorithm [START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF]. FMF has two versions implemented to perform the calculation of the correlations: one version written in C that can run on compute nodes typically comprised of several central processing units (CPUs), and another Nvidia CUDA C version that runs on computers with Nvidia graphics processing units (GPU) [START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF]. Because of its higher speed, for our calculations we use the version implemented for GPUs, and the High-Performance Computing infrastructures CIMENT of the Université Grenoble Alpes (UGA, https://gricad-doc.univ-grenoble-alpes.fr, accessed 21 November 2022). For this thesis we used approximately 1,400 hours of parallel GPU computations.

Post-processing

Once the daily CCs functions are calculated for each template, the next steps are:

Define a threshold to declare a new detection. Commonly this threshold is defined as NxMAD, where MAD is the median absolute deviation of the daily average CCs, and N is an empirically defined number. Several studies have shown that values of N in the range 8-12 provide reliable results (e.g., [START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF]Frank et al., 2017;[START_REF] Ross | Searching for hidden earthquakes in Southern California[END_REF]. The higher the N, the more conservative the threshold will be, but also the fewer events will be detected. Conversely, a very low N can lead to more detections, but also to more false positives. In Chapter 5, we propose a way to better constrain this value, based on the assessment of the number of false positives obtained from the scan with templates flipped in time.

Declustering of redundant detections. As previously mentioned in Sec. 2.1, this procedure is done to eliminate multiple detections of the same event from several templates. Thus, a time window usually of the same duration as the templates is defined and only the event with the highest average CC within the window is retained.

Magnitude estimation of the new detections. This is usually done by comparing waveform amplitudes between the templates and their respective detections. For example, [START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF] computes the magnitude of a detected event based on the median value of the maximum amplitude ratios for all channels between the template and the detected event, assuming that a tenfold increase in amplitude corresponds to one unit increase in magnitude.

Location of the new events. Usually, the initial location of a new detection is assumed to be the same location of the template that detected it. Although sometimes it is also possible to improve such locations using for example double difference algorithms such as HypoDD (Waldhauser & Ellsworth, 2000) or GrowClust [START_REF] Trugman | GrowClust: A hierarchical clustering algorithm for relative earthquake relocation, with application to the Spanish Springs and Sheldon, Nevada, earthquake sequences[END_REF].

In this thesis we will follow the procedure described above to improve local catalogs of the different sequences studied (more detail is given in each following chapter). By doing so, we obtain new high-resolution catalogs as the ones exemplified in Chapter 1 (Fig. 1.3, 1.4, 1.5). This type of catalogs enables a better monitoring of seismic dynamics to quantitatively study for example the interactions between earthquakes, changes in the seismicity rate, spatial migrations, among others. Thus, through the combination of our high-resolution seismological observations, together with other observations (e.g., GNSS) and modelling, we can achieve a better understanding of the physical processes that drive seismicity and the mechanics of the faults where it occurs.

Introduction

Rupture mechanism of intermediate-depth earthquakes (IDEs) is one of the biggest open questions in seismology [START_REF] Frohlich | Deep earthquakes[END_REF][START_REF] Houston | Deep earthquakes[END_REF]. They occur mostly at subduction zones in the so called double seismic zones (DSZ), where two parallel seismicity planes are observed. The upper seismicity plane (USP) is inferred to occur within the subducting oceanic crust and/or the upper mantle due to dehydration reactions [START_REF] Hacker | Subduction factory 2. Are earthquakes in subducting slabs linked to metamorphic dehydration reactions?[END_REF][START_REF] Kirby | Interslab earthquakes and phase changes in subducting lithosphere[END_REF]. In contrast, the lower seismicity plane (LSP) occurs in the lithospheric mantle, and its mechanism is still under debate [START_REF] Ferrand | Dehydration-driven stress transfer triggers intermediatedepth earthquakes[END_REF][START_REF] Ohuchi | Intermediate-depth earthquakes linked to localized heating in dunite and harzburgite[END_REF][START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle?[END_REF][START_REF] Reynard | Earthquakes and plastic deformation of anhydrous slab mantle in double zones[END_REF]. Hypothesis for this process include dehydration-embrittlement of antigorite [START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle?[END_REF], reactivation of pre-existing shear zones [START_REF] Reynard | Earthquakes and plastic deformation of anhydrous slab mantle in double zones[END_REF] and quasi-adiabatic shear-heating instabilities [START_REF] John | Generation of intermediate-depth earthquakes by self-localizing thermal runaway[END_REF]. Furthermore, laboratory experiments also suggest that faulting at intermediate depths may occur under dry conditions [START_REF] Ohuchi | Intermediate-depth earthquakes linked to localized heating in dunite and harzburgite[END_REF] or by dehydration-driven stress transfer under partially hydrated mantle conditions [START_REF] Ferrand | Dehydration-driven stress transfer triggers intermediatedepth earthquakes[END_REF].

The DSZ is usually characterized by compressional deep events within the oceanic lithospheric mantle in trench-outer rise regions prior to subduction, while the upper part of the oceanic lithosphere is affected by a tensional regime. This stress field in the trench-outer rise region is controlled by the plate bending caused by the interplay of slab pull forces and the rheology of the oceanic lithosphere (Contrereas-Reyes & Osses, 2010). Furthermore, based on a worldwide compilation of trench outer rise earthquakes, the position of the neutral plane of the stress field within the oceanic lithosphere have been estimated to be approximately coincident with the 400º-450ºC isotherm-depth [START_REF] Seno | Double seismic zones, compressional deep rise events, and superplumes[END_REF]. Likewise, Ruiz & Contreras-Reyes (2015) and Carrasco et al., (2019) found that tensional outer rise events are confined to the upper oceanic lithosphere to isotherm-depths shallower than 400ºC offshore central Chile. Off Japan, [START_REF] Kita | Existence of interplane earthquakes and neutral stress boundary between the upper and lower planes of the double seismic zone beneath Tohoku and Hokkaido, northeastern Japan[END_REF] reported that the transition from upper tensional to compressional intraplate earthquakes occur at 22 km beneath the base of the Eurasian plate beneath Tohoku. However, this transition occurs only ~10 km beneath the upper plate beneath Hokkaido. The authors attribute this difference to the complexity of the oblique subduction beneath Hokkaido, and they further pointed out that the stress regime is the combination of two stress fields: (i) the bending/unbending force, and (ii) the thermomechanical and petrological forces.

Although IDEs are characterized by shear slip on faults similar as inland crustal or interplate earthquakes, anomalous behavior is often observed for the first ones such as thermaldependence radiated seismic energies [START_REF] Wiens | Seismological constraints on the mechanism of deep earthquakes: Temperature dependence of deep earthquake source properties[END_REF], significant non-double-couple mechanisms [START_REF] Richardson | Low-frequency properties of intermediate-focus earthquakes[END_REF], and high stress drops [START_REF] Frohlich | Deep earthquakes[END_REF][START_REF] Poli | Global and variations of source duration and scaling for and earthquakes[END_REF][START_REF] Houston | Deep earthquakes[END_REF] among others. A remarkable characteristic of intermediate and deep earthquakes is their general lower aftershock productivity as compared with shallower events [START_REF] Frohlich | Aftershocks and temporal clustering of deep earthquakes[END_REF][START_REF] Wiens | A deep earthquake aftershock sequence and implications for the rupture mechanism of deep earthquakes[END_REF][START_REF] Zhan | Supershear rupture in a Mw 6.7 aftershock of the 2013 Sea of Okhotsk earthquake[END_REF][START_REF] Houston | Deep earthquakes[END_REF], although few cases with large activity exist [START_REF] Wiens | A deep earthquake aftershock sequence and implications for the rupture mechanism of deep earthquakes[END_REF]. Aftershocks activity has been related to fault properties as state of stress (Shebalien & Narteau, 2017), rupture processes (Poli et al., 2016), seismic coupling [START_REF] Hainzl | Linear relationship between aftershock productivity and seismic coupling in the Northern Chile subduction zone[END_REF] and the thermal state of the slab, with colder slabs showing a greater number of aftershocks [START_REF] Wiens | Effect of slab temperature on deep-earthquake aftershock productivity and magnitude frequency relations[END_REF]. A main issue is to understand if this general decrease in the number of aftershocks is real, or it is due to the lack of near-field data [START_REF] Li | Abundant aftershock sequence of the 2015 M w7. 5 Hindu Kush intermediate-depth earthquake[END_REF]. It is thus fundamental to tackle this issue, to better assess the aftershock activity [START_REF] Dascher Cousineau | What controls variations in aftershock productivity?[END_REF] and get new insights about the mechanism(s) and physical conditions controlling the occurrence of IDEs.

Since 2010, six intermediate-depth events with magnitudes (Mw) greater than 6 occurred in northern Chile between 17ºS and 21ºS (Figure 3.1 and Table A2.S1), where the Centro Sismológico Nacional de la Universidad de Chile (CSN; [START_REF] Barrientos | & National Seismological Center (CSN) Team[END_REF], Integrated Plate Boundary Observatory Chile (IPOC; GFZ, & CNRS-INSU, 2006) and Central Andean Uplift and the Geodynamics of the High Topography (CAUGHT; [START_REF] Beck | Central Andean uplift and the geodynamics of the high topography. Other/Seismic Network[END_REF] seismic networks have been deployed (Figure 3.1a). These events are of main interest for several reasons. First, they are registered by a dense network of broad-band and strong-motion sensors in the near-field providing an opportunity to avoid problems with detection levels (Poli et al., 2016;Shebalien & Narteau, 2017;[START_REF] Zhan | Supershear rupture in a Mw 6.7 aftershock of the 2013 Sea of Okhotsk earthquake[END_REF]. Second, the six events are located in a very limited volume, reducing issues related with attenuation and other propagation effects, which are possibly affecting the estimation and thus the comparison of source properties [START_REF] Prieto | Earthquake nests as natural laboratories for the study of intermediate-depth earthquake mechanics[END_REF]. Finally, their magnitudes (Mw ~ 6.3) and focal mechanisms are similar (Figure 3.1). We thus have a unique opportunity to compare source properties and aftershock activity, with limited bias due to propagation and rupture geometry. In addition, according to the CSN catalogue they exhibit different aftershock productivity. For the shallowest event (E1 in Figure 3.1c), there are more than 40 aftershocks reported, including some of them with magnitude equal or greater than Mw 4.0 (see Supplementary Information). For the deeper ones there is a decrease in the number of aftershocks with no aftershocks for the deepest event (E6 in Figure 3.1c).

To better characterize if differences in between these events exist, we study their aftershock activity and rupture properties, using template matching (TM) and kinematic inversions (KI). Furthermore, we developed a thermal model using the finite element method (FEM), with the aim of analyze the impact of slab temperature on seismic behavior. In the next sections, the tectonic setting and methodologies used in this study followed by their respective results are described, and finally the discussions and conclusions are presented. [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF]. Black and green dashed lines correspond to continental Moho [START_REF] Tassara | Anatomy of the Andean subduction zone: threedimensional density model upgraded and compared against global-scale models[END_REF] and oceanic Moho [START_REF] Contreras-Reyes | Abrupt change in the dip of the subducting plate beneath north Chile[END_REF], respectively.

Tectonic Setting

The northern Chilean convergent margin is characterized by the subduction of the oceanic Nazca plate beneath the continental South American plate (Figure 3.1) at a currently convergence rate of ~67 mm/a [START_REF] Khazaradze | effects of GPS measured crustal deformation rates along the south central Andes[END_REF], although the average convergence rate during the last 20 Ma is about 85 mm/a [START_REF] Demets | Geologically current plate motions[END_REF]. The Nazca plate in the study region (18º -21ºS) was formed in the East Pacific Rise 45-52 Ma (Mueller et al., 2008), whereas its crustal thickness is 6-7 km thick inferred by seismic constraints [START_REF] Patzwahl | velocity models of the Nazca plate subduction zone between 19.5° S and 25° S from seismic measurements during the CINCA95 project[END_REF]Ranero and Sallar s, 2004;[START_REF] Contreras-Reyes | Abrupt change in the dip of the subducting plate beneath north Chile[END_REF][START_REF] Maksymowicz | The Geometry of the Continental Wedge and Its Relation to the Rheology and Seismicity of the Chilean Interplate Boundary[END_REF]. The northern Chilean margin has a poorly sedimented trench (< 200 m;[START_REF] Maksymowicz | The Geometry of the Continental Wedge and Its Relation to the Rheology and Seismicity of the Chilean Interplate Boundary[END_REF] due to the reduced sediment supply from the Andes in the extremely arid region of the Atacama Desert. The continental slope is steep as the frontal part of the margin lacks a well-developed frontal prism (Contreras-Reyes 2018; [START_REF] Maksymowicz | The Geometry of the Continental Wedge and Its Relation to the Rheology and Seismicity of the Chilean Interplate Boundary[END_REF].

The northern Chilean margin has undergone subduction erosion since at least the Jurassic, which is inferred from the eastward migration of the volcanic arc (e.g. [START_REF] Rutland | Andean orogeny and ocean floor spreading[END_REF][START_REF] Stern | Subduction erosion: rates, mechanisms, and its role in arc magmatism and the evolution of the continental crust and mantle[END_REF] as well as the long-term arcward retreat of the trench, crustal thinning and subsidence of the outer forearc (e.g., [START_REF] Kukowski | Subduction erosion The mode of fore-arc material transfer along the Chilean margin?[END_REF][START_REF] Contreras-Reyes | Structure and tectonics of the Chilean convergent margin from wide-angle seismic studies: a review[END_REF]. Longterm subduction rates for northern Chile have been estimated with values of 40-45 km3Myr-1km-1 [START_REF] Kukowski | Subduction erosion The mode of fore-arc material transfer along the Chilean margin?[END_REF]. On the other hand, the oceanic Nazca plate presents well developed tensional faults seen in the bathymetric data [START_REF] Geersen | Active tectonics of the North Chilean marine forearc and adjacent oceanic Nazca Plate[END_REF], which were formed by bending of the oceanic plate [START_REF] Ranero | Relationship between at trenches and seismicity[END_REF]. These bend faults provide the pathways for seawater infiltration reaching the oceanic mantle according to seismic and geodynamic studies (see Contreras-Reyes & Osses, 2010 and references therein).

Aftershock productivity

Detection of aftershocks using template matching

The six studied IDEs have tensional focal mechanism and span a range of depth from 90 to 130 km (Figure 3.1), in a region where a double seismic zone (DSZ) has been previously reported [START_REF] Comte | A seismic zone in Arica, northern Chile[END_REF][START_REF] Dorbath | Double seismic zone of the Nazca plate in northern Chile: velocity structure, petrological implications, and thermomechanical modeling[END_REF][START_REF] Sippl | Seismicity structure of the northern Chile forearc from> 100,000 relocated hypocenters[END_REF][START_REF] Florez | Controlling factors of seismicity and geometry in double seismic zones[END_REF] and the seismic catalog has magnitude of completeness ~ Ml 3.5 [START_REF] Barrientos | & National Seismological Center (CSN) Team[END_REF]. It is thus likely that aftershocks are missing for some the studied events, which is improved using template matching (TM).

TM has been widely used to detect events in aftershocks sequences, where seismicity is known to repeat [START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF][START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF]Frank et al., 2017;[START_REF] Li | Abundant aftershock sequence of the 2015 M w7. 5 Hindu Kush intermediate-depth earthquake[END_REF]. In a similar way, we here apply TM to improve the detection of aftershocks following the analyzed events (Figure 3.1). We base our analysis on continuous three-component velocity seismological data from the CSN, IPOC, and CAUGHT seismic networks, using the five nearest stations to the source region of each event (Figure A2.S15). Data were continuously recorded at a sampling rate of 100 Hz for IPOC and CSN stations, and 40 Hz for CAUGHT stations. To detect non-reported aftershocks, we searched for events with similar waveforms to those reported as mainshock and aftershocks for the CSN (template candidates, see Supplementary Information). A space-time window needs to select reference events and build templates. A large window implies including more events, but also more background events as well, while a smaller window decreases this effect, but is more likely to miss some aftershocks. Since the six earthquakes have similar characteristics such as magnitude, focal mechanism, hypocenter and station coverage, we use the same spatial and temporal scale for all of them in order to preserve a homogeneous criterion that allows their comparison. After some tests with different windows size and considering some previous works (e.g., [START_REF] Persh | Strongly depth-dependent aftershock production in deep earthquakes[END_REF][START_REF] Dascher Cousineau | What controls variations in aftershock productivity?[END_REF] we define the spatial limit using a 3D radius of 25 km from each hypocenter, and 25 days of data after each mainshock, as this is the maximum number of days for which the five stations closest to each event were operating continuously. Both continuous data and template waveforms were bandpass filtered from 5 to 10 Hz because this frequency range exhibits better signal to noise ratios, and decimated to 25 Hz. To select high-quality template events, we follow Frank et al. (2017) and we estimate the signal-to-noise ratio (SNR) on the vertical component as the ratio between the RMS velocity during the first 25 s of the P-wave and the RMS velocity during a 25 s of noise before the P-wave arrival time. Only template candidates with SNR > 5 in at least 3 stations are used as templates. We define the resulting template waveforms as 30

-wave arrival at each station. Then we calculate the correlation coefficients between template waveforms and continuous data in a sliding window that preserves the seismic moveouts (Frank et al., 2017), to get a time series that represents the similarity of the continuous data to the matched-filter template. We search sample-by-sample considering a detection threshold that is 10 times the median absolute deviation (MAD) of the correlation sum to detect event significantly similar to the template. Events detected with this criterion, are then considered to occur at the same hypocenter (determined by the CSN) as their template, and we estimate their magnitude by computing the median amplitude ratio between the template event and the aftershock over the 5-station network, assuming that a tenfold increase in amplitude corresponds to one unit increase in magnitude. Figures 3.2b-c show an example for E2 and E5, where using the mainshocks waveforms as templates, two new detections with inferred magnitudes (Mw) of 1.8 and 2.7 occurred four days (for E2) and some hours (for E5) after mainshocks.

To avoid redundancy of events detected by multiple templates, we remove detections within 30 s of another detections, by keeping the event with the highest network correlation coefficient. Finally, we use the new catalogue (including CSN and new detected aftershocks) as templates candidates, and the process is repeated once.

Aftershock Activity

Figure 3.3a shows all new detections and catalog events, for 25 days after each mainshock. As is common for aftershock sequences, the number of events decreases over time, with a higher concentration close to the mainshock. On the other hand, the number of aftershocks shows a clear reduction as function of depth, down to E6 that is showing no aftershocks . The most productive event is E1, located just ~7 km from the top of slab (Figure 3.1). For this earthquake, we used 45 events as templates, and we found 2044 aftershocks, nearly 47 times more.

It is important to highlight that there are no aftershocks reported in the CSN catalog for E2, E5, and E6. For these events, we use the mainshock waveforms for detection. This strategy permits to find new events for E2 and E5, suggesting that the sole main event waveforms can improve the detection. Figure 3.2 shows an example of detections for E2 and E5 using mainshock waveforms as templates (see Table A2.S2 for more details about detections). According to the above, lack of detection seen for E6 is a robust feature of our analysis.

Thermal structure

Thermal-Fluid Numerical Model

To assess the thermal conditions at the depth of the analyzed events, we developed a 2D-thermal model using Finite Element Method (FEM) for northern Chile, along the profile (Figure 3.1). The model is constrained by the plate geometry of SLAB2.0 [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF], hypocenter data, and published thermal parameters in the area (see Table A2.S5).

Numerical Model

The thermal subduction model consists in a two-dimensional finite element system for the subducting lithospheric plate speed and asthenospheric fluxes. The heat equation is given by:

(3.1)
in the region of interest , where is the temperature and the first term corresponds to heat diffusion where is the thermal diffusivity depending on space (ratio of the thermal conductivity and the thermal capacity of the media), the second term corresponds to the heat transport due to lithospheric speed and asthenospheric fluxes , and the last term stands for the external heat sources such as radiogenic or frictional heat. No friction heat effects are considered in this study since in the case young lithospheres the effects of friction can be neglected (Vöelker et al., 2011). Thermal effects due to the subducted eroded material above the incoming/subducting lithosphere are also neglected for the sake of model simplicity.

Given a velocity field in the whole region , we obtain the steady state temperature as the limit of solution of the following variational problem discretized with Lagrange finite elements of type P2 (space ) and iterating for starting from until convergence, where is some fictitious numerical time step:

(3.2)
Given the units we use, we observe that it is sufficient to make the order of 15 iterations with a time step of 10 myrs corresponding to 150 myrs in total in order to reach a steady state.

In order to compute the heat transport field , we consider zero velocity in the continental crust and we impose a rigid and constant movement along the parallel shape lines of the incoming plate. We model the fluxes produced in the mantle wedge using a divergence free Stokes flux, considering both momentum and mass conservation:

(3.3) (3.4)
in the region occupied by the mantle wedge, where and are the velocity and pressure of the fluid flow and is the kinematic viscosity of the asthenosphere in the wedge mantle. Thus, only for the wedge mantle, the velocity field is computed from the ) and P1 for the pressure (space ) by solving:

(3.5)

Thermal Model Description (Boundary Conditions)

For the heat equation, we consider that the temperature on the edge condition of the incoming/subducting plate increases in deep by: (3.6) where =0 and =1350 ºC are the surface and mantle temperatures, where is the Gaussian error function, is the age of the plate (50 myrs) and is the thermal diffusivity of the incoming plate. This is a good approximation for deep temperature profiles for low aged plates that is the case for the entire Chilean trench (see [START_REF] Stein | A model for the global variation in oceanic depth and heat flow with lithospheric age[END_REF]Vöelker et al., 2011). For the other boundary conditions (Figure 3.4), we assume that varies quadratically from 1350º to 1400 ºC in interior of the right boundary of the mantle wedge, =0 at the surface, =1350 ºC at the bottom part of the subducting lithosphere and zero heat flux in the rest of the boundaries.

in the continental crust and in the subducting lithosphere where =67 mm/yr is the imposed constant speed of the incoming plate with the unit tangent vector to the subducting lithosphere shape lines depending on depth . The boundary conditions for the flow in the mantle wedge are on the upper part of the wedge, on the lower part of the wedge and free flow in the right edge of the wedge. See Text A2.S1 in the Appendix 2 for more details.

Thermal Model Description (Parameter Selection)

Topography and bathymetry are taken from published data [START_REF] Ryan | Global topography synthesis[END_REF]. The geometry of the subducted lithosphere is taken from Slab2.0 [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF]. In addition, to assess our model, we consider three different upper continental mantle depths of 60, 80 and 100 km. However, as they show roughly similar results for the location of the mainshocks, we consider the more realistic model of 60 km [START_REF] Wada | Common depth of decoupling: Reconciling diversity and uniformity of subduction zones[END_REF]. Figures A2.S17 and A2.S18 show models with 80 and 100 km, respectively). A detail list of the other parameters used to develop the two-dimensional model through finite element are listed in Table A2.S5 where values are taken from several authors [START_REF] Hyndman | Thermal constraints on the zone of major thrust earthquake failure: The Cascadia subduction zone[END_REF][START_REF] Herzberg | Temperatures in ambient mantle and plumes: Constraints from basalts, picrites, and komatiites[END_REF][START_REF] Wang | Thermal regime of the Southwest Japan subduction zone: effects of age history of the subducting plate[END_REF][START_REF] Oleskevich | The updip and downdip limits to great subduction earthquakes: Thermal and structural models of Cascadia, south Alaska, SW Japan, and Chile[END_REF]Vöelker et al., 2011[START_REF] Höink | Viscous coupling at the asthenosphere boundary[END_REF] for rheologic and geometric values. Please note that heat diffusivity takes different values in three different regions: the subduction lithosphere, the upper part of the continental crust and sediments and the rest of the domain. For the sources of temperature, we also consider different radiative sources in the upper and lower continental crust and sediments.

Thermal Conditions of the Slab

From our results of the slab thermal-conditions in northern Chile (Figure 3.3d), we observe that E1 and its aftershocks occurs near the 400 ºC isotherm-depth and does not exceed the 450 ºC isotherm, while the rest of the events occurs around or below this isotherm. Figures 3.3b and 3.3c illustrate the relationship between the temperature of the slab at which events occur, the distance from the top of the slab and number of aftershocks, respectively. Our analysis clearly reveals an inverse relationship between the number of aftershocks and temperature. It is important to highlight that E6, the deepest within the slab, occurs at the highest temperatures range between 700 800 ºC, and is the only event that does not exhibit aftershocks.

Seismic Source Properties

Seismic Source Kinematic Inversions

Having observed a clear difference of aftershock productivity as function of depth, we now attempt to resolve any difference in rupture properties. To obtain the slip distribution of each mainshock, we perform kinematic inversions using near-field stations on hard rock [START_REF] Leyton | Geophysical characterization of the Chilean seismological stations: First results[END_REF]) (Figure 3.1), including 19 strong-motion stations and one broad-band station (Figure A2.S16). Since it is not possible to distinguish the causative fault plane from aftershocks geometry, we consider the two possible fault planes, according to moment tensors reported by the USGS (Tables A2.S3 and A2.S4). We assume a finite-fault model with an apriori elliptical-patch slip distribution. This methodology has been successfully used by several authors [START_REF] Madariaga | Earthquake dynamics on circular faults: A review 1970 2015[END_REF][START_REF] Ruiz | The January 2019 (M w 6.7) Coquimbo Earthquake: Insights from a Seismic Sequence within the Nazca Plate[END_REF]. We inverted seven parameters. Five of them, correspond to geometrical characteristic: semi-axes and of the ellipse, the rotation angle of the ellipse and the location of its center into the fault plane, and the other two parameters are maximum slip and rupture velocity . From this inverted rupture model, we compute dynamic parameters using circular crack approximation obtained averaging both axes of the ellipse (Ruiz et al., 2019 and references there in). Strong motion and broad band records were corrected for the instrumental response and linear trend, filtered using a causal Butterworth filter with corner frequencies of 0.02 and 0.1 Hz, and then integrated to displacement. It is important to mention that the high corner frequency is controlled by our simple elliptical model and the 1-D velocity model [START_REF] Husen | Accurate hypocenter determination in the shallow part of the Nazca subduction zone in Northern Chile using a combined on-/offshore network[END_REF] used to simulate the wave propagation from the source to the receivers. The misfit between observed and synthetic records was computed using a L2 norm starting from origin time up to 75 s later to avoid including surface waves, and considering both possible rupture planes reported by the USGS (Tables A2.S3 and A2.S4).

Kinematic Parameters

Considering all six events, our results show similar geometries of rupture (Table 3.1), with an approximately circular rupture ranging from ~3.5 to ~5 km. Note: Green and yellow rows indicate parameters obtained using nodal plane 1 (NP1) and nodal plane 2 (NP2).

For more details see Tables A2.S4 and A2.S5.

Figure 3.5 shows an example of the kinematic inversion results, for event E1 considering NP1. The observed and synthetic displacement waveforms show a good agreement, with a maximum value about 0.4 cm for the East-West component in station A20P (waveforms for all events in Supplementary Figures A2.S1-A2.S12). In addition, although there are differences in the misfit values for each of the fault planes considered (Tables A2.S3 and A2.S4), these are small enough to give preference to one or the other. The estimated stress drop values vary between 7.5 MPa and 29.5 MPa (Table 3.1), which are high compared with typical values for shallow earthquakes [START_REF] Derode | Energy Budget of Earthquakes in Northern Chile: Comparison With Shallow Earthquakes and Implications of Rupture Velocity Models Used[END_REF], but similar to others observed for IDEs in Chile (Ruiz & Madariaga, 2011;[START_REF] Ruiz | The January 2019 (M w 6.7) Coquimbo Earthquake: Insights from a Seismic Sequence within the Nazca Plate[END_REF] and at global scale [START_REF] Poli | Global and variations of source duration and scaling for and earthquakes[END_REF][START_REF] Prieto | Earthquake nests as natural laboratories for the study of intermediate-depth earthquake mechanics[END_REF]. Despite a variability of rupture parameters exists (Table 3.1), we do not see any clear correlation between them and the depth of the events. We thus suggest that the analyzed events, while occurring under different thermal conditions, all have a similar rupture physics.

Discussion

While co-seismic rupture properties do not vary with distance from the top of the slab, we observe clear differences in the post-seismic activity pattern, with a decrease of the aftershocks as the distance from the top of the slab and temperature increase (700-800 ºC isotherm-depths).

Difference in the number of aftershocks between E1 and the rest of events could be related to the hypocenter-depth within the slab. E1 is located ~7 km from top of slab, close to the oceanic Moho (Figure 3.1b;[START_REF] Patzwahl | velocity models of the Nazca plate subduction zone between 19.5° S and 25° S from seismic measurements during the CINCA95 project[END_REF]Ranero and Sallar s, 2004). In addition, E1 aftershocks locations reported by the CSN delineate a potential pre-existing fault aligned and similar to those extensional faults located in the outer-rise region (Figure A2.S13). These tensional faults are likely caused by plate bending, and could provide a pathway for fluid infiltration, producing hydration into the crust and uppermost mantle [START_REF] Boneh | earthquakes controlled by incoming plate hydration along faults[END_REF]Contreras-Reyes et al., 2008;[START_REF] Iyer | Controls of faulting and reaction kinetics on serpentinization and double Benioff zones[END_REF]. However, how deeply water can infiltrate is still under debate, since the buoyancy of water (or equivalently, confining pressure) makes it difficult to bring water down even when faulting is deep [START_REF] Korenaga | On the extent of mantle hydration caused by plate bending[END_REF]. By contrast, numeric models suggest that stress variation induced by plate bending can produce subhydrostatic or even negative pressure gradients along normal faults, favoring thus downward pumping of fluids (Faccenda et al., 2009). On the other hand, previous studies in the DSZ of this region [START_REF] Dorbath | Double seismic zone of the Nazca plate in northern Chile: velocity structure, petrological implications, and thermomechanical modeling[END_REF] observed that the upper seismicity plane (USP) corresponding to oceanic crust, is characterized by intermediate Vp (~7.7 km/s) and low Vp/Vs (1.67) values, concluding that the USP is related to fluid releases associated with metamorphic reactions occurring within jadeite-lawsonite blueschists. Similarly, it has been observed that the USP has significantly larger b-values than the LSP [START_REF] Florez | Controlling factors of seismicity and geometry in double seismic zones[END_REF], suggesting an hydrated USP and a relatively dry lithospheric mantle, since previous observations have shown that high bvalues anomalies correlate well with regions where dehydration reactions are expected [START_REF] Wiemer | Mapping the anomaly at 100 km depth in the Alaska and New Zealand subduction zones[END_REF] or with hydrated fracture zones [START_REF] Schlaphorst | Water, oceanic fracture zones and the lubrication of subducting plate boundaries insights from seismicity[END_REF]. Thus, this hydrated zone could lead to a greater occurrence of aftershocks.

Furthermore, we observe a dominant tensional stress regime, where mainshocks with normal focal mechanism occur up to ~40 km from the top of the slab (i.e., 700 800 ºC isotherm). This isotherm-depth interval is much deeper than the expected depth for this type of focal mechanism compared to neutral plane of the stress field estimated by [START_REF] Seno | Double seismic zones, compressional deep rise events, and superplumes[END_REF], who show that the transition from tensional to compressional of outer rise events occur at the 400º-450 ºC isotherm-depth. If this is correct, our results show that the plate bending becomes stronger during subduction and is likely more predominant than the cyclic stresses associated to the seismic coupling at the thrust zone at depths where IDEs occur [START_REF] Seno | Double seismic zones, compressional deep rise events, and superplumes[END_REF]. Furthermore, seismic tomographic studies off south-central Chile (Contreras-Reyes et al., 2008) and off Nicaragua [START_REF] Lefeldt | Intraplate seismicity and related mantle hydration at the Nicaraguan trench outer rise[END_REF] show a reduction in the upper oceanic mantle velocities from the oceanic Moho up to 400º-450 ºC isotherm-depth. This lower bound has been interpreted as the maximum potential depth for hydro-alteration within the upper part of the oceanic lithosphere, where extensional stresses dominate, and fluids may not be able to penetrate any deeper than neutral plane (Contreras-Reyes et al., 2008;[START_REF] Lefeldt | Intraplate seismicity and related mantle hydration at the Nicaraguan trench outer rise[END_REF]Faccenda et al., 2009). Another proposed interpretation is that these seismic velocity anomalies could be related to small-crack porosities, which can be produced by a mixture between thermal cracking and bending-related faulting [START_REF] Korenaga | On the extent of mantle hydration caused by plate bending[END_REF]. As described above, our observations could suggest the existence of a deepening of the neutral plane from 400-450 ºC to 700-800 ºC isotherm-depths regardless of the presence of deep-water percolation.

A remaining question is, what could be controlling this deepening of the neutral plane along-dip? Northern Chile corresponds to the region where the oceanic Nazca plate is relatively old (Figure 3.1) and therefore colder and with a deep fragile-ductile system (Figure 3.6), with an age of 54 Ma near the trench axis [START_REF] Müller | Age, spreading rates, and spreading asymmetry of the world's ocean crust[END_REF] and 58 Ma estimated for the area where the mainshocks occur. A plausible mechanism for a deepening of the neutral plane along-dip in northern Chile would be the increase in bending stresses due to a greater slab-pull associated to the relative old, cold and heavier already subducted oceanic Nazca plate. In particular, in this segment of the Chilean margin, the largest known IDE has been registered for the slab-pull Tarapacá 2005 Mw 7.7 event [START_REF] Peyrat | Tarapacá earthquake (Mw 7.7, 2005, northern Chile): A event with horizontal fault plane constrained from seismologic and geodetic observations[END_REF]Legrand and Delouis, 2006;[START_REF] Kuge | Rupture characteristics of the 2005 Tarapaca, northern Chile, earthquake: Evidence for heterogeneous fluid distribution across the subducting oceanic plate?[END_REF][START_REF] Ruiz | Historical and recent large megathrust earthquakes in Chile[END_REF]. Moreover, the 1950 Ms 8.0 tensional event [START_REF] Kausel | The Ms= 8 tensional earthquake of 9 December 1950 of northern Chile and its relation to the seismic potential of the region[END_REF]) occurred a little further south, reflecting the tensional character of the stress field within the subducting slab. We also note how the moment tensors in northern Chile DSZ show a downdip extensive stress field throughout the seismically active area [START_REF] Sippl | Filling the gap in a double seismic zone: Intraslab seismicity in Northern Chile[END_REF], which may imply that plate bending, and unbending forces are small compared to slab-pull forces. Another mechanism previously proposed by [START_REF] Kita | Existence of interplane earthquakes and neutral stress boundary between the upper and lower planes of the double seismic zone beneath Tohoku and Hokkaido, northeastern Japan[END_REF] to explain the difference in the location of the neutral plane beneath Tohoku and Hokkaido in the northeastern Japan is related to differences in the thermal structure between the two regions. They concluded that the thermal age of the subducted lithosphere and the oblique component of convergence rate affect the density of the mantle wedge. This process triggers buoyancy force variations exerted in the slab affecting finally the stress field of the subducting slab. However, this mechanism should be more dominant along strike (thermal variations are more important in that direction) rather than in along dip. Yield stress is limited at the base of the plate by steady state creep, which depends on the cube of the stress and exponentially on temperature [START_REF] Kohlstedt | Strength of the lithosphere: Constraints imposed by laboratory experiments[END_REF], and computed considering steady state flow properties for olivine (green line; [START_REF] Goetze | The mechanisms of creep in olivine[END_REF], dry olivine (yellow line; [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF] and wet olivine (cyan line; [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF]. Isotherms were computed based on the cooling of a semi-infinite half-space mode (Turcotte & Schubert, 2002).

On the other hand, the occurrence of tensional events up to depths of 40 km as E6 is consistent with the brittle region predicted by a yield strength envelope (YSE) for an oceanic lithosphere under tensional stresses (Figure 3.6), which reaches depths of 40-42 km for olivine and dry olivine mineralogy, and it decreases for wet conditions. Likewise, using high-precision relocations, a DSZ south of Iquique (21.5 °S) is observed [START_REF] Rietbrock | A narrowly spaced zone in the subducting Nazca plate[END_REF], with predominantly extensional faulting in both seismic layers, and a clear down-dip tensional regime with T-axes oriented in slab-parallel directions [START_REF] Rietbrock | A narrowly spaced zone in the subducting Nazca plate[END_REF][START_REF] Sippl | Filling the gap in a double seismic zone: Intraslab seismicity in Northern Chile[END_REF]. The consistency between the strength envelope and the maximum depths of intraplate tensional events suggests that these earthquakes have ruptured through the entire brittle part of the oceanic lithosphere. Another interesting observation is that events E5 and E6 highlight for its isotherm-depth occurrence, about 700-800 °C, which are higher than the stability limit of antigorite at 600 650 °C commonly proposed as a lower limit for the LSP. As observed by [START_REF] Sippl | Seismicity structure of the northern Chile forearc from> 100,000 relocated hypocenters[END_REF]2019) these events are more common in the Northern Chile DSZ, and a runaway-type process could be responsible, under specific conditions in the lower plate (e.g. composition, hydration, stress state). Likewise, according to [START_REF] Mckenzie | Thermal structure of oceanic and continental lithosphere[END_REF] particularly large strain rates could also lead to earthquakes in regions below the 600 650 °C isotherm-depth. In our study region, this could be related to an increase in the slab pull. Finally, thermal cracking hypothesis [START_REF] Korenaga | On the extent of mantle hydration caused by plate bending[END_REF] could also play a role, since the deepest hydrated cracks could be more widely spaced than the most superficial cracks

Considering previous observations and our own results, we propose a conceptual model shown in Figure 3.7, where the processes that occur are as follows: (1) in the outer-rise region extensional bend-faulting occur and lead to partial hydration of the crust and upper mantle. However, fluids may not be able to penetrate any deeper than neutral plane (approximately 450 ºC isotherm-depth). Thus, the neutral plane separates a high-hydrated from a dry or poorly hydrated zones. ( 2) As the oceanic plate subducts, a deepening of the neutral plane occurs by the increase of the slab-pull forces, separating it from the 450 ºC isotherm, and increasing the region subject to extensional failure. ( 3) Normal (tensional) events occur at different distances from the top of the slab, but their behavior is mostly controlled by variable physical background conditions. In particular, in the high-hydrated region a greater number of aftershocks is observed, differently from what is observed in the dry deeper zone.

Conclusion

We presented a detailed study of seismic source and aftershock productivity for the six largest IDEs occurred in northern in Chile since 2010 and recorded by dense seismic arrays in the near-field. Our results show that, although all of these events are located at different depths and under different thermal conditions, all have similar rupture physics, with analogous geometries considering an elliptical-patch approach and stress drop values between 7 MPa and 30 MPa. On the other hand, a clear decrease of the number of aftershocks as the distance from the top of the slab and temperature increase is observed. We propose that this behavior could be controlled by the incoming plate hydration, where the 400 450 ºC isotherm-depths and neutral plane of the stress field acts as limits for hydration in the outer-rise region, which is deepened by the slab-pull as the slab subducts.

Chapter 4 -Lived Swarm in the Central

Introduction

In August 2020, a significant increase of seismic activity was reported in the Bransfield Basin, south of the South Shetland islands, Antarctica (Figure 4.1). Between 29 August 2020 and June 2021, the United State Geological Survey (USGS, 2020) reported 128 earthquakes with magnitudes larger than 4.0 (Figure 4.1). The seismicity was not characterized by any large mainshock (Figure A3.S1) that could potentially have triggered the prolific occurrence of earthquakes in the region. Similar observations were recently reported by [START_REF] Olivet | A seismic swarm at the Bransfield Rift, Antarctica[END_REF].

The central Bransfield Basin is a ridge separating the Antarctic plate to the South, from the South Shetland microplate [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF][START_REF] Olivet | A seismic swarm at the Bransfield Rift, Antarctica[END_REF]Taylor et al., 2008 and reference therein, Figure 4.1). The NW-SE extension of the ridge results from the combination of slab rollback from the past subduction of the Phoenix microplate under the South Shetland microplate, and transtensional motion between the Scotia and Antarctic plates [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF][START_REF] Gràcia | Morphostructure and evolution of the central and eastern Bransfield basins (NW Antarctic Peninsula)[END_REF][START_REF] Taylor | Kinematics and segmentation of the South Shetland islands-Bransfield Basin system, northern Antarctic Peninsula[END_REF]. The basin is also characterized by extensive volcanism [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF][START_REF] Taylor | Kinematics and segmentation of the South Shetland islands-Bransfield Basin system, northern Antarctic Peninsula[END_REF], occurring in several submarine structures, such as the Orca volcano, that is located 20 km southwest of the seismic swarm (Figure 4.1, Almendros et al., 2020;[START_REF] Olivet | A seismic swarm at the Bransfield Rift, Antarctica[END_REF]. The Orca Volcano consists of a large caldera surrounded by shallow magma reservoirs [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF]. The Orca volcano magmatism is mostly basaltic with mid-ocean ridge characteristics [START_REF] Barker | Rift propagation, detachment faulting, and associated magmatism in Bransfield Strait, Antarctic Peninsula[END_REF] and significant hydrothermal activity is also observed in the area [START_REF] Bohrmann | Hydrothermal activity at Hook Ridge in the central Bransfield Basin, Antarctica[END_REF]. The volcano s caldera and shallow magmatic bodies produce a positive magnetic anomaly extending along the axis ridge [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF]. The 2020 2021 swarm overlaps with the shallow magnetic anomaly located northeast of the Orca volcano [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF][START_REF] Olivet | A seismic swarm at the Bransfield Rift, Antarctica[END_REF]. Several other swarms connected with volcanic activity have also been observed along the Bransfield Ridge [START_REF] Almendros | Volcanotectonic activity at deception island volcano following a seismic swarm in the Bransfield Rift (2014 2015)[END_REF][START_REF] Dziak | Tectonomagmatic activity and ice dynamics in the Bransfield Strait back-arc basin, Antarctica[END_REF], but are mostly located in the western part of the ridge. No significant earthquakes have occurred in the central Bransfield basin since 1970 (Figure 4.1).

The Bransfield Ridge is of major geological interest as it represents a back-arc basin at the transition from rifting to ocean spreading [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF]. Detailed analysis of the 2020 2021 earthquake sequence therefore provides new insights about short term processes associated with the evolution from rifting to spreading, and attempts to distinguish if this swarm, and the related geodetic deformation, is the result of extensional tectonics or volcanic process [START_REF] Bergman | Earthquake swarms on the mid-Atlantic Ridge: Products of magmatism or extensional tectonics[END_REF][START_REF] Buck | Consequences of asthenospheric variability on continental rifting[END_REF][START_REF] Reiss | The impact of complex volcanic plumbing on the nature of seismicity in the developing magmatic natron rift, Tanzania[END_REF]. To address this question, we make use of the limited, but significant data available from the region (Figure 4.1). One seismic station (JUBA, network AI) located at the Carlini Base (King George Island, 20 km from the swarm, Figure 4.1) is used to improve the detection of earthquakes using template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF] and to characterize the waveforms associated which each detected event. We also analyzed data from the nearby GNSS stations (Figure 4.1), to assess any deformation associated with the seismic activity. Despite the limited data available, we were able to characterize the early phase and development of the largest swarm ever observed in the Bransfield Ridge area. Our observations to support the hypothesis of a volcanic origin of this earthquake sequence.

Data and methodology

Earthquakes Detection and Characterization

We use a seismic station located 20 km NW from the swarm centroid (Figure 4.1) to systematically detect and characterize seismic events in the study area. We first downloaded 1 year (1 June 2020, to 1 June 2021) of continuous three-component waveforms recorded at station JUBA (Istituto Nazionale di Oceanografia e di Geofisica Sperimentale, 1992) at 100 Hz sampling rate. The traces are first band-pass filtered between 1 and 9 Hz to enhance high frequency signals from the local seismicity and resampled at 50 Hz. Then, we extracted the events corresponding to the USGS catalog (Figure 4.1) within a 12s window (starting 2 s before P-wave arrival), ensuring the presence of both P-and S-waves associated with the swarm events. We visually controlled the quality of each seismogram and manually picked the P-and S-waves on good quality seismograms exhibiting phases well identifiable on the three components. This process yielded 114 earthquakes, the waveforms of which represent the templates to scan the continuous data. We then performed the detection of new events with template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF]) and a single station approach [START_REF] Bell | Dynamic earthquake triggering response tracks evolving unrest at Sierra Negra volcano, Galápagos Islands[END_REF][START_REF] Poli | Creep and slip: Seismic precursors to the Nuugaatsiaq landslide (Greenland)[END_REF][START_REF] Van Der Elst | Enhanced remote earthquake triggering at fluid-injection sites in the midwestern United States[END_REF]. Because the use of a single station can reduce the detection sensitivity and increase the presence of unwanted signals, we combined visual inspection and detection with fake templates (e.g., waveforms flipped in time, [START_REF] Cabrera | A clustered preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake[END_REF], to define an optimal detection threshold (3-components average correlation coefficient greater than 0.5). The correlation threshold was set after an extensive visual inspection of the detected waveforms, to exclude false events from our final catalog. At this stage several detections can results from multiple templates over the time window of a template (12 s). To remove multiple detections, we keep the event with the largest correlation coefficient as the final detection.

The template matching identified 36,241 earthquakes, 300 times more than the initial catalog (Figure 4.1). Figure A3.S2 shows examples of detected waveforms.

We estimated the magnitude of each newly detected event by computing the mean Swave amplitude ratio between the template events and our detections over the three components. Using the template event's catalog magnitude as a reference, the detection magnitude can then be determined assuming a ratio of 10 corresponds to a variation of oneunit magnitude (e.g., [START_REF] Kato | Accelerated nucleation of the 2014 Iquique, Chile Mw 8.2 earthquake[END_REF]Sánchez-Reyez et al., 2021). The average mean error between the estimated magnitudes and the USGS ones for our template is 0.2.

Figure 4.2a shows the temporal variations of the seismicity rate. Unlike what can be observed from the initial catalog, in which the first event occurs on 29 August 2020, the swarm begins on 7 August 2020, with an acceleration leading to the 1,200 events/day recorded on August 29th (Figure 4.2a). The rate of seismicity then decreases with a log-like behavior similar to an Omori law [START_REF] Omori | On the aftershocks of earthquake[END_REF], as is often observed for aftershock sequences although no mainshock was observed at the beginning of the sequence. Indeed, the largest event (M = 5.9) occurred on 6 November 2020 (Figure A3.S1). After this event, the seismicity dropped rapidly (Figures 4.2 and A3.S1). We study in more detail the time evolution of the sequence by calculating the distribution of waiting time (dt) in between events (e.g., Duverger et al., 2018, Figure 4.2b). We fit the waiting time with a gamma distribution:

(4.1)
where A is a constant [START_REF] Hainzl | Estimating background activity based on interevent-time distribution[END_REF]. The fitting results provides , a value which differs from mainshock-aftershocks style of seismicity ( ). Our result implies a clustered seismicity, but little interaction in between events [START_REF] Duverger | Dynamics of microseismicity and its rela-tionship with the active structures in the western Corinth Rift (Greece)[END_REF][START_REF] Hainzl | Estimating background activity based on interevent-time distribution[END_REF]. We can thus rule out static and dynamic events interaction as a driver for this sequence (e.g., events triggered by mainshocks), as discussed in previous works [START_REF] Duverger | Dynamics of microseismicity and its rela-tionship with the active structures in the western Corinth Rift (Greece)[END_REF][START_REF] Hainzl | Estimating background activity based on interevent-time distribution[END_REF]. The seismicity must therefore be driven by an external forcing (e.g., [START_REF] Bourouis | Evidence for coupled seismic and aseismic fault slip during water injection in the geothermal site of Soultz (France), and implications for seismogenic transients[END_REF][START_REF] Duverger | Dynamics of microseismicity and its rela-tionship with the active structures in the western Corinth Rift (Greece)[END_REF][START_REF] Hainzl | Estimating background activity based on interevent-time distribution[END_REF][START_REF] Perfettini | Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan[END_REF] likely to be a magmatic process occurring in the ridge. Similar observations (no clear mainshock, slow decay of the seismicity) were observed for the 2014 2015 Deception Island volcanic swarm, in the western part of the Bransfield Basin [START_REF] Almendros | Volcanotectonic activity at deception island volcano following a seismic swarm in the Bransfield Rift (2014 2015)[END_REF].

To gain more insights about the possible origins of this swarm we also assessed the bvalue [START_REF] Gutenberg | Seismicity of the Earth[END_REF]. For regular earthquake sequence, b-value is usually around 1 [START_REF] Frohlich | Teleseismic b values; or, much ado about 1.0[END_REF], and deviation from this average can provide information about stress and/or physical properties of the rock volume [START_REF] Farrell | Earthquake swarm and b-value characterization of the Yellowstone volcano-tectonic system[END_REF][START_REF] Mogi | Magnitude frequency relation for elastic shocks accompanying fractures of various materials and some related problems in earthquakes[END_REF][START_REF] Schorlemmer | Variations in earthquake-size distribution across different stress regimes[END_REF]. We use the magnitude frequency distribution of Figure 4.2c to estimate a b-value of 1.6 ± 0.1, using a least squares method and considering a magnitude completeness of 2.9, estimated from maximum curvature method [START_REF] Wiemer | A software package to analyze seismicity: ZMAP[END_REF]. We further study how the errors in magnitude estimation affect the estimation of the b-value, by simulating 1,000 catalogs with randomly perturbed magnitude up to 0.2 magnitude units. With this approach, we obtained a final average b-value of 1.6 with standard deviation of 0.06.

The estimated b-value is remarkably higher respect with values characterizing regular earthquake sequences [START_REF] Farrell | Earthquake swarm and b-value characterization of the Yellowstone volcano-tectonic system[END_REF][START_REF] Frohlich | Teleseismic b values; or, much ado about 1.0[END_REF]. It is also remarkably different than previous studies in this area [START_REF] Olivet | A seismic swarm at the Bransfield Rift, Antarctica[END_REF], which found a b-value of 1.2. Our detection of events with much smaller magnitudes allowed us to better resolve the b-value for this swarm. Such a high b-value coincides with other values observed in volcanic areas (e.g., [START_REF] Farrell | Earthquake swarm and b-value characterization of the Yellowstone volcano-tectonic system[END_REF][START_REF] Roberts | Are volcanic seismic b-values high, and if so when?[END_REF][START_REF] Wilks | Seismicity associated with magmatism, faulting and hydrothermal circulation at Aluto Volcano, Main Ethiopian Rift[END_REF].

No clear spatial migration of the seismicity can be observed from the USGS earthquake locations (Figure A3.S3). However, teleseismic events are usually characterized by significant uncertainties, which could bias this kind of observation. While we could not locate the microseismic events detected, due to small events being visible only at JUBA (Figure 4.1), we attempted to discern any spatiotemporal patterns from S-P time analysis, using a single station approach. We first analyzed if any general migration was visible, by plotting the S-P time picked at station JUBA (Figure 4.1) when selecting templates (Figure 4.3a). No clear global spatial migration was observed but rather some time-limited migration episodes and a shift of the seismicity farther from JUBA starting mid-October 2020 (Figure 4.3a). Despite the scarce information about the velocity structure in the area, we can get some insights about the extent of the swarm from the S-P time. We assume a Vp=6 km/s and . Using the mentioned Vp and Vs the spatial difference for the average S-P time of 2.6s is 21 km.

We also evaluated the relative P-and S-wave travel times for all events detected by a single template (a family of events). For some of the biggest families we used the templates as reference waveforms, and estimated the delay with respect to the reference events, by cross correlating 1 s of signal around P and S arrivals. For this analysis, the original sampling rate (100 Hz) at JUBA station was used.

Figures 4.3b and 4.3c shows the results for the two largest families, for P and S delays with correlation coefficients larger than 0.7. The seismicity is clustered within a small area, similar to previous volcanic swarms observed in this region (Almendros et al., 2018, Figure 4.3). The time delays show short times of coherent delays associated with bursts of seismic events. A possible diffuse migration toward the station can be observed at the beginning of the sequence (Figure 4.3b). From the analysis of the delay time, we also observed migration of the seismicity farther from JUBA station (Figure 4.3c), in late October and early November 2020, also inferred from S-P time differences (Figure 4.3a).

Geodetic Observations

We use GNSS observations from 2017 to June 2021, from a set of stations located around the area (Figure 4.1). All observations were processed in a network array including several IGS regional stations, applying the differential approach strategy with the Bernese GNSS Software V5.2 [START_REF] Báez | The Chilean GNSS network: Current status and progress towards early warning applications[END_REF]. We stack all daily solutions to generate time series of deformation with respect to the Antarctic plate (Figure 4.4). More details about processing of GNSS data can be found in Appendix 3, together with the estimated velocities (Table A3.S1 and Figures A3.S4, S5). 4.4) is observed from August 28th to June 2021 (when our analysis ends) on both the north and east components of station UYBA. The displacement vector is orthogonal to the Bransfield rift axis and is consistent with previous observations from the interseismic period [START_REF] Taylor | Kinematics and segmentation of the South Shetland islands-Bransfield Basin system, northern Antarctic Peninsula[END_REF]. This relative geometry agrees with rift opening dynamics [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF][START_REF] Taylor | Kinematics and segmentation of the South Shetland islands-Bransfield Basin system, northern Antarctic Peninsula[END_REF]. The evolution of GNSS displacement at station UYBA during the swarm (Figures 4.4c and 4.4d), closely follows the log-like decrease in cumulative events (Figure 4.3) for the first part of the sequence, suggesting that both seismicity and surface displacements are driven by the same process (either slow tectonic deformation or volcanic activity). In addition, after the occurrence of a Mw 5.9 event on 6 November 2020, we further observe that the seismicity rate starts to decay along with GNSS velocities (Figure 4.3).

Discussion and Conclusions

Despite a scarce geophysical instrumentation, we were able to document and characterize the largest earthquake swarm ever recorded in the Bransfield Basin [START_REF] Almendros | Volcanotectonic activity at deception island volcano following a seismic swarm in the Bransfield Rift (2014 2015)[END_REF][START_REF] Dziak | Tectonomagmatic activity and ice dynamics in the Bransfield Strait back-arc basin, Antarctica[END_REF]. We detected a long-lived swarm of more than 36,000 earthquakes that began on 7 August 2020, and accelerated on 29 August 2020, challenging previous findings from [START_REF] Olivet | A seismic swarm at the Bransfield Rift, Antarctica[END_REF]. The seismicity rate decreases following a log-like behavior, although no mainshock was observed at the beginning of the swarm (Figure 4.2). The detailed analysis of recurrence time (Figure 4.2b) allows us to rule out the mainshock-aftershock mechanism as the mechanism driving the occurrence of the large number of recorded events, and observed GNSS deformation (Figure 4.4). The cumulative number of earthquakes closely follows the deformation observed at the GNSS station on King George Island (Figures 4.1 and 4.4). Despite the noise in GNSS time series (Figures A3.S6 S13), probable from response to ice load changes, the orientation of the deformation is orthogonal to the ridge (Figure 4.1) and suggests that the ridge spreading is responsible for the 8 cm displacement (Figures 4.1 and 4.4). The deformation on King George Island could occur either in response to a dike intrusion [START_REF] Heimisson | Physically consistent modeling of Dike-induced deformation and seismicity: Application to the 2014 Bárðarbunga Dike, Iceland[END_REF] or could be associated with a large slow slip event with extensional geometry.

As previously mentioned, the increase in seismicity rate was not preceded by any clear mainshock (Figure A3.S1) and shows a very slow decay in time (Figures 4.2a and 4.2b). Both observations suggest a swarm-like sequence [START_REF] Mogi | Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (third paper)[END_REF], with limited interaction in between events (Figure 4.2b). This swarm is thus likely to be driven by external forcing, also responsible for the deformation observed in the GNSS data (Figure 4.4). The estimated b-value is 1.6, a significantly larger value than previous estimates based on a smaller seismicity catalog [START_REF] Olivet | A seismic swarm at the Bransfield Rift, Antarctica[END_REF]. We suggest that this difference arises mainly from our improved detection of events with template matching (36,241 events), which permits to outpace the number of events detected from a visual inspection approach, made by Olivet et al. (2021) (3,186 events). This large b-value can result from stress heterogeneity, significant thermal gradient and/or presence of magmatic fluids which has been observed in volcanic areas [START_REF] Farrell | Earthquake swarm and b-value characterization of the Yellowstone volcano-tectonic system[END_REF]. Beyond reporting a different b-value and more events than [START_REF] Olivet | A seismic swarm at the Bransfield Rift, Antarctica[END_REF] we also better characterize their time evolution (Figures 4.2a and 4.2b) and compare it with GNSS deformation recorded during the seismic activity (Figure 4.4) which allows us to better support our conclusions.

The spatiotemporal analysis of the seismicity (Figure 4.3) resolved with a single station approach, reveals that the swarm nucleated on a small region, with rapid migrations and quick activation of small seismic bursts (Figure 4.3). No clear large-scale spatially coherent migration of the seismicity, as the one related to dyke injection discussed by [START_REF] Roman | The origin of volcano-tectonic earthquake swarms[END_REF], is observed. This behavior can reflect a strong stress heterogeneity [START_REF] Farrell | Earthquake swarm and b-value characterization of the Yellowstone volcano-tectonic system[END_REF]) also suggested by the high b-value (Figure 4.2). We further observe that seismicity mostly lies in a region of high positive magnetic anomaly [START_REF] Almendros | BRAVOSEIS: Geophysical investigation of rifting and volcanism in the Bransfield strait, Antarctica[END_REF], interpreted as a shallow magmatic body. In addition, most of the large events show strike-slip mechanisms, similar to the swarm model of [START_REF] Hill | A model for earthquake swarms[END_REF] and [START_REF] Roman | The origin of volcano-tectonic earthquake swarms[END_REF]. The relatively small number of extensional earthquakes (Figure 4.1) suggests that deformation related to rifting [START_REF] Reiss | The impact of complex volcanic plumbing on the nature of seismicity in the developing magmatic natron rift, Tanzania[END_REF] is limited during this 2020 2021 long episode.

Taken together, our observations suggest a volcanic origin for the 2020 2021 Bransfield Ridge deformation episode [START_REF] Bergman | Earthquake swarms on the mid-Atlantic Ridge: Products of magmatism or extensional tectonics[END_REF]. We cannot however further discuss the detailed processes occurring during this swarm given the limited amount of data available in the region. We propose that the driving mechanism of the deformation can be either hydrothermal fluids [START_REF] Reiss | The impact of complex volcanic plumbing on the nature of seismicity in the developing magmatic natron rift, Tanzania[END_REF] or magma flows (e.g., shallow dike propagation [START_REF] Heimisson | Physically consistent modeling of Dike-induced deformation and seismicity: Application to the 2014 Bárðarbunga Dike, Iceland[END_REF] at crustal level favored by the presence of conjugate faults [START_REF] Hill | A model for earthquake swarms[END_REF]. The seismicity can also result from local increment of pore fluid pressure [START_REF] Sibson | Fluid involvement in normal faulting[END_REF]. This long-lasting volcanic activity is responsible for the significant deformation inferred from GNSS observations, while seismicity is a by-product of the magmatic activity [START_REF] Heimisson | Physically consistent modeling of Dike-induced deformation and seismicity: Application to the 2014 Bárðarbunga Dike, Iceland[END_REF], mainly occurring in limited areas with brittle characteristics and accumulation of stress [START_REF] Hill | A model for earthquake swarms[END_REF]. Our study also illustrates that, beyond the continuous 7 mm/yr extension [START_REF] Taylor | Kinematics and segmentation of the South Shetland islands-Bransfield Basin system, northern Antarctic Peninsula[END_REF] between the Antarctica Plate and the South Shetland microplate, rapid deformation episodes occurring at the ridge axial volcanic structures plays a main role in modulating the long-term extension. Finally, our study highlights the main role of magmatic structures in favoring the rifting process instead of tectonic deformation occurring in rifting bounding border faults [START_REF] Buck | Consequences of asthenospheric variability on continental rifting[END_REF][START_REF] Reiss | The impact of complex volcanic plumbing on the nature of seismicity in the developing magmatic natron rift, Tanzania[END_REF].

How faulting processes lead to a large earthquake is a fundamental question in seismology. To better constrain this pre-seismic stage, we create a dense seismic catalog via template matching to analyse the precursory phase of the Mw 6.3 earthquake that occurred in central Italy in 2009. We estimate several physical parameters in time, such as the coefficient of variation, the seismic moment release, the effective stress drop, and analyse spatio-temporal patterns to study the evolution of the sequence and the earthquake interactions. We observe that the precursory phase experiences multiple accelerations of the seismicity rate that we divide into two main sequences with different signatures and features: the first part exhibits weak earthquake interactions, quasi-continuous moment release, slow spatial migration patterns, and a lower effective stress drop, pointing to aseismic processes. The second sequence exhibits strong temporal clustering, fast seismicity expansion, and larger effective stress drop typical of a stress transfer process. We interpret the differences in seismicity behaviors between the two sequences as distinct physical mechanisms that are controlled by different physical properties of the fault system. We conclude that the earthquake is preceded by a complex preparation, made up of different physical processes over different time scales on faults with different physical properties.

Plain Language Summary

In this work we study the seismicity before the Mw 6.3 2009 L'Aquila earthquake. We first catalog nearly 5000 events from the continuous seismic record. We then analyse the spatiotemporal evolution of this sequence with several physical parameters. We observe that the sequence is divided into two main sequences. Our results indicate that several different physical mechanisms (e.g., aseismic deformation, stress transfer due to earthquake interactions) and potential heterogeneities in the fault system (e.g., distance between seismic regions) controlled how the earthquake sequence played out. Our observations show a complex spatiotemporal evolution during the precursory phase and challenge classic fault models that explain earthquake initiation as a process along a homogenous planar fault.

Introduction

The characterization of the physical processes occurring before major earthquakes is an essential step towards understanding when and where future earthquakes will nucleate. So far, physical models have been proposed to explain the processes that lead to large seismic events, including cascade, pre-slip, and progressive or migratory localization (Ellsworth & Beroza, 1995;[START_REF] Kato | The generation of large earthquakes[END_REF][START_REF] Mclaskey | Earthquake initiation from laboratory observations and implications for foreshocks[END_REF]. Which one of these mechanisms best represents the physics of the precursory phase of earthquakes is still under debate.

One of the most powerful observational tools to study the physical processes of earthquake nucleation are foreshocks, small earthquakes that precede some large mainshocks [START_REF] Bouchon | The long precursory phase of most large interplate earthquakes[END_REF][START_REF] Dodge | Detailed observations of California foreshock sequences: Implications for the earthquake initiation process[END_REF]. Foreshocks were first observed more than a century ago [START_REF] Omori | On the Fore-shocks of Earthquakes[END_REF]. Since then, many laboratory studies have focused on the precursory moment release [START_REF] Acosta | Can precursory moment release scale with earthquake magnitude? A view from the laboratory[END_REF], aseismic slip and stress changes [START_REF] Mclaskey | Foreshocks during the nucleation of stick slip instability[END_REF][START_REF] Mclaskey | Preslip and cascade processes initiating laboratory stick slip[END_REF], and other characteristics of the foreshocks during the initiation of laboratory earthquakes (McLaskey, 2019 and references therein). In addition, direct seismological observations in different seismotectonic settings such as strike-slip faults (e.g., [START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF][START_REF] Chen | California foreshock sequences suggest aseismic triggering process[END_REF][START_REF] Dodge | Detailed observations of California foreshock sequences: Implications for the earthquake initiation process[END_REF][START_REF] Dodge | Foreshock sequence of the 1992 Landers, California, earthquake and its implications for earthquake nucleation[END_REF][START_REF] Durand | A Two Scale Preparation Phase Preceded an M w 5.8 Earthquake in the Sea of Marmara Offshore Istanbul, Turkey[END_REF][START_REF] Ellsworth | Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks[END_REF][START_REF] Shelly | A high resolution seismic catalog for the initial 2019 Ridgecrest earthquake sequence: Foreshocks, aftershocks, and faulting complexity[END_REF][START_REF] Tape | Earthquake nucleation and fault slip complexity in the lower crust of central Alaska[END_REF][START_REF] Yoon | Foreshocks and mainshock nucleation of the 1999 M w 7.1 Hector Mine, California, Earthquake[END_REF], subduction zones [START_REF] Bouchon | The long precursory phase of most large interplate earthquakes[END_REF][START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake[END_REF][START_REF] Ruiz | Nucleation phase and dynamic inversion of the Mw 6.9 Valparaíso 2017 earthquake in Central Chile[END_REF], and extensional regimes [START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF] have been carried out to assess which model best explains the occurrence of foreshocks and the physical processes occurring during the precursory phase of large earthquakes. More recently, some studies have taken advantage of high-resolution detection methods, such as template matching and/or machine learning (e.g., [START_REF] Durand | A Two Scale Preparation Phase Preceded an M w 5.8 Earthquake in the Sea of Marmara Offshore Istanbul, Turkey[END_REF][START_REF] Gardonio | Seismic activity preceding the 2011 M w 9.0 Tohoku earthquake, Japan, analyzed with multidimensional template matching[END_REF][START_REF] Ross | Searching for hidden earthquakes in Southern California[END_REF][START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Shelly | A high resolution seismic catalog for the initial 2019 Ridgecrest earthquake sequence: Foreshocks, aftershocks, and faulting complexity[END_REF][START_REF] Yoon | Foreshocks and mainshock nucleation of the 1999 M w 7.1 Hector Mine, California, Earthquake[END_REF], and the availability of better field data (i.e., more stations near faults; e.g., [START_REF] Meng | Immediate foreshocks indicating cascading rupture developments for 527 M 0.9 to 5.4 Ridgecrest earthquakes[END_REF]; [START_REF] Savage | Possible precursory signals in damage zone foreshocks[END_REF]; [START_REF] Simon | Possible Precursory Slow Slip to two ML 3 Mainevents of the Diemtigen Microearthquake Sequence[END_REF]; [START_REF] Tape | Earthquake nucleation and fault slip complexity in the lower crust of central Alaska[END_REF]) to study foreshocks. These studies reveal an increased spatiotemporal complexity (i.e., fault interactions, volumetric processes, heterogeneous fault properties) of the processes taking place before large earthquakes. This complexity, mainly revealed by foreshocks patterns, is hard to reconcile with a single physical explanation of the precursory phase (cascade or preslip). In addition, the observed foreshocks patterns challenge the actual laboratory scale and theoretical models, which treat earthquake initiation as a process along a homogenous planar fault [START_REF] Dieterich | Earthquake nucleation on faults with rate-and state-dependent strength[END_REF][START_REF] Liu | Aseismic slip transients emerge spontaneously in three dimensional rate and state modeling of subduction earthquake sequences[END_REF][START_REF] Marone | The effect of loading rate on static friction and the rate of fault healing during the earthquake cycle[END_REF][START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF] or a combination of several planar fault segments [START_REF] Shimizu | Construction of fault geometry by finite-fault inversion of teleseismic data[END_REF], although some cases with non-planar fault geometry exist (e.g., [START_REF] Dutta | Simultaneous Bayesian Estimation of Non Planar Fault Geometry and Spatially Variable Slip[END_REF][START_REF] Zhang | Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics[END_REF].

To gain insight into the ongoing physical processes occurring near the nucleation region before a large earthquake, we study the Mw 6.1 2009 L'Aquila earthquake and its foreshock sequence. This event, which struck central Italy on 6 April 2009 (01:32 UTC) and caused damage and fatalities, was preceded by more than 500 small (M > 0.5) earthquakes since the beginning of January (Chiaraluce, Valoroso, et al., 2011). Based on the locations of the events, Chiaraluce, Valoroso, et al. (2011) reported that the sequence of foreshocks took place in two different faults: (a) the main fault, where the mainshock (Figure 5.1) occurs on 6 April 2009, that hosts most of the seismicity occurring from the beginning of January until 30 March and (b) an antithetic fault that is activated by an Mw 3.9 foreshock (hereafter F1, Figure 5.1) on 30 March 2009. On 5 April 2009 (5 h before the mainshock), the seismicity migrates back to the main fault after the occurrence of another Mw 3.9 foreshock (hereafter F2, Figure 5.1, Chiaraluce, Valoroso, et al., 2011). The co-seismic rupture took place in the Paganica fault [START_REF] Cheloni | Coseismic and initial post-seismic slip of the 2009 M w 6.3 L'Aquila earthquake, Italy, from GPS measurements[END_REF][START_REF] Falcucci | The Paganica fault and surface coseismic ruptures caused by the 6 April 2009 earthquake (L'Aquila, central Italy)[END_REF]; labeled PaF in Figure 5.1), generating exposed ground deformation [START_REF] Boncio | Coseismic ground deformation of the 6 April 2009 L'Aquila earthquake (central Italy, Mw6. 3)[END_REF][START_REF] Falcucci | The Paganica fault and surface coseismic ruptures caused by the 6 April 2009 earthquake (L'Aquila, central Italy)[END_REF] and maximum surface displacements of 8.1 and 16.5 cm in the vertical and horizontal directions, respectively [START_REF] Cheloni | Coseismic and initial post-seismic slip of the 2009 M w 6.3 L'Aquila earthquake, Italy, from GPS measurements[END_REF]. Joint inversion using GPS, strong motion, and Synthetic Aperture Radar (SAR) data indicate that the maximum slip on the fault is about 1.4 m [START_REF] Cirella | Complexity of the rupture process during the 2009 L'Aquila, Italy, earthquake[END_REF]. According to different rupture models (e.g., [START_REF] Cheloni | Coseismic and initial post-seismic slip of the 2009 M w 6.3 L'Aquila earthquake, Italy, from GPS measurements[END_REF][START_REF] Cirella | Rupture history of the 2009 L'Aquila (Italy) earthquake from non linear joint inversion of strong motion and GPS data[END_REF][START_REF] Cirella | Complexity of the rupture process during the 2009 L'Aquila, Italy, earthquake[END_REF][START_REF] Scognamiglio | Fast determination of moment tensors and rupture history: What has been learned from the 6 April 2009 L'Aquila earthquake sequence[END_REF], the slip was concentrated in two main asperities: a small patch updip from the hypocenter, and a second, larger asperity located to the southeast along strike. In this context, the foreshocks were located at the base of the activated fault plane in a region where almost no slip occurred during the mainshock rupture [START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF]. (Chiaraluce, Valoroso, et al., 2011;[START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF][START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF][START_REF] Vuan | Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake: Development of a Scalable Code Based on the Cross Correlation of Template Earthquakes Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake[END_REF] by estimating quantitative parameters of the spatiotemporal evolution of the foreshocks sequence. We focus on an area of 10 × 10 km surrounding the epicenter (Figure 5.1). We [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF] to scan 6 months of data before the main shock. We use a frequency band between 5 and 30 Hz. The inclusion of high frequencies (>20 Hz) compared to previous studies [START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF][START_REF] Vuan | Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake: Development of a Scalable Code Based on the Cross Correlation of Template Earthquakes Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake[END_REF] permitted us to detect more small events ( M < 1.0), which are best captured at high frequency. Our final catalog with coverage from 6 October 2008 6 April 2009 contains 4978 events; the first event occurring on 3 January 2009. No seismicity was detected from 6 October 2008 2 January 2009.

Using this new catalog, we analyze the seismic sequence of foreshocks by tracking the time evolution of temporal clustering (earthquake interactions), seismic moment release, and effective stress drop. We also study the spatio-temporal evolution of the events to better discuss the physical mechanisms that control the foreshock sequence, ultimately leading to the mainshock. Valoroso, et al. (2011); our 267 template events are drawn from this earthquake catalog and are shown in the zoomed map. The broadband stations we analyzed are shown by the red triangles. Black and cyan thin lines respectively represent traces of the active mapped faults and co-seismic surface ruptures, respectively [START_REF] Boncio | Coseismic ground deformation of the 6 April 2009 L'Aquila earthquake (central Italy, Mw6. 3)[END_REF]. PaF indicates the location of the Paganica Fault. Upper-right zoom: 267 events used as templates to scan continuous data. Their color and size are respectively coded to depth and magnitude. Beachballs (compressional quadrants in colors) represent source mechanisms (reported by INGV) for the mainshock (MS 6 April) and the two foreshocks Mw 3.9 1 week (F1 30 March) and 5 h (F2 5 April) before it. All of them correspond to normal (extensional) mechanisms.

Extending the Seismic Catalog

We apply template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF] to continuous seismic data collected by the Istituto Nazionale di Geofisica e Vulcanologia (INGV) from 6 October 2008 6 April 2009 (6 months). We use 10 broadband three-component stations (red triangles in Figure 5.1) from the Italian Seismic Network (INGV Seismological Data Center., 2006) and the Mediterranean Very Broadband Seismographic Network (MedNet Project Partner Institutions, 1990). Data were continuously recorded at a sampling rate of 100 Hz. Before using the data to study earthquakes, we performed a visual inspection of the spectrograms (Figure A4.S1) to find the frequency range that is the least affected by the strong anthropogenic noise in the Apennines [START_REF] Poli | Seasonal and Coseismic Velocity Variation in the Region of L'Aquila From Single Station Measurements and Implications for Crustal Rheology[END_REF]. We choose to filter the continuous data from 5 to 30 Hz after this analysis. The dataset was then organized into 24-hr continuous files with all gaps filled with zeros.

We consider 512 foreshocks reported by Chiaraluce, Valoroso, et al. (2011) as potential templates, which have relative horizontal and vertical location errors about 40 and 80 m, respectively (Chiaraluce, Valoroso, et al., 2011). We identify the highest-quality events by estimating the signal-to-noise ratio (SNR) of each event as the ratio between the RMS velocity during the first 3 s of the P and S waves (for vertical and horizontal components, respectively), and the RMS velocity during a 3 s of noise before the P and S wave arrival times (e.g., Cabrera et al., 2021;Frank et al., 2017). We estimate the arrivals using 1D velocity models for P [START_REF] Chiaraluce | The 2009 Aquila (central Italy) seismic sequence[END_REF] and S waves [START_REF] Herrmann | Regional moment tensors of the 2009 earthquake sequence[END_REF]. A signal is retained as a final template if it has a SNR for at least 12 components. We finally retained 267 template event waveforms (inset in Figure 5.1), defined as the 3.5 s time windows that start 0.5 s before the P-and S-wave arrivals at each station for the vertical and horizontal components, respectively, and filtered in an identical manner to the continuous data (bandpassed between 5 and 30 Hz). The template waveforms are then correlated against a sliding window of continuous data using a GPU-architecture and the Fast Matched Filter algorithm [START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF] to obtain daily correlation functions. We search sample-by-sample considering a detection threshold of 12 times the median absolute deviation (MAD) of the correlation function averaged over all stations and channels to detect events significantly similar to the template. We defined this detection threshold to minimize false detections by first scanning the continuous data using the templates flipped in time (see an example in Figure A4.S2). With this approach, the data are scanned using non-physical and acausal templates unlikely to detect anything, but with the same frequency content as the original templates. We test the number of detections using NxMAD with N in the range 9 12 (see Figure A4.S3), and we decided to use N=12 as this threshold provides only one false detection during the whole period of time (6 October 20086 April 2009). To remove double detections over the same time window, we merge consecutive detections with differential times less than 4 s; we keep the detection with the highest network-averaged correlation coefficient as the final detection.

We estimate the magnitude of each new event by computing the mean P-and S-wave amplitude ratio between the template event and the detection over the components with a detected event is determined, assuming that a ratio of 10 of the amplitude ratio corresponds to a variation of one unit of magnitude (e.g., Cabrera et al., 2021;Frank et al., 2017;[START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF].

We further attempt to relocate the newly detected seismicity with respect to the templates. For this scope, we use pair-wise cross-correlation (CC) between each template and its detections to measure differential delay times. For each event pair, we use waveform windows of 2 s starting 1 s before the phase arrival for both P-and S-waves. We then relocate each family of detections (a template and its detections) with GrowClust [START_REF] Trugman | GrowClust: A hierarchical clustering algorithm for relative earthquake relocation, with application to the Spanish Springs and Sheldon, Nevada, earthquake sequences[END_REF]. An event pair is only used if its cross-correlation coefficient (

) is 0.6 with a maximum source-receiver distance (

) of 80 km. We also considered a maximum rootmean-square differential time residual ( ) 0.2 for a proposed cluster merger to be allowed during relocation (see [START_REF] Trugman | GrowClust: A hierarchical clustering algorithm for relative earthquake relocation, with application to the Spanish Springs and Sheldon, Nevada, earthquake sequences[END_REF] for more details). This procedure resulted in 722 relocated events or 17% of the original catalog (Figure A4.S4). Although low, this percentage is not surprising given the configuration of the network. For example, [START_REF] Ross | Searching for hidden earthquakes in Southern California[END_REF] relocated 38.7% of events using a denser array of stations in California, and Simon et al. ( 2021) relocated 11.6% of their catalog in Switzerland, in both cases after template matching. This data reduction is because double difference relocations rely on high quality correlations at a single station, while template matching leverages an average correlation across the entire network to identify events that would otherwise go unnoticed. This means that some events that could be detected by template matching can have relatively low correlation coefficients that are not necessarily suitable for relocation. Although it is possible to increase the number of relocated events by relaxing some parameters like the minimum crosscorrelation coefficient ( ), we preferred to use parameter values similar to previous works (e.g., [START_REF] Ross | Searching for hidden earthquakes in Southern California[END_REF][START_REF] Trugman | GrowClust: A hierarchical clustering algorithm for relative earthquake relocation, with application to the Spanish Springs and Sheldon, Nevada, earthquake sequences[END_REF] to prevent degradation of the relocation quality.

As only a small number of events can be relocated with the approach described above, we consider the newly detected to have occurred at the same hypocenter (determined by Chiaraluce, Valoroso, et al. [2011]) as the template. We estimated the distance between the initial location of the detections and the relocated position of new detections. On average, horizontal and vertical distances between templates and new detections are in the order of 83 and 66 m, respectively (see Figures A4.S5-S6). These values are similar to other studies ( 100 200 m, [START_REF] Ross | Searching for hidden earthquakes in Southern California[END_REF]; [START_REF] Simon | Possible Precursory Slow Slip to two ML 3 Mainevents of the Diemtigen Microearthquake Sequence[END_REF]). 5.2). We estimate the magnitude of completeness (Mc) of our catalog, using the Lillefors test implemented by [START_REF] Herrmann | Inconsistencies and lurking pitfalls in the magnitude frequency distribution of high resolution earthquake catalogs[END_REF], which in general provides conservative values of the Mc (see examples in [START_REF] Herrmann | Inconsistencies and lurking pitfalls in the magnitude frequency distribution of high resolution earthquake catalogs[END_REF]) and allows us to ensure the s conservative in a statistical sense [START_REF] Clauset | Power-law distributions in empirical data[END_REF]. Therefore, we prefer Mc = 0.9, a more conservative value for the magnitude of completeness to show the stability of our later analysis (see Text A4.S1 for more details). Our catalog presents a decrease in the magnitude of completeness compared to the Mc = 1.8 derived for [START_REF] Vuan | Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake: Development of a Scalable Code Based on the Cross Correlation of Template Earthquakes Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake[END_REF] and the template's catalog, considering the same estimation described above.

This new catalog is the largest catalog for this precursory sequence to date [START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF] reported 3571 events and [START_REF] Vuan | Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake: Development of a Scalable Code Based on the Cross Correlation of Template Earthquakes Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake[END_REF] extended using one station up to 3786 events), and is created using many constraints to ensure high quality of the detections, such as the selection of the templates based on the SNR criteria for P and S waves, a higher frequency band, an N-value threshold selection determined from non-physical acausal templates, a relocation to measure the distance between templates and detections, and a longer period of time scanned.

Figure 5.3a shows that the seismicity starts on January 3 and lasts until April 6, when the mainshock occurs on the main fault. No seismicity is detected in the period between 6 October 2008 2 January 2009, so we consider the seismicity starting on the 3 of January as foreshocks of the 9 of April Mw 6.1 earthquake (Chiaraluce, Valoroso, et al., 2011;[START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF][START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF][START_REF] Vuan | Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake: Development of a Scalable Code Based on the Cross Correlation of Template Earthquakes Improving the Detection of Low Magnitude Seismicity Preceding the Mw 6.3 Earthquake[END_REF]. We observe that the rate of events strongly increases after a Mw 3.9 foreshock on 30 March (F1), which activates an antithetic fault (Chiaraluce, Valoroso, et al., 2011;[START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF]. This activation of the seismicity on the antithetic fault is evidenced in Figure 5.3b, which shows a summary of the vertical normalized waveforms for the AQU station (the closest one to the mainshock epicenter, see Figure 5.1) aligned on the P-wave arrival. A significant difference in the S-wave arrivals is observed after F1, at the same time as the spatial evolution reported by Chiaraluce, Valoroso, et al. (2011) and [START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF] (see also Movie S1). As a first analysis, we split the seismicity before and after F1 on 30 March into two sequences (hereafter S1 and S2, respectively). We observe that the respective cumulative event counts (Figures 5.3a We track the spatio-temporal evolution of several parameters that describe the style of the seismicity and provide hints about the physical processes active during the foreshock sequence. The mainshock is excluded from this analysis.

Tracking Foreshocks

We study and discuss the spatio-temporal evolution of the seismicity by tracking the time development of several parameters that characterize the style of seismicity. The parameters are estimated using moving windows of 100-events with a 99-events overlap, that is, the first estimate considers the first 100-events, and each subsequent estimate is just shifted by one event in time. This approach allows us to characterize the general evolution of the sequence rather than just focusing on specific time periods. At this point, it is necessary to consider the potential effects of the magnitude of completeness and the number of events used in each time window. To that scope, we performed tests considering only events with magnitudes larger than the magnitude of completeness and assessed the effect of varying numbers of events for windows-lengths and overlaps, to evaluate the stability of the results (see Figure A4.S7-S10). A jack-knife process was also carried out, removing 20% of the catalog in 100 realizations, to assess the uncertainties for each parameter (Figure A4.S7-S10).

Considering the robustness of the tests mentioned above, we present here the results for the entire catalog (Figures 5.4,5.5 and 5.6).

Temporal Clustering

Temporal clustering of seismicity, that is, how past events affect the occurrence of the future ones, is a key feature of seismicity and is thought to be principally related to static or dynamic stress transfer [START_REF] Freed | Earthquake triggering by static, dynamic, and postseismic stress transfer[END_REF]. Therefore, the study of temporal clustering probes the degree to which earthquake interactions drive the propagation of seismic sequences over external forcing or other physical processes [START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas[END_REF].

To quantify the level of time clustering of the seismicity, we estimate the coefficient of variation (COV) of the interevent times times ( ) (Figure 5.3a) as , where

is the standard deviation and is the average of the interevent times within the window [START_REF] Kagan | Long-term earthquake clustering[END_REF]. The COV is 0 for a periodic occurrence of seismicity, 1 for completely random Poisson occurrence, and larger than 1 for temporally clustered earthquakes; put plainly, the larger the COV is, the stronger the time clustering is [START_REF] Kagan | Long-term earthquake clustering[END_REF][START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas[END_REF]. Figure 5.4a shows the temporal evolution of the COV. During S1 we see slow oscillations of the COV, with generally low values (ranging from 1 to 2.5). We observe that decreases of the COV are often associated with accelerations of seismicity (Figure 5.4a). The lowest values (COV 1) for S1 are observed during an increase of the seismicity rate starting on 21 January (cyan dots in Figure 5.4a) and on 15 February. This observation suggests that the increment of seismicity rate is not due to interevent stress triggering (e.g., seismicity is not driven by a mainshock), and an external mechanism likely acts to increase the number of events [START_REF] Beaucé | Systematic detection of clustered seismicity beneath the Southwestern Alps[END_REF]. On the other hand, periods with increased seismicity rates within S2 exhibit episodic, strong temporal clustering followed by random seismicity akin to mainshockaftershock sequences [START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas[END_REF]. Summarizing, the COV evolution reflects changes in the seismicity style as a function of time, especially when moving from S1 to S2.

Evolution of Seismic Moment Release

The time evolution of the seismic moment ( ) release reflects the behavior of different types of seismic sequences and offers insights into the processes on activated faults [START_REF] Vidale | A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers[END_REF]. While a stable and gradual moment release by many earthquakes without a dominant large magnitude event is observed for swarm-type sequences [START_REF] Vidale | A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers[END_REF], most of the moment is released at once during mainshock-aftershocks sequences [START_REF] Mogi | Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (third paper)[END_REF].

To analyze the time evolution of the seismic moment release, we estimate for each 100event sliding window the ratio of the maximum seismic moment to the total seismic moment of the window (MaxMo/ , see Figure 5.4b). Values close to 1 indicate episodes when the largest event represents most of the Mo released, whereas values close to 0 are observed for windows without a dominant event in terms of moment release. We estimate the moment for each event using the definition of Mw of Kanamori (1977).

Our analysis shows a smooth evolution of the moment ratio during S1 (Figure 5.4b), with values ranging from 0.1 to 0.6. This implies that the seismic moment is released nearly uniformly within the window, rather than impulsively by some dominant event. From the beginning of the sequence until 15 February, the moment ratio is generally low, despite some periods with increased seismicity rates when larger events occur (Figure 5.4b). From 15 February until F1 we observe an increase of the moment ratio associated with larger events M2.7 and M 2.9 (17 February and 11 March, respectively; see Figures 5.2a and 5.4b). During S2, the evolution of moment release is more episodic, with large rapid releases of the moment, mainly associated with the largest foreshocks (e.g., F1 and F2). A comparison of Figures 5.4a and 5.4b shows that the COV and seismic moment release have similar patterns, especially during S2, with peak values associated with the largest foreshocks followed by a rapid decrease of COV and moment release. We note that for the moment ratio, the size of the selected window has an effect on the observed level of smoothing; the maximum magnitude (numerator) does not change with different window sizes while the total moment (denominator) does. We observe that smaller windows increase the relative values of small mainshock-aftershocks sequences (e.g., Figure A4.S9), which makes the evolution of the moment ratio more episodic and less smooth; this effect is diminished for larger windows (e.g., Figure A4.S10).

Spatial Evolution

To assess the spatial evolution of the seismicity as a function of time we track the alongstrike and along-dip position of the events. We first project all the seismicity onto the main fault plane for S1 (strike N133°E and Dip 50° according to [START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF]) and an orthogonal fault plane for S2 (when the antithetic fault is active), obtaining along-strike and along-dip distances measured from the position of the mainshock. We thus track the position of the seismicity centroid (estimated as the average of all event locations) for each time window and extract the along-strike and along-dip coordinates, which are plotted in Figure 5.4c. In addition, the 3D evolution is presented in the supplementary material (see Movie S1).

As previously observed by [START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF], the seismicity starts on the north-west segment of the fault (Figure 5.4c). We observe a slow but significant movement of the centroid along strike and dip until 13 February, mainly during increases of the seismicity rate. Observations along with the 3D spatial evolution presented in Movie S1 suggest seismicity reruptures the same fault segment. After the initial activity in the NW segment of the fault, a prominent along-strike migration occurs toward the south-east, observed in Movie S1 and also tracked by the large along-strike variation of the centroid (Figure 5.4c). This migration begins around 13 February, accelerates on 17 February covering 1.2 km in less than 24 hr, and fades on 18 February. At this point, the along-dip and along-strike position of the seismicity stabilizes until the end of S1. An example of the seismicity migratory behavior within S1 is shown in Figure 5.5a (group of cyan dots in Figure 5.4). On 21 January 2009, velocities of kilometers per day are required to reproduce the seismicity front (we discuss this in more detail in Section 5.4).

During the second sequence S2 the antithetic fault is activated after foreshock F1, and Figure 5.4c shows that the centroid of the seismicity is confined between 1.5 and 2.5 km from the hypocenter both along dip and strike. In addition, the 3D evolution of the seismicity during S2 (Movie S1) does not exhibit slow migrations as observed in S1, but rather a rapid spread of the seismicity across the antithetic fault. This latter behavior is exemplified by Figures 5.5b and 5.5c, with longer distances (kilometers) rapidly covered in seconds by seismicity after the occurrence of major events (e.g., F1, F2, and another event magnitude 3.2 in the middle of the sequence indicated by a yellow star). These distances are longer than the expected magnitudedependent rupture lengths of the major earthquakes which are about 400 800 m for a M3.9 event (e.g., [START_REF] Dascher Cousineau | What controls variations in aftershock productivity?[END_REF]Udias et al., 2014). The strong clustering that we observe in the COV (Figure 5.4a) at the moment of large foreshocks in S2 (e.g., F2) together with the rapid large scale spreading of the seismicity (Vel. = 3 km/s), suggests that stress triggering is the main mechanism driving the seismicity during the bursts in S2 [START_REF] Freed | Earthquake triggering by static, dynamic, and postseismic stress transfer[END_REF].

Event Index Evolution of the Sequence

The migration style of earthquake locations is generally considered to be an important characteristic to distinguish earthquake swarms from aftershock sequences [START_REF] Fischer | The Growth of Earthquake Clusters[END_REF]. While swarms typically show hypocenter migration that depends on the mechanism driving the swarm (e.g., pore-pressure diffusion, [START_REF] Shapiro | Estimating the crust permeability from fluidinjection-induced seismic emission at the KTB site[END_REF]; hydraulic fracture growth, [START_REF] Dahm | Bidirectional and unidirectional fracture growth during hydrofracturing: role of driving stress gradients[END_REF]; or slow slip, [START_REF] Schwartz | Slow slip events and seismic tremor at circumpacific subduction zones[END_REF]), aftershocks usually occur immediately across the entire fault plane and along the edges of the mainshock rupture due to stress transfer (e.g., [START_REF] Freed | Earthquake triggering by static, dynamic, and postseismic stress transfer[END_REF] although some slower migrations can be linked to afterslip (e.g., [START_REF] Perfettini | A model of aftershock migration driven by afterslip[END_REF]. Usually, the way to analyze such migration patterns is in the distance-time domain, where the independent variable is typically the time. However, as shown by [START_REF] Fischer | The Growth of Earthquake Clusters[END_REF] a complimentary analysis tool is to use the event order (e.g., event index) as the independent variable, which is also termed natural time [START_REF] Rundle | Natural time, nowcasting and the physics of earthquakes: estimation of seismic risk to global megacities[END_REF]. While the use of time as an independent variable allows resolving whether time controls the seismogenic process, the use of the event index indicates whether the seismogenic process itself controls the seismicity, that is, whether each rupture facilitates the nucleation of the next rupture [START_REF] Fischer | The Growth of Earthquake Clusters[END_REF]. [START_REF] Fischer | The Growth of Earthquake Clusters[END_REF] showed that an indexplot migration is linear or square-root for either external processes such as pore-pressure diffusion, hydraulic fracture, and slow slip or in case of an internal process, such as the creation of pore-space during ruptures. In contrast to the random (in space) occurrence of aftershock hypocenters along the mainshock fault plane.

Figure 5.6 shows the comparison between the time (Figure 5.6a) and the event-index (Figure 5.6b) plots for the along-strike position of the seismicity, centered in the mainshock. We observe a coherent spreading of seismicity during S1 (e.g., Figure 5.6a,b, cyan dots), suggesting that the active area is increasing due to the occurrence of an external seismicity mechanism. On the other hand, during S2 (Figure 5.6a,b, light orange dots) the event-index plot does not show any migration, even after removing the time dependence. Instead, we observe a continuous occurrence of events likely resulting from stress transfer, for which no migration patterns are expected (Figure 5.6b; [START_REF] Helmstetter | Diffusion of epicenters of earthquake aftershocks, law, and generalized continuous-time random walk models[END_REF]).

Effective Stress Drop

We further track the temporal evolution of the effective stress drop (

) measured by comparing the cumulative seismic moment and the areal extent of the sequence [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF][START_REF] Roland | Earthquake swarms on transform faults[END_REF]. The region enclosing the seismic events was measured using a Delaunay triangulation, after projecting all the seismicity onto the main fault plane for S1 (strike N133° and dip 50° according to [START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF]) and an orthogonal fault plane for S2 (where the antithetic fault is active). An example of this process is shown in Figure A4.S9. Following [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF], we impose a distance threshold between neighboring events to avoid outliers, with a maximum triangle leg length of 2.5 km according to the size of the hypocenter cloud. Similar to the previous analysis, the initial window contains 100 events to estimate the rupture area and the cumulative seismic moment (

). We then accumulate event by event, and for each window, we derive the effective stress drop as [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF]. is here the radius of an assumed circle with the same area as estimated from the triangulation. This procedure was carried out individually for each fault, and their respective results are plotted in Figure 5.7.

During S1, we observe a rapid increase in both the radius (the region enclosing seismicity) and the cumulative seismic moment (Figure 5.7a) until about 25 January. Then, both parameters become more stable until F1. The first part (S1) of the sequence releases a total seismic moment of 2.9x10 14 Nm ( Mw 3.6, without considering F1). The behavior of S2 is different : both and rapidly grow, reaching a radius and cumulative seismic moment greater than the values during S1 in a shorter time. Figure 5.7b shows the time evolution of the effective stress drop for both S1 and S2. The comparison between the effective stress drops for S1 and S2 highlights that during S1, the seismicity takes place over a wider area than the area of released seismic moment [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF]. This leads to a lower effective stress drop of 0.01 MPa for S1. Whereas in S2, the higher effective stress drop ( 0.1 MPa) indicates that most of the area enclosing the seismicity is seismically active. These values are of the order of effective stress drops estimated by [START_REF] Roland | Earthquake swarms on transform faults[END_REF] for seismic swarms along Southern California and the East Pacific Rise transform faults. In addition, the difference of almost one order of magnitude between S1 and S2 is similar to differences in the effective stress drop observed by [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF] for different seismic sequences such as injection-induced seismicity, natural earthquake swarms, and mainshock-aftershock sequences. In addition, we analyze the cumulative radius as a function of cumulative (Figure 5.7c). We observe that our measurements are characterized by a cubic scaling of the seismic moment with earthquake cluster radius (

) but for different constant stress drops. Such scaling is predicted in the case of fault models with brittle or mixed (brittle and ductile) rheology and homogeneous prestress [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF]. According to [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF], the scaling together with effective stress differences can be used to discriminate different physical processes driving a seismic sequence.

Discussions

ila earthquake reveals a sudden increase of earthquake activity in January 2009, 3 months prior to the mainshock (Figure 5.2a). From the beginning of the seismicity to the mainshock, almost 5000 foreshocks released a seismic moment of 3x10 15 Nm ( Mw 4.3, Figure 5.7a). Based on our estimated parameters (Section 5.3), we observe that the foreshocks sequence develops in two distinct phases and features a complex spatio-temporal evolution. The two stages of behavior that we report (mostly aseismic, S1, then mostly seismic, S2) has been observed in several other studies, in different tectonic settings (e.g., [START_REF] Durand | A Two Scale Preparation Phase Preceded an M w 5.8 Earthquake in the Sea of Marmara Offshore Istanbul, Turkey[END_REF][START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Ruiz | Nucleation phase and dynamic inversion of the Mw 6.9 Valparaíso 2017 earthquake in Central Chile[END_REF][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].

The first part of the sequence (S1) is characterized by a relatively low temporal interaction of the seismicity (Figure 5.4a), smooth moment release (Figure 5.4b), and a slow but significant movement of the centroid of the seismicity (Figure 5.4c). We also observe migrations lasting up to 7 days (Figures 5.5a and Figure 5.6). The linear velocity of these migrations ranges from 1 to 10 km/day (Figure 5.5a); these velocities are similar to those associated with seismic swarms driven by aseismic slip (e.g., De [START_REF] De Barros | Imbricated aseismic slip and fluid diffusion drive a seismic swarm in the Corinth Gulf, Greece[END_REF]. Finally, we observe migrations in time-space and event-index-space (Figure 5.6), which is indicative that an external seismogenic process controls the seismicity [START_REF] Fischer | The Growth of Earthquake Clusters[END_REF].

If seismicity is a byproduct of aseismic slip, its intermittent time evolution (Figure 5.3a) reflects a variable rate of aseismic slip during the first part of the sequence. Similar behavior is observed during slow slips in subductions zones, with bursts of aseismic slips mainly occurring in rapid episodes associated with bursts of tremors and/or low frequency earthquakes (e.g., [START_REF] Jolivet | The transient and intermittent nature of slow slip[END_REF][START_REF] Rousset | Weeks Long and Years Long Slow Slip and Tectonic Tremor Episodes on the South Central Alaska Megathrust[END_REF]. However, confirming the existence of aseismic slip using independent data as GNSS is difficult, as the expected surface displacement expected during the bursts of seismicity is smaller than the environmental signals often observed in GNSS data along the Apennines [START_REF] Amoruso | Surface temperature and precipitation affecting GPS signals before the 2009 L'Aquila earthquake (Central Italy)[END_REF]. For example, [START_REF] Borghi | Precursory slow-slip loaded the 2009 L'Aquila earthquake sequence[END_REF] proposed that a slow slip event started the 12 February and lasted for almost 2 weeks. However, [START_REF] Amoruso | Surface temperature and precipitation affecting GPS signals before the 2009 L'Aquila earthquake (Central Italy)[END_REF] showed that this signal was due to environmental noise likely caused by temperature and precipitation.

The observed migrations (Figure 5.5a) may also be explained by fluid diffusion [START_REF] Shapiro | Estimating the crust permeability from fluidinjection-induced seismic emission at the KTB site[END_REF]e.g., Ruhl et al., 2016), considering hydraulic diffusivities of 0.8 and 1.5 m2/s, which are within expected values for the crust [START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF][START_REF] Talwani | Pore pressure diffusion and the mechanism of reservoirinduced seismicity[END_REF]. If this was the case, it would be in agreement with the significant role of fluids reported in the region by several authors (e.g., [START_REF] Antonioli | Fluid flow and seismicity pattern: Evidence from the 1997 Umbria Marche (central Italy) seismic sequence[END_REF][START_REF] Lucente | Temporal variation of seismic velocity and anisotropy before the 2009 MW 6.3 L'Aquila earthquake, Italy[END_REF][START_REF] Poli | Seasonal and Coseismic Velocity Variation in the Region of L'Aquila From Single Station Measurements and Implications for Crustal Rheology[END_REF][START_REF] Savage | The role of fluids in earthquake generation in the 2009 Mw 6.3 L'Aquila, Italy, earthquake and its foreshocks[END_REF][START_REF] Terakawa | High-pressure fluid at hypocentral depths in the L'Aquila region inferred from earthquake focal mechanisms[END_REF].

The second part of the sequence (S2) starts with a magnitude 3.9 event (F1) on the 30 of March 2009, activating an antithetic fault (Chiaraluce, Valoroso, et al., 2011;[START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF] similarly to other recent normal fault earthquakes in the region (e.g., [START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF]. The activation of several faults highlights that the precursory process for this event is a complex volumetric process (Ben-Zion & Zaliapin, 2020; [START_REF] Savage | Possible precursory signals in damage zone foreshocks[END_REF], and is not limited to the fault plane. S2 is characterized by a high temporal clustering (Figure 5.4a) and large moment release (Figure 5.4b). These parameters suggest a strong interaction between seismic events, likely governed by stress triggering [START_REF] Freed | Earthquake triggering by static, dynamic, and postseismic stress transfer[END_REF]. No migration is inferred from the event-index analysis (Figure 5.6), and the speed at which seismicity spreads in time is significantly different from that observed during S1. Figures 5.5b and 5.5c show that after F1, the seismicity covers distances of kilometers in seconds to minutes and similar patterns are observed after the occurrence of another magnitude 3.2 events in the middle of S2 (yellow star in Figure 5.5b) and after F2. These velocities are not compatible with mechanisms such as fluid diffusion or aseismic slip, but rather are likely governed by static or dynamic stress transfer [START_REF] Freed | Earthquake triggering by static, dynamic, and postseismic stress transfer[END_REF].

The respective effective stress drops estimated for S1 and S2 are on the order of 0.01 and 0.1 MPa (Figure 5.7b). These values are in agreement with estimations in other seismotectonic contexts (e.g., [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF][START_REF] Roland | Earthquake swarms on transform faults[END_REF][START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas[END_REF], and the difference of between S1 and S2 (Figure 5.7b) provides new insights about the physical mechanisms that might take place during the precursory phase of the studied earthquake. [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF] estimated the effective stress drops for sequences such as hydraulic stimulations of geothermal reservoirs, seismic swarms, and mainshock aftershock-type are associated with effective stress drops from 0.1 to 3.0 MPa, dominating aseismic deformation (e.g., hydraulic fracturing). Considering the above classification, the low effective stress drop ( 0.01 MPa) of S1 (Figure 5.7b) suggests a dominant role of aseismic deformation during the first part of the sequence, with seismicity occurring over a large area with only a small fraction of the area occupied by asperities releasing seismic energy. In this model, aseismic slip is the main mechanism triggering the activation of distant asperities [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF]. Following the models proposed by es a fault with low asperity density embedded in a ductile matrix. On the other hand, the larger effective stress drops up to 0.1 MPa after F1 (Figure 5.7b) indicates that most of the area enclosing the seismicity is seismically active. In this case, the proximity of asperities favors the stress triggering as a mechanism for time clustering of events (Figure 5.4a) over short time scales (Figures 5.5b and 5.5c Both S1 and S2 show a similar cumulative moment versus radius scaling ( ). This scaling is observed either in the case of brittle fault rheology or in the mixed fault rheology models with homogeneous pre-stress, but with the different values of stress drop discussed above [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF]. However, in the case of a partly ductile fault with heterogeneous pre-stress, the seismic moment only scales with the square of the radius , which is not consistent with our observations (Fischer & Hainzl, 2017, Figure 5.7c). Considering that the mixed model is representative of S1, and the brittle model of S2 due to the variations of the effective stress drop (Figure 5.7c), we discuss possible differences between the fault rheologies in S1 and S2.

In the case of brittle asperities embedded in a ductile environment (mixed model during S1), numerical simulations suggest that two scenarios might occur. Either the asperities rupture simultaneously as a single earthquake or separately as individual events, depending on the distance between the asperities and the frictional strength of the ductile region [START_REF] Dublanchet | Interactions and triggering in a 3-D rateand-state asperity model[END_REF][START_REF] Kaneko | Towards inferring earth-quake patterns from geodetic observations of interseismic coupling[END_REF][START_REF] Yabe | Slip-behavior transitions of a heterogeneous linear fault[END_REF]. Thus, high density of the asperities and/or a small frictional parameter in the ductile region lead to simultaneous ruptures of the asperities, while a lower asperity density leads to isolated ruptures, producing a sequence of ruptures with diminished time interaction between each other [START_REF] Dublanchet | Interactions and triggering in a 3-D rateand-state asperity model[END_REF][START_REF] Kaneko | Towards inferring earth-quake patterns from geodetic observations of interseismic coupling[END_REF][START_REF] Yabe | Slip-behavior transitions of a heterogeneous linear fault[END_REF]. During S1, the low effective stress drop (Figure 5.7b) is resulting from void fault areas deformed aseismically among adjacent ruptures, which did not contribute to the seismic moment release [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF], 2021). In this scenario, the existence of large inter-asperities distances is also consistent with the low level of interaction between seismic events inferred from the COV values (Figure 5.4a)

For the case of brittle fault rheology (S2), the fault segment consists of densely distributed asperities that can rupture individually [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF]. For this, some mechanism that prevents the simultaneous rupture of the entire segment and leads to a piecewise rupturing of the fault segment by numerous small earthquakes is needed. Following [START_REF] Yamashita | Pore creation due to fault slip in a fluid-permeated fault zone and its effect on seismicity: Mechanism of earthquake swarm[END_REF] and [START_REF] Aki | Quantitative Seismology, Second Ed[END_REF], possible mechanisms might be the presence of barriers, non-stationary loading, or dilatancy due to pore creation, a process suggested by [START_REF] Lucente | Temporal variation of seismic velocity and anisotropy before the 2009 MW 6.3 L'Aquila earthquake, Italy[END_REF] after the activation of F1. In this model, due to the proximity between asperities, elastic stress plays an important role during the rupture process. This corresponds closely to what our observations indicate during S2: larger COV values (Figure 5.4a) imply episodic and rapid releases of the seismic moment (Figure 5.4b) and seismicity covering larger distances of kilometers in a short time from seconds to minutes (Figures 5.5b and 5.5c).

The observed cubic scaling between the accumulated seismic moment and radius is also indicative of re-rupturing for the two models mentioned above [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF]. The re-rupturing implies significant overlap between regions hosting subsequent seismic events. This behavior is observed during S1, as reactivation of earthquake families during multiple accelerations of seismicity (Figure 5.4c, Movie S1, Figure A4.S12).

The models of [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF] suggest that the rerupturing process is expected to continue until the stress is fully released within the whole fault segment. Interestingly, although the seismicity of S1 occurs on the fault plane that slipped during the mainshock (Chiaraluce, Valoroso, et al., 2011), there is no overlap between the coseismic slip and the foreshocks (Valoroso et al., 2013, Figure A4.S13). This suggests that this part of the fault released the full stress in an intermittent fashion through foreshocks (Figure 5.4b), as the localized fault properties prohibit the nucleation of a large slip episode. Similar behavior has been observed in modelling, where small events appear at the transition from the locked to creeping behavior toward the bottom of the seismogenic zone with decreasing values of the characteristic slip distance of the friction law [START_REF] Lapusta | Nucleation and early seismic propagation of small and large events in a crustal earthquake model[END_REF].

Conclusion

The analysis of our high-resolution seismic catalog highlights different physical mechanisms that each played a role during the precursory phase of the L'Aquila earthquake. Our results demonstrate how the faults involved in the sequence present quantitative differences in the earthquake activity they host. While the seismicity occurring on the main fault up to 1 week before the mainshock (S1) exhibits minimal time clustering, a smooth moment release, slow migrations, and a lower effective stress drop, the seismicity occurring on the antithetic fault after F1 (S2) shows strong episodic clustering and moment release, a rapid spreading of the seismicity and larger effective stress drop. Such differences in the seismicity behavior indicate that while an external process (aseismic or fluid diffusion, or likely a combination of both) is driving the seismicity in S1, stress transfer is the dominant mechanism during S2. A comparison of our observations with recent seismic swarm models [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF] indicates that during S1 a mixed rheology model of sparse brittle asperities embedded in a ductile environment is likely. On the other hand, brittle fault rheology with a dense population of asperities and small inter-asperity distances is more plausible for the antithetic fault during S2.

Our study shows a complex coalescence of different physical processes occurring during the precursory phase of a large earthquake. Moreover, we highlight how the quantitative analysis of spatio-temporal evolution of microseismicity can unveil complex precursory behaviors, which differ from nucleation models based on simple planar faults models [START_REF] Dieterich | Earthquake nucleation on faults with rate-and state-dependent strength[END_REF][START_REF] Liu | Aseismic slip transients emerge spontaneously in three dimensional rate and state modeling of subduction earthquake sequences[END_REF][START_REF] Marone | The effect of loading rate on static friction and the rate of fault healing during the earthquake cycle[END_REF][START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF] aiming for more complex scenarios (e.g., [START_REF] Dutta | Simultaneous Bayesian Estimation of Non Planar Fault Geometry and Spatially Variable Slip[END_REF][START_REF] Shimizu | Construction of fault geometry by finite-fault inversion of teleseismic data[END_REF][START_REF] Zhang | Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics[END_REF].

The catalog generated here is available at https://doi.org/10.5281/zenodo.4776701 (last accessed 20 May 2021).

Chapter 6

Leoncio Cabrera and Piero Poli Article submitted to Geophysical Research Letters (2022)

Introduction

Understanding the physical processes and conditions that lead to the initiation of an earthquake is one of the major challenges of seismology, with implications for earthquake prediction and risk assessment. Therefore, it is crucial to detect and study signals that allow us to relate the rupture of an earthquake to precursory physical processes, if any exist.

Nowadays, we know from laboratory experiments (e.g., [START_REF] Latour | Characterization of nucleation during laboratory earthquakes[END_REF][START_REF] Mclaskey | Earthquake initiation from laboratory observations and implications for foreshocks[END_REF][START_REF] Ohnaka | Scaling of the shear rupture process from nucleation to dynamic propagation: Implications of geometric irregularity of the rupturing surfaces[END_REF] and numerical models (e.g., [START_REF] Ampuero | Earthquake nucleation on rate and state faults Aging and slip laws[END_REF][START_REF] Dascalu | Fault finiteness and initiation of dynamic shear instability[END_REF][START_REF] Kaneko | The onset of laboratory earthquakes explained by nucleating rupture on a fault[END_REF][START_REF] Shibazaki | Transition process from nucleation to high-speed rupture propagation: scaling from stick-slip experiments to natural earthquakes[END_REF] that earthquakes are preceded by different phases: a stable quasi-static deformation phase (phase I), which evolves into an unstable acceleration phase (phase II), after which the large dynamic rupture occurs (phase III). However, direct measurements of the phases I and II in nature are hard, and the scientific community mostly relies on seismological observations, which are perhaps the most informative about the physical processes that precede large earthquakes.

Two of the main types of precursory seismological observations are foreshocks (e.g., [START_REF] Dodge | Detailed observations of California foreshock sequences: Implications for the earthquake initiation process[END_REF][START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF][START_REF] Ellsworth | Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks[END_REF][START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake[END_REF][START_REF] Ruiz | Nucleation phase and dynamic inversion of the Mw 6.9 Valparaíso 2017 earthquake in Central Chile[END_REF][START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF]Cabrera et al., 2022 and references therein) and the seismic signals related to the initial part (i.e., over a few seconds or less) of the mainshock waveform, often known as the seismic nucleation phase (e.g., [START_REF] Abercrombie | Local observations of the onset of a large earthquake: 28 June 1992[END_REF][START_REF] Abercrombie | The nucleation and rupture process of the 1981 Gulf of Corinth earthquakes from deconvolved broad-band data[END_REF]Beroza & Ellsworth, 1996;[START_REF] Colombelli | Evidence for a difference in rupture initiation between small and large earthquakes[END_REF]Ellsworth & Beroza, 1995[START_REF] Marone | The effect of loading rate on static friction and the rate of fault healing during the earthquake cycle[END_REF][START_REF] Kilb | The initial subevent of the 1994 Northridge, California, earthquake: Is earthquake size predictable?[END_REF][START_REF] Iio | Slow initial phase of the velocity pulse generated by microearthquakes[END_REF][START_REF] Mori | Initial rupture of earthquakes in the 1995 Ridgecrest, California sequence[END_REF]Poli et al., 2016;[START_REF] Tape | Earthquake nucleation and fault slip complexity in the lower crust of central Alaska[END_REF][START_REF] Umeda | High-amplitude seismic waves radiated from the bright spot of an earthquake[END_REF]. This initial phase is represented by a small emergent amplitude signal that is sometimes observed before the large impulsive amplitude onset of the P-waves (Fig. 6.1). Results from different studies focused on this initial phase have revealed a variety of new insights about the rupture initiation. [START_REF] Umeda | High-amplitude seismic waves radiated from the bright spot of an earthquake[END_REF] found that the duration of the initial phase increases proportionately with the earthquake size. [START_REF] Iio | Slow initial phase of the velocity pulse generated by microearthquakes[END_REF][START_REF] Iio | Observations of the slow initial phase generated by microearthquakes: Implications for earthquake nucleation and propagation[END_REF] and [START_REF] Iio | Slow initial phase generated by microearthquakes occurring in the western Nagano Prefecture[END_REF] showed that the emergent initial phase can be explained by models that predict slow slip velocities and/or rupture velocities immediately after rupture initiation, and that the initial slow phase is not a product of attenuation, but a source effect. Similarly, Shibazaki & Matsu'ura (1997) modeling far-field seismograms showed that the slow initial phase can be radiated in the theoretical acceleration phase (phase II) from the slow growth of a rupture. Using a larger dataset of earthquakes, Ellsworth & Beroza (1995) and Beroza & Ellsworth (1996) found that the size and duration of the initial phase scale with the final magnitude. In addition, selected group of earthquakes and seismic sequences have been studied in more detail such as the 1981 Gulf of Corinth earthquake [START_REF] Abercrombie | The nucleation and rupture process of the 1981 Gulf of Corinth earthquakes from deconvolved broad-band data[END_REF], the 1995 Ridgecrest events [START_REF] Mori | Initial rupture of earthquakes in the 1995 Ridgecrest, California sequence[END_REF][START_REF] Ellsworth | Observation of the seismic nucleation phase in the Ridgecrest, California, earthquake sequence[END_REF], the 1994 Northridge California earthquake [START_REF] Kilb | The initial subevent of the 1994 Northridge, California, earthquake: Is earthquake size predictable?[END_REF], intermediate depth events in the Hindu Kunsh nest (Poli et al., 2016) and events in central Alaska [START_REF] Tape | Earthquake nucleation and fault slip complexity in the lower crust of central Alaska[END_REF], among others.

On April 6, 2009 the Mw 6.1 earthquake struck central Italy causing damage and fatalities (Chiaraluce et al., 2011). It was preceded by thousands of foreshocks (Cabrera et al., 2022) and an emergent initial phase recorded by near-field accelerometers (Fig. 6.1a, [START_REF] Ellsworth | Supershear during nucleation of the 2009 M 6.3 Italy earthquake[END_REF][START_REF] Di Stefano | Fault zone properties affecting the rupture evolution of the 2009 (Mw 6.1) L'Aquila earthquake (central Italy): Insights from seismic tomography[END_REF]. In this work we study this initial phase to get new insights into the rupture initiation process of earthquakes. For this aim we take advantage of a dense azimuthal coverage of stations and we derive quantitative parameters [START_REF] Abercrombie | Earthquake directivity, orientation, and stress drop within the subducting plate at the Hikurangi margin, New Zealand[END_REF] of the rupture initiation (Fig. 6.1b). Then, the results are compared with theoretical models (e.g., [START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF] and previous observations (e.g., Beroza & Ellsworth, 1996). We do not analyze the large dynamic rupture after the initial phase, as previous studies have done so (e.g., [START_REF] Atzori | Finite fault inversion of DInSAR coseismic displacement of the 2009 L'Aquila earthquake (central Italy)[END_REF][START_REF] Cheloni | Coseismic and initial post-seismic slip of the 2009 M w 6.3 L'Aquila earthquake, Italy, from GPS measurements[END_REF][START_REF] Cirella | Complexity of the rupture process during the 2009 L'Aquila, Italy, earthquake[END_REF][START_REF] Scognamiglio | Fast determination of moment tensors and rupture history: What has been learned from the 6 April 2009 L'Aquila earthquake sequence[END_REF][START_REF] Trasatti | Finite element inversion of DInSAR data from the Mw 6.3 L'Aquila earthquake, 2009 (Italy)[END_REF]. Our results reveal how the earthquake struggled to begin due to the heterogeneity of the physical properties of the fault, which can be understood in the context of well-known nucleation models [START_REF] Rubin | Earthquake nucleation on (aging) rate and state faults[END_REF][START_REF] Dascalu | Fault finiteness and initiation of dynamic shear instability[END_REF][START_REF] Kaneko | The onset of laboratory earthquakes explained by nucleating rupture on a fault[END_REF][START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF] as a significantly large slip-weakening distance (Dc). 2011). Black and green thin lines represent traces of the active mapped faults and co-seismic surface ruptures, respectively [START_REF] Boncio | Coseismic ground deformation of the 6 April 2009 L'Aquila earthquake (central Italy, Mw6. 3)[END_REF]. The beach-ball (compressional quadrants in colors) represents the source mechanism (reported by the Istituto Nazionale di Geofisica e Vulcanologia, INGV) for the mainshock and EP (yellow diamond, Table 6.1). (b) Vertical component velocity seismograms [START_REF] Luzi | ORFEUS Working Group 5[END_REF] of the initial P-waves recorded for stations at different distances from the epicenter. Red and blue vertical lines indicate the arrival times of EP and IP, respectively.

Analysis of the Rupture Initiation

As illustrated in Figure 6.1b (see also Fig. A5.S1), the seismograms recorded at triggered regional accelerometers (Fig. 6.1a) show a short (~0.6-s) small emergent amplitude signal (Fig. 6.1b, EP), before the onset of large impulsive amplitude P-waves (Fig. 6.1b, IP). By locating the hypocenter of both phases (i.e., EP and IP), Di Stefano et al. (2011) interpreted EP as the initial rupture process of the earthquake, and defined it as a seismic nucleation phase (Ellsworth & Beroza, 1995;[START_REF] Ellsworth | Supershear during nucleation of the 2009 M 6.3 Italy earthquake[END_REF]. In this work we will call this initial phase the rupture initiation.

The rupture initiation keeps on growing from EP until the large amplitude onset (IP). Figure A5.S2 exemplifies this feature. This behavior is confirmed by comparing the waveforms of the rupture initiation and two Mw 3.9 foreshocks that occurred close to the hypocenter and were recorded at the closest broadband station located right on top of the hypocenters. We here assume that the rupture initiation is one event with duration from EP to IP. At the termination of the rupture, at time IP, the large slip identified by the large amplitude seismic waves onset (IP) is starting, in agreement with other studies [START_REF] Cirella | Complexity of the rupture process during the 2009 L'Aquila, Italy, earthquake[END_REF][START_REF] Colombelli | Evidence for a difference in rupture initiation between small and large earthquakes[END_REF][START_REF] Di Stefano | Fault zone properties affecting the rupture evolution of the 2009 (Mw 6.1) L'Aquila earthquake (central Italy): Insights from seismic tomography[END_REF]Ellsworth & Beroza, 1995[START_REF] Marone | The effect of loading rate on static friction and the rate of fault healing during the earthquake cycle[END_REF][START_REF] Iio | Slow initial phase of the velocity pulse generated by microearthquakes[END_REF]Poli et al., 2016).

We measure the time duration between EP and IP ( ) at 15 strong motion stations located in the near field (epicentral distance, <60 km) with good azimuthal coverage (Fig. 6.1; [START_REF] Luzi | ORFEUS Working Group 5[END_REF], for which both EP and IP are clearly registered and easily identified (see examples in Fig. 6.1b and for all of the stations in Fig. A5.S1). For this purpose, we use the raw vertical velocity seismograms from the Engineering Strong Motion Database (ESM) of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) recorded with a sampling rate of 200 Hz [START_REF] Luzi | ORFEUS Working Group 5[END_REF]. This high sampling rate permits to have a time resolution of 0.005s and allows us to ensure a correct recording of the rupture initiation (Beroza & Ellsworth, 1996). We do not apply any filter, and we manually pick the time arrival for EP and IP (see examples in Figs. 6.1, A5.S1). The time as a function of the azimuth, reveals a clear azimuthal variation (Fig. 6.2a), which indicates directivity due to the finiteness of the rupture process [START_REF] Abercrombie | Earthquake directivity, orientation, and stress drop within the subducting plate at the Hikurangi margin, New Zealand[END_REF]. Beyond the azimuthal directivity, we further observe how downgoing and up-going P-waves, that have similar azimuths, have different duration (e.g., see differences for the take-off angles of SUL-AQK and LSS-AQG-AQA-AQV), suggesting a significant vertical component for the rupture initiation [START_REF] Abercrombie | Earthquake directivity, orientation, and stress drop within the subducting plate at the Hikurangi margin, New Zealand[END_REF].

Before proceeding with further analysis, we assess whether the recorded signals and their azimuthal dependence are controlled by propagation effects. For this purpose, we model our observations as a Gaussian pulse of duration 0.6s (our smallest time measurement, Fig. A5.S3). We then apply an attenuation factor (Anderson & Hough, 1984) for the different travel times, for each station, and test a range of quality factors Q (see Text A5.S1 and Fig. A5.S4). The resulting waveforms (Fig. A5.S4) show a marginal effect of attenuation on the waveform shape and related duration, thus supporting the robustness of our azimuthal measures (Fig. 6.2a). The analysis of attenuation effects also reveals that our time measurements ( 0.6s) are systematically larger than potential attenuation effects, and can thus be used to infer rupture properties. This is in concordance with the work done by [START_REF] Iio | Slow initial phase generated by microearthquakes occurring in the western Nagano Prefecture[END_REF] who showed that the rupture initiation does not arise from the attenuation, but is a source effect. In addition, according to Beroza & Ellsworth (1996) for M>4.0 earthquakes, the duration of this rupture initiation should be greater than ~0.3s (considering their proposed scaling) and the attenuation should not have any role (Beroza & Ellsworth, 1996). Furthermore, we ruled out the existence of other seismic phases by modeling their arrival time using a local velocity model (Chiaraluce et al., 2011) and the ObsPy TauP Toolkit [START_REF] Beyreuther | ObsPy: A Python toolbox for seismology[END_REF][START_REF] Crotwell | The TauP Toolkit: Flexible seismic traveltime and ray-path utilities[END_REF]. Figure A5.S5 shows that the times observed for IP do not correspond to any other phase arrival.

We finally quantitatively study the rupture initiation in terms of its geometry (i.e., rupture plane orientation and rupture velocity), and estimate some rupture parameters.

Geometry of the Rupture Initiation

To retrieve the geometry of the rupture initiation (i.e., first ~0.6s of the signal), we use an approach that takes advantage of the azimuthal variation in the seismic radiation (i.e., directivity) and that has been used to analyze small earthquakes [START_REF] Abercrombie | Earthquake directivity, orientation, and stress drop within the subducting plate at the Hikurangi margin, New Zealand[END_REF].

Our main assumption is that the rupture can be represented by a line kinematic source model with uniform slip and constant rupture velocity. This model, although simplified, has shown good results in the study of small earthquakes in which kinematic inversions using heterogeneous models to account for high frequencies are not plausible (e.g., [START_REF] Abercrombie | Earthquake directivity, orientation, and stress drop within the subducting plate at the Hikurangi margin, New Zealand[END_REF]. In addition, considering a line source model facilitates comparison with theoretical parameters such as the critical rupture length (Lc, Ohnaka, 2000. See discussion section). Following [START_REF] Abercrombie | Earthquake directivity, orientation, and stress drop within the subducting plate at the Hikurangi margin, New Zealand[END_REF], the rupture duration ( ) at station given a line source of size is given by (6.1) where is the rupture velocity, V is the P-( ) or S-( ) wave velocity at the source, is the take-off angle, and and are the azimuth and dip of the rupture direction, respectively (the take-off angles and the dip are defined as 0°=vertically down, 90°=horizontal, and 180°=vertically up). We use Vp=6.5 km/s, based on a local three-dimensional tomographic model [START_REF] Di Stefano | Fault zone properties affecting the rupture evolution of the 2009 (Mw 6.1) L'Aquila earthquake (central Italy): Insights from seismic tomography[END_REF]. The take-off angles are calculated using the ObsPy TauP Toolkit [START_REF] Beyreuther | ObsPy: A Python toolbox for seismology[END_REF][START_REF] Crotwell | The TauP Toolkit: Flexible seismic traveltime and ray-path utilities[END_REF] in a local velocity model (Chiaraluce et al., 2011) considering the EP hypocenter previously estimated by Di Stefano et al. (2011). Then, assuming in (1), we use a grid search procedure to estimate the , , and that best explain the observations minimizing an L1 norm [START_REF] Abercrombie | Earthquake directivity, orientation, and stress drop within the subducting plate at the Hikurangi margin, New Zealand[END_REF]. In addition, the grid search is repeated 1000 times, randomly removing two stations each time. The result is stored every time, and the ensemble of the measures is used to evaluate the uncertainty of the estimates. Table 6.1 shows a summary of the results. The estimated strike ( ) and dip ( ) of the rupture initiation are, respectively, and (where dip is defined as 0°=vertically down, 90°=horizontal, and 180°=vertically up), which indicates a north-east initial rupture perpendicular to the fault system (N137°E, Chiaraluce et al., 2011) and going upwards, in agreement with the dip of the fault system of ~50° (Chiaraluce et al., 2011, Fig. 6.2a), kinematics models (e.g., [START_REF] Cirella | Complexity of the rupture process during the 2009 L'Aquila, Italy, earthquake[END_REF] and absolute locations of EP and IP (Fig. 6. 1, Di Stefano et al., 2011). Our estimated strike and dip of the rupture initiation are also consistent with the focal mechanism reported by the INGV for the mainshock (Fig. 6.1), which indicates that the geometry of the rupture initiation is similar to the geometry of the large dynamic rupture. We also corroborate the relative geometry by relocating EP and IP using the Growclust software [START_REF] Trugman | GrowClust: A hierarchical clustering algorithm for relative earthquake relocation, with application to the Spanish Springs and Sheldon, Nevada, earthquake sequences[END_REF]. For this purpose, we considered all the measured values at the strong motion stations (Fig. 6.1) as the differential delay times and a local velocity model (Chiaraluce et al., 2011). In addition, we performed two tests considering different initial locations for EP and IP, and another test considering the same initial location for EP and IP (see details in Text A5.S2). In both cases, the geometries retrieved from the relocation indicate a north-east and upwards propagating rupture (see Tables A5.S2-A5.S3), supporting the robustness of our estimations.

Our results also indicate a relatively slow rupture velocity ( ) of km/s and a linear rupture length ( ) of km for the rupture initiation. Assuming a constant rupture velocity, the duration of the rupture initiation is ~0.66 s (we discuss more in detail in the discussion section). Our slow rupture velocity is in agreement with an initial region of lower rupture velocity around EP reported by [START_REF] Cirella | Complexity of the rupture process during the 2009 L'Aquila, Italy, earthquake[END_REF]. Di Stefano et al. ( 2011) estimated a larger rupture velocity for EP (~2.2 km/s) from a 2km rupture between EP and IP using absolute locations. However, absolute locations are less sensitive to the relative geometry between EP and IP. Our relative relocation using GrowClust estimates a distance ~500m in agreement with L (see Tables A5.S2-A5.S3). ) for stations that have azimuths in common but different take-off angles (e.g., stations SUL-AQK and LSS-AQG-AQA-AQV) with a misfit of 0.015 s. This implies that we can reproduce the observed up-going and down-going P-waves, and that the estimate of the dip is therefore robust. (colored according to the takeoff angle) as a function of the azimuth compared with the best solution synthetics (orange dots). (b) Misfit between the dip of the rupture initiation (i.e., 0°=vertically down, 90°=horizontal, and 180°=vertically up) as a function of the strike of the rupture initiation. The black star contour indicates the minimum misfit and the red star indicates the dip and an angle perpendicular to the strike of the fault system (50° and N47°E, respectively according to [START_REF] Chiaraluce | The 2009 Aquila (central Italy) seismic sequence[END_REF] Misfit between the length of the linear source as a function of the rupture velocity. The black star contour indicates the minimum misfit.

Magnitude of the Rupture Initiation

To estimate the magnitude of the EP, we follow [START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF] and we compare the maximum amplitudes of the EP waveform with that of foreshock events within a threedimensional radius of 500 m, where the Mw is reported in a local seismic catalog (Chiaraluce et al., 2011). This procedure is more suitable than using the spectrum of the waveforms, as the time window of the signals is very short and they have a non-zero base-line (see [START_REF] Ruiz | Historical and recent large megathrust earthquakes in Chile[END_REF].

Since the strong motion records correspond to a trigger system that is activated when a certain acceleration threshold is exceeded, there are no foreshocks recorded for the stations shown in Fig. 6.1. Therefore, we use the four closest broadband stations (See Fig. A5.S7a;INGV Seismological Data Centre, 2006;MedNet Project Partner Institutions, 1990), for which the EP and IP identification is straightforward, and extract their waveform. Thus, we compute the magnitude of EP based on the mean value of the maximum amplitude ratios for all the vertical channels between the foreshocks and EP, assuming that a ratio of 10 in the amplitude ratios corresponds to a variation of one-unit of magnitude (e.g., Cabrera et al., 2021[START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF][START_REF] Essing | Spatiotemporal evolution of the Seismicity in the Alto Tiberina Fault System revealed by a Template Matching Catalog[END_REF]Frank et al., 2017;[START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF]. We perform the same procedure for 11 foreshocks (Fig. A5.S7a) considering a time window identical to at the corresponding station and starting at the P-wave arrival. The selected time window is long enough to capture the foreshocks P-waves, and at the same time short enough to avoid interfering with the S-wave signals (see for example Fig. A5.S2b). It is important to note that, although the broadband seismic records are clipped shortly after the onset of IP for the mainshock, they are not between EP-IP (see Fig. A5.S8). We obtain a magnitude of 3.99 0.18 for EP (see Fig. A5.S7a). The robustness of our estimation is observed in Figure A5.S2b, which shows the comparison between the waveforms of EP-IP and two Mw 3.9 foreshocks that occurred one week and 5 h before the mainshock. We note that although their amplitudes are similar, the amplitude of EP grows with time until the onset of IP (Figs. A5.S1, A5.S2). The seismic moment ( ) of the rupture initiation is 1.23 ×10 15 Nm, which corresponds to ~0.07% of the mainshock seismic moment (see more detail below).

Rupture Parameters of the Rupture Initiation

In this section we use the measures obtained above ( and ) to derive some rupture parameters of the rupture initiation. We assume a circular shear crack model with spatially constant stress drop, having a radius ( of 300 m (L/2) and a shear modulus ( of 30 GPa. Although a circular crack model may appear simplified, its adoption is reasonable given the small magnitude of the rupture initiation (sec. 6.2.2). Moreover, this model makes our results directly comparable with studies on conventional earthquakes (e.g., [START_REF] Venkataraman | Observational constraints on the fracture energy of subduction zone earthquakes[END_REF] and previous observations of seismic nucleation phases (Beroza & Ellsworth, 1996).

The average slip and the stress drop are estimated as (Beroza & Ellsworth, 1996) and [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], respectively. Our estimations indicate an average slip of 0.15 m which is in agreement with the slip estimated by [START_REF] Cirella | Complexity of the rupture process during the 2009 L'Aquila, Italy, earthquake[END_REF] for the region where the rupture begins and a stress drop of 19.9 MPa. Fig. A5.S9 shows that the estimated rupture parameters for the rupture initiation are in agreement with the scaling that relates seismic nucleation phases and their subsequent mainshocks proposed by Beroza & Ellsworth (1996). In addition, Fig. A5.S9 shows that the seismic moment of the rupture initiation is 0.01% to 10% of the seismic moment of the mainshock, in agreement with previous studies (Beroza & Ellsworth, 1996;Poli et al., 2016). Beyond the rupture parameters discussed above, we use the measured rupture velocity (Table 6.1) to infer the seismic efficiency assuming a crack model with constant rupture velocity (Husseini & Randall, 1976;Kanamori et al., 1998): (6.2) describes the proportion of the total energy radiated as seismic waves. If , the rupture is very efficient in radiating energy. If , the entire energy is dissipated mechanically and no energy is radiated (see [START_REF] Kanamori | The diversity of the physics of earthquakes[END_REF]. Thus, considering a shear wave velocity of according to a local velocity model (Chiaraluce et al., 2011), we obtain for the rupture initiation, which implies that nearly 80% of the energy was dissipated during the initial rupture.

We can now combine the stress drop ( ) with the seismic efficiency ( ) to estimate the total energy ( ), fracture energy ( ) and energy density ( ) using the expressions derived by [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] and Poli & Prieto (2016) , we obtain J and J. Finally, assuming that the fracture energy rate does not vary on the fault area we can estimate the fracture energy density ( as (6.5)

We estimate J/m 2 for the rupture initiation. This corresponds to the upper bound of previous estimated values of for events with an average slip of 0.15 m (e.g., [START_REF] Abercrombie | Can observations of earthquake scaling constrain slip weakening?[END_REF]; see a compilation in Cocco et al., 2022), which is explained by the higher percentage of the energy budget associated with fracturing energy during the initiation rupture.

Discussions

-s small amplitude rupture (~0.07% of the total moment), which preceded the large dynamic rupture (Fig. 6.1). From our geometrical analysis, we infer that the rupture initiation occurred in the fault plane where the large dynamic rupture takes place.

The inferred parameters for the rupture initiation (duration, , , ) follow the scaling observed in previous studies (Beroza & Ellsworth, 1996;[START_REF] Colombelli | Evidence for a difference in rupture initiation between small and large earthquakes[END_REF][START_REF] Iio | Slow initial phase of the velocity pulse generated by microearthquakes[END_REF]Poli et al., 2016;Fig. A5.S9), which suggests that scale-dependent quantities control the underlying physics of earthquake initiation [START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF].

Beyond the well-studied parameters for rupture initiation (e.g., Beroza & Ellsworth, 1996;[START_REF] Colombelli | Evidence for a difference in rupture initiation between small and large earthquakes[END_REF][START_REF] Iio | Slow initial phase of the velocity pulse generated by microearthquakes[END_REF]Poli et al., 2016, Fig. A5.S9), we also estimate the rupture velocity (Table 6.1). The inferred rupture velocity ( km/s) is slower than that of ordinary crustal earthquakes. For instance, considering (Chiaraluce et al., 2011), we obtain . This value contrasts with the usually observed for crustal earthquakes [START_REF] Kanamori | Frictional melting during the rupture of the 1994 Bolivian earthquake[END_REF]. However, is on the order of the estimated rupture velocities for tsunami earthquakes, which are about 1.0 km/s (e.g., Fuji & Satake, 2007;[START_REF] Bryant | The underrated hazard[END_REF]. These slow rupture velocity for tsunami earthquakes led to greater directivity [START_REF] Bryant | The underrated hazard[END_REF], as we here observe for the rupture initiation (Fig. 6.1, 6.2). Since tsunami earthquakes mainly occur at subduction zones, the slower rupture velocities have been associated to anomalous elastic properties of the uppermost part of the subduction zone, where a complex accretionary wedge is developed in the hanging wall of the megathrust (Sallarés & Ranero, 2019).

Furthermore, our estimated seismic efficiency ( ) indicates an ineffective rupture initiation in terms of radiated seismic energy. [START_REF] Venkataraman | Observational constraints on the fracture energy of subduction zone earthquakes[END_REF] compiled observations for several subduction zone earthquakes, and show that most earthquakes have seismic efficiencies in the range 0.25-1. By contrast, tsunami earthquakes and some deep earthquakes have very small seismic efficiencies (<0.25) and hence dissipate a large amount of energy during faulting. [START_REF] Venkataraman | Observational constraints on the fracture energy of subduction zone earthquakes[END_REF] proposed that low seismic efficiency for shallow megathrust earthquakes is a marker of complex rheology near the trench.

The low efficiency for the initial rupture implies a significant amount of fracture energy ( of total energy). This latter is indicative of a large slip-weakening distance (Dc, [START_REF] Beeler | Inferring earthquake source properties from laboratory observations and the scope of lab contributions to source physics[END_REF]Poli & Prieto, 2016). The large (and thus large Dc) is known to be the main parameter that controls the occurrence of small precursory phases in numerical models [START_REF] Lapusta | Nucleation and early seismic propagation of small and large events in a crustal earthquake model[END_REF]. The large Dc is also indicative of significant heterogeneity of fault physical properties [START_REF] Lapusta | Nucleation and early seismic propagation of small and large events in a crustal earthquake model[END_REF][START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF][START_REF] Scuderi | The role of fluid pressure in induced vs. triggered seismicity: Insights from rock deformation experiments on carbonates[END_REF].

The ensemble of our parameters indicates a high level of complexity during rupture initiation, which results in a slow rupture velocity and associated low seismic efficiency. One factor controlling this complexity could be the significant role of fluids reported in the region [START_REF] Antonioli | Fluid flow and seismicity pattern: Evidence from the 1997 Umbria Marche (central Italy) seismic sequence[END_REF]Pio et al., 2010;[START_REF] Poli | Seasonal and Coseismic Velocity Variation in the Region of L'Aquila From Single Station Measurements and Implications for Crustal Rheology[END_REF][START_REF] Savage | The role of fluids in earthquake generation in the 2009 Mw 6.3 L'Aquila, Italy, earthquake and its foreshocks[END_REF][START_REF] Terakawa | High-pressure fluid at hypocentral depths in the L'Aquila region inferred from earthquake focal mechanisms[END_REF]. This hypothesis is supported by laboratory experiments indicative of decreasing rupture velocity with increasing fluid pressure [START_REF] Passelègue | Initial effective stress controls the nature of earthquakes[END_REF]. Another factor could be the presence of a mixed rheology, with brittle asperities embedded in a ductile environment in the segment where the rupture starts, as proposed by [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF] from the study of the foreshocks preceding the mainshock. A third factor, somehow linked with our second hypothesis, that could play a significant role, is the effective contact surface on the fault. Although it is not possible to measure this parameter directly on the fault, sliding laboratory experiments using two brittle acrylic blocks separated by a rough interface (e.g., [START_REF] Rubinstein | Detachment fronts and the onset of dynamic friction[END_REF][START_REF] Gvirtzman | Nucleation fronts ignite the interface rupture that initiates frictional motion[END_REF] have shown that dynamic ruptures (i.e., laboratory earthquakes) are preceded by slow rupture fronts that gradually reduce the contact surface until sliding occurs. Considering our observations and the factors described above, our interpretation 6.3. Finally, if we assume that the initial rupture tracks the local physical properties of the initiation region, are our derived rupture parameters representative of the energy balance of the rupture initiation for slip-depend models [START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF]? Beyond the observed slow (Table 6.1), which is predicted for initial rupture stage [START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF], we note that the measured (~20 MPa) is in the order of the breakdown stress estimated by Ohnaka (2000) (1-100 Mpa). Furthermore, using the scaling relationships proposed by [START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF] for slipdepend models we can estimate the respective Dc and critical nucleation length (2Lc) associated with a Mw 6.1 mainshock (Note that Ohnaka calls 2Lc what we and other authors call Lc, i.e., the critical nucleation length, e.g., [START_REF] Mclaskey | Earthquake initiation from laboratory observations and implications for foreshocks[END_REF]. Namely, [START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF] proposes and , where represents the seismic moment of the mainshock. Thus, our estimated average slip and rupture size (Table 6.1) are relatively on the same order as and estimated from [START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF]. We thus propose that the derived rupture parameters from the detailed analysis of the initial rupture stage of the earthquake are closely related to scale-dependent quantities of slipdependent nucleation models.

Conclusion

Our full characterization of the rupture initiation including geometry, rupture velocity, and some other rupture parameters reveals important new insights about earthquake initiation, and helps to bridge the gap between laboratory experiments [START_REF] Latour | Characterization of nucleation during laboratory earthquakes[END_REF][START_REF] Mclaskey | Earthquake initiation from laboratory observations and implications for foreshocks[END_REF], numerical modelling [START_REF] Kaneko | The onset of laboratory earthquakes explained by nucleating rupture on a fault[END_REF], theoretical studies [START_REF] Ampuero | Earthquake nucleation on rate and state faults Aging and slip laws[END_REF][START_REF] Ohnaka | A physical scaling relation between the size of an earthquake and its nucleation zone size[END_REF] and observations of faults in nature.

We show that the earthquake started with a slow and seismically inefficient rupture initiation, akin to what is observed for tsunami earthquakes in complex media. The slow seismic efficiency implies a large fracture energy and thus large Dc, which is also indicative of significant heterogeneity in the fault. Based on the results listed above, we interpret the rupture initiation takes place in a complex environment, where factors such as the presence of fluids, a mixed rheology and the effective contact surface play a key role. In addition, our approach considering a line kinematic source model, although simplified, provides a new opportunity to compare the rupture length and average slip of the rupture initiation with theoretical parameters such as the critical rupture length (Lc) and slip-weakening distance (Dc), indicating that the formers are closely related to scale-dependent quantities of slip-dependent nucleation models.

Open Research

Data was downloaded from the Engineering Strong Motion Database (ESM) of the Istituto Nazionale di Geofisica e Vulcanologia [START_REF] Luzi | ORFEUS Working Group 5[END_REF]. Once in the webpage https://esmdb.eu/#/home (last accessed December 2, 2022), navigate to the section and choose the earthquake (ESM IDIT-1980-0012). Waveforms can be downloaded from the section. Computations were performed using the facilities of the University of Grenoble Alpes (UGA) High-Performance Computing infrastructures CIMENT (https://ciment.univ-grenoble-alpes.fr/wiki-pub/index.php/Welcome_to_the_CIMENT_site!). Access to the cluster requires registration by the UGA.

Article in preparation for submission

This chapter corresponds to work still in progress. In particular, we develop a methodology to classify the seismicity of a dense seismic catalog automatically, based on physical properties of the seismicity. To this aim, we first create a new seismic catalog for Central Italy (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) and then we classify the seismicity using hierarchical clustering. Although, the catalog generated in this chapter and the clustering methodology are finished, we are still working on the analysis and interpretation of the results in more detail. However, the main advances, results, conclusions and perspectives obtained so far are discussed.

Introduction

The analysis of earthquake catalogs is at the core of many studies ranging from scientific applications (e.g., [START_REF] Kato | The generation of large earthquakes[END_REF]) to risk assessment (e.g., Gerstenberger et al., 2016). These catalogs can provide valuable information on the type of sequences taking place (e.g., mainshock-aftershocks, swarm, foreshocks-mainshock-aftershocks, [START_REF] Mogi | Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (third paper)[END_REF], the mechanical properties of the faults where they occur, and the processes that drive them (e.g., stress transfer, aseismic slip, [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF]. Given the importance of seismic catalogs, seismology is currently experiencing a vertiginous increase in the generation of high-resolution catalogs covering years of records and increasingly large regions (e.g., [START_REF] Chamberlain | Illuminating the and Phases of the 2016 M7. 8 Earthquake With 10 Years of Seismicity[END_REF][START_REF] Ross | Searching for hidden earthquakes in Southern California[END_REF]. These advances have been made possible by three main drivers: denser networks of stations, improvements in computing, and the development and implementation of new techniques and data algorithms (See [START_REF] Arrowsmith | Big Data Seismology[END_REF] for a recent review). As a consequence, seismic catalogs are also becoming increasingly large and thus more complex to analyze.

A classic way to analyze a large seismic catalog is through the analysis of the spatiotemporal evolution of the seismicity (e.g., [START_REF] Fischer | The Growth of Earthquake Clusters[END_REF]. Recently some unsupervised learning strategy called clustering has been utilized to explore seismic catalogs and signals as well (e.g., Holtzman et al., 2018;Schoenball & Ellsworth, 2018;[START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF]Cesca et al., 2020 and references there in). In this approach, the events within the seismic catalog are automatically classified or grouped according to the similarity of their attributes (e.g., [START_REF] Cesca | Seiscloud, a tool for density-based seismicity clustering and visualization[END_REF]. In contrast to supervised learning strategies, clustering does not rely on a labeled training set and human expert knowledge [START_REF] Goodfellow | Deep learning[END_REF]. Thus, clustering algorithms can help to automatically classify earthquakes in order to extract information on seismic processes and faulting patterns out of large seismic datasets.

In this work, we show how clustering can be used to classify and study a dense seismic catalog using an unsupervised approach. To do so, we start by creating via template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF] a new and denser seismic catalog that connects both spatially and uila earthquake and the Amatrice-Visso-Norcia sequence (Mw 6.0, 5.9 and 6.5 respectively) in central Italy. Previous works have covered more limited regions, or shorter time windows, but none of them in their entirety (e.g., [START_REF] Chiaraluce | The 2016 central Italy seismic sequence: A first look at the mainshocks, aftershocks, and source models[END_REF][START_REF] Improta | Multi-segment rupture of the 2016 Amatrice-Visso-Norcia seismic sequence (central Italy) constrained by the first high-quality catalog of Early Aftershocks[END_REF][START_REF] Michele | structure of the 2016 2017 Central Italy seismic sequence from data recorded at the Italian National Network[END_REF][START_REF] Moschella | Minor earthquake sequences in the Amatrice-Norcia epicentral area (Central Italy)[END_REF][START_REF] Tan | earthquake catalog reveals how complex fault structures were activated during the 2016 2017 Central Italy sequence[END_REF]Vicic et al., 2020;[START_REF] Vuan | Loading rate variations along a midcrustal shear zone preceding the Mw6. 0 earthquake of 24 August 2016 in Central Italy[END_REF][START_REF] Waldhauser | Fault Planes, Fault Zone Structure and Detachment Fragmentation Resolved With Aftershock Locations of the 2016 2017 Central Italy Sequence[END_REF]. By using ~32,000 templates reported by INGV, 19 broadband stations and more than 31,000 hours of parallel GPU computations, we obtain a new and denser catalog of ~280,000 events (Section 7.3). We then explore this catalog using different physical parameters of the seismicity (features) and agglomerative clustering to classify seismicity into different groups (Section 7.4). Finally, we discuss what are the physical properties of each group, and what mechanism could be driving them (Sections 7.5 and 7.6).

Seismotectonic Setting

The Central Apennines have been the site of destructive seismic sequences in the last quarter century. In 1997 the Mw 6.0 Colfiorito earthquake ruptured the northwestern region [START_REF] Ripepe | Foreshock sequence of september 26th, 1997 Umbria-Marche earthquakes[END_REF]Fig. 7.1). Twelve years later, in 2009, the Mw 6.1 L'Aquila earthquake activated the southeastern section (Chiaraluce et al., 2011, Fig. 7.1). Finally, in 2016 and in between the two previously described regions, the Amatrice-Visso-Norcia sequence occurred, with moment magnitudes 6.0, 5.9 and 6.5, respectively (Chiaraluce et al., 2017, Fig. 7.1). Interestingly, while the Colfiorito and L'Aquila events were preceded by foreshocks (e.g., [START_REF] Ripepe | Foreshock sequence of september 26th, 1997 Umbria-Marche earthquakes[END_REF][START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF]Cabrera et al., 2022 and references therein), the 2016 sequence appears to have started without prior activity [START_REF] Chiaraluce | The 2016 central Italy seismic sequence: A first look at the mainshocks, aftershocks, and source models[END_REF][START_REF] Vuan | Loading rate variations along a midcrustal shear zone preceding the Mw6. 0 earthquake of 24 August 2016 in Central Italy[END_REF].

The region is dominated by a NW-SE Quaternary normal fault system (Pizzi & Galadini, 2009;Fig. 7.1), where an extension rate of 2-3 mm/yr perpendicular to the Apennines is . Here, the northern zone called Umbria-Marche domain contains the Monte Vettore-Norcia fault (Fig. 7.1). Further south in the Latium-Abruzzi domain, the Campotosto and Paganica faults are found (Fig. 7.1). The northern and southern zones are connected along the thrust Olevano-Antrodoco-Sibillini (OAS) tectonic alignment (Fig. 7.1).

Previous work has focused on the use of seismic catalogs to illuminate geological structures, although generally covering smaller areas. Chiaraluce et al. (2011) and Valoroso et occurred in 2009 revealing a complex anatomy of the fault segments with multiple antithetic and synthetic faults. After the 2016 sequence, [START_REF] Chiaraluce | The 2016 central Italy seismic sequence: A first look at the mainshocks, aftershocks, and source models[END_REF] and [START_REF] Improta | Multi-segment rupture of the 2016 Amatrice-Visso-Norcia seismic sequence (central Italy) constrained by the first high-quality catalog of Early Aftershocks[END_REF] studied the region around the Amatrice-Visso-Norcia epicenters using aftershocks. While the former used the reported aftershocks up to November 2016, the latter analyzed what they call "early aftershocks" (the seismicity of the first two days after each mainshock), but applying a more robust relocation scheme. [START_REF] Michele | structure of the 2016 2017 Central Italy seismic sequence from data recorded at the Italian National Network[END_REF] extended this analysis by using aftershocks until January 2018. Subsequently, Spallarossa et al. (2020) and [START_REF] Tan | earthquake catalog reveals how complex fault structures were activated during the 2016 2017 Central Italy sequence[END_REF] generated dense aftershock catalogs using automatic detection and machine learning schemes, respectively. Further, [START_REF] Waldhauser | Fault Planes, Fault Zone Structure and Detachment Fragmentation Resolved With Aftershock Locations of the 2016 2017 Central Italy Sequence[END_REF] relocated the catalog generated by Spallarossa et al. ( 2020) reaching a higher spatial resolution, and [START_REF] Moschella | Minor earthquake sequences in the Amatrice-Norcia epicentral area (Central Italy)[END_REF] studied some minor earthquake sequences in the Amatrice-Norcia epicentral area. All these works, focused in the Amatrice-Visso-Norcia region, revealed an intricate fault system, composed of two main segmented sub-parallel normal faults, NNW striking and west-dipping, and a fault complexity 2021). Finally, Vicic et al. (2020) studied the seismicity occurred between 2008 and 2016 for the L'Aquila and the 2016 sequence regions using template matching. However, within their templates they did not consider seismicity from 2010 to 2015, leading to a strong bias in the detection capability. One question that remains open is what the dynamics of seismicity is like in these different regions. Therefore, in this work we focus on different quantitative properties that seismicity exhibits in order to classify it using an unsupervised methodology. [START_REF] Gessel | The HIKE European Fault Database (EFDB) compiled in the framework of the GeoERA project HIKE[END_REF].

Extending the Seismic Catalog

We apply template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF] Institutions, 1990). These stations cover homogeneously the region between the epicenter of the L'Aquila and Amatrice-Visso-Norcia earthquakes (see Fig. 7.1). Moreover, they have a stable operability over time. Figure A6.S1 shows the number of available channels and the operability of each station over time. For each day, we only consider stations with less than two hours of data gaps (if there are any gaps). Data was continuously recorded at a sampling rate of 100 Hz, and we downsampled to 40 Hz to improve the computational efficiency. We also apply a 5-18 Hz filter to avoid the strong anthropogenic noise in the Apennines (see for example [START_REF] Poli | Seasonal and Coseismic Velocity Variation in the Region of L'Aquila From Single Station Measurements and Implications for Crustal Rheology[END_REF][START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF]. The dataset was then organized into 24-hour continuous files with all gaps filled with zeros.

We consider as templates the events reported by the INGV catalog in the same time period between 2009 to 2016, within the region defined by latitudes 41.75°N to 43.15°N, and longitude 12.5°E to 14.5° (~155 km x 220 km, Fig. 7.1), for which at least 15 phases (P and S) are manually picked by the professional operators of the INGV, of which at least 5 are P phase. Thus, we create an initial catalog of 31,988 seismic events, defined as the 1.5 s time windows that start 0.5 s before the manually picked P-and S-wave arrivals at each station for the vertical and horizontal components, respectively, and filtered in an identical manner to the continuous data (bandpassed between 5-18 Hz). The template waveforms are then correlated against a sliding window of continuous data using a GPU-architecture, more than 31,000 hours of parallel computing and the Fast Matched Filter algorithm [START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF] to obtain daily correlation functions. We search sample-by-sample considering a detection threshold of 10 times the daily median absolute deviation (MAD) of the correlation function averaged over all stations and channels to detect events significantly similar to the template. In addition, we impose a minimum number of 12 channels available to declare a detection. To remove double detections over the same time window, we merge consecutive detections with differential times less than 1.5 s; we keep the detection with the highest network-averaged correlation coefficient as the final detection. We estimate the magnitude of each new event by computing the median value of the maximum P-and S-wave amplitude ratio between the template event and the detection for all the channels. Using the template catalog magnitude as a reference, the magnitude of a detected event is determined, assuming that a ratio of 10 of the amplitude corresponds to a variation of one-unit of magnitude (e.g., [START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF]Frank et al., 2017;Cabrera et al., 2021). Finally, we consider the newly detected events to have occurred at the same hypocenter as the template reported by the INGV. Our new catalog contains 264,432 events with magnitude ranging from -2.0 to 6.5 (Fig. 7.2). This new catalog is the first and the largest catalog that connects both spatially and temporally the seismicity between the 2009 earthquake and the 2016 Amatrice-Visso-Norcia sequence (Fig. 7.3). A comparison of the time and space evolution of the new events is shown in Figure 7.3, together with the initial catalog. We observe that the new catalog contains almost nine times more events than the initial one (templates), while it preserves a similar time evolution for the cumulative number of events (Fig. 7.3a). This is seen as a first order quality estimate, as it reproduces the general evolution of the seismicity (e.g., [START_REF] Essing | Spatiotemporal evolution of the Seismicity in the Alto Tiberina Fault System revealed by a Template Matching Catalog[END_REF]. For example, the linear growth of the cumulative number of events between 2010 and the Amatrice event in 2016 (Fig. 7.3a) indicates that there is a component associated with background seismicity. On the other hand, the major seismicity rate increases follow the larger events in 2009 and 2016. In addition, our new catalog also captures in much greater detail the existence of different bursts of seismicity occurring at different locations in the fault system (Fig. 7.3b,c; note that strike 0 km corresponds to the epicenter of Amatrice, yellow star), which is also reflected in the daily count of events presented in Fig. 7.3a.

In the next section, we describe how we use the clustering algorithm to explore this new catalog in an automatic and data-driven way.

Clustering the Seismic Catalog

The general idea of clustering (also called unsupervised classification or exploratory data analysis) is to group objects (unlabeled data) according to the similarity between them (e.g., [START_REF] Xu | Clustering[END_REF]. The similarity between objects is measured based on a set of characteristics of each object called features. Thus, after grouping objects (based on some metric and method), objects within the same cluster will be more similar to each other than objects in separate clusters. In the following sections we describe the features used, and explain in detail the clustering algorithm employed.

Feature Preparation

For the preparation of the features, we make spatial groups of events ("seismic groups" hereinafter). For every template, we define a radius around it. Then, we define as a seismic group all seismicity contained within a sphere around the hypocenter of the template. Finally, for each seismic group we estimate eight physical parameters which have direct physical interpretation, which are used as features for the clustering algorithm. In this work we present the results using a radius of 2 km and minimum of 10 events within the sphere, although we performed tests with radii between 0.5-5 km showing a strong stability of the results.

In the following, we describe in detail the series of estimated features that can be related to seismicity properties. From the interevent times (i.e., the time between two consecutive events) we estimate the (1) mean, (2) median and (3) the coefficient of variation (COV, [START_REF] Kagan | Long-term earthquake clustering[END_REF]. While the interevent times provides information about the seismicity rate, the COV provides information about the temporal interaction of the seismicity, that is, how past events affect the occurrence of the future ones (e.g., [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF]Cochrane et al., 2018;[START_REF] Essing | Spatiotemporal evolution of the Seismicity in the Alto Tiberina Fault System revealed by a Template Matching Catalog[END_REF][START_REF] Kagan | Long-term earthquake clustering[END_REF][START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas[END_REF]. COV is estimated as the ratio of the standard deviation of the interevent times to the average interevent time. The COV is 0 for a periodic occurrence of seismicity, 1 for completely random Poisson occurrence, and larger than 1 for temporal interactive seismicity such as during mainshock-aftershock sequences. In simple words the larger the COV is, the stronger the interactions are (e.g., [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF]Cochrane et al., 2018;[START_REF] Essing | Spatiotemporal evolution of the Seismicity in the Alto Tiberina Fault System revealed by a Template Matching Catalog[END_REF][START_REF] Kagan | Long-term earthquake clustering[END_REF][START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas[END_REF]. We also estimate (4) the effective duration (difference between the 10th and 90th percentiles of the origin times, Cochrane et al., 2018), (5) the number of families (defined as the number of different templates within the sphere) and ( 6) the cumulative seismic moment (

). These parameters provide a measure of how long the seismicity lasts, the number of "asperities" that are active, and the total amount of seismic moment, respectively. Considering that we use a sphere of fixed radius, also gives an idea of the seismic moment released for a fixed volume (33.5 km 3 for a radius of 2 km), akin to effective stress drop [START_REF] Fischer | Effective stress drop of earthquake clusters[END_REF]. Furthermore, we analyze the time evolution of the seismic moment release by estimating (7) the ratio of the maximum seismic moment to the total seismic moment within the sphere ( , [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF]. While a stable and gradual moment release by many earthquakes without a dominant large magnitude event is observed for swarm-type sequences [START_REF] Vidale | A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers[END_REF], most of the moment is released at once during mainshockaftershocks sequences [START_REF] Mogi | Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (third paper)[END_REF]. We finally estimate (8) the time delay of the largest event normalized by the mean time delay of all the events ( , Chen et al., 2012). While values of indicate that the largest event occurs close to the beginning of the sequence (mainshock-aftershocks-type), large values of indicate that the main event occurs later in the seismic sequence (swarm-type, Chen et al., 2012). It is worth noting that within the features we do not provide spatial information such as latitude, longitude or depth.

Since all the parameters have different units, it is necessary to bring them to a similar scale. To do so, we normalize the feature matrix using the standard scaler method, which standardize features by removing the mean and scaling to unit variance [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. In the case of the accumulated seismic moment, we first take the logarithm of the accumulated seismic moment. Once the normalized feature matrix has been constructed, we proceed to clustering.

Hierarchical Clustering

Different algorithms are available to perform clustering (e.g., K-means, [START_REF] Lloyd | Least squares quantization in PCM[END_REF], density-based clustering, [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF], among others). In this work we use hierarchical clustering with a bottom-up approach, namely agglomerative clustering. In this approach, the similarity between all the objects in the data set (e.g., a specific distance in the feature space) is summarized in a similarity matrix [START_REF] Johnson | Hierarchical clustering schemes[END_REF]. At the beginning of the process, each object starts in a single cluster. Then, the clusters begin to merge based on the similarity matrix until all objects are unified into a single global cluster [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF]. The clustering process is summarized in a dendrogram (see Section 7.4.3), which reveals the hierarchical structure of the whole data set. One of the main advantages of this algorithm over other clustering algorithms, such as k-means, is that there is no need to establish a specific number of clusters beforehand. In addition, the nature of the seismological data is in accordance with the hierarchical structure [START_REF] Steinmann | Hierarchical exploration of continuous seismograms with unsupervised learning[END_REF].

To estimate the distance metric (similarity) and perform the hierarchical clustering, we use the normalized feature matrix obtained in the previous section and the Ward's minimum variance method [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF]. This method has been also used to cluster seismic events depending on their pairwise waveform similarity (e.g., Sanchez-Reyes et al., 2021) as well as continuous seismograms based on waveform features (e.g., [START_REF] Steinmann | Hierarchical exploration of continuous seismograms with unsupervised learning[END_REF]. An advantage of our approach using physical parameters as features is that, once the clustering is done, we can associate the parameters of the objects within each cluster and compare them with a physical sense.

Dendrogram

After the clustering is performed, the results are summarized in a dendrogram (e.g., Fig.

7.4a).

A dendrogram is a diagram of the distance matrix (similarity), and shows the hierarchical relationship between objects (e.g, [START_REF] Xu | Clustering[END_REF]. The key to interpret a dendrogram is to focus on the height at which any two objects are joined (e.g., vertical axis in Fig. 7.4a, [START_REF] Xu | Clustering[END_REF]. Simply put, the lower the link height, the greater the similarity between the objects. To extract different cluster solutions based on the distance matrix, the definition of a distance threshold is needed. While a larger distance threshold allows forming less clusters which are composed of more objects, a smaller distance threshold results in a larger number of clusters but composed of less objects. We show our results defining two different distance thresholds.

The first distance threshold is 250, which allows us to classify the seismicity into two large clusters labeled A and B (Fig. 7.4a). Cluster A corresponds to seismicity scattered over the entire study area (Fig. 7.4b), and contains 56.3% of the seismic groups (Fig. 7.4d). On the other hand, cluster B is more localized in space along the NW-SE fault system (Fig. 7.4b), and contains the remaining 43.7% of the seismic groups (Fig. 7.4d). We notice that although we do not provide spatial information to the clustering algorithm (e.g., latitude, longitude, depth), clusters A and B exhibit clear spatial patterns (Fig. 7.4b), which emerge from the physical properties represented in the features. When observing the average values of the different features (also called centroids) associated with each cluster (Fig. 7.4c), we notice clear differences between A and B. Cluster A has higher COV values than B (around 8 and 4, respectively), indicating a higher interaction of the events within cluster A. Cluster A also has lower interevent times and a shorter effective duration than B, indicating that in general the seismic groups in cluster A are active for a shorter time. In addition, it is observed that on average the seismic groups in cluster A have more families per volume (~200) than B (~50), and a higher cumulative seismic moment. The last two parameters and are quite similar for A and B, although with slightly higher values for A. The second distance threshold is 180, which gives us 6 clusters. Three subclusters for A labeled A1, A2 and A3, and three subclusters for B labeled B1, B2 and B3 (Fig. 7.5a). Considering now these six clusters, it is possible to see a greater degree of detail and variability among them in the centroid of their parameters (Fig. 7.5b), the distribution of their percentages (Fig. 7.5c), and their location (Fig. 7.6). First, we will focus on A1, A2 and A3. From the dendrogram (Fig. 7.5a) we can notice that A1 separates from A2 and A3 at a greater distance (more height in the dendrogram, ~215) than that at which A2 and A3 join each other (~190). This indicates that A2 and A3 are more similar to each other than they are to A1. In terms of percentages, A1 represents 3.8% of the total seismic groups and corresponds to the smallest cluster compared to A2 and A3 which represent 19.0% and 33.4%, respectively. If we look at the centroids of the different parameters for A1, A2 and A3 (Fig. 7.5b) we see that these three clusters always have either the highest values (e.g., COV, number of families per volume, and cumulative Mo) or the lowest values (median and mean interevent times as well as effective duration), except for and . In these two parameters, A2 stands out for higher and A1 for higher (Fig. 7.5b). From the map view presented in Figure 7.6 we can see that A2 and A3 present similar spatial distributions aligned along the projections of the normal fault system. On the other hand, A1 corresponds to events much more localized in space, and preferentially to the south-east of the Olevano-Antrodoco-Sibillini (OAS) thrust alignment in the Latium-Abruzzi domain [START_REF] Chiaraluce | The 2016 central Italy seismic sequence: A first look at the mainshocks, aftershocks, and source models[END_REF]Fig. 7.6). This localization in space of A1 is also evident from the along-strike temporal evolution of all the seismic groups that make up the cluster showed in Fig. 7.7. A1 corresponds mainly to seismicity activated after the occurrence of the L'Aquila earthquake in 2009, while A2 and A3 contain significant amount of seismicity which is already active before the earthquake. It is interesting to note that for A1, A2 and A3 the cumulative number of events shows different accelerations between 2009 and 2017. However, while clusters A2 and A3 contain a large number of earthquakes occurring after the 2016 mainshocks, A1 does not and captures mainly events occurring only after L'Aquila. Finally, if we plot some cross sections (Fig. 7.8) perpendicular to the alignment of the normal fault system we can see that the localization of A1, A2 and A3 also occurs at depth. In general, the three clusters do not exceed 12km depth, with A1 being the most concentrated. Focusing on clusters B1, B2 and B3 (Fig. 7.5a) we notice that B1 is a little further from B2 and B3, and corresponds to the cluster with the lowest number of events, representing only 0.3% (Fig. 7.5b). On the other hand, B2 and B3 represent 3.7% and 39.7% respectively, the latter being the biggest cluster (Fig. 7.5b). From the average values of the different parameters for B1, B2 and B3 (Fig. 7.5b) we see that from B1 to B3 the COV increases progressively, while the interevent times and effective duration decrease progressively. Moreover, in general B3 is the cluster that always differs from the other two (e.g., largest number of families per volume, highest cumulative seismic moment and lowest

). Figure 7.6 shows that clusters A1 and A2 correspond to the most scattered clusters in the study area without seismicity along the normal fault system. On the other hand, although A3 shows spread seismicity as well, it also contains seismicity along the normal fault system. Figure 7.7 shows that in general the three clusters cover homogeneously the distribution along strike (contrary to A1, A2 and A3), although with a much smaller number of events for B1 and B2. It is interesting to note that the cumulative seismicity curves (Fig. 7.7) show a linear increase for B1 and B2 over the entire time period, unaffected by the mainshocks. Something similar occurs for A3, although only between mid-2009 and until the occurrence of Amatrice. However, the behavior of this latter is much less episodic than for A1, A2 and A3. Finally, from Figure 7.8 we observe that in contrast to the concentration of A1, A2 and A3 at depth, clusters B1, B2 and B3 are spread over a much wider range of depths. In addition, they also cover a much larger area in the direction perpendicular to the strike.

Discussion

From the results presented above, we can observe that the seismicity is divided into two main clusters: on the one hand, the seismicity of cluster A, which is characterized by a higher temporal clustering (higher COV), shorter durations, higher number of families per volume and higher seismic moment (Fig. 7.4c). Looking at the spatial distribution of this seismicity, we see that it is mainly localized and aligned along the NW-SE normal fault system of the area (Fig. 7.4b). On the other hand, the seismicity of cluster B (Fig. 7.4a, 7.5a) is characterized by lower temporal interactions, longer durations, fewer families per volume and lower seismic moment (Fig. 7.4c). In addition, the spatial distribution of cluster B contrasts with that of cluster A, the former showing a high spatial spreading over the entire study region (Fig. 7.4b). This first result indicates that our hierarchical clustering approach using physical parameters of seismicity as features allows effectively classifying seismicity in an automatic and data-driven way. Furthermore, the classification of such seismicity shows that seismicity close to each other is characterized by similar properties, indicating that the different physical parameters chosen here as features contain somehow the spatial information not used during the clustering (e.g., latitude, longitude, depth).

Furthermore, the centroid of the parameters associated with each cluster allow a first interpretation of the type of seismicity of which the clusters are composed (Fig. 7 A for example, presents higher COV and a larger number of families per volume, indicating a higher level of interaction between seismicity (e.g., [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF][START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF][START_REF] Essing | Spatiotemporal evolution of the Seismicity in the Alto Tiberina Fault System revealed by a Template Matching Catalog[END_REF]Sanchez-Reyes et al., 2021;[START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas[END_REF]. In addition, cluster A shows larger than B, and shorter durations. This indicates that for cluster A, a larger amount of seismic moment is released in overall shorter duration sequences than B. This is congruent with a higher of A, indicative that seismic moment is released more episodically than in B (e.g., [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF][START_REF] Neves | A Earthquake Catalog for the 2004 Mw 6 Parkfield Earthquake Sequence Using a Matched Filter Technique[END_REF]. Considering all the aforementioned characteristics, we interpret cluster A as a composite of interactive seismicity, which points more to mainshock-aftershock type sequences and some swarm-like component. On the other hand, the lower COV and longer duration of cluster B indicate less interaction between seismicity (Fig. 7.4c). This is consistent with a considerably smaller number of families per volume as shown by the centroid of the number of families (Fig. 7.4c), indicating that the seismicity is more isolated as illustrated by the spatial distribution in Fig. 7.4b. Based on these characteristics, we interpret that contrary to what occurs in A, cluster B corresponds more with seismicity similar to background seismicity.

.4c). Cluster

As mentioned in Section 7.4.3, going deeper into the hierarchical structure we can analyze smaller clusters (in the number of objects that compose them, in this case seismic groups), but that present more similarities between them. Thus, the large seismicity dataset can be decomposed into six clusters (Fig. 7.5a), which are discussed below.

A1 is the most spatially localized cluster (Fig. 7.6). In general, it has similar characteristics with A2 and A3, but stands out for a smaller amount of accumulated seismic moment ( ) and the prominently largest value (Fig. 7.5). While the lower may indicate that there are no events of significant magnitude (mainshocks), the large value is indicative of a swarm-like sequence (e.g., Chen et al., 2012). In addition, the temporal evolution of A1 shows that the highest amount of seismicity occurs after the earthquake, followed by an episode starting a couple of months afterwards (Fig. 7.7). Based on this, we interpret that A1 corresponds to a combination of mainshock-aftershock seismicity but also with a strong swarm component.

A2 and A3 present very similar behaviors (high COV, high , low ). However, A2 stands out for a much higher than A3 (~0.65 and 0.25, respectively), indicating that in general there are larger dominant events in A2 [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF]. This is corroborated by looking at the magnitude distribution of each cluster, which shows that A2 contains the largest number of events of significant magnitude (M>4, Fig. A6.S2). Based on the above, we interpret that A2 corresponds to mainshock-aftershock type sequences, while A3 preferentially contains aftershocks.

B1, B2 and B3 correspond to the clusters with the lowest values of COV, number of families per volume, and the highest effective durations (Fig. 7.5b). Among the three, B1 and B2 represent the lowest percentages of seismic groups (Fig. 7.5c), and stand out for being the most scattered in the study area (Fig. 7.6) and with the lowest number of families per volume (Fig. 7.5). The latter indicates that in general they correspond to more spatially isolated seismicity. In addition, the cumulative number of events shows a linear trend over the entire time period, which is typical of background seismicity [START_REF] Bagh | Background seismicity in the Central Apennines of Italy: The Abruzzo region case study[END_REF].

On the other hand, B3 corresponds to a sort of transitional cluster. It shares the characteristics of B1 and B2, but among the three clusters, it is the one with the closest average parameters to the values of A1, A2 and A3. B3 differs from B1 and B2 in its higher COV (~4) and , as well as the lowest effective duration and among the three clusters. By observing the cross sections in Fig. 7.8 we notice that A3 also corresponds to the cluster that most extends both perpendicular to the strike, as well as in depth.

Based on the above, we interpret B1 and B2 as pure background seismicity. On the other hand, B3 corresponds to a combination of background seismicity, but probably also swarm and mainshock-aftershock (high COV) seismicity.

Conclusions

In this work we proposed a new way to explore a dense seismic catalog and to identify different types of seismicity in an automatic and data-driven way. Our approach is based on hierarchical clustering, which allows us to obtain different levels of classification according to the dendrogram, thus facilitating the exploration of the seismicity. The hierarchical clustering is applied to a set of features which are related to physical properties of the seismicity. A first advantage of this is that it facilitates the interpretation of the clusters obtained once the classification has been performed. To do so, it is only necessary to relate the classified clusters of seismicity with their respective features, and then interpret these features (physical parameters). Another advantage of our approach with respect to other classification methods such as K-means [START_REF] Lloyd | Least squares quantization in PCM[END_REF] is that there is no need to establish a specific number of clusters beforehand, which reduces the need for prior expertise in the process.

Our results show that although we do not provide spatial information to the clustering algorithm (e.g., latitude, longitude, depth), the spatial distribution of seismicity is somehow reflected in the different parameters used as features here. Thus, we can extract two main seismicity clusters. One represented by seismicity with low interaction levels (e.g., lower COV, number of families per volume and ), longer durations and scattered throughout the study area, that we interpret as mainly background seismicity. And on the other hand, a second group with high level of interaction (e.g., higher COV, number of families per volume and ), shorter durations and located along the main fault system, that we interpret mainly as mainshock-aftershock and swarm-type sequences. In addition, we show that within the two aforementioned clusters there are smaller clusters that differ from each other. We interpret these differences to be associated with the type of sequence that dominates the cluster.

Finally, our methodology is generalizable to the use of any other parameter of interest to perform the classification. Therefore, identifying which are the best physical parameters to extract a certain type of seismicity is fundamental. For this, the detailed study of how the different physical processes manifest themselves in seismicity presents an opportunity that should be explored extensively.

Chapter 8 Conclusion and Outlook

Conclusion

In this work we have used massive databases of continuous seismic records to create new high-resolution seismic catalogs and study different seismic sequences (e.g., Cabrera et al., 2021Cabrera et al., , 2022;;[START_REF] Herrera | Rupture properties of the 2020 Mw 6.8 Calama (northern Chile) intraslab earthquake. Comparison with similar intraslab events in the region[END_REF][START_REF] Poli | Volcanic Origin of a Swarm in the Central Bransfield Basin, Antarctica[END_REF]. Subsequently, we have quantitatively analyzed these seismic catalogs and complemented our analysis with other observations (e.g., GNSS) and modeling (e.g., seismic source, thermal model). All this has made it possible to elucidate the mechanisms that drive seismicity in different regions of the planet. Revealing the physical mechanisms that drive seismicity is important because it allows us to understand the physical processes that take place before, during and after a major earthquake, as well as the properties of sequences in which there is no major earthquake, such as a swarm sequence. The main conclusions reached in each study are presented below, followed by general conclusions.

The first chapter shows how mainshocks that are apparently very similar (e.g., magnitude, rupture geometry, stress drop), are followed by aftershock sequences which have different productivities (i.e., number of events). Although stress transfer seems to be the common mechanism for the six aftershock sequences studied, it is clear that the tectonic context in which they occur plays a role. We show that the six events are located at different depths within the slab, and thus at different thermal conditions. We also observe that there is a dramatic decrease in aftershock activity below the 400 coincides with the neutral plane in the outer-rise region. Based on the dominant tensional stress regime, where main shocks with normal focal mechanism occur up to 40 km from the top of the slab (i.e., 700 450 °C isotherm depth separates high-and low-hydrated zones, which define the aftershock productivity. Additionally, this hypothesis is supported by the study of three earthquakes (Mw~6.8) aligned latitudinally [START_REF] Herrera | Rupture properties of the 2020 Mw 6.8 Calama (northern Chile) intraslab earthquake. Comparison with similar intraslab events in the region[END_REF], but occurring at different depths in the Chilean subduction zone presented in the Appendix 1. In this latter, the combination of dynamic source modeling, template matching, and a 3-D tomographic model indicates that the level of hydration plays a role in the mainshock and the aftershock productivity. Thus, the study of sequences that can be compared with each other as shown in Chapter 3 and Appendix 1 opens an opportunity towards a more detailed understanding of the tectonic factors that manifest themselves through differences in seismicity patterns.

In the second chapter we have shown how to obtain information from a seismic swarm occurred in a remote and instrumentally limited region in Antarctica. Based on the distribution of the interevent times of more than 36,000 events, S-P time analysis, higher b-value (b=1.6), geodetic deformation and closeness to the submarine Orca Volcano, we conclude that the seismic swarm is related to a volcanic process, rather than a result of the extensional tectonics. This exemplifies how rapid deformation episodes occurring at the ridge axial volcanic structures plays a main role in modulating the long-term extension (7 mm/yr between the Antarctica Plate and the South Shetland microplate, [START_REF] Taylor | Kinematics and segmentation of the South Shetland islands-Bransfield Basin system, northern Antarctic Peninsula[END_REF]). In addition, we provide evidence of how the clustering of the seismicity (e.g., [START_REF] Beaucé | Systematic detection of clustered seismicity beneath the Southwestern Alps[END_REF][START_REF] Duverger | Dynamics of microseismicity and its rela-tionship with the active structures in the western Corinth Rift (Greece)[END_REF] can help us to differentiate between different physical processes such as tectonic extension or volcanism. However, due to the limited data in the region we cannot discuss in detail the volcanic driving mechanism of the deformation, although hydrothermal fluids (e.g., [START_REF] Reiss | The impact of complex volcanic plumbing on the nature of seismicity in the developing magmatic natron rift, Tanzania[END_REF] and/or magma flows [START_REF] Heimisson | Physically consistent modeling of Dike-induced deformation and seismicity: Application to the 2014 Bárðarbunga Dike, Iceland[END_REF]) can be hypothesized.

The third chapter shows how we can track a precursory phase from the statistical analysis of a dense seismic catalogue. By analyzing the time clustering, seismic moment release, centroid migrations, and the effective stress drop we show quantitative differences that point to different physical mechanisms (first aseismic slip and then stress transfer) and different rheologies (first ductile and brittle, and then brittle) of the faults involved. This indicates that the interaction between seismic and aseismic processes play an important role in the precursory phase of earthquakes, which differ from the classic nucleation models such as purely cascade or preslip (e.g., [START_REF] Ellsworth | Seismic evidence for an earthquake nucleation phase[END_REF]. Instead, an interplay of the two seems more plausible (e.g., [START_REF] Cattania | Precursory slow slip and foreshocks on rough faults[END_REF][START_REF] Mclaskey | Earthquake initiation from laboratory observations and implications for foreshocks[END_REF].

The fourth chapter presents how we use a line source approach to study the rupture initiati initiation (strike and dip) is in agreement with the geometry of the mainshock and the fault system. In addition, we estimate a slower rupture velocity (Vr= km/s), and thus a low seismic efficiency (

). We interpret these results as signs that the rupture is struggling to start, which makes about 80% of the energy budget going to fracture energy. We propose that these difficulties in initiating the rupture could be associated with the complexities of the medium. For instance, due to the presence of a mixed rheology, with brittle asperities embedded in a ductile environment as we proposed in chapter 5 [START_REF] Cabrera | Tracking the evolution of foreshocks preceding the Mw 6.1 2009 Earthquake[END_REF], or the significant role of fluids reported in the region [START_REF] Antonioli | Fluid flow and seismicity pattern: Evidence from the 1997 Umbria Marche (central Italy) seismic sequence[END_REF]Pio et al., 2010;Poli et al., 2020;[START_REF] Savage | The role of fluids in earthquake generation in the 2009 Mw 6.3 L'Aquila, Italy, earthquake and its foreshocks[END_REF][START_REF] Terakawa | High-pressure fluid at hypocentral depths in the L'Aquila region inferred from earthquake focal mechanisms[END_REF]. In addition, slow fronts have also been observed in laboratory experiments (e.g., [START_REF] Gounon | Rupture Nucleation on a Periodically Heterogeneous Interface[END_REF][START_REF] Gvirtzman | Nucleation fronts ignite the interface rupture that initiates frictional motion[END_REF][START_REF] Rubinstein | Detachment fronts and the onset of dynamic friction[END_REF]. These fronts gradually decrease the effective contact surface, leading to a complete rupture once the contact decreases sufficiently to no longer support the stress (e.g., [START_REF] Rubinstein | Detachment fronts and the onset of dynamic friction[END_REF]. However, to elucidate whether a similar process can take place in nature, future observations are required (see outlook section).

Finally, the fifth chapter shows an example of how unsupervised learning (clustering) can be used to explore a dense seismic catalogue. We first create the longest seismic catalog (~ and the 2016 Amatrice-Visso-Norcia sequence, and then apply agglomerative clustering. Our results indicate that the physical parameters used here as features are not randomly distributed over the study area. Instead, they contain information about how seismicity is distributed and behaves. It is important to mention that although we do not provide spatial information to the clustering algorithm (e.g., latitude, longitude, depth), our results show that the clusters do present spatial patterns. While the seismicity with low levels of interaction and longer duration is scattered throughout the study area, the seismicity with higher levels of interaction and shorter duration is aligned along the main fault system. We interpret this as one group primarily composed of background seismicity, and another composed mostly by mainshock-aftershocks and swarms. Our results show new opportunities to explore dense seismic catalogs from the classification of physical parameters that have direct seismological interpretation. Furthermore, they show the potential of using temporal clustering (note that we refer to the interactions between earthquakes [e.g., [START_REF] Beaucé | Systematic detection of clustered seismicity beneath the Southwestern Alps[END_REF][START_REF] Cattania | Complex earthquake sequences on simple faults[END_REF][START_REF] Duverger | Dynamics of microseismicity and its rela-tionship with the active structures in the western Corinth Rift (Greece)[END_REF], and not to the algorithm) as an important input and characteristic in the study of earthquake dynamics.

To summarize, we provide new observations that shed light on the physical mechanisms that control different seismic sequences and push further the borders of our understanding. In addition, we show how the medium in which these sequences occur plays a role in their signature, and how such information can be retrieved by combining different tools and analyses. A general characteristic that comes to light as a result of the use of dense and detailed seismic catalogs is the level of complexity involved during the progression of the seismic sequences in comparison with the current mechanical models that we have to explain such seismic sequences. Examples of sequences such as those studied in Chapter 3 and Appendix 1, as well as those in Chapter 5, show the variety of factors (e.g., temperature, hydration level, rheology) and processes (e.g., aseismic slip, stress transfer) that are manifested even on a small spatial scale. This highlights the need for better and more detailed mechanical models that allow, for example, to take into account seismic and aseismic processes as part of a whole.

Outlook and Perspectives

Collect more and new observations for comparison

The general methodology we have used in this work consists of generating a highresolution seismic catalog, and then analyzing it in detail with different tools, as well as being complemented with other observations. Although this allows us to obtain information on the processes and characteristics of the environment in which a sequence takes place, this has also shown us how different these processes and media can be. Therefore, it is only through a large number of observations that we will be able to obtain more robust and thus generalizable conclusions. Currently there are many efforts to study earthquake precursory phases in seismotectonic settings such as strike-slip faults (e.g., [START_REF] Bouchon | Extended nucleation of the 1999 Mw 7.6 Izmit earthquake[END_REF][START_REF] Chen | California foreshock sequences suggest aseismic triggering process[END_REF][START_REF] Dodge | Detailed observations of California foreshock sequences: Implications for the earthquake initiation process[END_REF][START_REF] Dodge | Foreshock sequence of the 1992 Landers, California, earthquake and its implications for earthquake nucleation[END_REF][START_REF] Durand | A Two Scale Preparation Phase Preceded an M w 5.8 Earthquake in the Sea of Marmara Offshore Istanbul, Turkey[END_REF][START_REF] Ellsworth | Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks[END_REF][START_REF] Shelly | A high resolution seismic catalog for the initial 2019 Ridgecrest earthquake sequence: Foreshocks, aftershocks, and faulting complexity[END_REF][START_REF] Tape | Earthquake nucleation and fault slip complexity in the lower crust of central Alaska[END_REF][START_REF] Yoon | Foreshocks and mainshock nucleation of the 1999 M w 7.1 Hector Mine, California, Earthquake[END_REF], subduction zones [START_REF] Bouchon | The long precursory phase of most large interplate earthquakes[END_REF][START_REF] Kato | Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake[END_REF][START_REF] Ruiz | Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake[END_REF][START_REF] Ruiz | Nucleation phase and dynamic inversion of the Mw 6.9 Valparaíso 2017 earthquake in Central Chile[END_REF], and extensional regimes [START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF]. However, no universal and generalizable results or patterns have been observed. In this sense, the study of seismic sequences of lower magnitude presents an opportunity (e.g., [START_REF] Meng | Immediate foreshocks indicating cascading rupture developments for 527 M 0.9 to 5.4 Ridgecrest earthquakes[END_REF][START_REF] Simon | Possible Precursory Slow Slip to two ML 3 Mainevents of the Diemtigen Microearthquake Sequence[END_REF], since there are many more events of magnitude 4-5 or smaller than over 6. In addition, this facilitates the study of more sequences in the same region, thus reducing the differences between sequences associated with a change in the tectonic regime. For this, proper instrumentation together with detection, location and analysis methodologies will play a fundamental role (see a recent review in [START_REF] Arrowsmith | Big Data Seismology[END_REF].

Testing Mechanical Models

As mentioned in the Cection 8.1 Conclusion, a common feature that emerged from our results is an increase in the level of complexity observed for the seismic sequences compared to current mechanical models. For example, the results of Chapter 5 suggest faults with different characteristics, and an evolution from a process preferentially driven by aseismic slip and later by seismic processes. This contrasts with the conceptual models of cascade and preslip currently under debate (Ellsworth & Beroza, 1995;[START_REF] Gomberg | Unsettled earthquake nucleation[END_REF], and points to the fact that this process could be a combination of the two, as indicated by laboratory experiments (e.g., McLasckey, 2019) and numerical modeling (Cattania & Segall, 2019). Moreover, this increase in the complexity of the process before an earthquake has also been observed in other seismic sequences (e.g., [START_REF] Durand | A Two Scale Preparation Phase Preceded an M w 5.8 Earthquake in the Sea of Marmara Offshore Istanbul, Turkey[END_REF][START_REF] Sánchez Reyes | The Imbricated Foreshock and Aftershock Activities of the Balsorano (Italy) Mw 4.4 Normal Fault Earthquake and Implications for Earthquake Initiation[END_REF][START_REF] Shelly | A high resolution seismic catalog for the initial 2019 Ridgecrest earthquake sequence: Foreshocks, aftershocks, and faulting complexity[END_REF][START_REF] Tape | Earthquake nucleation and fault slip complexity in the lower crust of central Alaska[END_REF], indicating that this could be a general feature that can be observed with highresolution seismic catalogs. Therefore, assessing for example whether the transition from aseismic to seismic behavior is part of the precursory process of the earthquake or not is of paramount importance.

To incorporate the additional complexity observed in the evolution of seismic sequences, at least two tasks with respect to modeling must be addressed: on the one hand, comparing our new observations with the current mechanical models, and assessing to what extent the models are representative of the observations. On the other hand, if our models do not explain the observations, new and improved models must be proposed and developed. For example, in Chapter 6 we modeled rupture initiation as a line source. Our results indicate that rupture initiation has a lower rupture velocity, which we interpret as a consequence of the complex environment in which the earthquake is initiated such as presence of fluids, variable rheology, and/or effective contact surface. A next step would be to incorporate such complexities into the source modeling, and assess the effect it has. To do so, several problems must be solved first. Carrying out more complex modeling (at higher frequencies for example) requires a better knowledge of the internal structure of the Earth (e.g., better velocity models). In addition, fault geometry also plays a role. Some authors have proposed that listric rather than planar geometries would be better suited to faults in central Italy (e.g., [START_REF] Porreca | Seismic reflection profiles and subsurface geology of the area interested by the 2016 2017 earthquake sequence (Central Italy)[END_REF]. The incorporation of these elements (among others) implies a high computational cost, so perhaps a first step would be to understand which factors (e.g., fluids, variable rheology, effective contact surface) are most relevant during rupture initiation, in order to prioritize them in the incorporation of modeling. An example of this direction is the work done by Catania & Segall (2019), who studied the effect of fault roughness on slip behavior, concluding that for their model he preseismic phase is characterized by feedback between creep and foreshocks. To understand which factors are more relevant in nature, more and new observations are also needed.

Bridging the gap between laboratory experiments, numerical modelling and nature

As we mention in Chapter 6, laboratory experiments (e.g., [START_REF] Latour | Characterization of nucleation during laboratory earthquakes[END_REF][START_REF] Mclaskey | Earthquake initiation from laboratory observations and implications for foreshocks[END_REF][START_REF] Ohnaka | Scaling of the shear rupture process from nucleation to dynamic propagation: Implications of geometric irregularity of the rupturing surfaces[END_REF] and numerical models (e.g., [START_REF] Ampuero | Earthquake nucleation on rate and state faults Aging and slip laws[END_REF]; that earthquakes are preceded by different phases: a stable quasi-static deformation phase, which evolves into an unstable acceleration phase, after which the large dynamic rupture occurs once a critical rupture length (Lc) has been reached. In addition, another type of experiment focused on the evolution of the contact surface before rupture (e.g., [START_REF] Gvirtzman | Nucleation fronts ignite the interface rupture that initiates frictional motion[END_REF][START_REF] Rubinstein | Detachment fronts and the onset of dynamic friction[END_REF] have shown that, slow fronts preceding the large rupture gradually decrease the effective contact surface, leading to a complete rupture once the contact decreases sufficiently to no longer support the stress. However, none of the above processes have been observed in nature yet. Two initial questions would be, does the same processes occur in nature? And if so, what capabilities do we need to observe them? To this end, using observations of the early onset of the rupture can bring us closer to how it all starts (e.g., [START_REF] Colombelli | Evidence for a difference in rupture initiation between small and large earthquakes[END_REF]Ide, 2019). Methodologies such as the one presented in Chapter 6, which assume a line source, although simplified, allow to retrieve information as long as a good azimuthal coverage and high sampling rate data are available. In addition, they directly provide a direct comparison of rupture lengths with theoretical parameters such as Lc. Therefore, considering these types of observations in well instrumented regions or in the planning of large-scale experiments which goal is to study a natural fault system in-situ such as the FEAR project (http://fear-earthquake-research.org/home/, last accessed November 15, 2022), would perhaps allow reaching previously unobserved resolutions.

Towards an automated and data-driven way to classify seismicity

In Chapter 7 we showed the potential for seismicity classification based on unsupervised learning (clustering). We obtained different seismicity clusters, characterized by different physical properties associated with the features used during clustering. Although the use of physical parameters as features facilitates the interpretation of the different clusters, this becomes more challenging as the number of clusters increases. In our first example with two clusters (A and B, Fig. 7.4), it is clear that cluster A corresponds to more interacting and shortduration seismicity associated with the fault system while cluster B is more similar to background seismicity. In a second example, we show that each cluster (A and B) can be decomposed into three smaller clusters (Fig. 7.5), which differ from each other by specific characteristics e.g., Cluster A1 has the highest (Fig. 7.5b) and is very spatially located in the south-east of the study region (Fig. 7.6). We interpret the A1 features to correspond to swarm-like sequences, given the higher values usually observed for swarms (e.g., Chen et al., 2012). This combination of "more specific" features requires a detailed knowledge of how the characteristics of every seismic sequence manifest themselves in the parameters with which we can quantify seismicity. Therefore, as future work, it would be interesting to explore which are the most relevant characteristics (features) to differentiate and extract different types of seismicity or sequences in order to study them separately and better understand, for example, their role and dynamics within a complex fault system (e.g., [START_REF] Collettini | Rock and fault rheology explain differences between on fault and distributed seismicity[END_REF] differentiated distributed microseismicity from on-fault seismicity using their frequencymagnitude distribution). One open question in seismology, for example is, are there any parameters that stand out for foreshocks? For this, understanding in detail "what makes a sequence special" is fundamental. Results such as those presented in Chapters 3, 4 and 5 are then essential ingredients for the algorithms presented in Chapter 7 and others, in order to make progress in this area.

Introduction

Seismicity within the subducting oceanic Nazca plate in the central Andes occurs at a wide range of depths and magnitudes. Intraslab earthquakes in this region can be as shallow as ~40 km depth, defining a deeper plane of seismicity aligned parallel to the plate interface in northern Chile (Bloch et al. 2014;[START_REF] Sippl | Seismicity structure of the northern Chile forearc from> 100,000 relocated hypocenters[END_REF]. At depths greater than 60 km, the lack of coupling on the plate interface results in a considerable decrease of thrust earthquakes, and only intraslab earthquakes occur, defining a double seismic layer within the Nazca plate that extends to ~140 km depth [START_REF] Comte | A seismic zone in Arica, northern Chile[END_REF][START_REF] Dorbath | Double seismic zone of the Nazca plate in northern Chile: velocity structure, petrological implications, and thermomechanical modeling[END_REF][START_REF] Sippl | Seismicity structure of the northern Chile forearc from> 100,000 relocated hypocenters[END_REF][START_REF] Florez | Controlling factors of seismicity and geometry in double seismic zones[END_REF]Lu et al. 2021). Beyond those depths, intraslab earthquakes are less frequent and more pervasively distributed within the subducting plate. Most of the recorded > 7.0 intraslab events in this subduction zone have been deep focus earthquakes within the 550 km 1994), the 1994 Bolivia earthquake (Kikuchi & Kanamori 1994), and the 2015 earthquake doublets in the Peru-Brazil border [START_REF] Ruiz | Nucleation phase and dynamic inversion of the Mw 6.9 Valparaíso 2017 earthquake in Central Chile[END_REF]. The shallower section of the Nazca plate in the central Andes has also ruptured with large intraslab earthquakes, such as the 1950 8.0 Antofagasta and 2005 7.8 Tarapacá earthquakes [START_REF] Kausel | The Ms= 8 tensional earthquake of 9 December 1950 of northern Chile and its relation to the seismic potential of the region[END_REF][START_REF] Peyrat | Tarapacá earthquake (Mw 7.7, 2005, northern Chile): A event with horizontal fault plane constrained from seismologic and geodetic observations[END_REF]. Additionally, starting in 2007 and within a period of eight years, two > 6.5 intraslab 2011; Herrera et al. 2017).

In this work, we study the rupture properties of a third intraslab event that occurred 6.8 Calama intraslab earthquake (Figure 1). Considering the peculiar spatiotemporal distribution of these three major intraslab earthquakes, we compare their mainshock properties and aftershock sequences, discussing them within the the central Andes. Our aim is to evaluate how these factors could control the mainshock and aftershock characteristics of these events. The Calama mainshock occurred within the subducting Nazca plate at 123 km depth on city of Calama in northern Chile, as reported by the Centro Sismológico Nacional (CSN) of the Universidad de Chile. The focal mechanism solution reported by the Global Centroid Moment Tensor (GCMT) catalog (Dziewonski et al. 1981;Ekström et al. 2012) shows that the rupture occurred on a normal fault (see Figure 1).

Several local seismic stations were operational at the time of the Calama earthquake. To carry out all the analyses shown in this work, we used strong motion and broadband waveforms from multiparametric stations of the Integrated Plate boundary Observatory Chile network (IPOC) (GFZ & CNRS-INSU 2006) and the CSN Network (Barrientos & National Seismological Center (CSN) Team 2018). Strong motion waveforms from the network of earthquake-triggered accelerometers of the CSN (Barrientos & National Seismological Center (CSN) Team 2018) were also used.

Earthquake detection and location

Earthquake detection using template matching

We used template matching [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF] to detect unreported earthquakes around the Calama mainshock. This was done by analyzing continuous broadband velocity waveforms of nine stations near the epicenter from the IPOC and CSN networks (Figure S1a). We used the three components of these stations and bandpass filtered the data from 5 to 30 Hz, because this frequency range exhibits better signal-to-noise ratios (Cabrera et al., 2021). The template events are earthquakes reported by the CSN that occurred within a defined space-time window around the Calama earthquake. When defining a space-time window, a large window might allow the inclusion of additional events, but also more background seismicity that may not be related to the target sequence. By contrast. a smaller window mitigates this effect, but it is more susceptible to miss some events (e.g., [START_REF] Dascher Cousineau | What controls variations in aftershock productivity?[END_REF]Cabrera et al. 2021). To determine the size of the region enclosing the seismicity of the Calama sequence, we followed the expression proposed by [START_REF] Dascher Cousineau | What controls variations in aftershock productivity?[END_REF] based on the source radius estimated by Wells & Coppersmith (1994), resulting in a radius of 21 km around the hypocenter. In terms of time, we scanned the waveforms from one month before to one month after the mainshock (between May 3, 2020, and July 3, 2020), since this is the maximum number of days for which the nine stations were operating continuously. This space-time window comprises the mainshock and other 25 earthquakes that occurred after (Dataset S1 in the Supplementary Material). The waveforms of each template event were extracted by cutting the continuous data 0.5 s before the P-wave arrival and 5 s after the S-wave arrival. Wave arrivals were estimated using a local 1-D velocity model [START_REF] Husen | Accurate hypocenter determination in the shallow part of the Nazca subduction zone in Northern Chile using a combined on-/offshore network[END_REF]. The length of templates was defined in this way due to the difficulty of estimating P-wave arrivals accurately, given the limitations of the 1-D velocity model (e.g., Frank et al. 2017;Cabrera et al. 2021). To avoid detection of distant events not related to the studied sequence, correlation coefficients between the template waveforms and the continuous data were calculated within a sliding window that preserves the seismic moveouts using the Fast Matched Algorithm [START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF]) and a GPU-architecture. This resulted in time series that represent the similarity of the continuous data with every single template. We used a daily detection threshold that is 12 times the median absolute deviation (MAD) of the correlation function, which was averaged over all stations and channels to define the detection of an earthquake significantly similar to the template. The events detected with this criterion are assumed to occur at the same hypocentral location as their template (determined by the CSN). Their magnitudes were estimated by computing the median amplitude ratio between the template event and the aftershock over the considered stations, assuming that a tenfold increase in amplitude corresponds to one unit increase in magnitude [START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF]. The resulting earthquake dataset of the Calama sequence now includes 108 events in the magnitude range of S1b shows the comparison of the frequencymagnitude diagrams between the initial catalog and the new catalog. A higher number of event detections is observed for M < ~3.5, which is the completeness magnitude of the CSN catalog (Barrientos & National Seismological Center (CSN) Team 2018). Figure S1c summarizes the normalized waveforms of all the events in the new catalog recorded at station AF01, which is the closest to the epicenter (see Figure S1a). No earthquakes were detected before the mainshock. The new catalog of the Calama sequence features only aftershocks (see Figure S1d).

Location of the mainshock and aftershocks

To obtain a better resolution of the possible fault plane, we carried out a location of the mainshock and its aftershocks from the dataset of 108 earthquakes that resulted from the template matching. The location was performed using the same stations that were used for template matching (Figure S1a). First, the arrival times of the P and S waves were manually picked using the SEISAN software (Havskov & Ottemöller 1999). Once the arrival times were determined, the location was performed using the LocIn software (Potin 2016) on a regional 3-D tomographic velocity model (Pastén-Araya et al. 2021;[START_REF] Cabrera | Northern Chile intermediate-depth earthquakes controlled by plate hydration[END_REF]) (Figure 2). Out of the 108 earthquakes, only the mainshock and 37 aftershocks could be reliably located (Dataset S3) due to high noise level in the waveforms and limitations on station coverage. Location results indicate that the hypocenter of the Calama mainshock occurred at 113 km depth. Aftershocks were located mostly updip from the hypocenter, between 100 km and 113 km deep, defining a subvertical rupture plane, consistent with the NE dipping fault plane 1 and 2). The same location method was applied for both the 2007 6.8 Michilla and 2015 6.7 Jujuy sequences, whose mainshock depths were previously reported at 43 km and 254 km, respectively (Ruiz & Madariaga 2011;Pastén-Araya et al. 2018;Herrera et al. 2017). Our location results show mainshock hypocentral depths of 43 km for Michilla and 228 km for Jujuy. Compared with the Michilla earthquake, location uncertainties are larger for the Calama and Jujuy events, since they occurred at greater depths and were located with a smaller number of available stations, with important azimuthal gaps (Table 1 and Figure S2a). The located aftershock sequences of the Michilla and Jujuy earthquakes exhibit nearly vertical spatial distributions, closely aligned with the orientations of the steeper east-dipping fault planes of their respective focal mechanisms (Figures 1 and2). These results are consistent with the previously reported aftershock distributions and selected fault planes for these two earthquakes. Additionally, we carried out a relocation of these events using a double-difference method. Although double-difference relocations tend to be slightly deeper (< 2 km) and slightly more clustered (Figure S2), results are similar to the trends obtained with the absolute location approach.

The Calama mainshock

Strong motion data were used to analyze both the ground shaking characteristics of the Calama mainshock and its rupture properties. This earthquake generated a maximum peak ground acceleration (PGA) of 0.13 g at the closest station (hypocentral distance of 132 km). In general, the observed ground shaking intensities are within the ranges predicted by current Due to the limitations of the 1-D velocity model, waveforms from a subset of eight stations around the epicenter were used for modeling (stations shown in Figure 3a). The Neighborhood Algorithm converged to a best dynamic rupture model that has a maximum coseismic slip of 1.59 m. The two axes of this elliptical model are 14.1 km and 24.4 km long (Figure 3a), with a rupture time of 5.6 s. Dynamically, the overall rupture had a MPa and MPa, nucleating within a circle of km with MPa inside. A distance m was required to nucleate the rupture. The model parameters started to converge towards these optimal values roughly after 10,000 sampled models (Figure S4). Figure S5 shows the distributions and optimal values of these parameters. Some model parameters (e.g., , , , , and ) are less Gaussian distributed than others, which could indicate trade-offs between them. In particular, the correlation is stronger between the stresses (Figure S6), since in the model formulation, and depend on . also has a tradeoff with in the friction formulation [START_REF] Madariaga | Earthquake dynamics on circular faults: A review 1970 2015[END_REF], see also Figure S6), which is a contributing factor for the observed trend of in Figure S4, in addition to the resolution limitations of the parameter search grid. If the full seismograms are considered, the overall misfit associated to the best dynamic model is 0.58. In this case, the high-amplitude SV waves in the radial and vertical components could not be properly modeled in some stations (Figure 3b), resulting in this large misfit. This is likely due to converted body and surface waves arriving behind the SV waves, which could be generated by structures that are not represented by a simple 1-D velocity model. A similar case was shown by Herrera et al. (2017) for the Jujuy earthquake that occurred further east. Following their formulation, if the misfit is calculated using only P and high-amplitude SH waves (highlighted seismogram sections in Figure 3b), its value is reduced to 0.24. This is the misfit formulation that was minimized in the inversion to obtain the described best dynamic model of the Calama earthquake.

Discussion

Comparing dynamic properties of mainshock ruptures

The Calama mainshock was modeled using a finite-fault model, where the rupture propagation is controlled by a slip-weakening friction law. The other two mainshocks at Michilla and Jujuy were previously modeled using the same dynamic rupture model and inversion method (Ruiz & Madariaga 2011;Herrera et al. 2017). This allows a comparison of the inferred dynamic parameters with no bias related to differences in methods. The dynamic rupture parameters are summarized in Table 2, including the similarity parameter (Madariaga & Olsen 2000), calculated assuming the characteristic rupture size as the average of the ellipse semi-axes. All dynamic parameters of these three earthquakes are rather similar. In particular, the stress drop does not seem to be correlated with depth, which has also been observed with global earthquake databases (Poli & Prieto 2016). Overall, the values of these three events fall within the empirically estimated ranges for intraslab earthquakes globally (e.g., [START_REF] Kanamori | Theoretical basis of some empirical relations in seismology[END_REF]Poli & Prieto 2016), and they are larger than the values of thrust earthquakes inferred with the same method in northern Chile (Otarola et al. 2021).

The Calama earthquake occurrence within the upper lithospheric mantle

Intraslab earthquakes mostly occur in double seismic zones (DSZ), which has been observed in different subduction zones (Brudzinski et al. 2007). This DSZ is characterized by an upper seismicity plane (USP) located in the oceanic crust and a lower seismicity plane (LSP) located 20-40 km below the USP in the upper lithospheric mantle. The subduction zone of northern Chile is not an exception, and this DSZ has also been recognized in that region [START_REF] Comte | A seismic zone in Arica, northern Chile[END_REF][START_REF] Rietbrock | A narrowly spaced zone in the subducting Nazca plate[END_REF], Bloch et al. 2014;[START_REF] Sippl | Seismicity structure of the northern Chile forearc from> 100,000 relocated hypocenters[END_REF][START_REF] Florez | Controlling factors of seismicity and geometry in double seismic zones[END_REF]Lu et al. 2021). For example, the mainshock and aftershocks of the Calama sequence located in this work are concentrated between 100 km and 113 km depth, indicating that the rupture occurred in the LSP within the oceanic lithospheric mantle, below the oceanic Moho (Figure 4). In contrast, the location results of the Michilla sequence indicate that its rupture occurred in the USP within the oceanic crust (Figure 4).

Different mechanisms have been proposed for the generation of intraslab seismicity (e.g., [START_REF] Frohlich | Deep earthquakes[END_REF][START_REF] Houston | Deep earthquakes[END_REF]. For the seismicity located in the USP within the oceanic crust, there is a consensus that it could be related to the presence of fluids linked to dehydration processes within the oceanic crust at different pressures and temperatures (e.g., [START_REF] Kirby | Interslab earthquakes and phase changes in subducting lithosphere[END_REF][START_REF] Hacker | Subduction factory 2. Are earthquakes in subducting slabs linked to metamorphic dehydration reactions?[END_REF]. Dehydration might cause the reduction of the effective normal stress, promoting brittle rupture of structures inherited from the faulting process in the outer-rise zone prior to subduction [START_REF] Ranero | Relationship between at trenches and seismicity[END_REF]Ruiz & Contreras-Reyes 2015;Pastén-Araya et al. 2018;Cabrera et al. 2021). However, the mechanisms that generate intraslab seismicity in the LSP are still a subject of debate (Duesterhoeft et al. 2014;[START_REF] Ferrand | Dehydration-driven stress transfer triggers intermediatedepth earthquakes[END_REF][START_REF] Ohuchi | Intermediate-depth earthquakes linked to localized heating in dunite and harzburgite[END_REF]Scambelluri et al. 2017). Mechanisms that point to a hydrated lithospheric oceanic mantle have been proposed to trigger seismicity in the LSP (Bloch et al. 2018;Cai et al. 2018). On the other hand, analysis of laboratory and field data suggests that faulting could be triggered in dry rocks within a partially hydrated oceanic lithospheric mantle [START_REF] Ferrand | Dehydration-driven stress transfer triggers intermediatedepth earthquakes[END_REF]Kita & Ferrand 2018). This process has been referred to as dehydration-driven stress transfer, which would not require the presence of a highly hydrated lithospheric mantle. Instead, a rupture could nucleate in a weakly hydrated portion of the lithosphere and propagate to dry regions of the lithosphere due to the stress transfer associated with volumetric change of the rock. Additionally, [START_REF] Florez | Controlling factors of seismicity and geometry in double seismic zones[END_REF] found that globally, LSP seismicity has consistently smaller b-values compared with the USP seismicity, which would also indicate a relatively dry environment in the LSP. According to hydrological and numerical models, dehydration of the subducted slab occurs mainly in three stages (Ulmer & Trommsdorff 1995;[START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle?[END_REF][START_REF] Hacker | Subduction factory 2. Are earthquakes in subducting slabs linked to metamorphic dehydration reactions?[END_REF]Rüpke et al. 2004). First, dewatering of subducting sediments leads to hydration of the mantle wedge at depths < 20 km (ANCORP Working Group 1999, Rüpke et al. 2004). Second, metamorphic dehydration reactions of the subducting oceanic crust increase pore pressure and decrease effective confining pressure, thereby promoting depth) [START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle?[END_REF][START_REF] Hacker | Subduction factory 2. Are earthquakes in subducting slabs linked to metamorphic dehydration reactions?[END_REF]. Third, at depths larger than 100 km, the subducting lithospheric mantle dehydrates (Rüpke et al. 2004) and triggers intraslab seismicity (Yuan et al. 2000;[START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle?[END_REF]) causing partial melting and leading to arc volcanism (Rüpke et al. 2004;[START_REF] Cabrera | Northern Chile intermediate-depth earthquakes controlled by plate hydration[END_REF]. In our study case, dehydration reactions of the upper lithospheric mantle are consistent with a zone of intermediate ratios (1 region of the Calama earthquake (Figure 2b mantle values of ~8.3 km/s (Figure 2a) at > 600 MPa, suggesting the presence of dry dunite/peridotite mantle rocks (Christensen 1996). In addition, the mantle wedge presents large indicates the occurrence of massive dehydration reactions from the subducting oceanic lithosphere (e.g., Rüpke et al., 2004).

In summary, our results indicate that the Calama earthquake is likely a good example of an event triggered by the dehydration-driven stress transfer mechanism in dryer conditions. By contrast, the Michilla earthquake occurred within the oceanic crust where (Figure 2b), suggesting that the presence of fluids and a reduction of the effective normal stress could favor earthquake occurrence. The oceanic crust cannot be resolved in the region of the Jujuy earthquake. Moreover, the tomographic model cannot resolve properly beyond 2b). Therefore, considering this and the location uncertainties of the Jujuy earthquake (see Table 1 andFigure S2), for now the available data shows that this event occurred somewhere within the uppermost oceanic lithosphere, likely at lithostatic pressures of about 7 GPa 4). At these P-T conditions, the uppermost oceanic/subducting lithosphere dehydrates, favoring brittle faulting (Rüpke et al., 2004).

Thermal conditions and aftershock rate

Several studies have suggested that temperature could be an important factor that controls the distribution of both thrust and intraslab seismicity [START_REF] Oleskevich | The updip and downdip limits to great subduction earthquakes: Thermal and structural models of Cascadia, south Alaska, SW Japan, and Chile[END_REF]Wang et al. 2015;Wei et al. 2017;Liu et al. 2021). To try to establish the degree of influence of temperature on the Calama sequence, as well as on the other two intraslab earthquakes, we used the thermal model of northern Chile proposed by Cabrera et al. (2021), which is well defined between the trench and the volcanic arc in the upper ~200 km (Figure 4). Clear common trends are observed between the isotherms and the seismicity distribution. While the thrust seismicity is concentrated along the 200°C isotherm, the intraslab seismicity defined by the DZS is mostly concentrated between the 300°C and 400°C isotherms. A decrease in the seismicity is observed at higher temperatures, which could indicate a transition from brittle to particularly in the zones of the Michilla and Calama earthquakes. The brittle/ductile transition 4). Cabrera et al. (2021) studied intermediate-depth seismicity in northern Chile between neutral surface and brittle/ductile transition zone becomes deeper within the subducting plate to be consistent with these findings (Figure 4).

The Calama sequence occurred between the 300°C and 400°C isotherms (Figure 4), and its aftershocks mostly occurred at shallower depths than the mainshock. Similar trends were observed for the aftershock distributions of the 2019 6.7 Coquimbo and 2018 7.1 Anchorage intraslab earthquakes, which also exhibited shallower aftershocks than the mainshock [START_REF] Ruiz | The January 2019 (M w 6.7) Coquimbo Earthquake: Insights from a Seismic Sequence within the Nazca Plate[END_REF]Ruppert et al. 2020;Liu et al. 2019). In particular, the Coquimbo mainshock occurred between the 600°C and 700°C isotherms within the subducting plate [START_REF] Ruiz | The January 2019 (M w 6.7) Coquimbo Earthquake: Insights from a Seismic Sequence within the Nazca Plate[END_REF]. However, its aftershocks mostly occurred at shallower (and colder) layers, at temperatures below 450°C. These examples indicate that temperature could play a significant role in the aftershock distribution of intermediate-depth intraslab earthquakes, which tend to occur in layers of lower temperatures.

The aftershock productivity of intraslab earthquakes is another aspect that is related to both the zone where they are triggered and the temperature. Cabrera et al. (2021) carried out an analysis of several intraslab earthquakes in northern Chile, finding that intraslab earthquakes that occur at greater depths below the 40 aftershocks, and would be associated with a dry environment. Conversely, those events that aftershocks and would be associated with a more hydrated environment. Our results corroborate this observation, particularly when comparing the cases of the 2007 Michilla and the 2020 Calama earthquakes, which occurred at depths where the thermal model is still well defined. The Michilla earthquake occurred within the oceanic crust between the 200°C to 300°C isotherms (Figure 4), producing a large number of aftershocks and a zone with persistent seismicity in time (Ruiz & Madariaga 2011;Fuenzalida et al. 2013;Pastén-Araya et al. 2018). Conversely, the Calama mainshock and its aftershocks occurred in the upper lithospheric mantle between the 300°C and 400°C isotherms. Within the first five days after the mainshock, the Calama earthquake produced a much smaller number o compared with the Michilla earthquake (313 events). Therefore, these observations, in combination with the observed differences of ratios between the Calama and Michilla earthquakes, suggest that the Calama earthquake occurred in a warmer and less hydrated environment than the region of the Michilla earthquake, and would be responsible for its lower aftershock productivity. This is consistent with observations obtained by Chu & Beroza (2022) in the subducting Pacific Plate in Japan. They found that the aftershock productivity is correlated with ratio, discussing that a high ratio can be a result of high fluid pressure and a larger number of faults and cracks that could be fluid-filled, or also oriented perpendicular to ray paths.

Conclusions

The 2020 6.8 Calama earthquake is an intraslab earthquake that occurred at in 6.8 Michilla and 2015 6.7 Jujuy intraslab events. It featured ground shaking intensities that are typical of Chilean intraslab earthquakes.

The hypocenter of the Calama earthquake was located at 113 km depth using a 3-D model. The same method was used to locate the hypocenters of the Michilla and Jujuy earthquakes, resulting in depths of 43 km and 228 km, respectively. At their located depths, we observed that the Michilla earthquake occurred within the oceanic crust, while the Calama earthquake occurred within the upper lithospheric mantle, below the oceanic crust. The resolution of our database does not allow exact interpretations of the Jujuy earthquake location within the uppermost oceanic lithosphere due to the larger uncertainties in earthquake, slab, and oceanic Moho locations at those depths.

The dynamic properties of the Calama earthquake were inferred through modeling of low-frequency waveforms, which is the same method that was previously used to model the Michilla and Jujuy earthquakes. Despite their different hypocentral depths and locations in different layers of the subducting oceanic plate, the dynamic properties of these three events are similar. Particularly, their stress drop values range between 10 MPa and 15 MPa, within the observed ranges of intraslab earthquakes, which are in general larger than stress drop values of thrust earthquakes.

Thermal and pressure conditions of the subducting plate likely control the spatial distri occurrence down to ~150 km depth. Additionally, the varying water content and thermal conditions of mantle rocks in the areas where intraslab earthquakes occur play an important role in their aftershock productivity. For instance, the Michilla earthquake occurred within the oceanic crust at temperatures between 200°C and 300°C, exhibiting a strong aftershock activity. The large ratio (> 1.8) at that location indicates a more hydrated environment that favors brittle rupture and an increase in aftershocks. On the other hand, the Calama earthquake occurred in the uppermost lithospheric mantle, where the ratio is smaller (between 1.72 and 1.76), and temperatures vary between 300°C and 400°C. This earthquake exhibited a smaller aftershock productivity, which is likely a result of a less hydrated environment, as suggested by the reduced ratios in this region. Our results show that even though the Michilla and Calama earthquakes occurred in regions of different thermal and compositional characteristics within the Nazca plate, curiously these factors do not significantly affect the dynamic characteristics of the mainshocks, which were found to be within the typical ranges of intraslab events. However, they do affect their aftershock productivity. Additional studies with a larger database of well-recorded earthquakes are necessary to confirm if this trend is observed in more events.

Text S1: Strong motion analysis

The processing of the acceleration waveforms included the removal of the instrument response, mean, and linear trend. Then, a fourth-order Butterworth bandpass filter between 0.1 and 35 Hz was applied. The peak ground acceleration (PGA) at each station was calculated from the geometric average of the maximum accelerations of the two horizontal components. Station AF01 (the closest to the epicenter) recorded the highest PGA of this event, which reached 0.13 g. Ground motion intensities decrease with increasing distance from the earthquake (Figure S3a). Additionally, we calculated spectral accelerations as a function of period (SA(T)) using the geometric average of the 5% damped response spectrum (Nigam & Jennings 1969) of the two horizontal components. PGA and SA(T) observations were compared with the predictions for intraslab earthquakes of two recent ground motion models (GMM) developed with Chilean data (Idini et al. 2017;Montalva et al. 2017). The site parameter required by both GMMs is VS30, which was obtained from the site database compiled by Herrera et al. (2020). Normalized total residuals were calculated for each station j:

where and are the observed and predicted ground motion intensities at station j for period T, respectively, and is the total standard deviation of the GMM for period T, usually provided in units. Residual results are shown in Figure S3b, where in general both GMMs perform well when predicting ground motion intensities at the selected periods, especially at distances within their calibration range, as shown by the nearly zero residuals. As expected, slightly larger residuals are observed at distances greater than this limit, but no systematic deviation from the zero trend as a function of distance is shown by any residual distribution. [START_REF] Kohlstedt | Strength of the lithosphere: Constraints imposed by laboratory experiments[END_REF], and computed considering steady state flow properties for olivine (green line; [START_REF] Goetze | The mechanisms of creep in olivine[END_REF], dry olivine (yellow line; [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF] and wet olivine (cyan line; [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF]. Isotherms were computed based on the cooling of a semi-infinite half-space mode (Turcotte & Schubert, 2002). [START_REF] Hyndman | Thermal constraints on the zone of major thrust earthquake failure: The Cascadia subduction zone[END_REF][START_REF] Wang | Thermal regime of the Southwest Japan subduction zone: effects of age history of the subducting plate[END_REF][START_REF] Oleskevich | The updip and downdip limits to great subduction earthquakes: Thermal and structural models of Cascadia, south Alaska, SW Japan, and Chile[END_REF]Vöelker et al., 2011[START_REF] Höink | Viscous coupling at the asthenosphere boundary[END_REF]. Text A3.S1. Detrend strategy for GPS velocities.

Parameter

Our strategy includes the next steps:

-We estimate a velocities field for several IGS stations around the Antarctic plate, from a period of time up to end 2020 (see examples in Fig. S6-S13). Then, we estimate the Euler pole and rate of the Antarctic plate (Table S1); -We use our Euler estimation to produce time-series with respect to Antarctic fix frame; -Then we detrend the time-series estimating residual velocities, heavy side and seasonal, which one are removed to get the residual deformation. Finally, we show detrended time-series in Figure 4 and Figure S5 (FREI station included). Text A4.S1. Estimating the magnitude of completeness using the Lillefors test.

We estimate the magnitude of completeness using the Lillefors test implemented by Herrmann est Mc for two significance level of and , obtaining Mc=0.8 and Mc=0.9, respectively. As indicated by Herrmann and Marzochi (2020), choosing is conservative in a statistical sense [START_REF] Clauset | Power-law distributions in empirical data[END_REF]. We therefore prefer the latter, more conservative value for the magnitude of completeness to show the stability of our analysis (see Fig. S8). In addition, the Lilliefors pvalues indicates that this catalog has different exponential distributions, in concordance with (Herrmann and Marzochi, 2020). This is likely due to the variability of the b-values in the al. ( 2012) and [START_REF] Sugan | The preparatory phase of the 2009 Mw 6.3 L'Aquila earthquake by improving the detection capability of low magnitude foreshocks[END_REF].

Fig. 1 . 1 6 Fig. 2 .

 1162 Fig. 1.1 Schematic diagram illustrating three types of earthquake sequences 2 Fig. 1.2 Example of a Mainshock-Aftershocks sequence 2 Fig. 1.3 Example of a Foreshocks-Mainshock-Aftershocks sequence 3 Fig. 1.4 Example of a seismic Swarm in the Corinth Gulf, Greece 5 Fig. 1.5 Example of the improved catalog of foreshocks before the 2011 Mw 9.1 Tohoku earthquake 6

Figure 1 . 1 .

 11 Figure 1.1. Schematic diagram illustrating three types of earthquake sequences (MS for mainshock): (a) Mainshock-Aftershock; (b) Foreshocks-Mainshock-Aftershocks; (c) Swarm. Modified from Scholz (2019).

Figure 1 . 2 .

 12 Figure 1.2. Example of a Mainshock-Aftershocks sequence. Spatial overlap between the calculated Coulomb stress increment (red zones) and the location of aftershocks (black dots) for the 1994 Mw=6.7 Northridge earthquake. The largest Coulomb stress change on optimally oriented thrust or strike-slip faults at depths of 3-10 km; the compressive axis of the regional stress is oriented N4°E, and . Locations of active surface faults, and M shocks during 3-6 months after the mainshock, are superimposed in black. Modified from Stein (1999).

Figure 1 . 3 .

 13 Figure 1.3. Example of a Foreshocks-Mainshock-Aftershocks sequence. The mainshock corresponds to the 24 April 2017 Mw 6.9 Valparaiso earthquake, Central Chile. The Time latitude evolution of the seismicity detected from 22 April to 12 May 2017 is presented. Foreshocks correspond to the events in the time range indicated by re repeating earthquakes. Modified from Ruiz et al. (2017).

Figure 1 . 4 .

 14 Figure 1.4. Example of a seismic Swarm in the Corinth Gulf, Greece. (a) Map of the Corinth Gulf, with the main outcropping faults (black lines), the seismic network (blue triangles), and the seismicity (orange dots, sized by magnitudes) detected and located by the CRL network between the Days 270 (27 September 2015) and 280 (7 October 2015). The exemplified swarm is the dense cluster below MALA station. The star MG indicates the local array used for detection. The triangles are the stations used for location. (b) Time distribution of the detected and odified from De Barros et al. (2020).

Figure 1 . 5 .

 15 Figure 1.5. Example of the improved catalog of foreshocks before the 2011 Mw 9.1 Tohoku earthquake.Earthquake migration toward the rupture initiation point of the mainshock. Space-time diagram of all detected events between 13 February and the mainshock origin time, with earthquake origin locations indicated in terms of the distance along the trench axis (blue circles scaled to magnitude). Red dashed lines, approximate locations of the fronts of earthquake migration; red stars, repeating earthquakes in the JMA catalog; green stars, newly detected events that were found to resemble those repeating events; black star, Mw 9.0 mainshock; yellow star,

Figure 2 . 1 .

 21 Figure 2.1. Illustrative summary of the template matching scanning process. (a) Continuous data are obtained with the corresponding preprocessing.In this example, the continuous traces of the vertical component are shown graphically for nine stations (the idea is generalizable to as many stations and components as available). (b) A template event is selected, and the window length to be used to extract the P (blue) and S (red) wave template is defined. (c) The waveforms with their respective moveouts are extracted. (d), (e) and (f) show the process of correlating the template with the continuous data and then shifting the window. For cases (d) and (e) the is not sufficient to declare a detection, while the case of (f) is (see details on how to define the threshold in the text).

Figure 2 . 2 .

 22 Figure 2.2. Example of template matching search over a single day. (a) Map of the region (Northern Chile) the (from equation 2.1) over one day. The threshold (cyan line) in this example is defined as , PB16.HHZ and 9 examples of detections on this same channel.

Figure 3 . 1 .

 31 Figure 3.1. Seismotectonic setting of the northern Chile region. a Orange circles indicate the mainshocks (labeled as E1, E2 and so on according to their distance from the top of the slab) and aftershocks epicenters (CSN catalogue), scaled and colored by magnitude and depth. Focal mechanisms are obtained from USGS, and also date and magnitude according to CSN are indicated. Black inverted triangles and cyan squares are the locations of the stations used for the kinematic inversions (KI) and template matching (TM), respectively. Yellow triangles represent active volcanoes. Red lines indicate faulting located on the outer-rise region. Isochrones for the Nazca plate (Müller et al., 2008) are indicated with magenta lines. (Inset) Global map with the specified study area. b, c indicated for each mainshock. The slab (black solid line) corresponds to the projection along cross-sections of the Nazca slab model SLAB2.0 (Hayes et al., 2018). Black and green dashed lines correspond to continental Moho (Tassara & Echaurren, 2012) and oceanic Moho (Contreras-Reyes et al., 2012), respectively.

Figure 3 . 2 .

 32 Figure 3.2. Example of new events detected using template matching. a Map of stations used for E2 (cyan inverted triangles) and E5 (green dots). Stars represent epicenter for each event (CSN), and their focal mechanism according to USGS are also plotted. (Inset) Global map with the specified study area. b Example of a new detection in the continuous data for E2. Cyan traces represent templates aligned considering their respective moveouts. The average correlation coefficient (Mean CC) is 0.10 and the inferred magnitude (Mw) 1.8. c Example of a new detection in the continuous data for E5. Green traces represent templates aligned considering their respective moveouts. The average correlation coefficient (Mean CC) is 0.12 and the inferred magnitude (Mw) for the new event is 2.7.

Figure 3 . 4 .

 34 Figure 3.4. FEM grid and boundary conditions used for thermal numerical modelling. Continental crust, mantle wedge and oceanic plate are indicated with green, black, and light brown colors, respectively. Capital letters circled represent where border conditions are imposed (see details in Text A2.S1).

Figure 3 . 5 .

 35 Figure 3.5. Example of kinematic inversion for E1. Ellipse shows the rupture area of the mainshock, with at maximum slip amplitude of 1.08 m. Yellow diamond is the National Seismological Center (CSN) mainshock epicenter. Inverted red triangles are near-field stations considered for inversion. Gray and blue lines are simulated and observed displacement waveforms respectively. Focal mechanism from USGS.

Figure 3 . 6 .

 36 Figure 3.6. Yield strength envelope for an oceanic lithosphere 58 Myr old. Stress differences or yield stress are limited at the top of the plate by frictional sliding rule according to Byerlee's frictional sliding rule (black curve).Yield stress is limited at the base of the plate by steady state creep, which depends on the cube of the stress and exponentially on temperature[START_REF] Kohlstedt | Strength of the lithosphere: Constraints imposed by laboratory experiments[END_REF], and computed considering steady state flow properties for olivine (green line;[START_REF] Goetze | The mechanisms of creep in olivine[END_REF], dry olivine (yellow line;[START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF] and wet olivine (cyan line;[START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF]. Isotherms were computed based on the cooling of a semi-infinite half-space mode(Turcotte & Schubert, 2002).

Figure 3 . 7 .

 37 Figure 3.7. Cartoon showing a conceptual model for the the evolution of subducting lithosphere in northern Chile. The topography of the plate in the outer-rise/trench region has been exaggerated to show better the deformation associated to plate bending. Scale is approximate everywhere else. Bubbles in oceanic lithosphere indicate hydration.

Figure 4 . 1 .

 41 Figure 4.1. Seismotectonic setting of Antarctica. Orange inverted triangles indicate GNSS stations and the green inverted triangle shows the location of the JUBA broad-band station. Pink stars represent historical seismicity with magnitudes equal or larger than 5.5 from the ISC (International Seismological Centre, 2021) occurred after the eruption of Deception Island (12 August 1970). Red arrows represent velocities of the stations (see text for details about GPS processing). Right inset: Focal mechanisms (colored compressional quadrants) reported by the USGS since 28 August 2020 and the location of the Orca volcano. Date dd/mm and magnitude are also indicated. The cyan focal mechanism is for the largest event (Mw5.9) also indicated in Figures 4.4 and A3.S4. The estimated NW displacement for UYBA with respect to the Antarctic plate is indicated with a black arrow (see text for details about GNSS processing).

Figure 4 . 2 .

 42 Figure 4.2. Statistics of the seismicity. (a) Number of events as function time with respect to 8 August 2020, with the vertical red dashed line highlighting the occurrence of the largest event in the sequence. (b) Distribution of inter-event waiting time (dt) for the detected events. The red line shows the fit of this data with equation (4.1) and = 0.6. (c) Cumulative magnitude-frequency distribution for the detected earthquakes. The red dashed line shows the regression for magnitudes ranging from 3 to 6. The corresponding b-value is 1.6.

Figure 4

 4 Figure 4.3. S-P times of earthquakes. S-P time from manual picking performed on 114 templates (a) and from delay estimated from cross-correlation of earthquakes detected (b c).

Figure 4 . 4 .

 44 Figure 4.4. GPS displacement at stations reported in Figure 4.1 for the east (a) and north (b) components. (c and d) are zooms during the swarm time for station UYBA (see Figure 4.1). (e) Cumulative number of events as function of time in red, and recurrence time of events in black dots.

  and 4.4b). However, no GNSS time series show a clear increase of velocity during the seismicity rate acceleration between 7 and 28 August 2020(Figures 4.2a

  , 4.3c 3e). The most important velocity change with a clear centimetric displacement (Figures 4.1

  and

Figure

  Figure 5.1.Valoroso, et al. (2011); our 267 template events are drawn from this earthquake catalog and are shown in the zoomed map. The broadband stations we analyzed are shown by the red triangles. Black and cyan thin lines respectively represent traces of the active mapped faults and co-seismic surface ruptures, respectively(Boncio et

  and 5.3c) of sequences S1 and S2 exhibit different time evolutions. The seismicity during S1 is characterized by a slow time evolution, with several accelerations occurring over a few days (Figures 5.2, 5.3a and 5.3c) and without any clear mainshock driving them (Figure5.2b). On the other hand, the cumulative number of events in S2 evolves with a log-like behavior similar to an aftershock Omori law[START_REF] Utsu | The centenary of the Omori formula for a decay law of aftershock activity[END_REF].

Figure 5 . 2 .

 52 Figure 5.2. (a) Magnitude-frequency distribution (0.1 bin) for events detected. (b) Estimated magnitudes (see

Figure 5 . 3 .

 53 Figure 5.3. Catalog generated using template matching. (a) Interevent times are plotted using black circles, defined as the elapsed time between consecutive events. Red line represents the cumulative number of events, and magenta, cyan, and blue vertical lines show the time of F1, F2, and mainshock (MS) events, respectively. Inset: a zoom of the time interval between F1 and the mainshock. (b) Normalized waveforms of the catalog for the vertical component of the AQU station, aligned 0.5 s before estimated P-wave arrival (black vertical dashed line). Event ID is chronologically ordered (i.e., the vertical axis is time-ordered). Time of occurrence of F1 and F2 are also indicated with black horizontal lines. (c) Comparison between the normalized cumulative event count for S1 and S2.

Figure 5 . 4 .

 54 Figure 5.4. Temporal evolution with a sliding 100-event window length and a 99-event overlap for: (a) Coefficient of variation of the interevent times (b) Ratio between the maximum value of Mo and its total amount within the window (c) Average along-strike and along-dip location of the seismicity measured relative to the mainshock (MS) and projected on the main fault. Time corresponds to the time of the last event within the 100-event window (see text). Interevent times plotted in Figure 5.3a is also indicated in Figures a, b and c (gray dots). Magenta, cyan, and blue vertical dashed lines show the time of F1, F2, and MS events, respectively. Cyan and light orange dots represent zoom in presented in Figure 5.5.

Figure 5

 5 Figure 5.5. (a) Example of radial spatial migration (measured from the first event of the sequence) for a burst during S1 (see Figure5.4). Yellow star represents the largest event within the sequence (M2.4 according toChiaraluce, Valoroso, et al., 2011). Red dashed lines represent best-fit fluid diffusion curves[START_REF] Shapiro | Estimating the crust permeability from fluidinjection-induced seismic emission at the KTB site[END_REF] for hydraulic diffusivity of 1.5 m2/s (all the seismicity in a) and 0.8 m2/s (95-percentile of the seismicity in (a). (b) Example of radial spatial migration for S2 (see Figure5.4). Magenta, cyan, and blue vertical stars show the time of F1, F2, and mainshock events, respectively, and yellow star represent another event magnitude 3.2 in the middle of the sequence. (c) Zoom for the first 60 min plotted in (b).

Figure 5 . 6 .

 56 Figure 5.6. Along-strike distribution of the seismicity centered in the mainshock for (a) time plot and (b) eventindex plot. The groups of seismicity are shown in Figure5.5a are highlighted in cyan (part of S1) and light orange (S2). Magenta, cyan, and blue vertical lines with stars show the time of F1, F2, and mainshock (MS) events, respectively. Note that the strike considered is the one of the main faults N133°E[START_REF] Valoroso | Radiography of a normal fault system by 64,000 high precision earthquake locations: The 2009 L'Aquila (central Italy) case study[END_REF].

Figure 5 . 7 .

 57 Figure 5.7. Cumulated moment, radius, and effective stress drop evolution. We use 100-events windows-length and 99 events overlapping for (a) cumulated radius (black line) and cumulated moment (red line). (b) Effective stress drop. Time corresponds to the time of the last event within the 100-events window (see text). (c) Scaling between the cumulated radius and cumulated moment for the first part of the sequence (S1, red dots) and the second part (S2, blue dots).

Figure 6 . 1 .

 61 Figure 6.1. Seismotectonic setting. (a) Location map for the earthquake, showing the precursory seismicity (black dots) reported by Chiaraluce et al. (2011), and the strong motion stations (inverted triangles) used in this study. The colors of the triangles indicate the duration between the initial emergent P-waves (EP, red star) and the large impulsive onset (IP, blue star) P-waves waves located by Di Stefano et al. (2011). Black and green thin lines represent traces of the active mapped faults and co-seismic surface ruptures, respectively[START_REF] Boncio | Coseismic ground deformation of the 6 April 2009 L'Aquila earthquake (central Italy, Mw6. 3)[END_REF]. The beach-ball (compressional quadrants in colors) represents the source mechanism (reported by the Istituto Nazionale di Geofisica e Vulcanologia, INGV) for the mainshock and EP (yellow diamond, Table6.1). (b) Vertical component velocity seismograms[START_REF] Luzi | ORFEUS Working Group 5[END_REF] of the initial P-waves recorded for stations at different distances from the epicenter. Red and blue vertical lines indicate the arrival times of EP and IP, respectively.

Figure 6 .

 6 Figure 6.2b, c shows the misfit between the dip of the rupture initiation as a function of the strike of the rupture initiation and the length of the line source as a function of the rupture velocity. Clear minima are observed in both plots, as well as in those that relate the rest of the parameters (Fig. A5.S6), which illustrates the robustness of our estimates. In addition, Figure 6.2a shows the comparison of the best solution synthetics and the observed rupture times as a function of the azimuth. Remarkably, our solution can reproduce rupture times () for stations that have azimuths in common but different take-off angles (e.g., stations SUL-AQK and LSS-AQG-AQA-AQV) with a misfit of 0.015 s. This implies that we can reproduce the observed up-going and down-going P-waves, and that the estimate of the dip is therefore robust.

Figure 6 . 2 .

 62 Figure 6.2. Observations, best model and misfit. (a) Observed rupture times (colored according to the takeoff angle) as a function of the azimuth compared with the best solution synthetics (orange dots). (b) Misfit between the dip of the rupture initiation (i.e., 0°=vertically down, 90°=horizontal, and 180°=vertically up) as a function of the strike of the rupture initiation. The black star contour indicates the minimum misfit and the red star indicates the dip and an angle perpendicular to the strike of the fault system (50° and N47°E, respectively according to

Figure 6 . 3 .

 63 Figure 6.3. Cartoon summarizing our observations and interpretation. (a) The rupture initiates in a complex zone, which hinders its upward propagation and reduces its rupture velocity. This leads to low seismic efficiency, indicating that most of the energy budget corresponds to fracture energy. Some complexities that may influence the zone are the presence of fluids, a mixed rheology (ductile and brittle) and/or variations in the effective contact surface (see the text for more detail). Once the rupture has left this zone (to the faster region), these complexities no longer play an important role. (b) Seismogram of station AQA (Fig. 6.1) showing the record of the rupture initiation and the large dynamic rupture.

Figure 7 . 1 .

 71 Figure 7.1. Location map. Black dots represent the templates from INGV used for template matching. The broadband stations we scanned are shown by the red triangles. The epicenter of the larger mainshocks are represented by stars linked to their respective focal mechanisms (compressional quadrants in colors, reported by INGV). Red and yellow lines represent traces of the active mapped normal and thrust faults, respectively[START_REF] Gessel | The HIKE European Fault Database (EFDB) compiled in the framework of the GeoERA project HIKE[END_REF].

  to continuous seismic data collected by the Istituto Nazionale di Geofisica e Vulcanologia (INGV) from 1 January 2009 to 31 December 2016 (8 years). We use 19 broadband three-component stations (red triangles in Fig. 7.1) from the Italian Seismic Network (INGV Seismological Data Center, 2006) and the Mediterranean Very Broadband Seismographic Network (MedNet Project Partner

Figure 7 . 2 .

 72 Figure 7.2. (a) Magnitude-frequency distribution (0.1 bin) for the new events detected (grey) and templates (red).

Figure 7

 7 Figure 7.3. (a) Comparison between the cumulative number of events for the initial catalog (templates, red line) and the new catalog (grey line). Same colors are used to indicate the daily counts. (b) Seismicity along strike of and Norcia earthquakes are indicated by the blue, yellow, orange and cyan stars, respectively.

Figure 7 . 4 .

 74 Figure 7.4. Dendrogram after the clustering is performed (example 1). An Euclidean distance of 250 is defined as threshold in order to get two clusters (labeled A and B) (b) Spatial distribution of the clusters. (c) Average value of the features for every cluster. (d) Relative size of the main clusters compared to the size of the entire data set.

Figure 7 . 5 .

 75 Figure 7.5. Dendrogram after the clustering is performed (example 2). An Euclidean distance of 180 is defined as threshold in order to get six clusters (labeled A1, A2, A3, B1, B2 and B3) (b) Average value of the features for every cluster. (c) Relative size of the main clusters compared to the size of the entire data set.

Figure 7 . 6 .

 76 Figure 7.6. Spatial distribution of the clusters presented in Fig. 7.5.

Figure 7 . 7 .

 77 Figure 7.7. Along Strike distribution of the clusters presented in Fig. 7.5. The cumulative number of events considering the templates and their respective detections are also indicated with grey lines.

Figure 7 . 8 .

 78 Figure 7.8. Cross sections for the clusters presented in Fig. 7.5. The cross sections are spaced every 15 km, and contain all events within 7.5 km of distance from the section.

  Strike (km) Distance Perpendicular to Strike (km) Distance Perpendicular to Strike (km) Dist ance Along St rike (km )

Figure 1 .

 1 Figure 1. Seismological context of the Calama earthquake. Stars show the epicenters of the Michilla, Calama, and Jujuy earthquakes. Their focal mechanisms from the Global Centroid Moment Tensor (GCMT) catalog (Dziewonski et al. 1981; Ekström et al. 2012) are also shown. The black dots indicate background seismicity reported by the Centro Sismológico Nacional (CSN) and relocated by Pastén.Araya et al. (2018). CC: Coastal Cordillera, ID: Intermediate Depression, DC: Domeyko Cordillera, SA: Salar of Atacama, VA: Volcanic arc. The red triangles correspond to the main active volcanoes. Cross section A-

Figure 2 .

 2 Figure 2. Seismicity within the 3-D tomographic model. (a) and (b) tomography models (Pastén-Araya et al. 2021) along cross section A-1. Hypocenters of the Calama mainshock and its aftershocks are shown with a yellow star and red circles, respectively. Red stars indicate the hypocenters of the Michilla and Jujuy earthquakes, and their aftershocks are shown with blue and green circles, respectively. The continental Moho was inferred by Tassara & Echaurren (2012). The plate interface as defined by Hayes et al. (2018) is shown by the continuous black line. The oceanic Moho defining the low limit of the oceanic crust was inferred by Contreras-Reyes et al. (2021). The oceanic crust is not accurately resolved below certain depths (segmented line extensions). The upper lithospheric mantle was defined based on the lower plane of seismicity reported by[START_REF] Sippl | Seismicity structure of the northern Chile forearc from> 100,000 relocated hypocenters[END_REF]. Red triangles show the main active volcanoes. The discolored areas of the tomographic model are regions of lower resolution.

Figure 4 .

 4 Figure 4. Cross section A-A'. Symbols of the Michilla, Calama, and Jujuy earthquakes and their aftershocks are as described in Figure 2. The black dots indicate the background seismicity reported by the CSN and by Pastén.Araya et al. (2018). The continental Moho was inferred by Tassara & Echaurren (2012). The orange isotherms correspond to the thermal model of northern Chile (Cabrera et al. 2021). The slab geometry is the Slab2.0 (Hayes et al. 2018). The oceanic crust and isotherms are not accurately resolved below certain depths (segmented line extensions). CC: Coastal Cordillera, ID: Intermediate Depression, DC: Domeyko Cordillera,

Figure S1 :

 S1 Figure S1: Template matching analysis to detect earthquakes of the Calama seismic sequence. (a) Map showing the broadband stations used, the Calama earthquake with its GCMT focal mechanism, and the aftershocks reported by the CSN (red dots). (b) Frequency-magnitude diagram of the original CSN catalog and the new catalog with events detected through template matching (TM). (c) Normalized waveforms of the new catalog for the vertical component of station AF01, aligned 0.5 s before the estimated P-wave arrival (grey dashed line). Events are sorted based on their occurrence time. (d) Daily number of events before and after the mainshock.

Figure S2 :

 S2 Figure S2: Locations of the three earthquake sequences. Mainshock locations are shown with stars and aftershocks with colored circles. Upper and lower boundaries of the oceanic crust are also shown, which are less resolved at greater depths (segmented lines). (a) Absolute locations. Error bars are also shown for the mainshocks. (b) Relocations obtained with a double difference method.

Figure S3 :

 S3 Figure S3: Strong motion analysis of the Calama earthquake. (a) Spatial distribution of the observed PGA at the analyzed stations. The earthquake moment tensor was obtained from GCMT and is located at the epicenter. The trench line was obtained from Bird (2003). (b) Zt residuals for PGA and SA(T) at three different periods, which are shown in different colors for each GMM. Hypocentral distance (Rhypo) is used by both GMMs for intraslab earthquakes. The maximum calibration distance of the GMMs is shown by the red lines.

Figure S4 :

 S4 Figure S4: Convergence of the 10 parameters of the dynamic rupture model and . All sampled models are shown by dots colored according to their misfit. The gray dashed line in each plot defines the start of the range where model parameters start to converge.

Figure S5 :

 S5 Figure S5: Histograms of the 10 parameters of the dynamic rupture model and , calculated within the range where model parameters start to converge (as defined in Figure S4). Cyan curves show the best-fitting Gaussian distributions of mean and standard deviation . The values of the best model ( ) are shown with red lines.

Figure S6: 2

 2 Figure S6: 2-D distribution between parameters of the dynamic rupture model. All sampled models are shown by dots colored according to their misfit. The white star in each plot shows the values associated to the best model.

Figure A2 .

 A2 Figure A2.S1. Observed (blue line) and simulated (red line) traces for E1 using nodal plane 1.

Figure A2 .

 A2 Figure A2.S2. Observed (blue line) and simulated (red line) traces for E2 using nodal plane 1.

Figure A2 .

 A2 Figure A2.S3. Observed (blue line) and simulated (red line) traces for E3 using nodal plane 1.

Figure A2 .

 A2 Figure A2.S4. Observed (blue line) and simulated (red line) traces for E4 using nodal plane 1.

Figure A2 .

 A2 Figure A2.S5. Observed (blue line) and simulated (red line) traces for E5 using nodal plane 1.

Figure A2 .

 A2 Figure A2.S6. Observed (blue line) and simulated (red line) traces for E6 using nodal plane 1.

Figure A2 .

 A2 Figure A2.S7. Observed (blue line) and simulated (red line) traces for E1 using nodal plane 2.

Figure A2 .

 A2 Figure A2.S8. Observed (blue line) and simulated (red line) traces for E2 using nodal plane 2.

Figure A2 .

 A2 Figure A2.S9. Observed (blue line) and simulated (red line) traces for E3 using nodal plane 2.

Figure A2 .

 A2 Figure A2.S10. Observed (blue line) and simulated (red line) traces for E4 using nodal plane 2.

Figure A2 .

 A2 Figure A2.S11. Observed (blue line) and simulated (red line) traces for E5 using nodal plane 2.

Figure A2 .

 A2 Figure A2.S12. Observed (blue line) and simulated (red line) traces for E6 using nodal plane 2.

Figure A2 .

 A2 Figure A2.S13. Potential pre-existing fault (yellow line) considering E1 mainshock and its aftershocks distribution (orange circles) reported by CSN. Red lines indicate faulting located on the outer-rise region.

Figure A2 .

 A2 Figure A2.S14. Yield strength envelope for an oceanic lithosphere 58 Myr old. Stress differences or yield stress are limited at the top of the plate by frictional sliding rule according to Byerlee's frictional sliding rule (black curve). Yield stress is limited at the base of the plate by steady state creep, which depends on the cube of the stress and exponentially on temperature[START_REF] Kohlstedt | Strength of the lithosphere: Constraints imposed by laboratory experiments[END_REF], and computed considering steady state flow properties for olivine (green line;[START_REF] Goetze | The mechanisms of creep in olivine[END_REF], dry olivine (yellow line;[START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF] and wet olivine (cyan line;[START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF]. Isotherms were computed based on the cooling of a semi-infinite half-space mode(Turcotte & Schubert, 2002).

Figure A2 .

 A2 Figure A2.S15. Near-field broad-band stations used for template matching (TM). Small stars represent epicenter for each mainshock and using the same color their respective focal mechanisms (from USGS). In addition, stations used for template matching of each event are plotted using the same color as their respective epicenter and focal mechanism: E1 inverted red triangles, E2 black squares, E3 green triangles, E4 cyan inverted triangles, E5 orange stars and E6 yellow circles.

Figure A3 .

 A3 Figure A3.S1. Time evolution of the magnitude for the 114 events in the USGS catalogue (red) and new detected events in black.

Figure A3 .

 A3 Figure A3.S2.Example of detected waveforms with single station template matching for three components. The y-axis are the template samples, P waves arrived around 50 samples, S waves at ~100 sampes. X-axis is the number of detected events in chronological order. Color is proportional to the amplitude of the seismic waves, each seismogram (row) is normalized by its maximum amplitude.

Figure A3 .

 A3 Figure A3.S3. Distance respect to the first event of the sequence as function of time, for events of the USGS catalog.

Figure A3 .

 A3 Figure A3.S4. Stations (red triangles) used to estimate Euler pole. Velocities and residues (see TableS1for more details) are shown with green arrows and red ellipses, respectively.

  Figure A3.S4. Stations (red triangles) used to estimate Euler pole. Velocities and residues (see TableS1for more details) are shown with green arrows and red ellipses, respectively.

Figure A3 .

 A3 Figure A3.S5. Detrended GPS (Fig. 1) displacement records (FREI included). a) and b) are the east and north time series for the full recording time, showing no deformation up to end of August 2020. c) and d) are zooms during the swarm time for the UYBA station. e) Cumulative number of events as function of time in red, and recurrence time of events in black dots.

Figure A3 .

 A3 Figure A3.S6. Velocity fit for UYBA station. See Text A3.S1 for more details.

Figure A3 .

 A3 Figure A3.S7. Velocity fit for OHI2 station. See Text A3.S1 for more details.

Figure A3 .

 A3 Figure A3.S8. Velocity fit for MBIO station. See Text A3.S1 for more details.

Figure A3 .

 A3 Figure A3.S9. Velocity fit for SPRZ station. See Text A3.S1 for more details.

Figure A3 .

 A3 Figure A3.S10. Velocity fit for SPGT station. See Text A3.S1 for more details.

Figure A3 .

 A3 Figure A3.S11. Velocity fit for FREI station. See Text A3.S1 for more details.

Figure A3 .

 A3 Figure A3.S12. Velocity fit for PALM station. See text Text A3.S1 for more details.

Figure A3 .

 A3 Figure A3.S13. Velocity fit for DUPT station. See Text A3.S1 for more details.

Figure A4 .

 A4 Figure A4.S1. Example of daily (30/04/2009) spectrogram for the vertical component of the AQU station.

Figure A4 .

 A4 Figure A4.S2. Example of flipped template waveform over time for AQU station. Letter P indicates the P-wave arrival.

Figure A4 .

 A4 Figure A4.S8. Example of parameters using events with magnitude 0.9, 100-events windows-length and 99 events overlapping for: (A) Coefficient of variation and interevent times plotted in Fig. 2A (black dots). (B) Ratio between the maximum value of Mo and its total amount within the window and interevent times plotted in Fig. 2A (black dots). (C) Average along-strike and along-dip location of the seismicity measured from the MS and projected on the main fault and interevent times plotted in Fig. 2A (black dots). (D) Cumulated radius (black line) and cumulated moment (red line). (E) Effective stress drop. Time corresponds to the time of the last event within the 100-events window (see text). Fuchsia, turquoise and blue vertical lines show the time of F1, F2 and MS events, respectively.

Figure A4 .

 A4 Figure A4.S9. Example of parameters using 50-events windows-length and 49 events overlapping for: (A) Coefficient of variation and interevent times plotted in Fig. 2A (black dots). (B) Ratio between the maximum value of Mo and its total amount within the window and interevent times plotted in Fig. 2A (black dots). (C) Average along-strike and along-dip location of the seismicity measured from the MS and projected on the main fault and interevent times plotted in Fig. 2A (black dots). (D) Cumulated radius (black line) and cumulated moment (red line). (E) Effective stress drop. Time corresponds to the time of the last event within the 50-events window (see text). Fuchsia, turquoise and blue vertical lines show the time of F1, F2 and MS events, respectively. Turquoise and light orange dots represent events plotted in Fig. 2 (d) and (e), respectively.

Figure A4 .

 A4 Figure A4.S10. Example of parameters using 200-events windows-length and 199 events overlapping for: (A) Coefficient of variation and interevent times plotted in Fig. 2A (black dots). (B) Ratio between the maximum value of Mo and its total amount within the window and interevent times plotted in Fig. 2A (black dots). (C) Average along-strike and along-dip location of the seismicity measured from the MS and projected on the main fault and interevent times plotted in Fig. 2A (black dots). (D) Cumulated radius (black line) and cumulated moment (red line). (E) Effective stress drop. Time corresponds to the time of the last event within the 200-events window (see text). Fuchsia, turquoise and blue vertical lines show the time of F1, F2 and MS events, respectively. Turquoise and light orange dots represent events plotted in Fig. 2 (d) and (e), respectively.

Figure A4 .

 A4 Figure A4.S13. Down-dip section (strike N133°E and dip 50° to th co-seismic slip from Cirella et al. (2012), the foreshocks (F1 and F2 are highlighted with fuchsia and turquoise stars, respectively) and the mainshock (blue star). Distances along-strike and along-dip are relatives to the mainshock.

Figure A5 .

 A5 Figure A5.S1. Vertical component velocity seismograms of the initial P waves recorded for stations at different distances from the epicenter (See Fig. 1). Red and blue vertical lines indicate the time arrival of EP and IP, respectively. Time relative to the event origin time.

Figure A5 .

 A5 Figure A5.S2. Comparison of the waveforms of the emergent phase and two foreshocks Mw 3.9. (a) Location of the mainshock (MS 6 April) and the two foreshocks Mw 3.9 1 week (30 March) and 5 h (5 April) before it. Beachballs (compressional quadrants in colors) represent source mechanisms (reported by INGV). Foreshocks previously reported (Chiaraluce et al., 2011) are also plotted with color and size coded according to depth and magnitude, respectively. (b) Vertical component velocity seismograms of the initial P wave recorded at station AQU (red triangle in a) for the emergent phase (EP, blue line), and the two foreshocks Mw 3.9 (orange and red lines. We use AQU station because the strong motion stations in Fig. 1 correspond to a trigger system, and they did not record any foreshock). We note that the amplitude of EP exhibits a growth with time until the onset of IP. (c) Cumulative energy of the signals in (a). The constant growth of EP is highlighted, while the foreshocks show a flat slope after the end of the P waves. (d) Envelope of the signals presented in (a) following Tape et al. (2018) but without filtering and a smoothing of N=5. The envelope for EP grows steadily, but the envelopes of the foreshocks show a decay after the P-wave.

Figure A5 .

 A5 Figure A5.S3. Gaussian pulse used to assess the effect of the attenuation. We consider a pulse, where is the time and the duration of the pulse is .

Figure A5 .

 A5 Figure A5.S4. Examples of the attenuation. The dashed lines represent the reference pulse indicated in Fig. S1 without attenuation, and the black continues lines represent the pulse after applying the attenuation factor (Anderson & Hough, 1984) considering its respective travel time and Q. Black dots with red lines indicate the time duration between EP-IP () measured for every station represented by its respective travel time.

Figure A5 .

 A5 Figure A5.S5. Examples of phase arrivals modeled using a local velocity model(Chiaraluce et al., 2011) and the ObsPy TauP Toolkit[START_REF] Beyreuther | ObsPy: A Python toolbox for seismology[END_REF][START_REF] Crotwell | The TauP Toolkit: Flexible seismic traveltime and ray-path utilities[END_REF]. The times are aligned at 0s according to the respective P-phase of each station. The red dots indicate all expected phases, and the blue stars indicate the times measured between EP and IP. The horizontal axis is cut off at 5s on purpose, although more phases can be modelled at longer times.

Figure A5 .

 A5 Figure A5.S6. Comparison of the misfits between the parameters retrieved and indicated in Table 6.1 (Strike , Dip , Rupture Velocity Vr and Length L). The upper (blue lines) and right (red lines) panels in each figure show a profile to highlight the location of the minimum.

Figure A5 .

 A5 Figure A5.S7. (a) Four closest broadband stations (red triangles) used to estimate the magnitude of EP. (b) Estimated magnitudes for EP based on the reference magnitudes of the 11 foreshocks.

Figure A5 .

 A5 Figure A5.S8. Waveforms of the mainshock recorded in the 4 closest broadband stations indicated in Fig. S5. It is possible to observe that although the stations clipped soon after IP, the record between EP-IP is not.

Figure A5 .

 A5 Figure A5.S9. Scaling. Comparisons between the measurements from the present study (red dots) and the scaling relationships proposed byBeroza & Ellsworth (1996) for the seismic nucleation phase.

Figure A6 .

 A6 Figure A6.S1. Availability of the stations used for template matching. Stations with data gaps equal or larger than 2 hours per day are removed.

Figure A6 .

 A6 Figure A6.S2. Magnitude distribution of the clusters presented in Fig. 7.5.

Table 3 .1. Resume

 3 of kinematic parameters and number of aftershocks obtained for each mainshock.

	Event			Stress drop (MPa)	Number of aftershocks	Distance from the top of the slab (km)
	E1	3.49 5.13 1.08 3.67 5.68 0.94	1.07 1.19	18.2 14.7	2044	7
	E2	4.38 6.99 0.65 4.12 6.58 0.73	1.59 1.51	8.5 9.9	123	12
	E3	5.05 6.16 1.12 4.14 6.44 1.24	0.67 0.82	14.3 16.7	30	14
	E4	4.70 3.90 0.73 4.56 5.60 0.52	1.22 1.91	12.8 7.5	46	20
	E5	5.89 4.25 0.71 4.88 6.11 0.60	1.68 2.20	10.0 8.1	16	38
	E6	3.12 4.42 1.61 3.43 6.40 0.95	0.68 0.80	29.5 14.1	0	41

Table 6 . 1 .

 61 Parameters retrieved for the rupture initiation considering a line source model.

	Parameter	Value	Searching space	Grid search step
				1
				1
				0.1
				0.1

  for a circular shear crack model of constant stress drop. is especially useful because it allows us to compare our results with values of as a function of estimated from laboratory experiments, geological observations, modeling and seismological observations (for a review see Cocco et al., 2022):

		(6.3)
		(6.4)
	Assuming	and

Table 1 .

 1 Location of the Michilla, Calama, and Jujuy earthquakes. These hypocenters were inferred in this work based on a 3-D velocity model. Their absolute errors were estimated based on the 90% confidence level.

	Event	Origin time	Lon. [°]	Lat. [°]	Depth [km]	RMS [s]	N° of stations	Azimuthal gap [°]	Horizontal error [km]	Vertical error [km]
	Michilla	2007-12-16 08:09:17.28	70.1828	22.9962	43.5	0.05	24	81	2.5	3.0
	Calama	2020-06-03 07:35:34.82	68.5173	23.2502	113.4	0.73	9	174	8.0	12.0
	Jujuy	2015-02-11 18:57:20.3	66.8584	23.0936	228.7	0.8	13	193	9.0	14.0

Table 2 .

 2 Comparison of the best dynamic models of the Michilla, Calama, and Jujuy earthquakes. For the Calama earthquake, values of their posterior mean and standard deviation are also shown in parenthesis.

	Parameter	Michilla	Calama	Jujuy
	Semi axis [km]	4.0	7.08 (7.05±0.04)	7.94
	Semi axis [km]	10.12	12.21 (12.24±0.05)	4.87
	Center [km]	0.85	12.83 (12.95±0.1)	12.71
	Center [km]	-2.0	14.56 (14.82±0.26)	11.63
	Rotation angle [°]	85.9	159.8 (162±1.8)	203.4
	[MPa]	14.97	10.05 (9.86±0.14)	11.87
	[MPa]	19.18	11.87 (11.81±0.12)	14.37
	[MPa]	23.65	15.35 (15.11±0.25)	16.1
	[km]	0.98	1.46 (1.47±0.01)	1.09
	[m]	0.65	0.7 (0.67±0.02)	0.41
		2.08	1.5 (1.53±0.02)	1.97

Table A2 .

 A2 S1. Hypocenters reported by the Chilean National Seismological Center

	Event Date UTC	Latitude (º)	Longitude	Depth (km)	Magnitude
	E1	10/10/2017 06:32:20	-69.739	-18.489	99.7	6.3
	E2	06/03/2011	-69.440	-18.309	112.4	6.2
		12:31:57
	E3	23/03/2015 04:51:38	-69.275	-18.415	122.2	6.3
	E4	01/11/2018 22:19:52	-69.409	-19.649	101.4	6.2
	E5	14/05/2012	-70.237	-18.115	119.6	6.4
		10:00:40
	E6	21/01/2018 01:06:40	-69.579	-18.861	129.3	6.3

Table A2 .

 A2 S2. Detailed results for template matching

	Event	Number of	Number or	Distance from the top
		templates	aftershocks	of the slab [km]
	E1	45	2044	7
	E2	1	123	12
	E3	3	30	14
	E4	11	46	20
	E5	1	16	38
	E6	1	0	41

Table A2 .

 A2 S3. Results for kinematic inversions considering nodal plane 1 (NP1) according to USGS moment Tensor. and are semi-axis for elliptical patch, maximum slip and velocity of rupture.

	Event Strike	Dip	Rak				Stress	Misf
		(º)	(º)	e				drop	it
				(º)				(MPa
								)
	E1	180	18	-70	3.49 5.13	2.14 6.72 -0.70 1.08 1.07	18.2 0.43
	E2	210	26	-37	4.38 6.99	4.90 6.16 2.52 0.65 1.59	8.5 0.36
	E3	215	16	-49	5.05 6.16 -5.04 -4.90 0.21 1.12 0.67	14.3 0.37
	E4	176	24	-82	4.70 3.90	0.98 0.84 0.73 0.73 1.22	12.8 0.48
	E5	206	22	-22	5.89 4.25	5.21 6.86 -0.14 0.71 1.68	10.0 0.47
	E6	189	26	-75	3.12 4.42	2.01 6.44 -3.08 1.61 0.68	29.5 0.32
		Minimum Value		3.00 3.00	0 -7.00 -7.00	0.5 0.1
		Maximum Value		7.00 7.00	6.28 7.00 7.00	4.0 5.0

Table A2 .

 A2 S4. Results for kinematic inversions considering nodal plane 2 (NP2) according to USGS moment Tensor. and are semi-axis for elliptical patch, maximum slip and velocity of rupture.

	Event Strike	Dip	Rake							Stress	Misfit
		(º)	(º)	(º)							drop
											(MPa)
	E1	339 73	-96	3.67	5.68 -3.08	4.06	0.50 0.94 1.19	14.7	0.44
	E2	335 75 -111	4.12	6.58 -2.80	5.18	0.23 0.73 1.51	9.9	0.38
	E3	353 78 -101	4.14	6.44	5.60 -0.70 -0.28 1.24 0.82	16.7	0.37
	E4	348 66	-93	4.56	5.60	3.21	6.86 -0.70 0.52 1.91	7.5	0.47
	E5	316 82 -111	4.88	6.11 -6.72	6.86	0.29 0.60 2.20	8.1	0.50
	E6	352 65	-97	3.43	6.40	6.72	6.30	0.01 0.95 0.80	14.1	0.32
		Minimum Value		3.00	3.00	0 -7.00 -7.00	0.5	0.1
		Maximum Value		7.00	7.00	6.28	7.00 7.00	4.0	5.0

Table A2 .

 A2 S5. Parameters used to develop the two-dimensional model through finite element

  The lower continental mantle has a conductivity equal to the upper case, but allowing local flows in the area, driven by the movement of the incoming lithosphere.

		Value
	Oceanic lithosphere speed	66 mm/yr
	Plate age	50 myrs
	Oceanic mantle temperature	1450 ºC
	Upper temperature (surface)	0 ºC
	Continental mantle edge temperature	900 1400 ºC
	Upper crust heat source	1.9 µW/m 2
	Lower crust heat source	0.4 µW/m 2
	Sediment heat source	1.9 µW/m 2
	Upper crust depth	17 km
	Lower crust depth	35 km
	Maximum model length	400 km
	Extension to the sea before subduction	100 km
	Upper continental mantle depth	50, 85 or 120 km
	Incoming lithosphere thermal conductivity	2.9 W/(m K)
	Incoming lithosphere heat capacity	3.3 MJ/(m 3 K)
	Continental crust thermal conductivity	2.5 W/(m K)
	Continental crust heat capacity	2.5 MJ/(m 3 K)
	Asthenosphere dynamical viscosity	10 20 Pa s

Note:

Table A5 .

 A5 S1. Absolute location parameters for EP and IP estimated by Di Stefano et al. (2011).

		Origin Time	Lat. [°]	Lon. [°]	Dep. [km]
	EP 2009-04-06 01:32:40.74 42.3548 13.3813	8.55
	IP	2009-04-06 01:32:41.61 42.3577 13.3822	6.64

Table A5 .

 A5 S2. Relocations retrieved using GrowClust and different initial locations for EP and IP (Test 1, see Text A5.S2). Origin Time, Latitude (Lat), Longitude (Lon) and Depth (Dep) are indicated, and their respective time (et), horizontal (eh) and vertical (ez) errors. Estimated strike and dip from EP to IP are also indicated (Dip is defined as 0°=vertically down, 90°=horizontal, and 180°=vertically up).

	Origin Time	et (s)	Lat (°) eh (m) Lon (°) eh (m) Dep (km) ez (m)
	EP 2009-04-06 01:32:39.67 0.001 42.3554 0.053 13.3803 0.053	7.76	0.018
	IP 2009-04-06 01:32:41.33 0.001 42.3571 0.053 13.3832 0.053	7.43	0.018
	Parameter	Value
	Strike EP-IP	N58°E
	Dip EP-IP	143°
	Dh EP-IP	462m

Table A5 .

 A5 S3. Relocations retrieved using GrowClust and the same initial location for EP and IP (Test 2, see Text A5.S2). Origin Time, Latitude (Lat), Longitude (Lon) and Depth (Dep) are indicated, and their respective time (et), horizontal (eh) and vertical (ez) errors. Estimated strike and dip from EP to IP are also indicated (Dip is defined as 0°=vertically down, 90°=horizontal, and 180°=vertically up).

	Origin

Time et (s) Lat. (°) eh (m) Lon (°) eh (m) Dep (km) ez (m)

	EP 2009-04-06 01:32:39.67 0.002 42.3570 0.047 13.3808 0.047	6.82	0.012
	IP 2009-04-06 01:32:41.33 0.002 42.3584 0.047 13.3836 0.047	6.46	0.012
	Parameter	Value
	Strike EP-IP	N50°E
	Dip EP-IP	138°
	Dh EP-IP	452m

Acknowledgements Acknowledgements vii

Acknowledgments

We thank the Incorporated Research Institutions for Seismology Data Management Center and Centro Sismológico Nacional who provided the raw data used in this work. LC thanks the scholarship PFCHA/Mag sterNacional/2018-22189566, granted by CONICYT. SR and ECR thank the support of Chilean Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT), Grants 1170430 and 1170009, respectively. PP and LC were supported by the European Union Horizon 2020 Research and Innovation Programme (grant agreements, 802777-MONIFAULTS). AO acknowledges FONDECYT Grants 1191903 and 1201311, Basal CMM ANIF PIA AFB170001 and Fondap CR2 Grant 15110009. We also acknowledge the support of CONICYT PIA/Anillo de Investigación en Ciencia y Tecnología ACT172002. We also thank the editor Eiichi Fukuyama, Jun Korenaga and another anonymous reviewer for their comments, which improved the manuscript. We used Generic Mapping Tools 6 software for preparing some plots.

Acknowledgments

PP and LC received funding from the European Research Council (ERC) under the European Union Horizon 2020 Research and Innovation Programme (grant agreements, 802777-MONI-FAULTS). JCB was also supported by ANID PIA (ACT192169). SR and JCB were supported by Fondecyt project (N° 1200779, ANID, Chile). JBA is supported by Fondecyt project (N° 3200633, ANID, Chile). SR thanks to Programa de Riesgo Sísmico (PRS) from University of Chile, Chile.

Acknowledgments

LC and PP were supported by the European Union Horizon 2020 Research and Innovation Programme (grant agreements, 802777-MONIFAULTS). WBF acknowledges support from the National Science Foundation under award EAR-2103408. LC, PP, and WBF were also supported by the Thomas Jefferson Fund award #093. We thank the editor Rachel Abercrombie, the associate editor, and Alessandro Vuan and an anonymous reviewer for their comments that improved our manuscript.

Acknowledgments

LC and PP were supported by the European Union Horizon 2020 Research and Innovation Programme (grant agreements, 802777-MONIFAULTS).

Acknowledgments

The author C. Herrera designed and wrote most of the manuscript. He also carried out the strong motion analysis and the dynamic modeling of the Calama mainshock. F. Pastén-Araya designed and contributed with the discussion that relates the tomography and thermal models with the analyzed earthquake sequences. F. Pastén-Araya and R. Madariaga thank the Programa de Riesgo Sísmico (PRS) of the Universidad de Chile. L. Cabrera carried out the template matching using the University of Grenoble Alpes (UGA) High-Performance Computing infrastructures CIMENT. He was supported by the European Union Horizon 2020 Research and Innovation Programme (grant agreements, 802777-MONIFAULTS). B. Potin and E. Rivera contributed with the location and relocation analyses of the three studied earthquakes and their aftershocks. Part of the computations were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. S. Ruiz, R. Madariaga, and E. Contreras-Reyes contributed with editing assistance and review of the manuscript during preparation. E. Contreras-Reyes also acknowledges the support of PIA/FONDEYT grant 1210101. We finally thank Jörg Renner, Frederik Tilmann, and another anonymous reviewer for their thorough reviews that improved this work.

Data Availability Statement

Seismological data are available through the IRIS Data Management Center (IRISDMC) at http://service.iris. edu/fdsnws/dataselect/1/ and can be obtained using the IRIS DMC FDSNWS web service. The corrected GPS time series are at https://www.csn.uchile.cl/red-sismologicanacional/red-gps/.

Article published in Journal of Geophysical Research: Solid Earth (2022) https://doi.org/10.1029/2021JB023888

Data Availability Statement

Data was downloaded from the Istituto Nazionale di Geofisica e Vulcanologia (INGV, 2006) using obspyDMT (https://github.com/kasra-hosseini/obspyDMT, [START_REF] Hosseini | Obspydmt: A python toolbox for retrieving and processing large seismological data sets[END_REF]. The fast matched filter [START_REF] Beaucé | Fast matched filter (FMF): An efficient seismic matched filter search for both CPU and GPU architectures[END_REF] used in this study can be found at https://github.com/beridel/fast_matched_filter. Computations were performed using the University of Grenoble Alpes (UGA) High-Performance Computing infrastructures CIMENT (https://ciment.univ-grenoble-alpes.fr/wiki-pub/index.php/Welcome_to_the_CIMENT_site!).

Data availability statement

Waveform data from multiparametric stations were downloaded from the International Federation of Digital Seismograph Networks (FDSN) web services using the ObsPy toolkit [START_REF] Beyreuther | ObsPy: A Python toolbox for seismology[END_REF]. Waveforms from the earthquake-triggered network of accelerometers of the CSN can be accessed from their website (evtdb.csn.uchile.cl/). The earthquake catalogs used in this study can be accessed from their respective websites: CSN catalog (www.sismologia.cl), GCMT catalog (www.globalcmt.org). Maps were created using Generic Mapping Tools (Wessel et al. 2013).

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.

Chapter 5 -Temporal Evolution of Foreshocks Leoncio Cabrera, Piero Poli and William Frank Abstract Understanding under which physical conditions large earthquakes begin, is a key question in Earth science. Laboratory experiments and numerical models have shown that earthquake nucleation has distinct phases: a quasi-static and an acceleration stage, followed by dynamic propagation. However, obtaining observations of such or similar processes in nature is complex. Here, we report on the rupture initiation of the Mw 6.1 2009 L'Aquila earthquake. From the detailed analysis of seismic waves recorded at several stations, we identify an ~0.6-s signal preceding the large dynamic rupture. From the geometrical characterization and rupture parameters of this initial phase, we infer that the rupture struggled to initiate exhibiting a slow rupture velocity ( km/s) and low seismic efficiency ( ) due to a complex environment in the region where the rupture starts. We also show that the parameters of the rupture initiation are representative of scale-dependent quantities for slip-dependent nucleation models.

Plain Language Summary

Understanding the process leading to a large earthquake is a key question in Earth science with implications for earthquake prediction and risk assessment. Although results from laboratory experiments and numerical simulations show that earthquake nucleation is composed of several preliminary stages, these stages and the associated processes are very difficult to observe in natural earthquakes. In this work we study the Mw 6.1 2009 L'Aquila earthquake, and show that it started with a slow rupture, where most of the energy was used to propagate the rupture. Our results show that the rupture struggled to initiate, and it did it slowly probably due to the conditions in which it starts. In addition, our results show similarities with previous work done for earthquakes in different regions of the world and also using theoretical models. Corresponding author: Carlos Herrera (carlos@onurseemann.com)

Summary

We study the 2020 6.8 Calama earthquake sequence that occurred within the subducting oceanic Nazca plate. The mainshock is modeled via waveform inversion using a dynamic rupture model, while detection and location techniques are used to better characterize its aftershock sequence. We analyze the local seismotectonic and thermal context of the subducting Nazca plate to understand the trigger mechanism of this earthquake and how it compares with other significant earthquakes in the vicinity. The stress drop and the related dynamic rupture parameters of the Calama mainshock are similar to those of the nearby 2007 6.8 Michilla and 2015 6.7 Jujuy intraslab earthquakes, which occurred to the west (trenchwards) and to the east (under the back-arc) of the Calama earthquake, respectively. The sequences of these three events were located using a 3-D tomographic velocity model. While the Michilla earthquake sequence occurred within the oceanic crust at temperatures of ~250°C, the Calama sequence occurred within the upper lithospheric mantle at ~350°C and exhibited a smaller aftershock productivity than Michilla. Additionally, the 3-D tomographic model shows intermediate less hydrated environment that could be responsible for the smaller aftershock productivity of the Calama earthquake. ground motion models for Chilean intraslab earthquakes (see Text S1 and Figure S3 in the Supplementary Material). The low frequency rupture properties of the mainshock were inferred via inversion using a finite-fault model. Following the method used to model the Michilla and Jujuy earthquakes (Ruiz & Madariaga 2011;Herrera et al. 2017), the rupture model used in this work assumes an elliptical coseismic slip distribution with semi-axes and , centered at ( , ) within the fault plane. This ellipse is also allowed to rotate around its center. The rupture nucleates at the hypocenter within a circular area. The overall rupture propagation in this model is controlled by a slip-weakening friction law (Ida 1972). This allows the determination of dynamic rupture parameters, such as: stress drop ( ), yield stress ( ), slipweakening distance ( ), and a nucleation of radius with a stress acting inside it [START_REF] Madariaga | Earthquake dynamics on circular faults: A review 1970 2015[END_REF]). The finite fault was centered at the hypocenter and was oriented using the strike, dip, and rake of the NE-dipping plane of the focal mechanism reported by GCMT, as suggested by the spatial distribution of the located aftershocks. Prior to inversion, the baseline-corrected acceleration waveforms were integrated to velocity and filtered between 0.02 and 0.1 Hz with a Butterworth bandpass filter. Finally, the horizontal channels were rotated into radial and transverse components. To create the modeled waveforms, the wave propagation was simulated with the AXITRA code (Bouchon 1981;Coutant 1989) based on a 1-D velocity model [START_REF] Husen | Accurate hypocenter determination in the shallow part of the Nazca subduction zone in Northern Chile using a combined on-/offshore network[END_REF]. The inversion was performed using the Neighborhood Algorithm (Sambridge 1999), which in this case minimizes the misfit ( ) to find the best fitting model: which runs over the samples of the observed (

) and predicted ( ) waveforms. The three components (radial, transverse, and vertical) were used in the inversion. , where is time and the duration of the pulse (Fig. S1). For the attenuation operator we consider an exponential attenuation described by 1 (Anderson & Hough, 1984), where is the frequency, the travel time for each single station and the quality factor. Since we do not have a Q model for the study area, we model the effect of different Q in the range 100-600, where 600 represents almost no-attenuation and 100 high attenuations. The results are presented in Fig. S3.

Text A5.S2. We relocate EP and IP using the GrowClust software [START_REF] Trugman | GrowClust: A hierarchical clustering algorithm for relative earthquake relocation, with application to the Spanish Springs and Sheldon, Nevada, earthquake sequences[END_REF], a local velocity model (Chiaraluce et al., 2011) and all the measured values at the strong motion stations (Fig. 1) as the differential delay times. As described in the main text, we performed two tests: In the first one (test 1), we used as initial locations for ED and IP those estimated by Di Stefano et al. (2011) listed in Table S1. The results presented in Table S2 indicate a northeast rupture (strike of 58°) and up-dip (143° dip). In the second test (test 2), we assigned for both ED and IP the same initial location of IP (mainshock) estimated by Di Stefano et al. (2011). The results presented in Table S3 likewise indicate a northeast (strike 50°) and up-dip (dip 150°) rupture, in agreement with the results presented in the main text.

Article in preparation for Submission

This chapter corresponds to work still in progress. In particular, we develop a methodology to classify the seismicity of a dense seismic catalog automatically, based on physical properties of the seismicity. To this aim, we first create a new seismic catalog for Central Italy (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) and then we classify the seismicity using hierarchical clustering. Although, the catalog generated in this chapter and the clustering methodology are finished, we are still working on the analysis and interpretation of the results in more detail. However, the main advances, results, conclusions and perspectives obtained so far are discussed.