Temporal and spatial evolution of the 2021 eruption in the Tajogaite volcano (Cumbre Vieja rift zone, La Palma, Canary Islands) from geophysical and geodetic parameter analyses

The eruption of the Tajogaite volcano in 2021, in the Cumbre Vieja rift zone (La Palma, Canary Islands), is the most recent but also the most damaging and longest-lasting eruption on the island since historical records began. We analyzed geophysical parameters (seismicity, surface deformation, thermal anomalies, and volcanic activity) with the aim of reconstructing the temporal development of the magma plumbing system. The seismic sequence (> 9000 events, mbLg ≤ 5.0) can be grouped into nine distinct phases, which correspond to well-separated spatial clusters and distinct earthquake regimes. A strong premonitory thermal anomaly was calculated on July 2021 in the South of La Palma. The source of the magma was within the upper mantle (40–20 km depth). A second magma budget (< 15 km depth) indicated a short-term stagnation of alkaline magma within the lower crust, near the Moho. The outline of these earthquakes delineates the magma pathway. Just before and immediately after the eruption (phases 2 and 3), the high-frequency seismicity in the mantle disappeared completely, reappearing ten days after the eruption started. The pattern of deformation began before the eruption with a breakpoint around September 11, after which, deformation accelerated sharply, followed by stabilization in phase 3 and a decrease in the following phases. The level of deformation suggested that a shallow magma intrusion was established. Based on the historical eruptive pattern, it is possible to infer two possible hazard-forecast scenarios: (a) the volcano could enter a prolonged phase of tranquility for more than one century or, (b) a new eruption could occur between ca. 2045 and 2065, with a slightly shorter duration than this recent one. Therefore, if the more pessimistic scenario is assumed, appropriate territorial planning is required.


Introduction
La Palma (28° 36′ 46″ N-17° 51′ 58″ W) and the neighboring island of El Hierro are the westernmost islands of the Canarian archipelago, which comprises seven larger islands and three smaller ones, forming a roughly east-west oriented chain from 97 to 403-516 km west of the passive margin of the African continent (Fig. 1). Both islands are active volcanoes, with the last eruption observed at El Hierro in 2011 and frequent historic eruptions documented on La Palma Klügel et al. 1999Klügel et al. , 2005 including the recent one in 2021. La Palma is the second youngest and historically the most active island in the Canarian archipelago. It accounts for most of the historic eruptions that have occurred in the archipelago over the last 600 years, and all the eruptive events have taken place along the Cumbre Vieja rift to the south of the island. The previous eruption on La Palma occurred in 1971 and is known as the Teneguía eruption, after a local aboriginal place name. The 2021 Tajogaite eruption (in allusion to the locality where the first vent was formed) is important in a local context as it provides a link between the current magma plumbing system and the systems that supplied the preceding 1971, 1949 and  Marzol 2000). C Simplified volcanostructural features of La Palma (based on Ancochea et al. 1994;Carracedo et al. 1999;Acosta et al. 2003;Camacho et al. 2009;Fernández et al. 2021), D Eruptive cones opened in the Tajogaite eruption, 2021 with the dates earlier eruptions, thus allowing us to trace the temporal seismic development of the magma plumbing system beneath La Palma.
The 2021 eruption on La Palma is the most damaging and long-lasting eruption since historical records began. According to the Copernicus (European satellite system) disaster monitoring programme (Copernicus Emergency Management Service, EMS), lava flows extended over an area of 1219 ha, affecting, or damaging 3216 buildings, destroying at least 1676 buildings registered in the Spanish cadastre, forcing the evacuation of about 7500 residents (2329 directly affected), while also severely affecting or destroying irrigation systems on several farms, 73.8 km of roads and a 370 ha. area of cropland (mostly bananas but also some vineyards and avocado plantations). The airport at La Palma operated 74% of the scheduled flights (2130 out of 2889 flights). About 34% of the cancellations occurred as a result of operational closure of the airport due to ash, and the other 66% were cancelled when airlines detected ash in the airspace. Initial estimates put losses at more than 1025 million US dollars. The government has assigned 70 million US dollars in direct aid for the local farming and fishing industries in the impacted area.
On 14 November, a 72-year-old man died from injuries sustained when the roof of his house collapsed while he was clearing ash off it. His death was the only one indirectly associated with the eruption.
The eruption of the Tajogaite volcano at La Palma attracted the attention of the world due to its wide media coverage and its effects on the population. At least 528 accredited scientists visited the volcano during the eruption to contribute their experience, learn lessons for future events and observe the successful functioning of the emergency management in order to transmit these observations to their countries.
The eruptive event at Cumbre Vieja has been the subject of various research investigations using InSAR techniques. The approaches have been highly diverse, ranging from characterizing the seismicity that preceded and accompanied the eruption (De Luca et al. 2022), characterizing the underground dynamics (Fernández et al. 2021) or determining the volume and rate of effusion during the eruption to control lava flows (Plank et al. 2023).
The main novelty of this work lies in the collection and analysis of geodetic and geophysical parameter data from a multidisciplinary perspective, leading to results with important implications for future disaster management planning. We specifically present seismological data, geothermal anomalies and recent DINSAR deformation results, using the combined findings to propose a preliminary model for the lithospheric architecture and volcanic processes at the 2021 Tajogaite eruption. The seismicity analysis led us to identify 9 seismic phases where different seismicity patterns were recognized. These phases correlated well with other volcanological (eruptive column, emission of sulfur dioxide, new vents), geodetic and thermal data, including the existence of clear premonitory anomalies two months before the eruption started. Finally, based on the recent eruptive history, two possible hazard scenarios for future eruptions are proposed. The results obtained in this work could help to establish improvements in the emergency measurement monitoring and forecasting of future eruptions.

Geological setting
The Canary archipelago comprises seven islands and several islets that form a roughly E-W oriented chain west of the African continent, emplaced since the Lower Miocene on a Jurassic oceanic crust on the African plate (Roest et al. 1992;Schmincke and Sumita 2010). A pronounced decrease in age from the highly eroded eastern centers  to the youthful westernmost islands of La Palma and El Hierro (1.7 and 1.1 Ma, respectively) reflects migration of the basaltic shield-stage volcanism that has dominated all the Canary Islands at a rate of ca. 1.6-1.9 cm/a (Carracedo 1999;Schmincke and Sumita 2010). However, it is not the case that there has been a simple age progression of volcanism in the Canary Island chain comparable to that displayed in the Hawaiian-Emperor chain, as Holocene (i.e., Gran Canaria, Fuerteventura) and historical eruptions have taken place at Lanzarote (1730-1736 and 1824), Tenerife (1704Tenerife ( -1706Tenerife ( , 1798Tenerife ( , 1909, El Hierro (2011) and La Palma (since 1585). Thus, the entire Canarian Archipelago can be considered a volcanically active area (Carracedo et al. 1998;Carracedo 1999;Schmincke and Sumita 2010). The island of La Palma in the Canary Archipelago is positioned on the slowmoving African plate, approximately 100 km north of the assumed present-day emergence of the Canary mantle plume (Hoernle et al. 1991), although a NW African suture zone has also been proposed. A consensus regarding plume versus suture zone models for the Canaries has not yet been reached, and a combination of both is likely, although the plume model is the most widely accepted (Carracedo et al. 1998;Schmincke and Sumita 2010).
La Palma is situated slightly to the east of the M25 oceanic magnetic anomalies, suggesting the underlying oceanic crust to be ca. 156 Ma in age (Watts 1994). Several aspects (xenoliths, plutonic rocks, geochemical signature) provide convincing evidence that greater thickness, lower temperature, and higher rigidity of old oceanic lithosphere underlay all the Canary Islands. The depth to the Moho beneath La Palma is thought to be around 14-15 km (Bosshand and Macfarlane 1970;Ranero et al. 1995;Llanes 2006; Schmincke and Sumita 2010).

Brief volcanological history of La Palma
La Palma is a salient example of a steep-sided composite volcano. Its volcanic edifice rests on the oceanic crust at a sea-floor depth of 3000 m, and the subaerial topographic elevations reach around 2426 m a.s.l. (meters above sea level). It is a triangular-shaped island measuring 45.5 km from north to south and 27.5 km from east to west at its widest part (Fig. 1C). It has an area of 708.4 km 2 . The geology of La Palma has been reviewed by Carracedo et al. (2001).
Volcanic activity at La Palma began with the formation of a submarine complex of seamounts and associated plutonics between 3 and 4 Ma. The oldest Pliocene rocks on the island, belonging to the uplifted submarine seamount series complex (Fig. 1), are exposed on the lower walls and floor of the Caldera de Taburiente (Staudigel and Schmincke 1984;Staudigel et al. 1986;Carracedo et al. 1999).
The submarine evolution (old volcanic series) of La Palma was followed by subaerial volcanic complexes that rest unconformably on the uplifted seamount and were emplaced during three volcanic cycles, represented by the volcanoes of the Taburiente, the smaller Bejenado volcano and the elongated ridge to the south known as the Cumbre Vieja rift, from which all prehistoric to historic rocks were issued. The younger Bejenado complex grew in the southwest-oriented landslide scar of the Taburiente volcano (Carracedo 1999;Carracedo et al. 2001).
Because the eruption centers on La Palma moved successively southwards, the volcanic edifices are aligned with decreasing age from north to south. The Taburiente caldera, which is cut by the Barranco de las Angustias canyon, is a type of flank collapse structure arising about 0.56 Ma ago, generating submarine debris avalanche deposits and sculpted 1 3 by subsequent erosion. Three older debris avalanche deposits (0.8-1 Ma) might be related to the three evolutionary stages of Taburiente volcano (Urgeles et al. 1999;Guillou et al. 2001). Following the decline of activity at the Taburiente and Bejenado volcanoes, the volcanic activity propagated southwards to form the Cumbre Vieja rift zone from ca. 125 ka onwards. The north-south elongated Cumbre Vieja volcano (25 km N-S, 19 km E-W) forms the southern half of La Palma and marks the region from which all historic eruptions have originated . García and Jones (2022) confirmed previous studies that suggested the existence of a western flank lying over collapse debris material and hyaloclastites.

Structural framework
Cumbre Vieja is the last volcanic expression of a succession of growing and overlapping volcanic centers in La Palma (Ancochea et al. 1994;Carracedo 1999;Urgeles et al. 1999).
For most of its history, from before 125 ka ago to around 20 ka, the Cumbre Vieja volcano was characterized by a triple ''Mercedes Star'' volcanic rift zone geometry, a very steep slope, and high-aspect-ratio rift zone (Fig. 1C). The three rift zones were unequally developed, with a highly productive south rift zone and weaker NE and NW rift zones: the disparity in activity was probably due to topographic gravitational stresses associated with the west facing Cumbre Nueva collapse structure underneath the western flank of the Cumbre Vieja (Carracedo 1994). From 20 ka to about 7 ka, activity on the NW volcanic rift zone diminished and the intersection of the rift zones migrated slightly to the north. More recently, the triple rift geometry has been replaced at the surface by a N-S-trending rift zone which transects the volcano, and by E-W-trending en-echelon fissure arrays on the western flank of the volcano. The NE rift zone has become completely inactive. This structural reconfiguration indicates weakening of the western flank of the volcano. The 1949 eruption near the summit of the Cumbre Vieja was accompanied by the development of a west facing normal fault system along the crest of the volcano. The geometry of this fault system and the timing of its formation in relation to episodes of vent opening, during the eruption, indicate that it is a surface rupture along a developing zone of deformation. Reconfiguration of the volcanic rift zones and underlying dyke swarms of the volcano was used to determine the present degree of incipient stage flank instability of the western flank of the Cumbre Vieja with a seaward movement ).

Recent eruptive history of the Cumbre Vieja rift zone
The Cumbre Vieja rift zone eruptions commenced with pre-historic events dating back some 125 ka. The eruptions have created cinder cones and craters and produced fissurefed lava flows that have reached the sea several times . Eruptions recorded since the fifteenth century have produced mild explosive activity and lava flows that damaged populated areas. There have been 8 previous volcanic eruptions in the Cumbre Vieja region since the fifteenth century, namely the events at Tacande (1430-1440), Montaña Quemada (1470-1492), Tahuya (1585), Martín-Tigalate (1646), the Fuencaliente (a small vent resting upon the northern flank of the San Antonio volcano) (1677-1678)), El Charco (1712), San Juan/Nambroque or Llano del Banco/Hoyo Negro/Duraznero (1949) and Teneguía (1971), making La Palma the most active of the Canary Islands (Araña and López 1974;Araña and Ortiz 1984;Carracedo et al. 2001). The geological map shows the distribution of the products of historic eruptions along the Cumbre Vieja rift zone Table 1 Historical eruptions at Cumbre Vieja volcano and brief details of relevant aspects (based on Araña and López 1974;Araña and Ortiz 1984;Romero 1989;Klügel et al. 1999Klügel et al. , 2000Carracedo et al. 2001;Johansen et al. 2005;Pankhurst et al. 2022) The duration was taken in accordance with the beginning and ending moment of the eruption rather than the calendar days Eruption  Characteristics of the historic eruptions from the Cumbre Vieja rift zone include the connection to phonolitic domes , probably due to new magma exploiting the intense fracture systems associated with the domes that act as easy pathways (e.g., Roque Teneguía in the case of the 1971 event). Another common feature is the eruption from multiple vents, generally aligned along several-hundred-meter to kilometer long fissures, usually parallel but occasionally oblique to the rift axis (Klügel et al. 2000;Carracedo et al. 2001). Notably, the higher-altitude vents are predominantly characterized by explosive phreatomagmatic eruptions, whereas the lower-altitude vents usually produce larger scoria deposits and lava flows. The historic eruptions of the Cumbre Vieja rift, including the 1971 Teneguía eruption, produced basanitic to phonolitic lava flows and associated scoria deposits Klügel et al. 1999;Barker et al. 2015;Pankhurst et al. 2022). Table 1 summarizes the historical eruptions at Cumbre Vieja volcano and brief details of relevant aspects. However, in the literature (even in the same reference) there are certain contradictions as regards the date, duration, and lava field area of the eruptions. The data on eruptions of the XV century are also supported by radiocarbon dating, hence, the wide range of estimated duration. Figure 2 shows the frequency, duration and the area covered by lava flows in the different eruptions since the XV century.
From the historical record of eruptions, several interesting aspects of the Cumbre Vieja volcano emerge (Table 1 and Fig. 2): (a) eruptive activity over the last 6 centuries has been concentrated in the southern sector of the island. (b) The volcano had two important recent eruptive cycles from 1430 to 1712, then a prolonged eruptive pause of almost 237 years, before reactivating in the XX and XXI centuries. (c) Excluding the eruptions of the XV century (lacking historical details), the periods of inactivity between eruptions prior to the XIX century became increasingly shorter, while in the XX and XXI century they have been longer. (d) Between 1585 and 1971, the duration of the eruptions became increasingly shorter (and the affected area smaller), a pattern that was broken by the recent eruption.

The 2021 Tajogaite eruption
Several seismic crises have been observed since 2017 along with background surface deformation on the island identified during the 2006-2010 period using geodetic techniques such as Global Navigation Satellite System (GNSS), Advanced Differential Synthetic Aperture Radar Interferometry (A-DInSAR) and microgravimetry. These findings reveal that, although there were no significant associated variations in gravity, there was clear surface deformation, which was spatially and temporally variable (Escayo et al. 2020). An increase in the 3 He/ 4 He ratio was also observed toward the southern part of the Cumbre Vieja volcano. It was interpreted as the result of residual degassing of volatiles from magma bodies stored at lithospheric levels beneath the southern part of the volcano (Padrón et al. 2012).
The two seismic swarms were registered in October 2017 and February 2018 with b-values of 1.6 ± 0.1 and 2.3 ± 0.2, respectively, suggesting an intense magmatic fluid contribution, gas and/or magma. Both swarms were linked to changes in gas emissions. This dataset was compatible with a stalled magmatic intrusion at a depth of ca. 25 km, with an estimated volume of between 5.5 × 10 −4 and 3 × 10 −2 km 3 (Torres- González et al. 2019). The geodetic techniques employed only a few months before the eruption allowed the fluid migration induced by magma injection at depth to be tracked and the existence of dislocation sources below Cumbre Vieja volcano to be identified (Fernández, et al. 2021).
The 2021 eruption was preceded by 14 seismic crises since 2017, focused in the area where the eruption eventually took place. After a multi-year period of elevated seismicity, the eruption began at 13:10 (local time, UTC) on September 19, 2021, in a forested area called Cabeza de Vaca of Montaña Rajada, in Las Manchas section of the municipality of El Paso, on the SW flank of the island, about 20 km NW of the 1971 eruption and at an elevation of about 920 m a.s.l. Two fissures opened (200 m apart) and multiple vents produced lava fountains, ash plumes, and lava flows that traveled to the W. By November 10, the highest point of the main cone was 1130 m a.s.l. Lava continued to flow west through pre-existing lava channels and tubes, over older flows, and occasionally formed new branches. Occasional short-lived overflows of lava ponds occurred at the main crater vents. The eruption continued during November 11-16, characterized by Strombolian explosions and lava fountaining from multiple vents, with advancing and sometimes branching lava flows and daily ash emissions. On October 1, a new vent opened and then, on October 4, the main cone of the volcano collapsed, increasing the lava flow. On October 9, part of the north face of the volcano collapsed again, causing the lava within the crater to spill out. A strong, 6 km high column was observed on December 12, and small landslides were reported at the main active cone. After 12 weeks (85 days and 8 h) of Strombolian activity and emission of lava, the eruption ceased rather abruptly on December 13, 2021, at 22:21 local time. At the end of the event, at least a dozen vents of composite cinder cones were preserved.
In general, the eruption was characterized by Strombolian eruptive activity to lava fountain (500 m height) from a group of at least nine major volcanic vents aligned along a 600-m-long fissure and dominated by emitting tephra, producing a welded spatter blanket around the vents, distal tephra deposits, and mainly clastogenic lavas and lava flows throughout the duration of the eruption. The eruption was also characterized by rhythmic explosions that produced ash, lapilli, and bombs in the vicinity of the vents. Part of the cone collapsed several times on October 4, 10-11; on October 10, large accretionary lava balls of about 10 m in diameter were generated and rolled downslope. In some cases, the eruptive column rose to altitudes as high as 5 to 7 km. According to the 'Instituto de Vulcanología de Canarias (Involcan)' and 'Instituto Geográfico Nacional of Spain (IGN)', the total tephra volume was greater than 15 × 10 6 m 3 , corresponding to a Volcanic Explosive Index of 3.
The lavas reached the sea on September 28, November 10 and November 15 and formed a lavafall on the cliff and a coastal composite lava delta, 550 m wide by 500 m from the coast, called Isla Baja. The lava flow field was made up of overlapping flows that were emitted mainly as lava flows, with several kipukas (free terrain surrounded by lava) in the lava field (about 3.5 km at the widest point to about 6.2 km long); in total covering an existing land area of 1241.1 ha, and a new ~ 69 ha piece of land which extends out to sea. The average thickness of the lava flow is estimated at around 12 m, with a maximum thickness of up to 70 m (PEVOLCA bull., December 25, 2021), while the main vents reached a height of around 200 m above the pre-eruptive land surface and 1100 m a.s.l. with about a dozen preserved craters. The total volume of emitted materials is estimated to be ≥ 205 × 10 6 m 3 , of which the lava flows are estimated to be about 190 × 10 6 m 3 . The evolution of the lava flow velocity was estimated using Sentinel-1 imagery along with an offset tracking technique by Kyriou et al. (2022). The maximum lava flow velocity was calculated to be 2 m/day.

Methodology and data
A multiparametric analysis was developed to identify behavior patterns of different geophysical and geodesic parameters in the pre-eruptive and eruptive period. The parameters analyzed and the subsequent methods of analysis are classified into four blocks: (1) seismicity, (2) deformation of the surface, (3) thermal anomalies and (4) other characteristic parameters of the eruption. A flowchart of the applied analysis methodology is shown in Fig. 3. In the following subsections, the methodology and data used in each block of analysis are described.
The data used in our analysis are publicly available, published by different agencies and institutions. One limitation in this regard is that certain data that would have been of interest to analyze are not yet publicly available, such as relocalised hypocentral locations, although the IGN has shared preliminary results. Future research will attempt to address this information gap.
On the one hand, we try to establish possible prognostic indicators for future eruptions and, on the other, infer the magma pathway from its upper mantle source to the surface, including the possible magma storages during the eruption. The specific methodology for each block is described below.

Seismicity
The seismological data have been extracted from the IGN seismic catalog, published on its website (https:// www. ign. es/ web/ sis-area-sismi cidad). The complete sequence is composed of 9088 events, covering an interval of 193 days. A first descriptive analysis of the distribution of epicenters and hypocenters of the complete series was carried out, with the corresponding histograms of distribution of the number of events by magnitude and depth. Next, the temporal evolution of the seismicity was analyzed, based on the visual inspection of the number of events per day, identifying inflection points in the corresponding histograms. As regards the changes in frequency of events, 9 phases were identified, and each phase was analyzed independently to identify the distribution of magnitudes (mbLg), depths, intensities, and recurrence laws. To characterize the recurrence of each series, the well-known Gutenberg-Richter law was estimated, which establishes a linear proportion between the number of earthquakes that exceed each degree of magnitude and the magnitude itself. The seismicity of each series was then fitted to the expression: where N is the cumulative number of earthquakes of magnitude bigger than M, "a" represents the earthquake rate, and "b" is the slope of the relationship, this being a critical parameter describing the ratio between the number of earthquakes of smaller and larger magnitude in the series (Gutenberg and Richter 1944) (Fig. 4).

Surface deformation
The surface deformation was monitored during the pre-eruptive and eruptive phases identified in the seismic analysis, through multitemporal InSAR monitoring (MT-InSAR). Taking into account the nature of the terrain and the slope of the western side of La Palma, the SBAS (small-baseline) approach was selected. The degree of uncertainty using SBAS methods is similar to that of Persistent Scatterers Interferometry (PSI techniques), in the order of ~ 1 mm/year (Lanari et al. 2004).
LICSBAS ) is an open-source InSAR time series analysis package that integrates with LiCSAR (Gonzalez et al. 2016;Lazecky et al. 2020), developed by the Centre for Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET) at the University of Leeds. LICSBAS pre-processing allows users to save time and disk space while obtaining the results of InSAR time-series analysis in the frame of JASMIN, the UK's collaborative data analysis environment (Lawrence et al. 2013). During the LiCS-BAS processing, interferograms with many unwrapping errors are automatically identified by loop closure and removed. Reliable time series and velocities are derived with the aid of masking using several noise indices. The straightforward implementation of atmospheric corrections to reduce noise is achieved with the Generic Atmospheric Correction Online Service for InSAR (GACOS).
The reference points for the analysis were the GNSS stations located in the area, namely LP02, LP03 and LP04 and their displacements published online (IGN 2021) ( Table 2). Copernicus Emergency Management Services provided updated vector files of the eruptive vents and graded extension of lava flows (Copernicus EMS 2021). For the analysis of deformation of La Palma Island, we worked with the frames provided by COMET (2021) using the available data from Sentinel-1 A/B. The frames 060A_00001_030604 and 169D_00001_020800 include 1222 interferograms obtained between January 2019 and December 2021, 522 from Sentinel-1 images on the ascending track 60 and 700 from the descending track 169, respectively (COMET 2021).
A first processing (P1) was carried out "in the middle of the eruption," to check the feasibility of quick MT-InSAR processing during emergencies. This processing included available COMET interferograms from September 1, 2020 to October 14, 2021. They were processed using LICSBAS on a 100*100 m grid, taking the GNSS station LP02 (IGN 2021) located in the northern block of the island as reference point. The stability of LP02 allowed the velocities and deformations to be compared with the GNSS stations. Ascending and descending geometries were considered and processed, performing vectorial decomposition to determine net vertical movements as proposed by García et al. (2018) based on the developments by Fialko et al. (2001) and Diao et al. (2016).
A second processing (P2) was carried out with the updated interferograms released by COMET on December 2021, providing a longer record on the ascending Sentinel-1 track 60 and descending Sentinel-1 track 169, up to the 1222 interferograms calculated from January 2019 to December 2021 (COMET 2021).
The "area of maximum deformation" was delineated by selecting the points with velocities over 20 mm/year for the monitored periods, resulting in 2838 points from the 100*100 m grid in P1 processing (September 1, 2020 to October 14, 2021). We calculated the main statistics of deformation (mm) for each date: average, maximum, minimum, standard deviation, 95% percentile and 5% percentile. Velocities were calculated for the phases identified from the seismic activity. The deformation estimated with MT-InSAR was validated with the time series of displacements of available GNSS stations LP03 and LP04, taken from official plots published at the IGN website (2021) using WebPlotDigitizer 4.5 (Rohatgi 2021).

Thermal anomalies
The thermal activity in a previous period of the volcano eruption in La Palma was studied. For this purpose, spectral anomalies in the thermal range were computed using satellite images from the ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) sensor and remote sensing techniques. The method applied addresses the unsupervised technique of target detection, also called anomaly detection. Since this technique assumes no prior knowledge about the target or the statistical characteristics of the data, the only available option is to look for objects that can be differentiated from the background. The goal of this technique is to find a reduction in background information that allows the detection of atypical points or outliers. In our case, we search for pixels in the multispectral images corresponding to surface areas that have very different thermal responses compared to their natural surroundings. In this regard, several strategies for anomaly detection based on image spectrometry (Koltunov and Ustin 2007) have been developed in the thermal range as an anomaly detection technique for early detection.
There are several, different methodologies and algorithms for anomaly detection that automatically identify these atypical points, allowing a determination to be made regarding their consideration as either interesting or non-interesting information. Expert-based methods, in contrast, determine the information through previous knowledge of the research subject. Points identified as atypical in anomaly detection should be evaluated based on the background information to determine their suitability for inclusion in further analyses. The RX algorithm (Reed and Xiaoli 1996), adopted as standard by the scientific community, for example, essentially calculates the Mahalanobis distance from the test pixel to the mean of the background.
Remote sensing pre-processing techniques have been used to extract information from satellite data and to generate ASTER (Hook et al. 2001) mosaics so accurate spatial coregistration between the different scenes becomes crucial for subsequent analysis of thermal anomalies. Pre-processes included geometric correction, radiometric correction, water masking and a number of arithmetic and statistical image analyses, as well as a MNF algorithm (Minimum Noise Fraction) to reduce noise from single ASTER images.
The single ASTER scenes are 1750 × 1483 pixels, with a spatial resolution of 15 m for VNIR (Visible and Near Infrared) bands, 30 m for SWIR (Short Wave Infrared) bands to 90 m for TIR (Thermal Infrared) bands, and radiometric resolution of 8-bit, 8-bit to 12-bit, respectively. The 15 ASTER image channels are grouped in a port of 4 channels in the VNIR (0.52-0.86 µm), a port of 6 channels in the SWIR (1.60-2.43 µm) and 5 channels between 8125 and 11,650 µm in the MIR and TIR. For the ASTER images used for this analysis, the SWIR port was not available. We used four ASTER scenes acquired through the USGS Earth Explorer (http:// earth explo rer. usgs. gov/).
A variation of the Projection Pursuit (PP) algorithm (Malpica et al. 2008), widely proven to be an efficient spectral anomaly detection method in different multi-and hyperspectral sensors, which uses a thermal subspace, was applied to all ASTER mosaic data. The method is implemented by first defining a TIR spectral subspace, then reducing the dimensionality and extracting end members that are projected on the thermal background, and finally searching for the best projection by maximizing the thermal index THIN (Rejas et al. 2014) (2).
where conv = image convolution median filter. Lµi = radiance at sensor for wavelength.
In this method, the first variables of the Principal Components Analysis (PCA) particularized for thermal spectrum wavelengths define the subspace of the thermal spectral background. The Principal Components (PC) are then projected in the thermal subspace before maximizing THI, adopting the RX standard in its global variant. For ASTER images, we calculated a thermal index THI, profiting from the separability between covers in the emissive spectrum of the sensor.
We generated image convolutions using a median filter. These were subsequently used to establish a ratio between 8.63 µm and 11.31 µm wavelengths, weighted/corrected by the ratio between each channel gain. We established thresholds for the resultant variable to highlight detected pixels as possible anomalies.

Other parameters
Finally, the temporal evolution of other variables was also analyzed, taking account of the values per day for: number of earthquakes with felt intensities, maximum microseismic intensity; maximum magnitude, maximum and minimum depth, SO 2 emissions; height of eruptive column. The opening of eruptive vents was also considered. A final representation of all of these variables, normalized between 0 and 1, is produced in order to examine possible correlations. All of this latter information came from the website of the IGN (https:// ign-esp. maps. arcgis. com).

Analysis of geophysical and geodetic parameters
The evolution of different geophysical and geodetic parameters was analyzed to identify behavior patterns before and during the eruption, which may provide some information to understand the dynamic processes that mark the eruptive activity. Possible correlations between parameters that could provide indicators for future forecasts and that could contribute to volcanic monitoring and surveillance were also examined.
In a first step, seismicity was analyzed, covering a period prior to the beginning of the seismic crisis that triggered the eruption, up until the eruption was officially declared ended. The deformation for these same phases was then monitored using DINSAR. There was also a follow-up of other characteristic parameters of the eruption, such as the height of the eruptive column, emission of sulfur dioxide, volume of lava and opening of eruptive vents. The analysis was completed with the detection of thermal anomalies two months before the beginning of the eruption.

Evolution of seismicity
The evolution of seismicity was analyzed before and during the eruption, covering the period from June 11, about 3 months before the beginning of the pre-eruptive phase that began on September 11, to December 20, when the activity had begun to decrease, entering the final phase of the eruption. The magnitudes of the seismic events are in the range mbLg [1.3-5.0] and depths z [0-46.9 km]. This number of events allows a descriptive analysis to be performed of the temporal distribution of different seismic parameters. Figure 4 shows the epicentral and hypocentral distribution of the events for the total sequence analyzed; the histograms of magnitude and distances of the events are also shown.
Several important facts are observed in this figure: (1) Seismicity spread through the southern half of La Palma, with three depth ranges in the recorded events: superficial crust (z ≤ 5 km), deep crust (5 < z ≤ 15 km) and mantle (z > 15 km). From here on, we refer to these three ranges as SC, DC, and M, respectively; (2) The magnitude histogram shows a practically normal distribution, with a mean value around 3 and a deviation σ ± 2; 3) In the depth histograms, three modes can be seen: a first environment at z [10-15 km] with a dominant frequency, with about 6500 events, and another two with less frequency in the ranges of z ≤ 5 km and z [30-38 km], with approximately 1000 events in each of these ranges.
In an initial processing, the histograms of the number of earthquakes per day were represented together with the temporal evolution of the maximum and minimum magnitudes and the depths associated with these magnitudes, considering only earthquakes with mbLg > 1. The results are shown in Fig. 5, where the cumulative seismic moment is also represented.
Through a visual analysis of the graphs, it was possible to identify 9 phases with different seismic activities, the first phase being prior to the triggering of the pre-eruptive swarm, which illustrates the background seismicity. Phase 2 can be considered a precursor, which covers the period from September 11 to 19, the day on which the eruption began. Phase 3 covers the ten days following the beginning of the eruption and is characterized by a sharp decrease in seismic activity, which begins to rebound again in phase 4, when two series clearly appear at average depths of 10 and 35 km. Both series are maintained until the end of the period covered in the analysis, although in this final period, five more series are identified based on the variation in the number of earthquakes per day.
The representation of the epicenters in the identified phases and the corresponding foci or hypocenters are shown in Fig. 6. In this figure, the coexistence of the two series at different depths is clearly perceived from phase 4. It is also possible to appreciate how the patterns of seismicity (in plan, epicenters) in phases 2 and 3 present a different orientation from those of the rest of the phases, particularly noteworthy is the fact that the pattern of phase 1, prior to the triggering of the pre-eruptive swarm, resembles that of phases 4 onwards, when deep seismicity has started to appear. It is also noteworthy that in phases 2 and 3, just before and after the beginning of the eruption, the most superficial volcano tectonic seismicity appears (z < 5 km) and the seismicity of the mantle disappears.
In Fig. 6 and in the rest of the figures presenting the temporal evolution of the variables, the same color code is used to identify the phases.
In order to explain the distribution of the depths found, the crust model for La Palma is considered in accordance with the depth to the Moho beneath the island, which is thought to be around 14-15 km (Bosshand and Macfarlane 1970;Ranero et al. 1995;Llanes 2006;Schmincke and Sumita 2010). Figure 7 shows the distribution of depths and magnitudes for the complete period analyzed, from June 2021, together with the model of the crust and upper mantle. The foci are then located in the different phases, visualizing the seismicity associated with the shallow crust (z ≤ 5 km), deep crust (5 < z ≤ 15 km) and mantle (z > 15 km). It is noteworthy that in phase 1 foci also appear in the mantle, as in phases 4 onwards. However, in phases 2 and 3, just before and immediately after the beginning of the eruption, the volcano tectonic seismicity of the mantle disappears completely, reappearing ten days after the beginning of the eruption. In both phases, the most superficial activity is registered in the upper crust, something that is obviously explained by the opening of magma-emitting vents. The appearance of surface activity coincides temporarily with the disappearance of deep volcano tectonic seismicity.
A descriptive analysis of the seismicity of each phase was then carried out, separating the intermediate and deep seismicity in the series of phases 1 and 4 onwards. The modeling of the series is beyond the scope of this work and will be the subject of future research. It is then simply a question of following the evolution of the magnitudes against the depth of each series and of estimating the Gutenberg-Richer law that characterizes the recurrence in each case (Gutenberg and Richter 1944).
For this purpose, the seismicity of each series was adjusted, representing the number of earthquakes that exceed each interval of magnitude (Eq. 1), considering intervals of 0.1. The model parameters "a" and "b" were determined in the fit, the latter of these two parameters being a critical parameter that describes the size distribution of events.
A high/low b-value indicates a relatively larger/smaller proportion of small events. In particular, the b-value can give an indication of the stress conditions in the crust (Scholz 1968;Wyss 1973) and identify volumes of active magma bodies (Wiemer et al. 1998) or roots of regional volcanism (Monterroso and Kulhanek 2003). The b-value is also used to distinguish purely tectonic (b ≅ 1.5-1) from volcano-tectonic earthquakes (b > 1.5), which are mainly characterized by hydraulic fracturing of the host rock induced by overpressurized magma. The determination of the b value for volcano seismicity is an important task for characterizing different phases of seismicity and thus, changes in the dynamic processes. Once the parameters a and b of the different series identified have been estimated, their temporal variations will be analyzed in the successive phases of the eruption. The results are shown in Sect. 4.1.2.

Analysis of each seismic phase
The results of the analysis of the series in the 9 phases identified so far are presented below, indicating the number of events in each range of magnitude and depth and the values of the parameters a and b of the Gutenberg-Richter law that characterize the recurrence in each case.
Phase 1 (June 17-September 10) Background seismicity In this phase, seismic activity was only registered in the mantle, with 49 events in a range of low magnitudes, mbLg [1.4-2.5] and depths of [20-37 km]. A trend of increasing magnitude with depth was observed during this phase.
Phase 2 (September 11-19) Pre-eruptive-eruption onset Within this phase, two subseries were identified, with events in the shallow crust (0 to 5 km) and deep crust (5 to 18 km), respectively, with the seismicity of the mantle disappearing. Figure 8 shows the best fit of this series and of all the others analyzed according to the Gutenberg-Richer law, differentiating the most superficial seismicity (z ≤ 5 km) from the intermediate depth (5 > z ≤ 15 km) and, where appropriate, the seismicity of the mantle (z > 15 km). In the 9-day interval that this phase lasts, 295 events were recorded in a range of magnitudes mbLg [1.3-3.9] and depths z of [0-5 km] and 1118 events in a range of mbLg [1.3-3.6] and z from [5.1-17 km]. The best fit according to the Gutenberg-Richter law for each subseries is Log N = 5.02-1.3 M for the shallowest and Log N = 6.63-1.5 M for the deepest.
Phase 3 (September 20-30). Initial eruption and decrease in seismicity This phase covers a period of ten days, just after the beginning of the eruption, perceiving a clear decrease in seismicity, although the pattern observed in the previous phase was maintained, with two subseries in shallow and deep crust, respectively, while the mantle remains dormant. This phase lasted ten days, one more than the pre-eruptive phase, during which 31 events were recorded in a range of magnitudes mbLg [2.2-3.9] and depths z of [0-5 km] and 80 events in a range of mbLg [2.3-3.6] and z of [6.8-15.7 km]. The best fit according to the Gutenberg-Richter law for each subseries are Log N = 3.9-1.1 M for the shallowest and Log N = 4.9-1.2 M for the deepest. The a and b parameters of both substrings are slightly lower than those of the previous phase.
Phase 4 (October 1-18). Reactivation of seismic activity In this phase, the activity was reactivated, clearly changing the pattern of seismicity with respect to the two previous phases, with two events being noteworthy: (1) disappearance of shallow activity (z ≤ 5 km) and (2) appearance of deep volcano tectonic seismicity of the mantle, with z > 15 km. In addition, the activity in the deepest crust was maintained, between 5 and 18 km, this activity being much more intense than in the previous phase. This phase lasted 18 days, almost double that of the previous phase, with 1223 events recorded in a range of magnitudes mbLg [2.3-4.0] and depths z (in km) of [6.7-17.5 km] and 225 events in a range of mbLg [2.4-4.5] and z (in km) of [19-42 km]. There were a total of almost 1500 events, compared to 110 in the previous phase. The best fit according to the Gutenberg-Richter law for each subseries are Log N = 9.9-2.3 M for the intermediate depth and Log N = 4.8-0.9 M for the deepest. The parameters a and b of both subseries are clearly different from those of the previous phase, with almost double the value of b for the crust seismicity, which is noteworthy, and a low value of less than 1, for the seismicity of the mantle, where the highest magnitudes were recorded up to 4.5 (mbLg).
Phase 5 (October 19-November 4). Intensification of seismic activity In this phase, the activity intensified even more, both in number of events and magnitude, although the two subseries of deep crust and mantle remained. The duration of this phase was 17 days, practically the same as the previous reactivation phase, with 2339 events being recorded in a range of magnitudes mbLg [2.4-4.1] and depths z (in km) of [8-17.5 km] and 139 events in a range of mbLg [2.7-4.9] and z of [19-41 km]. It is noteworthy that the number of intermediate depth events almost doubled compared to the previous phase, but the mantle events reduced by almost half, although the maximum magnitude of these increased to 4.9. In total, almost 1000 more events were recorded than in the previous phase. The best fit according to the Gutenberg-Richter law for each subseries is Log N = 9.2-2.0 M for the intermediate depth and Log N = 4.06-0.67 M for the deepest. The parameters a and b of both subseries are clearly different from those of the previous phase, and the fact that the values of a and b are quite similar to those of the previous phase is noteworthy, b being somewhat lower for the seismicity of the mantle, which was by then below 0.7. This is explained by the increase in the maximum magnitude, which already reached 5 (mbLg).
Phase 6 (November 5-13). Decrease in seismic activity In this phase, the activity decreased again, which led us to believe that this could have been the beginning of the end of the eruption. The seismicity pattern was maintained, although with lower intensity with respect to the previous phase, maintaining the series in the deep crust and the mantle. It lasted for 9 days, practically a third of the duration of the two previous reactivation phases. In total, 394 events were recorded in a range of magnitudes mbLg [2.4-4.0] and depths z Phase 9 (December 3-20). Decrease in activity and end of eruption The activity began to decrease from December 3, which can be considered the beginning of the end of the eruption. Our analysis continued until December 21, when the activity was practically null, although the Pevolca Emergency Plan decided to wait until a few more days of inactivity had passed before establishing the official end to the eruption. Finally, the announcement of the end of the eruption was made on December 25. As regards the seismic activity in this final phase, it is worth highlighting that 33 events were recorded in a range of mag- If the period from September 11 to December 20, 2021, is analyzed, that is, from phase 2, just before the start of the eruption, a common cyclical sequence of activity is observed: triggering-decrease-reactivation-intensification. Table 3 shows a summary of the seismic parameters found in the different phases between June 11 and December 20, 2021.

Variation of the b parameter
The Gutenberg-Richer laws for the different phases are shown in Fig. 8, differentiating the shallowest seismicity (z ≤ 5 km) from that of intermediate depth (5 < z ≤ 15 km) and deep seismicity (z > 15 km). Important variations were found in the three depth ranges. The most superficial series (z < 5 km), which occur in phases 2, 3 (before and after onset of eruption) and 9 (final), present values of b between 1 and 1.5. Those with events in the deepest crust that appear in all phases except the first (background seismicity) present values of between 1.5 and 2.3, except for the phase that started just after the beginning of the eruption, the b value of which declines to 1.17. The lowest values are found in the mantle series, below 1.29, and less than 1 in four of the series. It is in these series where the highest magnitudes were recorded, reaching mbLg equal to 5 in six 6 cases.
Regarding the values of the b parameter, a substantial difference can be seen in the series of the shallow crust, deep crust, and mantle. In the case of the former, they range between 1 and 1.3, in the deep crust the values increase in the range of 1.5-2.3, and decrease notably in the mantle, remaining below 1, except for series 1 and 8, which are precisely of lesser and greater activity, respectively. The average value of b for shallow crust earthquakes in phases 2 and 3 is 1.2, which is typical of tectonic earthquakes.
The seismicity of the deep crust, with an average depth of 12 km, presents notably higher values of b, with an average of 1.8, and it appears in all phases, except the  initial background phase. These earthquakes reflect ruptures which are being generated at this level during magma migration. The area where magma initially intruded into the lower crust is characterized by high fracturing, thus favoring anomalously high b-values, greater than 1.5 (Ibañez et al. 2012). The seismicity of the mantle, at a depth of around 35 km, presents drastically lower values of b, with an average of 0.8, reflecting an increase in the maximum magnitudes. In the eruption of the Tajogaite, six earthquakes of mbLg 5 were recorded, a magnitude significantly higher than that recorded in the shallow and deep crust.

Deformation analysis
Ground deformation was analyzed in a first "fast response" processing (P1) for the period September 1, 2020, to October 14, 2021. This processing highlighted an area of maximum deformation (ca. 31 km 2 ) in the SW slope of the Cumbre Vieja complex delimited by coordinates 17. 9196 W,28.5213 N and 17.8576 W,28.6123 N (WGS 84), with uplift velocities over 20 mm/year. The area is framed by the former Todoque village (North) and the GNSS stations LP03 (East) and LP04 (South) (Fig. 9A).
The deformation pattern was analyzed during the pre-eruptive and eruptive phases, emphasizing the behavior of the area of maximum deformation with higher uplift velocities than the rest of the island during the monitored period (September 2020-October 2021) (Fig. 9A). The deformation trend in this area was an uplift that began to increase shortly before the eruption (mid-August) with a breaking point around September 11, after which, deformation accelerated sharply reaching a global maximum by the beginning of October (Fig. 9B). A vectorial decomposition of displacements was performed by Detektia Earth Surface Monitoring, integrating the results of Sentinel-1 ascending and descending tracks. The decomposition shows a net uplift movement of the points located at the identified deformation area (Fig. 9A). The vectors of displacement of GNSS stations were drawn in Fig. 9A based on data presented in Table 2 (IGN 2021).
The average vertical deformation time series were calculated for all the points included in the surrounding 500 m buffer of the GNSS stations and were compared with the recorded displacements at these stations. LICSBAS processing provided few points with deformation time series around LP03 station due to decorrelations, with the nearest point at 427 m (ellipsoidal length). Figure 10 shows the comparison of the deformation series of LP03 and LP04 stations with the average vertical deformation of points located within the 500 m buffer (7 points for LP03 and 70 for LP04).
The second LICSBAS processing (P2) allowed the pattern of deformations during the identified seismic phases up until December 19, 2021, to be analyzed (Fig. 11A). Note the sharp increase in deformation parallel to the peak of seismicity around September 11, 2021, (Phase 2), eight days before the eruption (September 19). The last phases showed fluctuations up until the extinction of the volcanic activity.
In order to better characterize the pattern of deformation, the velocities of deformation in LOS direction were calculated for each of the identified phases. This analysis allows the characterization of deformation during the pre-eruptive and eruptive periods in relative terms (Fig. 11B). The results show a net increase in velocities in phase 2, with relatively stable phases 3 and 4, followed by a reactivation in phase 5. Phases 7 and 8 showed fluctuations up until the end of the eruption.

Thermal anomalies
Data from satellite sensors were analyzed in a global, multiscale spatial and temporal strategy to detect thermal activity. We used radiance images from the ASTER sensor, acquired on the 17 April 2017 and the 17 July 2021 on the La Palma Island to compute thermal anomalies.
The results are consistent with previous events in the Turrialba volcano in Costa Rica (Rejas et al. 2014) as regards the fact that thermal spectral anomalies detected previously at Turrialba and recently in La Palma show a strong spatial and temporal correlation with the volcanological activity in both volcanoes. We have applied a color table to the grayscale thermal anomaly image (see Fig. 12) to emphasize the main thermal anomalies computed. In this case, we can see the strong thermal anomalies for July 2021 in green, orange, and yellow, respectively, in descending order. Blue, black and white anomalies are not significant in this color table. These strong thermal anomalies in green and orange must be understood as very different thermal responses, from a statistical perspective, to the rest of the entire island of La Palma. We can observe that the strong thermal anomalies calculated for July 17, 2021, correspond to the South of the dividing line of La Palma island and are congruent with the seismic activity corresponding to the first phase in June-July 2021 (see Fig. 12A) of the period analyzed. In contrast, the results for April 17, 2017, show no strong thermal anomalies (see Fig. 12C). In Fig. 12D, strong differences between thermal spectral anomalies between July 17 (2021) and April 4 (2017) are shown in red and could reflect their relationship with the beginning of activity in 2021. The main thermal spectral anomaly between these dates according to data from ASTER (in relative units, it should be noted) is in the North-South profile, shown in yellow, at the site of the volcano eruption, where a strong thermal anomaly difference can be observed for the volcanic cone of Tajogaite. This finding is congruent with the Phase 1 seismic temporal analysis (Fig. 12B). These results could prove the suitability of the anomaly-detection strategy for monitoring states previous to volcanic activity. Fig. 11 A Temporal evolution of the pattern of deformation in LOS direction (mm): average, 95% and 5% percentiles for the deformation area processed using LICSBAS with automatic reference point. B Average and standard deviation of the deformation velocity (mm/day) in LOS direction for the seismic phases 1 3

Other parameters
The evolution of other parameters which are characteristic of the eruption has also been analyzed in order to explore possible correlations. Figure 13 shows the histograms of the seismic events felt by the population, together with the temporal evolution of maximum intensity, height of the eruptive column, and SO 2 emission, highlighting differences among the phases identified in the seismicity analysis. Sulfur dioxide emissions are first registered two days after the start of the eruption, first showing a sharp rise and then falling off during phase 3 as seismic activity declines. They then rise again during reactivation phase 4, remain stable during intensification phase 5, and fall again when the decline in seismicity begins. From November 10 onwards, these emissions were low. However, the height of the eruptive column shows no correlation with seismicity.
From the graphs above, it is clear that: (a) few earthquakes were felt during seismic phases 3 and 6 (b) the maximum seismic intensity per day increased slightly over the time Fig. 12 A Thermal anomalies in relative units from ASTER July 17 (2021) 07/17/2021 (left) with ENVI color table "beache", detail in red frame for the volcanic cone in Cumbre Vieja area (B), strong thermal anomalies in green, orange and yellow, respectively. C Thermal anomalies from ASTER April, 17 (2017) (right) with the same color table, no strong thermal anomalies detected. D ASTER thermal anomaly differences (in relative units) between July, 17 (2021) and April, 4 (2017), strong thermal anomalies in red (right), and North-South profile in yellow, strong thermal anomaly difference for the volcanic cone of Tajogaite in red frame analyzed, (c) the height of the eruptive column presents a general pattern of decrease with a few peaks, (d) SO 2 emissions reached a peak in phase 5, corresponding to a notable peak in the number of earthquakes felt per day. Note that the representation of all the parameters begins in phase 2, just prior to the eruption, although in the case of height of eruptive column and SO2 emission there are no data until phase 3, from the beginning of the eruption

Discussion
From the multiparametric analysis presented in this work, interesting results are derived that shed light on the dynamics of the eruptive process that took place and highlight some conclusions for future forecasts. In addition to seismicity, thermal anomalies and deformation, other variables were analyzed in this work, such as the height of the eruptive column and the emission of sulfur dioxide. The temporal evolution of all the variables analyzed is shown in Fig. 14, in which all of them are represented with values normalized between 0 and 1.
The analysis of the seismicity covering a period of 6 months, of which three are prior to the beginning of the eruption, made it possible to identify 9 phases with different activity patterns that can be associated with different processes of the eruption dynamics, such as magma movement and stress due to volcano loading and flexure.
The first phase illustrates a background seismicity from June to September 11, 2021, when the seismic swarm that triggered the eruption, a week later, was released (phase 2). Phase 3 was characterized by a drastic decrease in seismicity just after the onset of the eruption and from then on, phases of activation, intensification and weakening alternated until the end of the eruption, both in frequency and magnitude (Fig. 13). The next sequence of seismic activity was observed from phase 1 to the end of the eruption: triggering (September 11)-decrease (September 19, beginning of eruption)-reactivation (October 1)-intensification (October 19)-decrease (November 5)-reactivation (November 14)intensification (November 26)-decrease and final phase (December 3).
The final decline of the activity began on December 3, accompanied by a decrease in sulfur dioxide and magma, which indicated the end of the eruption. The phases of decreased activity (phases 3 and 6) lasted approximately three times less than those of activation and intensification, which meant that in the three months that the eruption lasted, seismic activity was practically constant, with almost 850 earthquakes felt by the population.
Three clusters of volcano tectonic earthquakes at different depth ranges were identified, located in the superficial crust, at z ≤ 5 km, deep crust at z between 5 and 15 and upper mantle at z between 30 and 40 km. The former appear only in phases 2 and 3, before and after the onset of the eruption, and the other two deeper clusters occur in all phases. Additionally, another cluster at z of 20 km appears in phase 8. It is noteworthy that during the phases in which deep seismicity appears at more than 35 km depth, the superficial seismicity of less than 5 km disappears. This happens specifically in the pre-eruptive phase and just after the onset of the eruption (phases 2 and 3), suggesting a relaxation in the volcanic system. It should also be noted that the opening of eruptive vents coincided with specific decreases in seismic activity, which can be explained by a relaxation of the pressure exerted by the magma. The pattern of seismicity on the ground in phases 2 and 3 was very different from that of the other series, including the first one prior to the swarm that triggered the eruption. The epicenters show a clear trend toward the southwest in this case. However, from phase 4 onwards, when the reactivation of the seismic activity began, the deep series of more than 35 km appeared, producing the strongest earthquakes in the entire sequence, with magnitudes greater than 3. Therefore, the biggest events were also the deepest. The earthquakes with depths of less than 10 km were, in general, of magnitudes less than 3.
Our results partially agree with those of D' Auria et al. (2022), who identified three distinct structures: a shallow localized region with less hydrothermal alteration at depths of under 3 km, a spatially extensive, consolidated, oceanic crust extending to ~ 10 km depth and a large sub-crustal magma-filled rock volume intrusion extending from ~ 7 to 25 km depth. According to our observations, this intrusion could reach a depth of ~ 30 to 35 km.
The distribution of the b values with depth may be related to the stress state in the shallow crust, deep crust and mantle. In the series where superficial seismicity appears (2, 3 and 9), the value of b in shallow crust is lower than that in deep crust, which indicates greater fracturing of the rocks at less depth in the same time intervals. The average value of b in surface crust in phases 1 and 2, just before and after the eruption, is 1.2., which is typical of tectonic earthquakes. This value can be explained by the highly fractured rock, generated by the migration of the magma toward the surface, where the pore pressure increases, thus decreasing the effective stress, which generates many shallow earthquakes of low magnitude, with a value of b appreciably greater than 1. The decrease in the value of b from the deep to the superficial crust has also been observed by Mezcua and Rueda (2023), who interpret this decline as a sign of magma cooling during its ascent.
The b value in all the series identified throughout the eruption in the deep crust (15 km) fluctuates between 1.17 and 2.31 and is considerably higher than the values presented by the simultaneous series of the mantle, which vary between 0.67 and 1.29. The values obtained here are consistent with those estimated by Mezcua and Rueda (2023), who found variations in the value of b in the interval from 9 to 15 km depth, with an average of 1.5. The authors associate this characteristic with a semi-rigid framework or dike network instead of a magmatic chamber.
This deep seismicity of the mantle presents certain similarities with that observed in the Kilauea volcano (Hawaii), where deep earthquakes with high magnitudes have been recorded, this seismicity being consistent with stresses from volcanic loading and flexure acting on preexisting faults. According to Wolfe et al. (2004), the stresses from magma movement on preexisting faults in regions with high differential ambient stresses may additionally help trigger mantle earthquakes. At Kilauea, it has also been suggested (Eaton and Murata 1960;Klein et al. 1987;Ryan 1988;Tilling and Dvorak 1993) that earthquakes define a pipe-like magma conduit that extends through the lithosphere to depths as great as 60 km, thereby delineating the magma pathway from its source to the surface (Wright and Klein 2006;Matoza et al. 2020). One of the main features of the relocated seismicity at Kilauea is that a cluster of earthquakes near a 30 km depth collapses to reveal a horizontally aligned fault zone (Wolfe et al. 2003(Wolfe et al. , 2004, while in other regions of Hawaii there are also many clusters of earthquakes parallel to the southeast coast at depths of 25-50 km, indicating the existence of several tectonic fault zones in the mantle throughout this region, which are consistent with stresses from volcano loading and flexure acting on preexisting faults. This could also explain the deeper earthquakes of greater magnitude that occurred in the mantle in the eruption of Tajogaite. In this discussion section, it is worth highlighting certain differences from the eruption which took place ten years ago on the island of Hierro, also in the archipelago of the Canaries. The Hierro eruption occurred undersea and began on October 10, 2011. This episode was linked to almost 10,000 located events. During the syn-eruptive phase and after the end of the eruption in March 2012, the seismic activity continued for several years (2012-2014) with variable intensity. Over the whole period more than 22,000 events were recorded and located by the volcano monitoring network of the Spanish Instituto Geográfico Nacional deployed in El Hierro. Three seismic phases were identified in this case: (1) Phase 1, from July 1, 2011, to September 20, 2011, includes events located directly underneath El Hierro at depths of about 8-13 km, which have the typical appearance of volcanic earthquakes; (2) phase 2, between September 20, 2011, and October 15, 2011 when the epicentres migrated south of the El Golfo and (3) Phase 3, between October 15, 2011, and January 1, 2012, stands out almost as an independent sequence. The values of parameter b were 2.25 in the first phase, 1.34 in the second and 1.12 in the third (Ibañez et al. 2012).
Finally with regard to seismicity, it is worth mentioning a possible correlation detected between the occurrence of earthquakes and the confining tidal stress, according to Miguelanz et al. (2023). The authors find that tidal stress variations contribute to the onset of the eruption and that certain explosive events, typical in Strombolian type volcanoes, seem to occur preferentially when the tidal stress rate is high. Their analysis supports the hypothesis that tides can modulate earthquake activity in volcanic areas, especially during magma migration at shallow depth. It is possible that their conceptual model could be applied in the Canary Island archipelago.
Furthermore, the characteristics of multispectral and thermal satellite data, both spatial and spectral, for La Palma Island, were studied using different remote detection methods with ASTER data in 2017 and 2021. Two different methods were tested in La Palma to compute spectra-anomaly detection with similar spectral background, with no prior knowledge of the thermal spectrum of the volcanic cone in Cumbre Vieja. It was found that a high-quality, higher spectral resolution in the thermal spectrum improves accuracy in the monitoring and detection of thermal anomalies. The spectral mixing and sources of image error, such as clouds, associated directly with the spatial resolution in the ASTER satellite data used, has a significant impact on the characterization of the spectral background, and thus on the spectral anomalies calculated.
A spatial and temporal correlation between seismic activity and thermal anomalies calculated from ASTER satellite scenes was detected in phase 1 of the La Palma Island eruption. However, this relationship is strongly dependent on the spectral background calculated as well as on effects associated with the volcanic activity. Spectral anomaly detection is used in a wide variety of applications and is important because it sometimes constitutes the first step in multispectral satellite image analysis. Any significant deviations from the spectral background of the image can cause alerts that attract the attention of human operators to the anomalies and drive them to further investigate the peculiarities in the image in relation with the volcanic activity. The thermal THIN index was chosen because of previous research experience reporting its usefulness in detecting outliers.
There was a notably strong deformation pattern at the beginning of September, which increased abruptly after day 11 (phase 2) when the seismicity also increased. Both phenomena indicated the proximity of the eruption, associated with the rapid ascent of magma. After the eruption, there was a stable phase, after which the deformation decreased slightly, although without reaching the level of phase 1, indicating a potential emplacement of magma at shallow levels. This agrees with observations by Kyriou et al. (2022), who stated that: (a) the island was uplifted before the eruption due to magma influx; and (b) the eruption was followed by a relatively rapid subsidence, arising from the removal of magma and gas. This study was based on the interpolation of the multitemporal interferograms, as well as the pre-and post-eruption LOS displacements. These patterns also coincide with the pre-eruptive and co-eruptive deformation processes studied by De Luca et al. (2022).
In addition to seismicity, thermal anomalies and deformation, other variables were analyzed in this work, such as the height of the eruptive column and the emission of sulfur dioxide. The temporal evolution of all the variables analyzed (starting in phase 2) is shown in Fig. 14, in which all of them are represented, normalized between 0 and 1.
The normalization of each variable x i was performed according to the following equation: where x i represents the corresponding variable; i.e., maximum magnitude, deformation of the ground, and so on. x min is the minimum value of the variable in the total sequence. x max is the maximum value of the variable in the total sequence. This Figure allows a direct comparison of the indicators and the behavior of the system during the eruptive period. It should be noted that the opening of the vents marks the changes of phase between Phases 2-3 (first eruption) and 3-4. The rest of the vents (C7, C8 and C9) were opened during Phase 5 (October 19-November 4), at the same time as an intensification of seismic activity in magnitude (D), intensity (E) and number of events (Fig. 14).
Phase 2 (11-19 September) presents changes in most parameters, with increasing seismicity and deformation, leading to the beginning of the eruption. The pattern of emissions and deformation seems to be related. Following the maximum deformation, there is a period of increasing seismic depths and magnitudes, coinciding with the reactivation of deep series in the crust and mantle (Phases 7 and 8).

Final remarks and future lines of research
Finally, several noteworthy facts and future lines of research can be established. The source of the magmatic plumbing system beneath La Palma was originally established at between 45 and 15 km depth (Barker et al. 2015, and references therein), which is the case of the Tajogaite eruption, where most of the deep seismicity was located between 40 and 20 km, with prolonged magmatic storage within the rigid upper mantle. Seismicity and surface deformation data indicate that a batch of alkaline magma was emplaced at shallower depth, creating a new reservoir at a depth of about 10-15 km, where magma evolved until the initiation of the eruption, suggesting a short-term stagnation of the magmatic system within the lower crust near the Moho. These earthquakes outline a pipe-like magma conduit that delineates the magma pathway from its upper Mantle source to the surface. These seismicity clusters are similar to those reported during 2011-2012 submarine eruption off the coast of El Hierro (Ibáñez et al. 2012;Martí et al. 2013).
According to the literature, an incipient stage of flank instability may exist on the western flank of the Cumbre Vieja with a seaward movement . It is possible Fig. 15 Lithospheric E-W section of La Palma and the multiple stages of magma transport (magma plumbing system) based on 2021 seismicity. Several magma storages and bodies can be delineated through a zigzagging pathway from the upper Mantle source to the surface. The numbers in the cross section indicate the different depths of major magma and tectonic (faults) features. The source of the magma at 40 km is beyond the E-W cross section that dike emplacement and rifts correspond to the relaxation and accommodation that occurs due to slip on deep faults at the contact between the volcanic edifice and the seafloor sediments or at even deeper levels. The lithospheric deformation induced by the loading of La Palma volcano, plus the listric fault systems at the contact between the volcanic edifice and the seafloor sediments, could explain the persistence of volcanism along the N-S Cumbre Vieja rift zone. Similar persistence of the volcanic process was proposed for Hawaiian volcanic rift zones (Dieterich 1988;Thurber and Gripp 1988). The focal mechanisms analyzed by Del Fresno et al. (2022) were mostly of the strike-slip type, with a clear thrust component in some deeper events. These solutions could shed light on the tectonic model in this case.
A possible interpretation of the results found in the seismicity analysis leads to the establishment of a preliminary model for the process of magma ascent from the mantle to culminate with the eruption at the Tajogaite in 2021, shown in Fig. 15, in which, several magma storages and bodies can be inferred.
Recalculated hypocentral solution imagery released by the IGN shows a clear constraining of the hypocenter clouds to well-delimitated ramps and structures. Modeled in Fig. 15 and shown to scale as a W-E cut and a S-N cut, mantle earthquakes between depths of 40 km to 29 km delimitate a steeply inclined ramp structure of about 50° inclined toward the East, shown at (A), and active during most phases. An intermediate magma ponding level (B) appearing during phase 8 can be inferred as a planar structure inclined southward at about 40° between 24 and 21 km depth. The underplating body at the Moho horizon appears as two distinct but connected horizons and phases. A first horizon occurs with rapid progression through the crust, identified by a steep and constrained plane that trends northeast from the moho level at about a 45° angle to a depth of about 4 km during phase 2 (D). Further underplating at the moho level appears at phase 3 and onwards as a tubular structure with vertical development between 11 and 14 km depth, modeled at (C). Finally, a cloud of shallow hypocenters at about 3 km depth is inferred to be a surface holding chamber (E), from which a near-vertical dike transports magma to the eruption site.
Our model is consistent with the one proposed by Mezcua and Rueda (2023), who report an ascending path of magma to the surface along two paths, namely, one in an almost vertical direction and the other ascending with a 45° dip, reaching the location where the eruption occurred. According these authors, the changing path to the surface was interpreted as a dike event, and the two zones of simultaneous seismicity in the mantle and deep crust are attributed to the opening of vertical conduits for rapid magma ascent toward the subcrustal reservoir thought existing dikes. Magma ascent from the deeper layer in the mantle to the shallow reservoir may have occurred though dikes opened in the previous activity.
This model also agrees with the plumbing system proposed by Dayton et al. (2023), who interpret a scenario where older melt ascended through dikes that cut through the rigid lithospheric mantle (causing earthquakes at a depth of ~ 30 km), from a deeper storage zone.
Another compatible interpretation is given by Luca et al. (2022), who first identified a sill-like source in the pre-eruption phase, relevant to the seismicity interval from 8 to 16 September. This source is interpreted as the effect of the temporary accumulation of magma during its transport toward the surface. The two main swarms from September 12 to 15 and from September 18 to 19 accompanied the propagation of the sill-like source and dikes that formed the magma conduit of the eruption.
A strong thermal anomaly was detected just two months before the beginning of the eruption, which coincided temporarily with the appearance of seismicity in the mantle, between 20 and 37 km deep, in phase 1. After that, the seismic swarm was triggered, when the activity of the mantle disappeared and the most superficial activity appeared, coinciding with a strong increase in deformation. This pattern could be interpreted as a possible precursor of future eruptions, a strict monitoring of the three variables analyzed being recommendable.
Last but not least, as stated in the section on recent eruptive history, a probability exists that the volcano will reactivate again around the middle of this century. Thus, based on Fig. 2 and Table 1, two possible scenarios can be inferred: (a) that the volcano could eventually enter in a prolonged phase of eruptive tranquility, as occurred between 1712 and 1949, or (b) that a new eruptive phase may occur at some time between ca. 2045 and 2065, but with a slightly shorter duration than the recent eruption. Hence, it is better to contemplate the more pessimistic of the scenarios, in which adequate territorial planning is of the essence. It could be of particular use as a tool to re-evaluate the previous long-term study of the expected volcanic hazard at La Palma (Rodriguez 2015), and to improve forecasts and degree of precision with regard to areas which may be affected by the next eruption. Thus, it may provide important information to be taken into consideration in the territorial planification process.
Because the 2021 eruption was preceded by background surface deformation on the island according to data obtained during the 2006-2010 period, with 14 seismic crises since 2017, along with changes in the gas emissions (Padrón et al. 2012;Torres-González et al. 2019;Fernández, et al. 2021) focused in the area where the eruption eventually took place, it is clear that volcano monitoring for the next eruption must be undertaken in a way that ensures an alert can be activated sufficiently in advance, updating the hazard management and population preparedness plans based on the present experience.
are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.