The 2021 volcanic eruption at La Palma, Canary Islands, marked the island's most voluminous eruption in history, exceeding 0.2 km³ of extruded magma and necessitating the evacuation of approximately 7000 residents. Prior to this event, understanding of Cumbre Vieja's feeding system was limited, relying mainly on petrological and gravity studies. The eruption, however, generated an unprecedented amount of seismic data, providing a unique opportunity to map magma pathways and reservoirs. The Canary Islands, situated on the passive NW margin of the African plate, are of volcanic origin. La Palma, the second youngest island, has a history of volcanic activity marked by a north-to-south migration, culminating in the Cumbre Vieja volcano. Previous studies suggested pre-eruptive magma storage at depths of ~15-26 km, possibly extending to 50 km, with magma ascending to an accumulation zone at 7-15 km during eruptions. Long-term seismic precursors before the 1949 and 2021 eruptions hinted at progressive magma accumulation. The 2021 eruption, uniquely monitored by an improved seismic network, presented an opportunity to investigate these processes with unprecedented detail. The Instituto Geográfico Nacional (IGN) seismic network, enhanced in recent years, recorded anomalous activity beneath Cumbre Vieja beginning in October 2017, with multiple seismic swarms occurring at depths of 20–35 km. These swarms, accompanied by changes in gas emissions, pointed towards magmatic activity and unrest at depth. The pre-eruptive period culminated in a more intense swarm at shallower depths in September 2021, leading to the eruption on September 19, 2021. The eruption lasted for 85 days and 8 hours, exhibiting variable explosive behavior and significant lava flow.

Previous research on La Palma's volcanism primarily relied on petrological and gravity studies to infer the characteristics of the magma plumbing system. Petrological studies of lavas and xenoliths suggested pre-eruptive magma storage in upper mantle reservoirs at ~15-26 km depth, and possibly down to 50 km. These studies also indicated that during eruptive episodes, magma ascends forming dikes and sills, temporarily stagnating at an accumulation zone at 7-15 km depth. Gravity studies provided additional constraints on the subsurface structure, but lacked the detailed resolution to directly map magma pathways. Long-term seismicity before previous eruptions suggested progressive magma accumulation in pre-eruptive reservoirs months to years before eruptions. However, the 2021 eruption offered the first opportunity for comprehensive local monitoring, allowing for a more detailed investigation into the magma plumbing system and its dynamics during the eruption.

The study utilized seismological methods to analyze the volcanic reactivation, improving upon the preliminary IGN catalogue. A relative location method based on waveform cross-correlation was used to relocate 8488 earthquakes associated with the 2021 unrest. This improved the accuracy of hypocenter locations and spatial resolution. Furthermore, a moment tensor (MT) catalogue was developed for 156 events, providing information on the source mechanisms and stress field. The relative relocation was performed using the HypoDD algorithm, minimizing residual travel-time differences between earthquake pairs. Different subsets of seismic stations were selected for relocation based on the temporal evolution of the network and the characteristics of the earthquakes. Waveform correlation on P and S phases was performed using time windows of 2 and 3 seconds respectively, with earthquakes pairs weighted according to their waveform correlation. For moment tensor inversion, the Grond software was used, fitting P and S waveforms in the time domain. Data were manually reviewed to ensure quality. Synthetic seismograms were computed using a regional velocity model. The inversion aimed to resolve centroid depth, time, and the six independent moment tensor components. Stress inversion was performed using a least squares method to determine the direction and relative magnitude of principal stresses. A density-based clustering algorithm (DBSCAN) was employed to classify relocated seismicity based on hypocentral locations and moment tensor similarity. The eruption tremor was analyzed using the Real-time Seismic Amplitude Measurement (RSAM) from continuous seismic data, while GNSS data were processed to obtain daily time series of deformation.

The analysis revealed two main seismicity clusters during the co-eruptive phase: a shallow cluster (9.9-12.7 km depth) and a deep cluster (32.8-38.1 km depth). The shallow cluster displayed a roughly circular shape and contained multiple seismogenic volumes and subclusters. The moment tensor solutions for the shallow cluster revealed four families, primarily characterized by oblique to strike-slip faulting with positive compensated linear vector dipole (CLVD) and negative isotropic components. The deep cluster showed more spatial heterogeneity with several small seismogenic regions. Moment tensor solutions for the deep cluster also indicated four families, mostly with strike-slip to oblique mechanisms, with some exhibiting thrust or normal components. The spatial clustering analysis and moment tensor inversion results highlighted the heterogeneity of the seismicity, suggesting a complex interplay between magma movement and stress changes. Stress inversion results indicated stress rotations in both clusters, suggesting the influence of magma reservoirs and their drainage on the local stress field. Pre-eruptive seismicity, from 2017 to August 2021, revealed multiple swarms at 18–32 km depth, coinciding with the pre-eruptive magma storage zone, suggesting progressive formation of a mushy reservoir. The co-eruptive seismicity showed a distinct lack of earthquakes in this intermediate depth range, indicating an open magma pathway with efficient magma transfer. The temporal evolution of seismicity suggests four stages: 1) eruption onset with decreased seismicity and high-amplitude tremor; 2) shallow cluster activation with subsidence and lateral deformation; 3) deep cluster activation, suggesting delayed triggering by shallow reservoir drainage; and 4) eruption end with seismicity ceasing first at depth, then at shallow levels. The study also observed earthquake doublets, interpreted as evidence for pulse-like magma flow. Focal mechanism reversals, observed in both shallow and deep clusters, were attributed to stress heterogeneities introduced by the depleting reservoirs. The overall findings suggest the presence of three distinct magma storage levels: a shallow reservoir (9-13 km), an intermediate mushy reservoir (18-32 km), and a deep, more distributed reservoir (33-38 km).

The findings provide a comprehensive model of the complex magma plumbing system beneath Cumbre Vieja. The integration of seismological data with geodetic and petrological data supports the existence of three interconnected reservoirs, contrasting with previous models proposing a single reservoir. The study reveals a dynamic interplay between magma accumulation, withdrawal, pressure changes, and seismicity. The observed temporal changes in seismicity, including the initiation and progression of the shallow and deep clusters, reveal the dynamic response of the magmatic system to pressure changes and drainage. The absence of seismicity in the intermediate reservoir during the co-eruptive phase suggests an efficient, high-temperature pathway for magma transfer. The stress rotations observed in the moment tensor solutions further highlight the influence of magma reservoir dynamics on the local stress field, influencing dike propagation and fault activation. The detailed spatiotemporal analysis sheds light on the eruption triggering mechanisms and the complex interaction between different parts of the magmatic system. The model presented contributes to a more comprehensive understanding of volcanic behavior, aiding in improved early warning and monitoring systems for oceanic islands.

This study presents a detailed analysis of the 2021 La Palma eruption, revealing a complex three-reservoir magma plumbing system. The integrated analysis of seismicity, moment tensors, and geodetic data provides a dynamic model of magma storage, transfer, and eruption triggering. The findings contribute significantly to our understanding of volcanic processes and improve eruption monitoring strategies for oceanic islands. Future research could focus on refining the reservoir geometry and connectivity through advanced tomographic imaging, investigating the role of volatiles in triggering the eruption and influencing reservoir dynamics, and further developing coupled models of magma flow and stress evolution in volcanic systems.

The study's interpretation is based on the available seismic data and may be affected by uncertainties in velocity models and earthquake location accuracy. The stress inversions provide only approximate estimates of stress orientations due to the highly heterogeneous stress field near the reservoirs. The resolution of the imaging is limited by the spatial distribution of the seismic stations, potentially affecting the precise mapping of reservoir boundaries. Furthermore, the model relies on a simplified representation of the reservoir geometry and magma rheology, which might not fully capture the complexities of the natural system.