Journal Pre-proof Seismicity during the recent activity (2009–2020) of Turrialba volcano, Costa Rica Leonardo van der Laat, Mauricio M. Mora, Javier Fco. Pacheco, Philippe Lesage, Esteban Meneses PII: S0377-0273(22)00182-2 DOI: https://doi.org/10.1016/j.jvolgeores.2022.107651 Reference: VOLGEO 107651 To appear in: Journal of Volcanology and Geothermal Research Please cite this article as: L. van der Laat, M.M. Mora, J.F. Pacheco, et al., Seismicity during the recent activity (2009–2020) of Turrialba volcano, Costa Rica, Journal of Volcanology and Geothermal Research (2022), https://doi.org/10.1016/ j.jvolgeores.2022.107651 This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2022 The Author(s). Published by Elsevier B.V. https://doi.org/10.1016/j.jvolgeores.2022.107651 https://doi.org/10.1016/j.jvolgeores.2022.107651 https://doi.org/10.1016/j.jvolgeores.2022.107651 ● Turrialba volcano slowly transitioned from a closed- (2010) to an open-system (2017) ● Seismic tremor decreases precedes eruptive phases at Turrialba volcano ● Pre-eruptive tremor abatement is accompanied with SO2 flux decrease ● Chromatic signals are common seismic precursors (e.g. tornillos, harmonic tremor) ● Magnitude 5.5 earthquake correlates with the transition to an open-vent system Jo ur na l P re -p ro of Journal Pre-proof Declaration of interests 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. The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Jo ur na l P re -p ro of Journal Pre-proof Seismicity during the recent activity (2009-2020) of Turrialba volcano, Costa Rica - Author Statement August 15, 2022 Author contributions Leonardo van der Laat: Conceptualization, Methodology, Software, For- mal Analysis, Investigation, Writing - Original Draft, Visualization; Mauricio Mora: Supervision, Conceptualization, Methodology, Software, Formal Analy- sis, Investigation, Writing - Original Draft, Resources, Funding acquisition, Vi- sualization Javier Fco. Pacheco: Supervision, Conceptualization, Resources, Formal Analysis, Investigation, Funding acquisition, Philippe Lesage: Formal Analysis, Writing - Review & Editing; Esteban Meneses: Writing - Review & Editing, Funding acquisition. 1 Jo ur na l P re -p ro of Journal Pre-proof Jo ur na l P re -p ro of Journal Pre-proof Seismicity during the recent activity (2009-2020) of Turrialba volcano, Costa1 Rica2 Leonardo van der Laata,b,c,d,1,∗, Mauricio M. Morab,c, Javier Fco. Pachecod, Philippe Lesagee, Esteban Menesesf 3 aDepartment of Earth and Environmental Sciences, University of Michigan, 1100 North University, Avenue, 4518 North University Building, Ann4 Arbor, MI 48109-10055 bEscuela Centroamericana de Geologı́a, Universidad de Costa Rica, Apdo. 11501-2060, San José, Costa Rica6 cRed Sismológica Nacional, Universidad de Costa Rica, Apdo. 11501-2060, San José, Costa Rica7 dObservatorio Vulcanológico y Sismológico de Costa Rica, Universidad Nacional, Apdo. 2386-3000, Heredia, Costa Rica8 eUniversité Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, Université Gustave Eiffel, ISTerre, 38000 Grenoble, France9 fColaboratorio Nacional de Computación Avanzada, Centro Nacional de Alta Tecnologı́a, Apdo. 1174-1200, San José, Costa Rica10 Abstract11 Turrialba is a stratovolcano located at the easternmost part of the Costa Rican volcanic front. After remaining qui-12 escent for more than a century, in 1996 it started to show signs of unrest, until a first phreatomagmatic explosion13 occurred on January, 2010. Since then, the activity evolved from phreatic to magmatic, in a series of distinct eruptive14 phases. In this paper, we investigate the seismic records that span the whole eruptive process (2010-present), in order15 to identify precursory signals and characterize the volcanic evolution. A long-term analysis was carried out based on16 the continuous records, as well as seismic catalogs (volcano-tectonic seismicity, harmonic tremor, etc.). In addition,17 the gradual character of the evolution of this eruption allowed for the analysis of independent precursory stages. Thus,18 we inspected in detail the most important of those periods, particularly, prior to the first 2010 phreatomagmatic erup-19 tion, and prior to the 2016 transition to an open vent system. Temporary tremor amplitude decreases were found to20 precede most of the eruptive phases. In total, 9 pre-eruptive tremor abatement periods were identified spanning several21 days (5-44), which often concurred with a decrease in the SO2 flux. The analysis of the volcano-tectonic seismicity22 highlights the migration of magma from a deep (6-10 km) reservoir beneath the neighboring Irazú volcano towards23 Turrialba volcano, especially between the years 2015 and 2016. This activity peaked on December 2016 when a Mw24 5.5 earthquake took place between both volcanoes. Harmonic tremor episodes thrived in the later phase when the25 system finally opened (2017-2018). In the short-term, compounded tonal seismic signals were identified as precursor26 events, such as long-period events followed by harmonic tremor or by a multichromatic coda similar to tornillo-type27 events. The co-occurrence of tremor amplitude decreases and tonal seismic signals is interpreted to be caused by a28 sealing of the hydrothermal system, which blocked the circulation of fluids and permitted the resonances in the inner29 cavities. This process leaded to pressure accumulation and the consequent eruptions. Thus, trough a series of cycles of30 sealing and rupture the system of conduits gradually opened. The seismic characterization of this eruption constitutes31 insightful knowledge useful for monitoring and risk assessment purposes.32 Keywords: tremor, precursor, volcano, seismology, LP, tornillo33 1. Introduction34 Turrialba volcano is a (3,340 m a.s.l.) basaltic–andesitic stratovolcano located at the eastern edge of the Central35 Volcanic Range (CVR) of Costa Rica, 35 km east-northeast of San José, the capital and socioeconomic center of the36 ∗Corresponding author Email addresses: laat@umich.edu (Leonardo van der Laat), mauricio.mora@ucr.ac.cr (Mauricio M. Mora), javier.pacheco.alvarado@una.ac.cr (Javier Fco. Pacheco), lesage@univ-smb.fr (Philippe Lesage), emeneses@cenat.ac.cr (Esteban Meneses) 1Present address: Department of Earth and Environmental Sciences, University of Michigan, 1100 North University, Avenue, 4518 North University Building, Ann Arbor, MI 48109-1005 Preprint submitted to Elsevier August 15, 2022 Jo ur na l P re -p ro of Journal Pre-proof country. It is one of the five active volcanoes in Costa Rica and it has erupted with significance at least 10 times during37 the last 7,000 years (Reagan et al., 2006; Alvarado et al., 2020). The last major eruption took place between 1864 and38 1866 and the volcano then remained quiescent for more than a century. In 1996, the volcano started to show signs of39 unrest that turned conspicuous in 2007 evidenced by geochemical changes in fumarole gas composition and increased40 seismicity (Hilton et al., 2010; Martini et al., 2010; Vaselli et al., 2010). The first phreatic eruptive phase occurred in41 January 5, 2010. Since then, the activity evolved from phreatic to magmatic, in a slow process of opening of conduits.42 Understanding precursory seismicity at short and long terms is globally one of the major tasks in order to reduce43 volcanic risk. At Turrialba volcano, we were able to record part of the initial unrest before the first eruption in 201044 and all of the eruptive evolution towards the open-vent phase, this offers an important opportunity to understand45 seismic source processes at different stages of the volcanic activity of Turrialba.46 In this work we focus on the precursory activity of the main eruptive stages at Turrialba volcano since 2010 up to47 2020 by means of close inspection of seismic records with different time-frequency analysis, in order to understand48 the evolution of the eruption dynamics. We also carry out a detailed description and classification of Turrialba volcano49 seismicity. This is the first time that overall seismicity of Turrialba volcano is presented together with an interpretative50 model of the evolution of the system over the 11 years of activity.51 1.1. Geological setting52 Turrialba volcano, Costa Rica, is located in a complex geotectonic background where the Cocos, Caribbean, Nazca53 plates and the Panama microplate interact (Figure 1). The Cocos plate subducts beneath the Caribbean plate and the54 Panama microplate at an approximate rate of 90 mm/year (LaFemina et al., 2009). Subduction beneath Costa Rica55 is highly complex due to variations on the Cocos plate slab morphology and its petrological characteristics, but also56 because of its double interaction with the Caribbean plate and Panama microplate (Von Huene et al., 1995; Gazel57 et al., 2021). The southeastern portion of the Central American Volcanic Front (CAVF) in Costa Rica and Panama,58 developed on top of oceanic crust terranes, including the western edge of the Caribbean Large Igneous Province59 (CLIP) (Gazel et al., 2021). The oceanic crust subducting along northwestern Costa Rica was produced at the East60 Pacific Rise, while that beneath central Costa Rica and Panama was formed at the Cocos-Nazca spreading center and61 includes the 13.0–14.5 Ma Galapagos hotspot track (Von Huene et al., 1995; Werner et al., 1999; O’Connor et al.,62 2007). The subduction of the Cocos plate under the Caribbean plate generates profuse volcanism, resulting in the63 modern Costa Rica volcanic front (CRVF), that is oriented NW-SE and partitioned in three segments: at the NW the64 Guanacaste volcanic range (GVR), the Tilarán volcanic range (TVR) and to the SE the Central volcanic range (CVR).65 From NW to SE, the CVR comprises Platanar, Poás, Barva and Irazú-Turrialba volcanoes. Turrialba, in particular, is66 the only volcano that is shifted 10 km towards the NE from the CVR main axis.67 Turrialba volcano lies on sedimentary deposits of the Limón Basin and its origin dates back to the Middle to68 Upper Pleistocene. Ruiz et al. (2010) and Soto (2012) consider three construction phases: 1) Proto-Turrialba that69 began around 1 Ma; 2) Paleo-Turrialba whose last phase corresponds to the Finca Liebres volcano (250 ka); and70 3) Neo-Turrialba, the current phase. The lava and pyroclast compositions vary from basalt to dacite, being andesite71 and dacite the most important compositions (Ruiz et al., 2010; Soto, 2012). DeVitre et al. (2019) proposed that72 a multi-stage magma mixing occurred at Turrialba, starting at deep mid-crust reservoirs with a mingling between73 high-Nb decompression melt induced by corner-flow from the backarc and a low-Nb volcanic front component. A74 second stage resulted from the chaotic mixing at shallow reservoirs between the resulting basaltic andesite and the75 rhyolitic end-member, delivering the presumably differentiated residual magma from previous events such as that of76 1864–1866.77 Geophysical evidence suggests that the feeding system beneath Turrialba volcano is made up of two magma78 reservoirs. Conde et al. (2014) deduced two levels of magma storage from the clustering of localised seismic events,79 one probably at 4–6 km depth and the shallow one at roughly 1 km below the summit. de Moor et al. (2016), Müller80 (2018) and Badilla and Taylor (2019) proposed a mid-crust reservoir located in the range between 5-10 km, which81 is located beneath Irazú volcano according to Lücke et al. (2010) and Müller (2018). Di Piazza et al. (2019) also82 interpreted a deep magma reservoir (∼13 km deep) for the andesitic magma involved in the ∼2 ka sub-Plinian El83 Retiro eruption. A second, shallower and smaller, reservoir could be in the range 1-2 km (de Moor et al., 2016;84 Badilla and Taylor, 2019).85 Tectonic faults cutting Turrialba volcanic structure are characterized by left (NE-SW) and right-lateral (SW-NE)86 strike-slip systems (Soto, 1988; Montero et al., 2013; Calvo et al., 2019) (Figure 2). Turrialba volcano has three87 2 Jo ur na l P re -p ro of Journal Pre-proof Figure 1: Simplified tectonic setting of Costa Rica and location of Turrialba volcano. Abbreviations are: Guanacaste volcanic range (GVR), Tilarán volcanic range (TVR), Central volcanic range (CVR), the Central Costa Rica Deformed Belt (CCRDB) and North Panama Deformed Belt (NPDB), Panama Fracture Zone (PFZ) and Nazca plate (Np) (Denyer et al., 2014). Black continuous and dashed lines represent faults. Triangles denote reverse faulting and point towards the overriding plate. Bathymetry and topography from Ryan et al. (2009). Active volcanoes are denoted by red triangles: Rincón de la Vieja (RV), Arenal (A), Poás (P). Main cities are indicated with white letters: SJ: San José, Alajuela (A), Heredia (H) and Cartago (C). summit craters: the Southwest, which is the active one, the Central and the Northeast. They follow a NE-SW direction88 together with the pyroclastic cones Armado and Tiendilla to SW and the avalanche amphitheater to NE, which is the89 result of an structural control by the Elia and Ariete faults which are part of the NE system (Soto, 1988; Linkimer,90 2003; Calvo et al., 2019) (Figure 3). On the other hand, the NW fault system, contains the Rı́o Sucio fault and the91 Liebres faults. The Rı́o Sucio fault is presumed to be the source of a Ms 5.9 historical earthquake that occurred92 on December 30, 1952, 7 km NE from the crater of Irazú volcano, known as the Patillos earthquake (Montero and93 Alvarado, 1995). More recently, another important Mw 5.5 earthquake took place on December 1st, 2016 at 00:2594 UTC, 7 km SE from the crater of Irazú volcano. This event, known as the Capellades earthquake, might have ruptured95 part of the Liebres fault (Linkimer et al., 2018) (Figure 2).96 1.2. The 2010-present eruption97 Tephrostratigraphic studies evidence at least 10 eruptions at Turrialba volcano during the last 7,000 years (Reagan98 et al., 2006; Alvarado et al., 2020). Since 1450 B.C. eruptions occurred with a period of around 230 years (Alvarado99 et al., 2020). The last and only known historical eruption, prior to recent activity, took place from September 1864 to100 February 1866, producing ash plumes that spread west, passing through most of the central valley and reaching the101 Pacific coast (Reagan et al., 2006; Alvarado et al., 2021).102 The first signs of unrest related to the most recent eruption could possibly be traced back to the year 1982 when103 two seismic swarms occurred in the Irazú-Turrialba region (Güendel, 1984). Those were interpreted as tectonic by104 Alvarado et al. (1986), Barquero and Alvarado (1989) and Fernández et al. (1998). In 1987 there was still no sign of105 activity (Figure 4A). However during the second half of the 1990s Turrialba volcano showed physical, geochemical106 and geophysical signs of unrest (Tassi et al., 2004; Martini et al., 2010). The seismic activity increased slowly and107 progressively from several isolated swarms during the 1990s and early 2000s to a more abundant and varied seismicity108 since 2005 (Martini et al., 2010). The increase of volcanic activity became more visible in 2007 (Figure 4B), when109 the composition of fumarolic gases changed from a hydrothermal to a magmatic signature and seismicity increased110 3 Jo ur na l P re -p ro of Journal Pre-proof 1000 1000 3000 3000 Longitude [°] L a ti tu d e [ °] Irazú volcano 2016 M W 5 5 Capellades 1 52 M S 5 Patillos 5 km N L ie b re s fa u lt Capella des fa ult R ío S ucio fault Turrialba volcanoTurrialba volcano E lia fa ul t RSN Seismic stations Fault trace Epicenter OVSICORI Infrasound sensor Figure 2: Irazú-Turrialba volcanic complex. Black lines are the main active faults Montero and Alvarado (1995); Montero et al. (2013); Linkimer et al. (2018). The 1952 Patillos earthquake 2016 Capellades earthquake (Mw 5.5) are shown with stars. The RSN and OVSICORI seismic stations are indicated with blue and yellow triangles, respectively. The rectangle indicate the zoom area showed in Figure 3. Topographic contour levels are shown every 100 and 500 m with thin and thick lines, respectively. 4 Jo ur na l P re -p ro of Journal Pre-proof Southwest crater (active) Central crater Northeast crater B 1 0 .0 4 °N 1 0 .0 2 °N 9 .9 9 °N 83.79°W 83.77°W 83.75°W + - + - Tiendilla Armado Turrialba volcano CVTR VTUN Southwest (active) Central Northeast 0 0.75 1.5 km A 1 0 .0 4 °N 1 0 .0 2 °N 9 .9 9 °N 83.79°W 83.77°W 83.75°W VTCE CVTC VTCG CVTV CVTQ VTRT E lia fa ul t E lia fa ul t Arie te fa ult +- Craters Pyroclastic cone Avalanche/erosion scarp Fault trace Fault movement Seismic stations: RSN OVSICORI Infrasound sensor Figure 3: (A): Digital elevation model of Turrialba volcano with volcano-tectonic geomorphic features and seismic stations (triangles). Summit craters and Tiendilla and Armado pyroclastic cones, Elia and Ariete faults show a NE-SW trend. (B): Panoramic view of Turrialba volcano summit taken by Mauricio M. Mora on November 26, 2015, showing the Southwest (active), Central and Northeast craters. 5 Jo ur na l P re -p ro of Journal Pre-proof (Hilton et al., 2010; Vaselli et al., 2010; Martini et al., 2010).111 The first phreatomagmatic eruption occurred on January 5, 2010 (10:57 local time). It opened two vents in the SW112 inner wall of the active crater that became one on January 8, 2010 (Alvarado et al., 2016a), as it is observed in figures113 4F and 5A. Minor ash emissions occurred on June, 12 and July, 13 later that same year (Soto and Mora, 2012).114 On January 14, 2011, residents of a neighboring town (2.3 km from the active crater) informed about an ash fall115 and the presence of a strong smell of sulfur. However, no seismic activity was found that could indicate the occurrence116 of an explosion or volcanic eruption (Martı́nez and Pacheco, 2011). Later that year, on July, a new vent was observed117 on the W wall of the active crater. Another vent opened on January 12, 2012, on the SE external flank of the West118 crater. This event was accompanied by ash emission for a few hours. It was followed by a second ash emission from119 the same vent on January 18, 2012 (Figures 4F and 5C). On May 21, 2013 occurred a new ash emission that lasted120 2 hours (Avard et al., 2012, 2013; Alvarado et al., 2016a). The ash came out simultaneously through the 2010 and121 2012 vents (Figure 5D). For the 2012 and 2013 eruptions it was concluded, based on geophysical and geochemichal122 records, that the events were caused by an excessive accumulation of gases at shallow depth (Avard et al., 2012, 2013).123 On October 29, 2014, the volcano emitted ash for almost 13 hours in a sustained manner and ended with a 25-124 minute high-energy explosion that destroyed most of the eastern wall of the West Crater (Alvarado et al., 2016a)125 (Figures 4C and 5D, left). The explosion was followed by a practically continuous emission of ash that lasted about126 3 days. The right image on Figure 5D shows the deposits scattered by the eruption towards the east. This eruption127 established a breaking point in the eruptive dynamics of Turrialba volcano, which shifted from annual isolated eruptive128 events to well-defined eruptive phases (Figure 5 right panel).129 The first eruptive cycle lasted from October to December 2014, the second from March to June 2015 (an eruption130 on March 13, 2015 at 11:07 local time is shown in Figure 5E) and the third from October, 16 to October, 31, 2015.131 The last two cycles were separated by an eruptive episode on August 15, 2015. After the last cycle in 2015, some132 eruptions occurred between November and December 2015 in a context of low seismicity and expanded the area of the133 active crater as it is shown in Figure 4D. Subsequently, there were a few ash emissions and small isolated explosions134 between January and early April 2016.135 A new eruptive stage began on 2016, with intense seismicity since April 27 and followed by a long-lasting ash136 emission on April 30. Finally, the explosive eruptive activity started on May 1, 2016 continuing until August 1 (Mora,137 2016). After a short period of relative quiescence from August 1 until November 2016 a new explosive eruptive phase138 took place, with eruption columns reaching up to 4 km high (Mora and Soto, 2016). An example of this activity can139 be observed on Figure 5F.140 On December 2016 a new open-vent eruptive phase started and waned on the first semester of 2019. It was141 characterized by frequent phreatomagmatic explosions, sustained ash emissions and Strombolian activity. Between142 July 2017 and May 2018 a lava pool was observed at the bottom of the active crater (Ruiz et al., 2017; Alvarado143 et al., 2020). The activity decreased during the first semester of 2019. A minor eruptive phase occurred between June144 and August 2020. Since then and until 2021 some infrequent eruptions have been observed, some of them explosive,145 that reached from 500 to 1000 m height. During this last eruptive period the active crater reached its current shape146 (Figure 4E and g). During 2021 infrequent phreatic explosions have occurred.147 Based on the relative significance of each of the eruptive phases described above and their seismological charac-148 teristics that we will develop later, we selected four specific pre-eruptive stages (PES) to be analysed in more detail,149 in order to identify and characterize the precursory seismicity. The selected pre-eruptive stages will be referred from150 now on as PES and a roman number as follows:151 • PES-I: from December, 2009 to January, 2010, which preceded the eruption of January 5, 2010.152 • PES-II May, 2013, that preceded the May 21, 2013 eruption.153 • PES-III from March to April, 2016, that preceded the open vent transition starting on April 30.154 • PES-IV September, 2016, that preceded the definitive step towards an open vent system.155 1.3. Previous seismological observations156 Here we provide a brief overview of the seismological research conducted at Turrialba volcano. Mora et al. (2001)157 carried out a short campaign on March 12-14, 2001. They recorded 30 LP events with impulsive onsets and dominant158 6 Jo ur na l P re -p ro of Journal Pre-proof March 2007 1987 November 2014 November 2015 May 2018 83°46.0'W 83°45.9'W 83°45.8'W 83°45.7'W 1 0 °1 .0 'N 1 0 °1 .1 'N 1 0 °1 .3 'N 1 0 °1 .2 'N 83°46.0'W 83°45.9'W 83°45.8'W 83°45.7'W 1 0 °1 .0 'N 1 0 °1 .1 'N 1 0 °1 .3 'N 1 0 °1 .2 'N 0 50 100 m 0 50 100 m 2001 2021 N Active SW crater: Central crater: Outer rim Feb. 2021 Mar. 2015 Nov. 2001 2010 2012 vents rim A B C D E F G rim Feb. 2021 Nov. 2001 bottom level 2011 Figure 4: Temporal evolution of the active crater at Turrialba volcano from 1987 to 2021. (A-E): Photographs of the SW active crater and the central crater between 1987 and 2018. (F-G): Google Earth satellite images provided by Maxar Technologies and CNES/Airbus, respectively. The growing rim of the active crater is shown in (B) and (F). The progressive filling of the Central crater with pyroclastic material in indicated by the bottom level shown in (F). Photos in (A) and (B) are courtesy of Gerardo J. Soto. Photos from (C) to (E) were taken by Mauricio M. Mora. 7 Jo ur na l P re -p ro of Journal Pre-proof 20 10 20 11 20 12 20 13 20 14 20 15 20 16 20 17 20 18 20 19 20 20 20 21 ~ a n n u a l is o la te d e ru p tiv e e v e n ts e ru p tiv e c y c le s o p e n v e n t: q u a s i-c o n tin u o u s e ru p tiv e a c tiv ity D e c lin e 2 0 2 0 2 0 1 2 2 0 1 4 2 0 1 8 2 0 1 6 2 0 1 5 2 0 1 3 2 0 1 0 A B C D F G E H Figure 5: Left panel: Photograph examples of the activity. (A) Left, new vents formed by the eruption on January 4, 2010 and right: the vents joined. (Photographs by Mauricio M. Mora on April 13, 2011). (B) New vent formed in January 2012. (Left photo by Mauricio M. Mora and right by Javier F. Pacheco). (C) Left: Turrialba volcano degassing on January 31, 2013 (Photograph by Mauricio M. Mora) and right: eruption of May 21, 09:02 local time, taken by OVSICORI webcam. (D) Images obtained on November 1, 2014 after the explosion occurred on October 29, 2014. (Photographs by: Mauricio M. Mora). (E) Eruption on March 13, 2015 at 11:07 local time obtained by OVSICORI webcams at Turrialba (left) and Irazú (right) volcanoes. (F) Sequence of photographs from the RSN webcam, during an eruption on September 27, 2016 at 07:22 local time. (G) Sequence of photographs from the RSN webcam of an eruption on October 5, 2018 at 19:47 local time. (H) Left: Eruption on August 1 at 07:47 local time by OVSICORI webcam. Right: Image of the active crater degassing obtained on May 13, 2020 (Photo by Mauricio M. Mora). Right panel: Phases of the recent eruption (2010-2021) at Turrialba volcano. Solid red lines and pink areas show the start and span of each eruptive phase respectively. 8 Jo ur na l P re -p ro of Journal Pre-proof frequency in the range 2-4 Hz, that were grouped in 3 families, denoting non destructive sources acting at that time.159 Tassi et al. (2004) provided a first classification and description of seismicity at Turrialba. They observed hybrid160 events with dominant frequency in the range 2.1-3.5 Hz and LP events with dominant frequency <1.8 Hz that were161 located 4-6 km beneath a 4 km elongated area comprising the summit craters and the northern flank. Martini et al.162 (2010) observed LP (0.8–1.8 Hz), tremor (2.4–3.1 Hz) and tornillo-type events. They also described two main VT163 swarms in April-May 2007 and September 2007.164 A temporary installation of a broadband seismic network of 13 stations was carried out between March and165 September 2013 (Eyre et al., 2013). LP events recorded during that time period have been analysed in several papers.166 Eyre et al. (2013) and Zecevic et al. (2016) located those events through joint location-moment tensor inversion. They167 both found an elongated epicentral distribution oriented NE-SW at the summit, located southeastward of the West and168 Central craters. Thun et al. (2015) observed step-like ground deformation near the summit in LP event recordings.169 Different models for the source mechanism of these events have been proposed. While Bean et al. (2014) suggested170 that LP events are caused by slow, quasi-brittle failure, Chouet and Dawson (2016) favored a fluid-driven crack model.171 In 2018 Jiwani-Brown et al. (2020) installed 20 temporary broadband seismic stations around the Irazú-Turrialba172 volcanic complex. Their network was supplemented with 45 permanent stations from the regional networks (OVSI-173 CORI and RSN) and was utilized to produce an ambient noise tomography.174 Since 2009 the RSN and OVSICORI observatories started to develop a network of broadband instruments at175 Turrialba volcano for monitoring purposes. In this paper we provide an overview of the seismological observations of176 the recordings of this network.177 2. Data and methods178 2.1. Seismic monitoring systems179 The seismic surveillance at Turrialba volcano began in 1990, when the Observatorio Vulcanológico y Sismológico180 de Costa Rica (OVSICORI-UNA) installed a 1 Hz Ranger SS1 seismometer in the summit. In 1996, three short-period181 vertical L4C Mark Products seismometers were added. The first permanent broadband station, CIMA, was installed182 by RSN on September 2009, and it was equipped with a 30 s CMG-6TD. On September, 2011 it was replaced with a183 Guralp CMG-3T 120 s seismometer and the code changed to CVTR (Figure 3). This recording site became the RSN’s184 reference station until the present. On the other side, OVSICORI set up VTUN as its reference station since April185 2010 (Figure 3). It was equipped with a 120 s Nanometrics Trillium Compact and a Nanometrics Taurus recorder. In186 the following years, both observatories consolidated a broadband seismic network which is currently made up of more187 than a dozen stations (Figure 2). The seismic observations are complemented with visual monitoring from permanent188 webcams deployed by OVSICORI in 2010 and RSN since 2014 and infrasound instruments (Chaparral Physics 25189 Vx) deployed in March, 2017 by OVSICORI.190 2.2. Seismic events at Turrialba volcano191 The seismicity at Turrialba volcano exhibit a wide variety of signals, the most significant ones can be described as192 follows (Figure 6):193 • Volcano-tectonic (VT): typical high frequency events (1-30 Hz) generated by brittle fracture and movement194 along faults within the volcano. Those are divided in distal (dVT) and proximal (pVT), according to the distance195 to the active crater, and in shallow (<2 km) and deep (>2 km). dVT and deep events usually exhibit clear P and196 S phases, while pVT are shallow and have emergent onsets, non-detectable S phase and lower frequency (<5197 Hz) contents than dVT.198 • Long-period (LP): These signals have an emergent onset and a clear central pulse. The frequency content is199 generally restricted in a band between 0.5 and 4 Hz and the duration is of a few seconds (<5s).200 • LP and tremor: This family is characterized by an LP event followed by a short tremor (a few minutes long).201 The tremor can have a broad-spectrum (1-20 Hz), hence it is called LP-T or a harmonic one (LP-HT) (van der202 Laat and Mora, 2017; van der Laat et al., 2018). The tremor may have also high-energy pulses (spasmodic203 character). The harmonic tremor coda exhibit a characteristic frequency gliding, very similar to the whooshes204 events observed at Arenal volcano (Benoit and McNutt, 1997; Mora, 2003; Lesage et al., 2006).205 9 Jo ur na l P re -p ro of Journal Pre-proof • Very long period (VLP): Pulses in the range from 2 and 100 seconds. A high frequency phase is usually206 superimposed.207 • Tornillo (TOR): This is a variety of seismic-volcanic signals that are distinguished by a tonal or multitonal coda208 whose frequencies are stationary along the duration of the event and whose amplitude decays exponentially209 (van der Laat, 2020).210 • Compounded tornillo (cTOR): They are characterized by an initial broad-spectrum event that often includes a211 VLP phase, followed by a multitonal coda typical of tornillo events (van der Laat, 2020). The frequency peaks212 of the coda are usually very numerous (up to 25) compared to typical tornillos. The first phase is impulsive and213 more energetic than the coda.214 • Double phase (DP): These signals are characterized by an emergent first low-amplitude phase followed by215 a high amplitude one, variable duration, a characteristic frequency range of 2-12 Hz with peak frequencies216 between 5 and 9 Hz, and are followed by an acoustic phase (Pacheco, 2018). These events started in 2015 and217 often accompanied visually confirmed eruptive activity. For the most part, the associated eruptive activity were218 ash emissions (Figure S1), although some strombolian explosions were accompanied by a similar two-phase219 signal (Figure S2). Pacheco (2018) found very low values of volcano acoustic-seismic ratio (VASR) indicating220 a greater seismic than acoustic coupling of the explosions.221 • Drumbeats (DB): Repetitive LP events with characteristic frequencies between 2-14 Hz and inter-event intervals222 in the range of 6 - 10s. Sequences of drumbeat events were observed in 2010, 2012, 2013, 2017 and 2018.223 Lesage et al. (2018) analyzed the sequence recorded in 2017. According to them, the clearest episode occurred224 on 27 January, 2017 from 03:00 to 05:45 UTC and it was composed of about 600 small earthquakes with225 remarkably similar waveforms but with varying amplitudes.226 • Tremor: It is a ground vibration that lasts tens of minutes to several hours. Its spectrum can be broad band227 (1-5 Hz) and non-harmonic (T) or more narrow-band and be harmonic, with one and sometimes two sets of228 evenly spaced dominant peaks. Mora and Pacheco (2018) and Mora et al. (2019) performed a detailed spectral229 analysis from a catalog of 2044 episodes of harmonic tremor that covers the period between 2009 and 2019.230 In the long-term, they observed a progressive increase in the number of episodes of harmonic tremor with time231 as well as significant variations in the fundamental frequency (range of 0.4 and 3.5 Hz) and in the number of232 overtones.233 2.3. Analysis of continuous recordings234 For the long-term analysis we use the real-time seismic energy measurement (RSEM) (De la Cruz-Reyna and235 Reyes-Dávila, 2001) computed from seismic records obtained at CVTR (previously CIMA) and complemented with236 those obtained at VTUN to fill the gaps (∼5%). Records were first corrected for instrument response and then band237 pass filtered between 1 and 24 Hz. The RSEM was calculated using 1-minute segments. The results were decimated238 by taking the daily median value. Finally, the data were displayed using a 5-day moving median. To correct the239 RSEM amplitudes due to different stations used, a correction factor was obtained to compensate for site and trajectory240 differences.241 We complemented the seismic observations with SO2 measurements obtained by Conde et al. (2014) and fur-242 ther completed by OVSICORI-UNA using a scanning dual-beam Mini-DOAS (Conde et al., 2014). We also used243 petrographic ash analysis measurements obtained from OVSICORI databases (Supplementary Figure S12), partially244 published in Alvarado et al. (2016a).245 2.4. Analysis of discrete events246 In order to seek for significant precursory events of the eruptive phases, a close inspection of the time domain247 records and detailed time-frequency analysis were carried out. We also used the REDpy algorithm (Hotovec-Ellis and248 Jeffries, 2016) to group the events and obtain specific catalogs.249 We focused on TOR, LP-HT, DP, HT and VT events since each of those types of signals carries significant informa-250 tion about the volcanic system. We implemented specific analysis for each type to better constrain their characteristics:251 10 Jo ur na l P re -p ro of Journal Pre-proof dVT pVT LP VLP TOR cTOR DP DB T HT LP-T LP-HT 5 s 5 s 10 s 50 s 20 s 400 s 400 s 20 s 100 s 20 s 10 s 50 s50 s 30 20 10 20 10 20 10 30 20 10 30 20 10 10 1 0.1 6 4 2 15 10 5 15 10 5 10 5 10 5 20 10 Bandpass filter (10-100 s) A B C D E F HG I J LK 20 Figure 6: Typical seismicity of Turrialba volcano. For each type of event the velocity trace (above) and the corresponding spectrogram (below) is depicted. For the cTOR event in (F) an additional bandpass-filtered trace (middle) is provided. For the double phase event (DP) both seismic (above) and infrasound (middle) traces are shown in (G). All events, except the DP event, were recorded by RSN’s CVTR station. The DP’s seismic and infrasonic records come from OVSICORI’s VTCE station. Abbreviations are as follows: dVT: distal volcano-tectonic, pVT: proximal volcano-tectonic, LP: long period, DB: drumbeat; VLP: very long period, DP: double phase; TOR: tornillo; cTOR: long period event with tonal coda; T: tremor; HT: harmonic tremor; LP-T: long period event with tremor; LP-HT: long period event with harmonic tremor. 11 Jo ur na l P re -p ro of Journal Pre-proof • TOR and cTOR events analysis: As these signals share a multitonal coda with stationary non-regularly spaced252 spectral peaks, we analyzed each peak independently. For each event, the signal was first bandpass filtered in253 the 1-30 Hz range. To extract the peaks from the Fourier spectrum, a threshold was obtained by smoothing254 the spectrum with a Savitzky-Golay filter (Savitzky and Golay, 1964) of order 1 and a 3 Hz window. Only the255 peaks whose amplitude exceeded twice the corresponding value in the smoothed spectrum were extracted. In256 addition, the search for peaks was restricted to a minimum distance of 0.3 Hz between peaks and a minimum257 amplitude of 0.04 (value normalized by the maximum amplitude).258 • Location of VLP phases of cTOR events: Because very low frequency waves are not affected by site effects,259 particle motion analysis was used to locate the VLP phase of cTOR events (Caudron et al., 2018). The signal260 was first bandpass filtered between 8 s and 16 s and then we used a 20 s window around the pulse. The261 azimuth and angle of incidence of wave propagation were obtained from the orthogonal regression adjustment262 of the movement in the horizontal plane and in the radial vertical plane, using the algorithm implemented in263 the ObsPy package (Beyreuther et al., 2010). The final solution was retrieved from the intersection of the264 propagation vectors obtained at several stations. The mean square error in each dimension is reported.265 • DP events analysis: We also carried out an specific overview of these events in order to determine which were266 associated with eruptive activity. We used the images from the permanent webcams. The events were classified267 in three types: 1) uncertain, due to adverse atmospheric conditions or the presence of a constant gas and/or ash268 plume from the active crater; 2) not eruptive; and 3) eruptive. Types 2 and 3 are certain because clear images269 were available.270 • LP-HT events: A semi-automatic processing based on the cepstrum (Noll, 1967) was implemented to seek for271 the fundamental frequency f1 of the harmonic tremor associated to these events. f1 corresponds to the dominant272 cepstrum quefrency. The analysis is performed in a moving window of 5 s with a step of 1 s. For each event we273 recorded the initial, final, minimum, maximum, average and standard deviation of f1. In addition, the number274 of harmonic overtones were manually counted.275 • Harmonic tremor: We first built a catalog with the extracted traces by manual inspection. We used the spec-276 trogram to obtain the fundamental frequency f1 at each window in order to account for the gliding. We then277 obtained the average, median and mode of the fundamental frequency, and the number of overtones.278 • VT seismicity was obtained from a joint catalogue compiled from OVSICORI and RSN data from 2010 on-279 wards. For this period the locations are good enough to be analyzed. From 2010 backwards we do not have280 good station coverage as the seismic monitoring systems at Irazú and Turrialba were incipient as we mentioned281 in subsection 2.1. The VT catalogue comprise the locations, depth, magnitude (Ml) and uncertainties coming282 from routine processing. The velocity model used for the location comes from Quintero and Kissling (2001)283 and was modified to include a layer between 0 and 4 km of elevation above sea level with a P and S-wave284 velocities of 3.2 km/s and km/s, respectively.285 3. Results286 This research was developed using two time scales. We first analyzed the seismicity in a long-term basis to287 characterize the evolution of the activity between 2009 and 2020 (Section 3.1). Secondly, we carried out a detailed288 analysis focusing on the precursory seismicity of the main eruptive stages (Section 3.2).289 3.1. Long-term analysis290 3.1.1. RSEM and pre-eruptive tremor abatement (PETA)291 To obtain a long-term overview of the seismic activity at Turrialba volcano we used the RSEM (Section 2.3)292 and compared it with the SO2 measurements (Figure 7). In general, the highest RSEM values occurred prior and293 immediately after the first eruption on January 2010 and during the period between May 2016-May 2017. In between294 those periods (years 2011-2015) the RSEM remained in lower levels. After May 2017 the RSEM decreased gradually.295 12 Jo ur na l P re -p ro of Journal Pre-proof Figure 7: RSEM and SO2 flux between 2009 and 2020. Vertical red lines and bands indicate the eruptive phases onset and time-span, respectively. Vertical blue bands indicate pre-eruptive tremor abatement (PETA). Vertical gray bands indicate periods of tremor and SO2 flux abatement that do not precede an eruptive event or phase. The horizontal line in the lower part of the plot indicates the station used for the RSEM computation as in the legend. The roman numbered labels indicate the pre-eruptive stages (PES) analysed in detail in this work. The Mw 5.5 Capellades earthquake and its April foreshock swarm are indicated by a large and a small red star, respectively. 13 Jo ur na l P re -p ro of Journal Pre-proof Figure 8: Evolution of the characteristics of the PETA periods. The abatement magnitude was computed as the difference between the maximum RSEM amplitude during the 7 days prior the PETA period and minimum amplitude during the PETA period. The Mw 5.5 Capellades earthquake is indicated by a red star. Roman number labels, vertical red dot lines and color bands are the same as in Figure 7. A relevant precursory feature identified in this analysis is a decrease in the seismic energy prior to most of the296 eruptive phases. For clarity, we will refer to this feature as a pre-eruptive tremor abatement, hereinafter abbreviated as297 PETA. Each PETA period is indicated with a vertical blue band in Figure 7 and following time series representations.298 This behaviour was observed in 9 occasions, being more conspicuous for the January, 2010 and April, 2016 eruptive299 phases. In general, we observe a decrease in the abatement magnitude, i.e. the amplitude difference between the300 maximum RSEM amplitude during the 7 days prior the PETA period and the minimum amplitude during the PETA301 period (Figure 8). On the contrary, the duration of the PETA increased up until the April, 2016 period, and drastically302 dropped later that year (Figure 8).303 In some PETA period cases the decrease was drastic and the low amplitudes were sustained prior to the eruptive304 phase onset (inverse plateau), while in others, the decrease in amplitude was progressive and constant (Figure S3 in305 Supplementary Material). Two of the PETA periods identified occurred with a considerable temporal separation (∼1306 month) before the corresponding eruptions in early 2011 and early 2012 Figure 7A.307 Since a correlation between the seismic activity and the magmatic gas flux has been established in other studies308 (Conde et al., 2014; Chiodini et al., 2017) we broadly compare the RSEM to the SO2 flux (in red in Figure 7). For309 example, for the 2014 eruptive phase, de Moor et al. (2016) reported the concomitant decrease of SO2 flux and seismic310 amplitude. In this work, we observe this behavior in all PETA periods occurring during the period with available data.311 Furthermore, we observe two other periods when RSEM and SO2 decrease, in September, 2012 and September,312 2013, but that were not followed by an eruption. These periods are indicated by the vertical gray bands Figure 7.313 Although, by the end of these periods no eruptive activity took place, some interesting seismic activity is noted. For314 the September, 2012 period, the RSEM increase was associated with a drumbeat sequence of LP events (Figure S4 in315 Supplementary Material). For the September, 2013 period, a sharp increase of TOR events was observed (Figure 10).316 3.1.2. Volcano-tectonic seismicity317 The VT seismicity beneath the Irazú-Turrialba volcanic complex was distributed in two zones: 1) a broad (∼12318 km) zone beneath the Irazú edifice between 0 and 4 km below sea level (dVT); and 2) a tighter cluster at Turrialba319 volcano between 0 and 3 km above sea level (pVT, Figure 9). Additionally, we observe that from 2015 to 2017 the VT320 seismicity progressively moved from Irazú (dVT) to Turrialba volcano (pVT) (Figure S5 in Supplementary Material).321 The temporal distribution shows two periods of major activity. The first occurred from 2012 to 2014 with a sharp322 increase in October, 2012, after the Mw 7.6 Nicoya earthquake (September 5, 2012) (Lupi et al., 2014). The second323 occurred from mid 2014 to 2017. During this period occurred the Capellades seismic sequence (main shock depicted324 as a red star in figures 2 and 9), the climax of volcano-tectonic seismicity between 2009 and 2021 (Figure 9C and D).325 14 Jo ur na l P re -p ro of Journal Pre-proof This sequence started with some foreshocks on April 24, 2016 (Linkimer et al., 2018), during a PETA prior to the326 May 2016 eruptive cycle.327 We used the relationship between the cumulative seismic moment of the VT events and the geodetic estimation328 of the intruded magma volume found by White and McCausland (2016). Based on the VT seismicity, we estimated329 a volume equal to 0.016 km3. This result is consistent with the estimation based on geodetic data carried out by330 Müller (2018) who modeled a magma reservoir beneath the Irazú volcano at 8 km depth where an intrusion of 16331 Mm3/year occurred during 2015 and 2016 (Battaglia et al., 2019). If we accumulate this rate of intrusion in the course332 of two years, it yields a total volume change of 0.032 km3, which is in the same order of magnitude of our results.333 Furthermore, this match is in the variability range of the data collected by White and McCausland (2016) and Meyer334 et al. (2021) (Figure S6 in Supplementary Material). As those authors explain, the cumulative seismic moment is335 dominated by the few larger events. Thus, the Capellades earthquake must be related to the magma intrusion in order336 to account for the geodetic estimation by Müller (2018).337 3.1.3. Long-term evolution of discrete TOR, DP and HT events338 TOR events were first observed in 2012 and their number had a significant increase during the second half of 2013,339 after the May eruption (Figure 10A). This was coincident with the RSEM and SO2 decrease period that occurred in340 September 2013 (gray vertical band in Figure 10). In the following 2014 and 2015 years the number of events341 decreased, although an important increase was observed in October 2015 (Figure 10A). Between 2012 and 2013, the342 frequency of the first spectral peak decreased from 10 to 2 Hz (Figure 10B). Then, in 2014 this frequency went up343 suddenly up to 12-18 Hz. Finally during 2015 it diminished again down to 2 Hz. In 2016 cTOR became dominant344 showing lower frequencies down to 0.6 Hz by September 2016. Preliminary estimations of the coda attenuation345 quality factor of representative TOR and cTOR events (Lesage, 2008, 2009) indicate that TOR events could have346 beeen related to the resonance of a suspension of ash in gas (Q ∼600), while cTOR events could have been related to347 bubbly water or bubbly basalt (Q ∼60-80), following Kumagai and Chouet (2000) (Supplementary Tables S1 and S2).348 DP events started in 2015 (Pacheco, 2018) and reached a maximum number in November 2015 (Figure 10C).349 Afterwards, a progressive but fluctuating decrease in the number of events can be observed (Figure 10C). Several DP350 events accompanied eruptions, becoming more usual since 2017 (Figure 10C).351 Harmonic tremor (HT) has been observed in Turrialba volcano since 2007 (Martini et al., 2010). However, our352 catalog extends from 2009 to 2018. Between 2009 and 2013 HT episodes were rare, but as of 2014 they became more353 conspicuous, particularly in 2018 when they reached a maximum (Figure 10D). In parallel, the number of overtones354 also showed a slight increase (Figure 10F). The fundamental frequency ( f1), however, fluctuated between 0.5 to 3 Hz.355 The lowest frequencies were reached during 2015 (Figure 10E).356 3.2. Short-term analysis of the studied pre-eruptive stages (PES)357 In section 1.2 we defined 4 pre-eruptive stages to be studied in detail from the seismological point of view. Each358 PES was dominated by a different type of seismic-volcanic event, but some families were present in several stages359 (Table 1). LP-HT events were observed in PES I-III, but mostly in PES-I. LP-T events were only observed in PES-II.360 The cTOR events were observed mainly in PES-III and PES-IV. Finally, DP events were only recorded on PES-III and361 PES-IV. We highlight that the most important family in number is the LP-HT, which has the highest average rate and362 is present in all stages. Also, in all of the stages a swarm of pVT events was observed towards the end of the stage.363 Table 1: Average rate (events per day) of type of event at each pre-eruptive stage PES Start date End date Duration [days] LP-HT LP-T TOR/cTOR DP PES-I 2009-12-19 2010-01-03 16 66 0 0 0 PES-II 2013-04-21 2013-05-21 29 3 42 0.2 0 PES-III 2016-03-14 2016-04-27 44 9 0 21 0.6 PES-IV 2016-09-02 2016-09-12 11 0 0 4 0.6 LP-HT events were more frequent during the PES-I. In this 16-day period a total of 1054 events occurred, with364 a maximum rate of 8 events per hour. On the other hand, PES-II and III were longer (29 and 44 days, respectively),365 but presented fewer LP-HT events (71 and 429, respectively). No LP-HT were detected on PES-IV. We statistically366 compared their characteristics per PES (Figure 11). The amplitude and the mean fundamental frequency (Mean f1)367 15 Jo ur na l P re -p ro of Journal Pre-proof Figure 9: Volcano-tectonic (VT) seismicity in the Irazú and Turrialba volcanic complex. (A): Map: Black lines are the main active faults Montero and Alvarado (1995); Montero et al. (2013); Linkimer et al. (2018). The XX’ dashed line indicates the trend of the profile in (B); (B): XX’ profile; (C): magnitude time series and cumulative seismic moment of VT events; (D): Number of VT events per month, the number of events in August 2015 and December 2016 is indicated next to each bar, in order to provide a vertical scale adjusted to the minor variations. The color of each point is dependent on the order of occurrence as shown in the colorbar. The red star in (A) to (C) indicates the Mw 5.5 Capellades main shock. The Liebres fault is indicated with a dashed line in (A) and (C). Color coding for vertical lines and bands as in Figure 7. We use a minimum magnitude cut-off of 1.2 based on the analysis of the Gutemberg-Richter Law. 16 Jo ur na l P re -p ro of Journal Pre-proof Figure 10: Evolution of the tornillo-type (TOR), double-phase (DP) and harmonic tremor (HT) events. (A): Number of TOR and cTOR events per month. (B): Frequency of the first spectral peak of TOR and cTOR events. (C): Number of DP events per month. (D): HT duration per month. (E): HT fundamental frequency ( f1). (F); HT number of harmonic overtones. Roman number labels, vertical red dot lines and color bands are the same as in Figure 7. The Mw 5.5 Capellades earthquake is indicated by a red star. 17 Jo ur na l P re -p ro of Journal Pre-proof showed a progressive decrease along each PES, while the number of harmonic overtones remained the same for PES-I368 and II, but increased in PES-III. In two instances of LP-HT events we observe intermittence of the odd harmonics369 (Supplementary Figure S13) which can be understood as period doubling (Julian, 1994; Rust et al., 2008; Hagerty370 et al., 2000).371 LP-T events occurred only during PES-II. LP-TH and LP-T share some similarities on the LP part, with a dominant372 frequency of 1.5 Hz and typical duration of 20 s. Remarkably, LP-HT occurred in the second half of the PES, when373 the RSEM started to increase prior to the May 21, 2013 eruption (Figure S7 in Supplementary Material).374 cTOR events were important precursors during PES-III and IV, when 943 and 245 events occurred, respectively.375 For PES-III, cTOR events were grouped by cross-correlation, using the REDPy algorithm (Hotovec-Ellis and Jeffries,376 2016). We first filtered the signals using a Butterworth bandpass filter between 1-10 Hz and order 4. Then, we set377 minimum correlation coefficient of 0.7 to group the events. From the 943 manually detected events, 72% (678 events)378 were grouped into 66 families. Of the grouped events, 537 events (80%) belong to six families made up of more than379 14 events each. Using a stacked waveform per cluster, the particle motion analysis yields towards a consistent source380 location beneath the central crater ∼500 m deep (Figure S8 in Supplementary material).381 We also extracted each spectral peak for the cTOR events for a detail time-frequency analysis (Figure 12). In382 contrast to the typical proportional shifting of all spectral peaks (gliding) of the HT spectral peaks, we can observe383 that those of cTOR and TOR events varied independently. Due to the narrow frequency range of variation of the384 spectral peaks (<1 Hz), the scale in Figure 12 was adjusted for each peak. The same data is displayed in a common385 frequency and amplitude scale in Supplementary Figure S9 for direct comparison between peaks. During PES-III, a386 general pattern of frequency variation can be established: a prolonged decline followed by two cycles of rise and fall.387 Mora and Alvarado (2016) reported three VT swarms that occurred on April 15, April 24, and April 26 2016,388 repectively (vertical blue lines in Figure 12). Those were concomitant with each minimum frequency of the cTOR389 peaks. For each VT swarm, RSN located at least 1 event. Although, the rest of the events were not located due390 to a poor signal-to-noise ratio (SNR) and/or network coverage, the catalog was completed using the REDPy routine391 (van der Laat and Linkimer, 2019) (Figure S10 in Supplementary material). The April 15 dVT swarm was located in392 the north flank of Turrialba volcano. The April 24 dVT swarm, was located between Irazú and Turrialba volcanoes393 and it is a foreshock swarm of the Mw 5.5 Capellades (Linkimer et al., 2018). Finally, the April 26 pVT swarm, was394 located near the active crater of Turrialba volcano at shallow depths (<1 km).395 During PES-IV, the frequency of cTOR spectral peaks continued to decrease (>0.75 Hz), especially during the last396 24h before the September 13, 2016 eruption onset (Figure S11 in Supplementary material).397 4. Discussion398 4.1. Interpretation of the pre-eruptive tremor abatement (PETA) periods399 Although not termed PETA, similar behaviour (i.e. a decrease in seismic amplitude previous to eruptive activity)400 has been reported in other volcanoes and in a variety of time-scales. For example, at Redoubt volcano, a short-lived401 (seconds to minutes) seismic pause was observed prior to some Vulcanian explosions following gliding harmonic402 tremor (Dmitrieva et al., 2013; Hotovec et al., 2013). Similarly, at the Suwanosejima volcano, some Vulcanian explo-403 sions have been preceded by a sudden cessation of tremor on a short-term scale (minutes) (Nishimura et al., 2013). In404 Telica volcano, Roman et al. (2016) showed that episodes of seismic inactivity on a medium-term scale (minutes to405 hours) systematically preceded the explosions in May 2011. In the same volcano, Geirsson et al. (2014) and Rodgers406 et al. (2015) reported a decrease in tremor amplitude and in the number of discrete seismo-volcanic events over a407 long-term period (weeks), between March and April prior to the May 2011 eruption. This last example is similar to408 Turrialba volcano, in which the PETA periods last for several days to over a month.409 In other volcanoes this tremor quiescence has been interpreted as a pause in the circulation of fluids and/or ash410 in the shallow part of the volcanic edifice due to the sealing or blockage in the conduits. The resulting overpressure411 would eventually provoke rock failure and, consequently, the eruptive activity (Nishimura et al., 2013; Geirsson et al.,412 2014; Rodgers et al., 2015; Roman et al., 2016). Particularly, for Turrialba volcano, the concomitant decrease in413 SO2 flux with PETA periods suggests the blockage in the system. The cause of this sealing can be related to the414 precipitation of minerals from rising magmatic gases (de Moor et al., 2016; Stix and de Moor, 2018). For Turrialba415 volcano, this hypothesis is supported by the collection of blocks of highly indurated hydrothermal breccia that were416 18 Jo ur na l P re -p ro of Journal Pre-proof Figure 11: LP-HT events characteristics at each pre-eruptive stage (PES). (A) Amplitude of the LP event part; (B): mean fundamental frequency f1 of the harmonic part; (C): number of overtones of the harmonic tremor part. In each case, the box extends from the lower to upper quartile values of the data, with an orange line at the median. The whiskers extend from the box to show the range of the data. 19 Jo ur na l P re -p ro of Journal Pre-proof Figure 12: Evolution of the cTOR events spectral peaks during PES III. (A): Number of events per day. (B-F): frequency variation of the first 5 spectral peaks of the coda. Each point represents a peak of each event. The spectral amplitudes of the peaks are represented by a color scale. The solid blue lines indicate VT swarms. The first two are distal swarms (April 15 and 24), while the third one is a proximal swarm (April 27). The vertical red band indicates the eruption phase that began in April 29. 20 Jo ur na l P re -p ro of Journal Pre-proof ejected during the violent explosion on October 29, 2014 (de Moor et al., 2016), which was preceded by a PETA417 period of 34 days.418 The observations here provided suggest that the sealing of the hydrothermal system was a recurrent process that419 leaded to most of the eruptive phases prior to the open-conduit stage. We observe that PETA periods became in-420 creasingly longer in duration (Figure 8) which could reflect the progressive opening of the system, since a more open421 system would require longer duration to over-pressurize. Also, we observe that the magnitude of the PETA tended422 to decrease, which could mean that required pressure difference to generate the eruption was decreasing due to the423 progressive opening of the system.424 4.2. Precursory seismo-volcanic events: cTOR, LP-HT and LP-T425 For the generation of the multitonal resonance typical of cTOR events, a closed cavity is required where a high426 pressure gradient is generated due to a high flow rate (Konstantinou, 2015; Fazio et al., 2019). According to Kumagai427 and Chouet (2000), the estimated quality factor of 60–80 for a representative cTOR event would correspond to bubbly428 water or bubbly basalt. We propose that the VLP pulse of cTOR is generated by a slow fluid movement in the shallow429 magma body. The downward polarity of this pulse could be associated with the rise of a gas pocket in the magmatic430 column, which causes a downward acceleration of the surrounding fluid. If the fluid is massive enough, it can generate431 a force that couples to the rock walls during a momentary deflation (Imai, 1983; Waite, 2015). The escaping mixture432 of magma and gas from the magma body then produces pressure transients that excite the oscillation of the cavities.433 cTOR frequency decrease was the general trend in both PES-III and PES-IV, as observed in other volcanoes434 (Torres et al., 1996; Molina et al., 2004). This could be due to a variety of factors, such as increase in the fluid435 density and/or a growth of the cavities or the increase of gas fraction in magma (Kumagai and Chouet, 2001; Molina436 et al., 2004; Konstantinou, 2015). Despite this general trend, during the two final weeks of PES-III we observe two437 cycles of frequency increase and decrease. We have already seen that the beginning of each cycle coincides with a438 VT swarm. In particular, the VT swarm occurring in April 26, 2016 is proximal and no more cTOR were detected439 after that swarm. Hence, we propose that the first two dVT swarms (April 15 and April 24, 2016, Figure 12) may440 indicate magma injection events, which in turn could increase degassing and gas temperature in the shallow magma441 column and, therefore, the acoustic velocity and the frequency of oscillation would also increase. Conversely, the last442 pVT swarm might be related to the rupture of the hydrothermal system seal, which previously conditioned the cTOR443 resonances.444 Unlike for cTOR events, it was not possible to obtain a source location for LP-HT. However, high attenuation445 observed at distant stations (VTLA and VTCV) indicate a shallow source. A possibility is that the LP phase is446 the result of an internal explosion due to the accumulation of gases in the sealed zone of the hydrothermal system447 (Figure 13B). The consequent HT could be generated by one of multiple processes: acoustic resonance of a fluid-448 filled fracture (Chouet, 1988); thermoacoustic oscillations (Busse et al., 2005; Montegrossi et al., 2019); Dirac comb449 effect by sustained repetition of a pressure pulse (Lesage et al., 2006); flow-induced vibrations to the walls of a conduit450 (Julian, 1994); the excitation of standing waves in a reservoir by flow in a coupled channel (analogous to the clarinet)451 (Rust et al., 2008); gas flow through a permeable medium (Girona et al., 2019), among others. The observation452 of period doubling in some LP-HT events (Supplementary Figure S13) points to the flow-induced vibration model453 (Julian, 1994). Notice that these observations cannot be explained by the shift from matched end (open-open or closed-454 closed) to mismatched end (open-closed) conditions, or vice versa, in a 1-D resonator since only odd overtones would455 be observed in the mismatched case. Rust et al. (2008) observed period doubling in an laboratory experiment where456 air flows up a vertical slit in a gelatin block. Thus, we prefer the interpretation in which the HT phase results from the457 high-speed flow of the gases in a narrow conduit after being released during the internal explosion (LP phase).458 During PES-II (2013) the most important precursor class was the LP-T, which share similarities with LP-HT. We459 suggest that these two type of event could share the same source but under different conditions. During this stage,460 the mean RMS amplitude of the tremor phase was 0.19 µm/s for the LP-T and 0.33 µm/s for the LP-HT events. The461 fact that the amplitudes of the LP-HT were greater than those of the LP-T could indicate that there was an increase462 in pressure and/or flow velocity that generated the tremor phase, reaching the critical pressure necessary to generate463 the harmonic signal (Rust et al., 2008). Additionally, the RSEM increased at the same time that the LP-HT events464 occurred, this could further support the idea of a general pressure increase in the system.465 21 Jo ur na l P re -p ro of Journal Pre-proof 4.3. Interpretation the evolution of the 2010-present eruption466 We propose a conceptual model for the 2010-present eruption in Figure 13. This model was derived from the467 results of this work and the geological and vulcanological information known up to this moment. We will address468 each stage of the eruption according to this model in next subsections.469 4.3.1. Slow reactivation from 1996 to 2009470 As we previously mentioned in section 1.2, Turrialba volcano featured signs of unrest possibly dating back to471 1980s, but surely since 1996. Inflation of the volcano and VT swarms reflect the early magma intrusion processes472 (Martini et al., 2010; Soto, 2012). Furthermore, significant increases in low-frequency seismicity were observed473 during this period (Tassi et al., 2004; Mora et al., 2001; Martini et al., 2010)). Finally, the surface conditions also474 underwent changes such as the increase in the flow and decrease in pH of the magmatic gases (SO2), as well as the475 opening of fissures in the walls of the West Crater (Martini et al., 2010) (Figure 13A).476 4.3.2. Annual eruptive behaviour from 2010 to 2013477 The first five years of eruptive activity at Turrialba volcano comprise 4 main eruptions that occurred at almost478 annual intervals (Figure 5).479 Our seismic records start only a few months prior to the first eruption on January 5, 2010. In this period we observe480 relatively high RSEM amplitudes (Figure 7A). Among other explanations, the first PETA could reflect the formation481 of a seal in the hydrothermal system due to mineral precipitation. In this context, a high number of LP-HT signals482 were generated (Table 1). When the pressure in the seal exceeds the resistance of the rock, an internal explosion483 occurs (LP phase) and the gas accumulated is released at high speed into a narrow conduit or fissure (HT phase). The484 late stage of this PETA period was accompanied by pVT events that reflect brittle rupture in the hydrothermal system485 (Figure 13B). The process culminated with the opening of the 2010 vent. During this first eruption, the volcano486 emitted a low content of fresh ash (∼5%) (Alvarado et al., 2016a,b), indicating the destruction of pre-existing rock to487 form the new conduit and the interaction of the shallow magmatic body with the hydrothermal system.488 Other minor ash emissions occurred in June, July 2010 and January 2011 (Section 1.2). Before the last one we can489 observe a PETA (Figure 7A) between November and December 2010. After the first eruption we observe a gradual490 decrease of the RSEM amplitude, probably related to a relative stabilization of the system after that eruption.491 The January 2012 eruption also took place after a PETA (Figure 7A). After this eruption we observe the first few492 TOR events (Figure 10) indicating the formation of a cavity that can host the oscillation of a fluid at high pressure493 (Figure 13C). Also, after this eruption we note an small increase in the number of low-magnitude VT events prior,494 and thus not induced by the Mw 7.6 Nicoya earthquake (September 5, 2012, Figure 9D).495 The PETA periods observed previous to the small eruptions of January 2011 and January 2012 show ∼1 month496 temporal separations with their corresponding eruptions Figure 7. This could be interpreted as a delay between the497 actual seal rupture at depth and the superficial eruptive process.498 Additionally, in September 2012 we observe a joint RSEM and SO2 flux pattern similar to the PETA but that did499 not result in an eruption (Figure 7B). Instead, at the end of this process, when the amplitude of the RSEM increased in500 October 2012, we note the earliest observed LP drumbeat sequence (Supplementary Figure S5) which could suggest501 the repetitive movement of an intruding magma body (Figure 13C) as has been interpreted in the literature (Matoza502 and Chouet, 2010; Kendrick et al., 2014; Bell et al., 2017).503 The May 2013 eruption precursory seismicity was characterized by: a) an overall context of low seismicity (Fig-504 ure 7B); b) a second drumbeat sequence in March 2013 (Lesage et al., 2018); c) a PETA longer in duration and505 smaller in magnitude than the prior (Figure 8), and; d) the dominance of LP-T events and only a small proportion of506 LP-HT. These observations suggest that, among different possible processes, in this case the seal in the hydrothermal507 system was not very effective, resulting in internal explosions that weren’t energetic enough to pressurize the gas that508 circulates the shallow conduits and was therefore unable to generate the harmonic resonance.509 In general, from 2010 to 2013, we can highlight the following facts: a) the RSEM progressively decreased after510 the first eruption and remained low after 2012; b) the PETA magnitude decreased with each new eruption while its511 duration became longer (Figure 8); and c) the VT seismic moment decreased until the first half of 2013. These512 observations reflect, during the first 4 years of activity, how the system progressively decompressed while losing the513 ability to build up a significant amount of pressure each time it opened a new vent. At the same time, the small magma514 22 Jo ur na l P re -p ro of Journal Pre-proof 0 100 200 300 400 Distance [m] 0 100 200 300 400 Distance [m] 0 100 200 300 400 Distance [m] pre-vent 2010 fractures gas ascent through cracks 2700 2800 2900 3000 3100 3200 3300 A ltitu d e [m a .s .l.] A ltitu d e [m a .s .l.] 2011-20132009-2010 2016 new vents D CA 0 100 200 300 400 Distance [m] 2700 2800 2900 3000 3100 3200 3300 2017 - 2018 F B E 2014-2015 Deep explosions (DP) and unblocking of conduits Juvenile ash content increase (< 20 %) 1996-2009 Reactivation Annual phreatic eruptions 2010 - 2013 Eruptive cycles Quasi-continuous activity LP seismicity Magmatic degassing (SO 2 ), pH decrease fractures opening Magmatic intrusion (dVT, inflation) Mar. 2016 June 2017 Renewed magmatic intrusion (dVT, inflation) Crater bottom levelHydrothermal systemSW crater inner wall MagmaRock LP-HT pVT pVT DP TOR cTOR pVT HT LP-T/HT LP-HT (VLP) 2010 2012 Figure 13: Conceptual model of the evolution of the Turrialba volcano shallow system. The model was conceptualized based on the observations made in this research, together with contributions from other authors mentioned in the text. Each sketch represents a SW-NE cross-section of the southwest active crater. The positions of the crater bottom (horizontal dashed lines in E and F) were obtained by Ruiz et al. (2017) and the date of each measurement is indicated. The positions of both the hydrothermal system and the magma body are approximate. Note that only for pVT and the VLP phase of cTOR events the source locations could be obtained by seismic methods. For LP-T/HT and DP events we speculate about their position based on their characteristics discussed in the text. Acronyms: DP: double-phase event; HT: harmonic tremor; LP-HT: long-period event with harmonic tremor; LP-T: long-period event with tremor; TOR: tornillo-type event; cTOR: compounded tornillo-event; pVT: proximal volcano-tectonic event. 23 Jo ur na l P re -p ro of Journal Pre-proof body continued to decompress and rise as reflected by the drumbeat swarm observed in March 2013. Additionally, the515 decrease of VT seismic moment could reflect the end of the first intrusion that reactivated Turrialba volcano (Figure 9C516 and D).517 4.3.3. Eruptive cycles behaviour from 2014 to 2016518 After the 2013 eruption the RSEM level remained low until the end of 2014 (Figure 7B). VT seismicity began to519 increase progressively since late 2014 (Figure 9D, Figure 13D). Concomitantly, the frequencies of the TOR events520 increased up to ∼10.5 Hz (Figure 10B). During October 2014 a subtle PETA is observed together with a notable521 decrease on SO2 flux and a few low amplitude HT episodes (Figure 10). These observations can be interpreted as a522 greater blocking or sealing of the system, which finally led to the violent eruption on October 29, 2014 as mentioned523 in section 1.2 (Figure 5D). These observations are in agreement with de Moor et al. (2016) who collected breccia-type524 ejecta bearing a matrix of hydrothermal minerals (e.g. silica, clays, zeolite) that support the hypothesis of a seal in the525 hydrothermal system leading to explosive eruptions. The eruptive activity lasted until mid-December, followed by a526 significant but short increase of VT seismicity (number and magnitude), HT and SO2 (Figure 7B and Figure 10D-F).527 This suggests a more open system at the same time that the injection of magma continued.528 After the last eruptions in December 2014, a relative calm of around 3 months followed during which RSEM529 remained low (Figure 7C). Eruptive activity resumed on March 8, 2015, right after the VT events and HT episodes530 increased (Figure 9D and Figure 10D). The fundamental frequencies (f1) of the HT episodes were below 1 Hz while531 TOR events frequencies dropped from 11 to 5 Hz (Figure 10A, B and E). DP events appeared, some of them, associ-532 ated to eruptions (Figure 13D). The next cycle started on October 16, 2015. In contrast to the previous cycle, it was533 preceded by a PETA and a notable decrease of the SO2 emission (Figure 7C). HT episodes had fundamental frequen-534 cies (f1) below 1 Hz (Figure 10D-F). Between the two 2015 cycles seismicity remained relatively low, with only few535 TOR events and HT episodes. Furthermore, an isolated eruptive event took place on August 15, 2015, preceded by a536 notable peak of VT seismicity (Figure 9D).537 In October occurred the highest number of DP and an increase in the number of TOR events, which had ceased538 since May of that year (Figure 10). For the latter, the frequencies decreased down to 2 Hz, (Figure 10B) which could539 indicate, among other explanations, the growth of the resonant cavities and opening of conduits due to the intense540 interaction between the magmatic gases and the hydrothermal system during the 2014 and 2015 eruptive phases.541 We interpret that DP events, which commonly accompanied discrete eruptions, were part of this process of opening542 of the conduits which began at the margin of the rising magma body in interaction with the hydrothermal system543 (Figure 13E).544 Subsequently, during the March to April 2016 PES III, there was an increase in DP events, a high number of545 cTOR and LP-HT events, and some episodes of harmonic tremor, but in less quantity than in 2015. In this case, the546 amplitudes, fundamental frequencies and frequency gliding of LP-HT events were lower than in 2009. During PES-I547 the mean f1 was 3.2 Hz, while during PES-III it was 1.6 Hz. Assuming the same resonator, this difference could be548 due to changes in various factors (fluid velocity, resonator geometry, velocity changes in the cavity, etc.). In 2010549 the first vent was barely opened (Alvarado et al., 2016a; Elizondo et al., 2019), while in 2016 there were already550 multiple widened vents. Thus, the frequency decrease could be explained by an increase of the length of the resonator.551 Similarly, the much lower frequencies of the cTOR (1.6 Hz) than the preceding TOR events (>5 Hz) (Figure 10)552 could reflect the formation of a larger cavity (Figure 13E). The PETA period associated with this eruption reached the553 longest duration (Figure 8), suggesting that the system needed a longer time to pressurize because it was almost open.554 Once the pressure exceeded the resistance of the rocks in the hydrothermal system due to the accumulation of gases,555 the seal broke releasing the accumulated material (April 27 pVT swarm). The eruptive phase began on April, 29, and556 continued intermittently until August, 1.557 After a pause, a short PETA is observed between September 1 and 13. A new sequence of cTOR occurred558 during this stage, with frequencies lower than before (∼1 Hz, Figure 10) which continued to decrease down to the559 observed overall minimum (∼0.8 Hz) prior to the September 13 eruptive phase onset (Supplementary Figure S11).560 During the eruptive phase, the RSEM showed relatively high values, comparable to those in 2010 (Figure 7). This561 activity continued until the end of November when a new short pause occurred. It was precisely during this relative562 quiescence that the M 5.5 Capellades mainshock occured, on December, 1st 2016 (Linkimer et al., 2018). This563 sequence represents the climax of the VT activity, both in number and in magnitude (Figure 9).564 24 Jo ur na l P re -p ro of Journal Pre-proof The moment tensor inversion of the Capellades main shock indicates that the dextral strike-slip double couple565 component represents 84.5% of the total seismic moment (Linkimer et al., 2018). Although, this observation makes it566 possible to rule out the significant influence of fluids on the rupture mechanism (Linkimer et al., 2018), fluids can still567 trigger the rupture via changes in pore pressure or local static stress change. Here we provide several observations568 reflecting a relationship between the sequence and the volcanic activity that support the idea that this sequence was569 triggered by the magmatic or hydrothermal fluids:570 1. The volume change in the deep magmatic reservoir modeled from deformation data (Müller, 2018; Battaglia571 et al., 2019) coincides with the one estimated from the seismic moment (Supplementary Figure S6) based on572 the relationships found by White and McCausland (2016) and Meyer et al. (2021).573 2. The position of the deep magma reservoir (Müller, 2018; Battaglia et al., 2019) is aligned in a SW-NE trend574 with the centroid of the Liebres fault activity (Linkimer et al., 2018) and the active crater of Turrialba volcano.575 3. The Liebres fault (Figure 2) (Linkimer et al., 2018) is located between the deep reservoir (5-10 km deep)576 (Müller, 2018) and the intermediate reservoir (1.65-3.95 km deep) (Badilla and Taylor, 2019), at depths similar577 to that of the intermediate reservoir.578 4. A migration of distal seismicity in the SW-NE direction occurred between 2015 and 2017 (Supplementary579 Figure S5).580 5. An increase in cTOR frequencies just after the Capellades foreshock swarm on April, 24.581 6. Changes in eruptive activity after the main shock:582 (a) An abrupt increase in fresh content in erupted ash (Supplementary Figure S12).583 (b) Quasi-continuous character of the eruptive activity.584 (c) Frequent strombolian activity (Alvarado et al., 2020).585 7. Change in tremor amplitude trend (RSEM, Figure 7). The climatic point occurred during the preceding pre-586 eruptive period (September to November, 2016). On the other hand, the Capellades main sequence coincides587 with a PETA, after which a downward trend was sustained until the second half of 2019.588 8. An increase in the amount of harmonic tremor, a pronounced increase in its frequencies and sudden decrease in589 the number of harmonics just after the Capellades mainshock (Figure 10).590 These observations support the hypothesis that the Capellades sequence (Linkimer et al., 2018) occurred, at least591 partially, in response to the stress applied by the magmatic injection. However, the geometry of the Liebres fault might592 be conditioned by the NNE regional compressive stress in accordance with the dextral fault systems proposed in the593 literature (Montero and Alvarado, 1995; Montero et al., 2013; Calvo et al., 2019).594 4.3.4. Open system and quasi-continuous eruption: 2017 - 2019595 The changes detailed above, that occurred in relation to the Capellades sequence, indicate the complete opening596 of the system. This was confirmed by the observation of an undermining of ∼120 m of the active crater, an open597 conduit evidenced by a magma body at the bottom of the crater from July 2017 to probably May 2018 (Ruiz et al.,598 2017; Alvarado et al., 2020) and by strombolian activity. The seismo-acoustic signal of the strombolian explosions599 were DP events (Figure S2), although with lower frequencies the typical values (5-9 Hz), and shorter delays between600 the seismic and acoustic signal indicating a shallower magma source.601 Another VT swarm occured in December 2017 (Figure 9) and was located SW of the Turrialba cone at sea level602 depth. This swarm culminated with the migration of dVT seismicity that was observed since 2015 with a SW-NE trend603 (Supplementary Figure S5). At the same time, the amount of harmonic tremor attained a minimum (Figure 10D), after604 which it increased until reaching its maximum between May and June 2018.605 Fundamental frequency changes during this phase could be associated to the magma column level or with vari-606 ations of the level of exsolution in the magma column (Figure 13F). During 2017, the HT displays its highest fun-607 damental frequency values (Figure 10). This could mean a decrease in the length of the conduit. Then, in April608 2018 started a decrease in the fundamental frequency. Inversely, this could indicate that as the magma level descends,609 the resonator grows and therefore the generation of harmonic tremor increases and its frequency decreases. Also, a610 change in the gas fraction could be the cause of the frequency changes observed.611 25 Jo ur na l P re -p ro of Journal Pre-proof 4.3.5. Activity decline612 During the first semester of 2019, the eruptive activity began to wane and in 2020 only isolated small eruptions613 occurred between June and July (OVSICORI-UNA, 2021). In addition, between 2019 and the date of writing this614 paper, a deflation of the Irazú-Turrialba volcanic complex was observed (OVSICORI-UNA, 2021).615 5. Conclusions616 In this paper we provided the first seismological overview of the recent activity of Turrialba volcano, as early as a617 few months prior to the first eruptive phase and until its decline in recent years (2019-2020). We analysed more than618 11 years of records in order to provide a general description of the opening process of the system after ∼150 years of619 dormancy. The main focus of this work was the identification of precursory signals that could enhance the anticipation620 of eruptive processes in a surveillance context, as well as the understanding of the geological processes involved.621 Our analysis of VT seismicity on the Irazú-Turrialba volcanic complex supports the idea that the main magma622 reservoir associated with the Turrialba activity is located beneath Irazú volcano. Furthermore, we conclude that,623 although the Liebres fault geometry might be conditioned by regional tectonic stresses, its triggering resulting in the624 largest earthquake (M 5.5) in the volcanic complex may have been due to magma dynamics.625 The recent eruption (2010-present) of Turrialba volcano has been a long-term process that evolved through dif-626 ferent eruptive phases. The pauses between those phases allowed for in-depth observation of different precursory627 processes for each phase. In the long term, volcanic tremor amplitude decrease was identified as an eruption precur-628 sory signal. In total, 9 pre-eruptive tremor abatement (PETA) periods were identified, ranging from 5 to 44 days in629 duration. Often, PETA periods concurred with a decrease in the SO2 flux. We interpret these observations a result630 of the sealing of the hydrothermal system, most probably due to the formation of hydrothermal minerals. Once the631 system was over-pressurized the seal was breached and the eruption occurred. This process accompanied most of the632 eruptions previous to the open-conduit phase, and the variations of the PETA features (duration and magnitude) reflect633 the progressive opening of the system.634 In order to better understand the process ongoing during these periods and to identify short-term precursor signals,635 we analysed in detail four pre-eruptive stages (PES I-IV). In the short term, the main precursor seismic classes were636 the LP and tremor events, both harmonic and broadband spectrum (LP-HT and LP-T), compounded tornillo-type637 events (cTOR) and distal and proximal volcano-tectonic events (dVT and pVT). LP-HT events were observed in three638 PES (I-III) and the overall evolution of their characteristics reflect the opening of the system. cTOR events were639 dominant in PES-III and IV.640 DP events played an important role in the process of opening the conduits through internal explosions. At first,641 these explosions showed minor surface expression, probably because the magma was still at depth (e.g. Tameguri642 et al. (2002) determined initial sources for explosions at Sakurajima volcano to be 2 km deep). Once the system was643 definitely open in 2017, these explosions showed a Strombolian character. It was during this open-conduit stage that644 harmonic tremor thrived due to the formation of a resonant conduit fully developed and influenced by the shallow645 magma body.646 Several lines of investigation can be taken in the future to expand our understanding of Turrialba volcano. For647 example, VT seismicity can be further relocated using double-difference techniques. Quality factor attenuation for648 TOR and cTOR events must be obtained systematically for each event in the catalog in order to track variations in the649 fluid involved.650 To conclude, we consider that this study contributes to understand the general process of the recent eruption (2010-651 2021) at Turrialba volcano. Our decadal analysis informs about the slow transition from a closed- to an open-conduit652 system. Furthermore, our in-depth analysis of pre-eruptive stages provides insights into eruption precursor processes.653 This new knowledge constitutes valuable information for the volcano monitoring community.654 Author contributions655 Leonardo van der Laat: Conceptualization, Methodology, Software, Formal Analysis, Investigation, Writing656 - Original Draft, Visualization; Mauricio Mora: Supervision, Conceptualization, Methodology, Software, Formal657 Analysis, Investigation, Writing - Original Draft, Resources, Funding acquisition, Visualization Javier Fco. Pacheco:658 26 Jo ur na l P re -p ro of Journal Pre-proof Supervision, Conceptualization, Resources, Formal Analysis, Investigation, Funding acquisition, Philippe Lesage:659 Formal Analysis, Writing - Review & Editing; Esteban Meneses: Writing - Review & Editing, Funding acquisition.660 Acknowledgments661 This work was funded by the research project Red en Sismologı́a Computacional para el Estudio de los Volcanes662 Activos en Costa Rica (113-B8-767) and partially from the research project Implementación de un sistema de mon-663 itoreo automático en tiempo real para los volcanes activos de Costa Rica (113-B9-664) and Vigilancia Sı́smica de664 Costa Rica (113-B5-704), all subscribed at the Vicerrectorı́a de Investigación, Universidad de Costa Rica. Part of665 the work was funded by the research project Vigilancia de los volcanes de Costa Rica por medio de su actividad666 sı́smica (SIA 0208-16). We also thank Geoffroy Avard for providing SO2 and ash data. Seismic instrumentation of667 the RSN and the OVSICORI is financed by the National Emergency Law Nº 8488. We thank the technical staff from668 RSN (Luis Fernando Brenes and Jean Paul Calvo) and OVSICORI (Christian Garita, Antonio Mata, Daniel Rojas and669 Hairo Villalobos) for the instrumentation maintenance. L. van der Laat was partially funded by Centro Nacional de670 Alta Tecnologı́a (CeNAT) Scholarship.671 References672 Alvarado, G.E., Barquero, R., Boschini, I., Chiesa, S., Carr, M.J., 1986. Relación entre la neotectónica y el vulcanismo en Costa Rica. CIAF 11,673 246–264.674 Alvarado, G.E., Brenes-André, J., Avard, G., Pereira, R., Galve, J.P., Campos-Durán, D., de Moor, J.M., Sánchez, R., 2021. La actividad eruptiva675 del volcán Turrialba (Costa Rica) en el siglo XIX: reinterpretación de los documentos históricos y de los depósitos. Revista Geológica de676 América Central , 1–41.677 Alvarado, G.E., Brenes-André, J., Barrantes, M., Vega, E., De Moor, J.M., Avard, G., Dellino, P., Mele, D., Devitre, C., Piazza, A., Rizzo, A.,678 2016a. La actividad explosiva del volcán Turrialba (Costa Rica) en el perı́odo 2010 - 2016. Revista Geológica de América Central 55, 7–60.679 doi:10.15517/rgac.v55i0.26965.680 Alvarado, G.E., Esquivel, L., Sánchez, B., Matamoros, G., 2020. Peligro volcánico de Turrialba, Costa Rica. Technical Report. Comisión Nacional681 de Prevención de Riesgos y Atención de Emergencias (CNE). San José, Costa Rica.682 Alvarado, G.E., Mele, D., Dellino, P., de Moor, J.M., Avard, G., 2016b. Are the ashes from the latest eruptions (2010–2016) at Turrialba volcano683 (Costa Rica) related to phreatic or phreatomagmatic events? Journal of Volcanology and Geothermal Research 327, 407–415. doi:10.1016/684 j.jvolgeores.2016.09.003.685 Avard, G., Pacheco, J.F., Fernández, E., Martı́nez, M., Menjı́var, E., Brenes, J., van der Laat, R., Duarte, E., Sáenz, W., 2012. Turrialba volcano:686 Opening of a new fumarolic vent on the southeast flank of the West Crater on January 12th, 2012, as consequence of shallow decompression.687 Technical Report. Observatorio Vulcanológico y Sismológico de Costa Rica, Universidad Nacional. Heredia, Costa Rica.688 Avard, G., Pacheco, J.F., Martı́nez, M., Sáenz, W., 2013. Reporte técnico: Emisión de cenizas al volcán Turrialba el 21 de mayo del 2013. Technical689 Report. Observatorio Vulcanológico y Sismológico de Costa Rica Universidad Nacional (OVSICORI-UNA).690 Badilla, D., Taylor, W., 2019. Perfil de magnetotelúrica en el volcán Turrialba, in: 3er Congreso Geológico, Universidad de Costa Rica, San José,691 Costa Rica. p. 35.692 Barquero, R., Alvarado, G.E., 1989. Los enjambres de temblores en el arco volcánico de Costa Rica. Technical Report.693 Battaglia, M., Alpala, J., Alpala, R., Angarita, M., Arcos, D., Euillades, L., Euillades, P., Muller, C., Narváez, L., 2019. Monitoring Volcanic694 Deformation, in: Reference Module in Earth Systems and Environmental Sciences, pp. 1–42. doi:10.1016/b978-0-12-409548-9.10902-9.695 Bean, C.J., De Barros, L., Lokmer, I., Métaxian, J.P., O’Brien, G., Murphy, S., 2014. Long-period seismicity in the shallow volcanic edifice formed696 from slow-rupture earthquakes. Nature Geoscience 7, 71–75. doi:10.1038/ngeo2027.697 Bell, A.F., Hernandez, S., Gaunt, H.E., Mothes, P., Ruiz, M., Sierra, D., Aguaiza, S., 2017. The rise and fall of periodic ‘drumbeat’ seismicity at698 Tungurahua volcano, Ecuador. Earth and Planetary Science Letters 475, 58–70. doi:10.1016/j.epsl.2017.07.030.699 Benoit, J.P., McNutt, S.R., 1997. New constraints on source processes of volcanic tremor at Arenal Volcano, Costa Rica, using broadband seismic700 data. Geophysical Research Letters 24, 449–452. doi:10.1029/97GL00179.701 Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y., Wassermann, J., 2010. ObsPy: A python toolbox for seismology. Seismological702 Research Letters 81. doi:10.1785/gssrl.81.3.530.703 Busse, F.H., Monkewitz, P.A., Hellweg, M., 2005. Are harmonic tremors self-excited thermoacoustic oscillations? Earth and Planetary Science704 Letters 240, 302–306. doi:10.1016/j.epsl.2005.09.046.705 Calvo, C., Madrigal, K., Merayo, F., Salazar, M., Fallas, C., Alvarado, G.E., Sánchez, B., Sánchez, R., 2019. Modelo volcanotectónico del graben706 cuspidal complejo del Turrialba (Costa Rica) y su relación con los colapsos sectoriales bajo un régimen transpresivo y transtensivo. Revista707 Geológica de América Central 61, 57–77.708 Caudron, C., Taisne, B., Neuberg, J., Jolly, A.D., Christenson, B., Lecocq, T., Suparjan, Syahbana, D., Suantika, G., 2018. Anatomy of phreatic709 eruptions. Earth, Planets and Space doi:10.1186/s40623-018-0938-x.710 Chiodini, G., Giudicepietro, F., Vandemeulebrouck, J., Aiuppa, A., Caliro, S., De Cesare, W., Tamburello, G., Avino, R., Orazi, M., D’Auria, L.,711 2017. Fumarolic tremor and geochemical signals during a volcanic unrest. Geology 45. doi:10.1130/G39447.1.712 Chouet, B.A., 1988. Resonance of a fluid-driven crack: radiation properties and implications for the source of long-period events and harmonic713 tremor. Journal of Geophysical Research 93, 4375–4400. doi:10.1029/JB093iB05p04375.714 27 Jo ur na l P re -p ro of Journal Pre-proof Chouet, B.A., Dawson, P.B., 2016. Origin of the pulse-like signature of shallow long-period volcano seismicity. Journal of Geophysical Research:715 Solid Earth 121, 5931–5941. doi:10.1002/2016JB013152.716 Conde, V., Robidoux, P., Avard, G., Galle, B., Aiuppa, A., Muñoz, A., Giudice, G., 2014. Measurements of volcanic SO 2 and CO 2 fluxes by717 combined DOAS, Multi-GAS and FTIR observations: a case study from Turrialba and Telica volcanoes. International Journal of Earth Sciences718 103, 2335–2347. doi:10.1007/s00531-014-1040-7.719 De la Cruz-Reyna, S., Reyes-Dávila, G.A., 2001. A model to describe precursory material-failure phenomena: Applications to short-term forecast-720 ing at Colima volcano, Mexico. Bulletin of Volcanology 63, 297–308. doi:10.1007/s004450100152.721 Denyer, P., Aguilar, T., Montero, W., 2014. Cartografı́a geológica de la penı́nsula de Nicoya, Costa Rica: estratigrafı́a y tectónica. Editorial UCR,722 San José, Costa Rica.723 DeVitre, C.L., Gazel, E., Allison, C.M., Soto, G., Madrigal, P., Alvarado, G.E., Lücke, O.H., 2019. Multi-stage chaotic magma mixing at Turrialba724 volcano. Journal of Volcanology and Geothermal Research 381, 330–346. doi:10.1016/j.jvolgeores.2019.06.011.725 Di Piazza, A., Vona, A., Mollo, S., De Astis, G., Soto, G.J., Romano, C., 2019. Unsteady magma discharge during the “El Retiro” subplinian726 eruption (Turrialba volcano, Costa Rica): Insights from textural and petrological analyses. Journal of Volcanology and Geothermal Research727 371, 101–115.728 Dmitrieva, K., Hotovec-Ellis, A.J., Prejean, S., Dunham, E.M., 2013. Frictional-faulting model for harmonic tremor before Redoubt Volcano729 eruptions. Nature Geoscience 6, 652–656. doi:10.1038/ngeo1879.730 Elizondo, V., Alvarado, G.E., Soto, D., 2019. Evolución espacio-temporal de las bocas eruptivas de los volcanes Irazú, Arenal, Turrialba y Poás en731 el tiempo histórico (Costa Rica). Revista Geológica de América Central 61, 35–55.732 Eyre, T.S., Bean, C.J., De Barros, L., O’Brien, G.S., Martini, F., Lokmer, I., Mora, M.M., Pacheco, J.F., Soto, G.J., 2013. Moment tensor inversion733 for the source location and mechanism of long period (LP) seismic events from 2009 at Turrialba volcano, Costa Rica. Journal of Volcanology734 and Geothermal Research 258, 215–223. doi:10.1016/j.jvolgeores.2013.04.016.735 Fazio, M., Alparone, S., Benson, P.M., Cannata, A., Vinciguerra, S., 2019. Genesis and mechanisms controlling tornillo seismo-volcanic events in736 volcanic areas. Scientific Reports 9. doi:10.1038/s41598-019-43842-y.737 Fernández, M., Mora, M.M., Barquero, R., 1998. Los procesos sı́smicos en el volcán Irazú (Costa Rica). Revista Geológica de América Central ,738 47–59.739 Gazel, E., Flores, K.E., Carr, M.J., 2021. Architectural and tectonic control on the segmentation of the central american volcanic arc. doi:10.740 1146/annurev-earth-082420-055108.741 Geirsson, H., Rodgers, M., LaFemina, P., Witter, M., Roman, D., Muñoz, A., Tenorio, V., Alvarez, J., Jacobo, V.C., Nilsson, D., Galle, B.,742 Feineman, M.D., Furman, T., Morales, A., 2014. Multidisciplinary observations of the 2011 explosive eruption of Telica volcano, Nicaragua:743 Implications for the dynamics of low-explosivity ash eruptions. Journal of Volcanology and Geothermal Research 271, 55–69. doi:10.1016/744 j.jvolgeores.2013.11.009.745 Girona, T., Caudron, C., Huber, C., 2019. Origin of Shallow Volcanic Tremor: The Dynamics of Gas Pockets Trapped Beneath Thin Permeable746 Media. Journal of Geophysical Research: Solid Earth doi:10.1029/2019JB017482.747 Güendel, F., 1984. Enjambres sı́smicos en el volcán Irazú. Technical Report. OVSICORI-UNA. Heredia, Costa Rica.748 Hagerty, M., Schwartz, S.Y., Garcés, M.A., Protti, M., 2000. Analysis of seismic and acoustic observations at Arenal Volcano, Costa Rica,749 1995–1997. Journal of Volcanology and Geothermal Research 101, 27–65. doi:https://doi.org/10.1016/S0377-0273(00)00162-1.750 Hilton, D.R., Ramı́rez, C.J., Mora-Amador, R., Fischer, T.P., Füri, E., Barry, P.H., Shaw, A.M., 2010. Monitoring of temporal and spatial variations751 in fumarole helium and carbon dioxide characteristics at Poás and Turrialba Volcanoes, Costa Rica (2001-2009). Geochemical Journal 44,752 431–440. doi:10.2343/geochemj.1.0085.753 Hotovec, A.J., Prejean, S.G., Vidale, J.E., Gomberg, J., 2013. Strongly gliding harmonic tremor during the 2009 eruption of Redoubt Volcano.754 Journal of Volcanology and Geothermal Research 259, 89–99. doi:10.1016/j.jvolgeores.2012.01.001.755 Hotovec-Ellis, A.J., Jeffries, C., 2016. Near Real-time Detection, Clustering, and Analysis of Repeating Earthquakes: Application to Mount St.756 Helens and Redoubt Volcanoes, in: Seismological Society of America Annual Meeting, Seismological Society of America Annual Meeting.757 Seismological Society of America Annual Meeting, Reno, Nevada.758 Imai, H., 1983. Experimental study on mechanism of implosion earthquakes during the 1973 successive eruptive activity of Asama volcano.759 Bulletin of the Earthquake Research Institute, University of Tokyo 58.760 Jiwani-Brown, E.A., Planes, T., Pacheco, J.F., Mora, M.M., Lupi, M., 2020. Using ambient noise tomography to image tectonic and mag-761 matic features of the Irazú-Turrialba volcanic complex at regional and local scales, in: EGU General Assembly 2020. doi:10.5194/762 egusphere-egu2020-17795.763 Julian, B.R., 1994. Volcanic tremor: nonlinear excitation by fluid flow. Journal of Geophysical Research 99. doi:10.1029/93jb03129.764 Kendrick, J.E., Lavallée, Y., Hirose, T., Di Toro, G., Hornby, A.J., De Angelis, S., Dingwell, D.B., 2014. Volcanic drumbeat seismicity caused by765 stick-slip motion and magmatic frictional melting. Nature Geoscience 7, 438–442. doi:10.1038/ngeo2146.766 Konstantinou, K.I., 2015. Tornillos modeled as self-oscillations of fluid filling a cavity: Application to the 1992-1993 activity at Galeras volcano,767 Colombia. Physics of the Earth and Planetary Interiors doi:10.1016/j.pepi.2014.10.014.768 Kumagai, H., Chouet, B.A., 2000. Acoustic properties of a crack containing magmatic or hydrothermal fluids. Journal of Geophysical Research:769 Solid Earth 105, 25493–25512. doi:10.1029/2000jb900273.770 Kumagai, H., Chouet, B.A., 2001. The dependence of acoustic properties of a crack on the resonance mode and geometry. Geophysical Research771 Letters 28, 3325–3328. doi:10.1029/2001GL013025.772 van der Laat, L., 2020. Análisis de los precursores sı́smicos de los periodos eruptivos de 2010, 2013 y 2016 en el volcán Turrialba, Costa Rica,773 Bachelor’s thesis. Universidad de Costa Rica , 1–155.774 van der Laat, L., Linkimer, L., 2019. La sismicidad tectónica en el edificio del volcán Turrialba durante el 2016, in: 3er Congreso Geológico, Uni-775 versidad de Costa Rica, San José, Costa Rica. URL: http://geologia.ucr.ac.cr/descarga-de-documentos.html?view=download&776 id=38.777 van der Laat, L., Mora, M.M., 2017. Evolución de las señ