Magmatic plumbing and dynamic evolution of the 2021 La Palma eruption

The 2021 La Palma eruption began on 19 September and lasted over 85 days, forming a new edifice on the western flank of Cumbre Vieja. It was the longest and most voluminous historical eruption on La Palma, with an extruded volume exceeding 0.2 km3 and the evacuation of approximately 7000 residents. Prior to 2021, knowledge of the Cumbre Vieja feeding system was limited to petrological and gravity studies. The Canary Islands are an oceanic volcanic archipelago; La Palma is the second youngest island, with a history of alternating eruptive episodes and flank collapses, and a north-to-south migration of volcanism culminating in the Cumbre Vieja rift. Petrological barometry suggests pre-eruptive magma storage in upper mantle reservoirs at roughly 15–26 km depth (possibly to 50 km), and an accumulation or underplating zone at 7–15 km depth where more evolved magmas may form. Seismic unrest in 2017–2021, including multiple swarms, indicated progressive magma accumulation at depth. The 2021 eruption is the first locally monitored on La Palma with an expanded IGN network, providing an unprecedented dataset to map magma pockets and pathways and to resolve the dynamics of the magmatic plumbing system.

Background work on La Palma and Cumbre Vieja includes petrological and geophysical evidence for magma storage at upper mantle depths (approximately 15–26 km and possibly down to 50 km) and transient magma accumulation in the 7–15 km depth range during eruptive episodes. Historical studies document the migration of eruptive activity southward on La Palma and the structural control of the Cumbre Vieja rift. Prior to 2021, long-term seismicity precursors were noted before the 1949 and 2021 eruptions, suggesting progressive magma accumulation months to years in advance. Between 2017 and June 2021, seven weak swarms (ML ≤ 1.9) occurred at 20–35 km depth with durations of 2–10 days; changes in gas emissions (hydrogen, helium isotopic ratios, thoron) before or after these swarms evidenced magmatic unrest. Petrological and deformation studies had hypothesised shallow storage near ~12 km depth and deeper reservoirs, but a comprehensive seismological moment tensor catalogue and high-resolution relocations for La Palma were lacking before this work.

The study integrates multiple seismological and geodetic methods. - Hypocentral relocation: A relative relocation of the IGN catalogue (2017–2021) was performed using the double-difference method (HypoDD). Waveform cross-correlation was applied on P and S phases with 2 s and 3 s windows, respectively, considering correlations ≥ 0.75 and event pairs within 1 km; phase-pick time differences considered event pairs with >8 common phases and <3 km separation. The relocation used the IGN velocity model and included approximately 16 million time differences (~15 million from waveform correlation and ~1 million from picks), averaging 20 phase data and 2000 time delays per event. Station sets varied with time and depth, including temporary stations and sites on neighbouring islands for deeper, larger events. - Moment tensor inversion: Full moment tensor inversion for events with ML ≥ 2.5 was performed using Grond, fitting P and S bodywaves in the time domain. Data from 13 stations on La Palma and 6 on neighbouring islands (up to ~130 km) were used, mainly broadband sensors. Instrument responses were deconvolved and empirically validated. Manual data curation removed noisy or saturated traces; P and S phases were picked manually to align with synthetics. Synthetics were computed for a regional velocity model using QSEIS. The inversion fit 1 s windows on P vertical and S transverse components with a 1–5 Hz bandpass. Optimization included 100,000 iterations (10,000 random global, 90,000 focused). Bootstrap ensembles (100 chains) provided uncertainties. Robust solutions were obtained for 156 earthquakes (73 shallow, 83 deep). Centroid depths and times and six independent MT components were estimated; epicentral coordinates were fixed to the catalogue values. - Stress inversion and analytical modelling: Least-squares optimization estimated principal stress orientations and relative magnitudes in spatially defined subvolumes, assuming homogeneous stress in each. Due to heterogeneous stress near reservoirs, results are approximate. A simplified analytical model of a tectonically loaded block with a pressurized circular reservoir was used to interpret stress rotations, dike paths, and seismicity clustering. - Seismicity and MT clustering: Density-based clustering was applied separately to shallow and deep volumes. Spatial clustering thresholds: shallow Nmin=150, ε=0.0002 (200 m radius); deep Nmin=50, ε=0.0004. MT clustering used Kagan-angle similarity: shallow Nmin=2, ε=0.2 (k ≤ 24°), deep Nmin=2, ε=0.15. Waveform-based clustering used 8 s windows, 6–16 Hz bandpass, and station-specific networks; families were defined at average correlation ≥ 0.7 across stations. - Eruption tremor: RSAM was computed at station CENR using 10-min windows, normalized by median and standard deviation, to track tremor amplitude. - GNSS deformation: Daily GNSS time series were processed with Bernese v5.2 on a regional network (>30 stations), ITRF2014, FES2004 ocean loading, IGS antenna models and orbits.

- Seismicity architecture and volumes: Relocation resolved two dominant co-eruptive clusters separated by an aseismic gap: a shallow cluster at 9.9–12.7 km depth (90% interval), roughly circular with ~3 km diameter and two main seismogenic volumes (6145 events), and a deep cluster at 32.8–38.1 km depth extending ~10 km NE (1054 events) with several small active subregions. A total of 8,488 earthquakes were relocated for the 2021 unrest. - Pre- and co-eruptive evolution: Seven swarms (2017–June 2021) at ~18–32 km depth indicated progressive pressurization and intrusion within an intermediate mushy reservoir. Immediate pre-eruptive seismicity began on 11 September 2021, with upward and lateral migration along a curved dike path from ~11 km depth to the first vent on 19 September. - Moment tensors and stress: 156 MTs (73 shallow, 83 deep) cluster into four families in each depth range. Shallow mechanisms are mostly oblique to strike-slip with positive CLVD and overall negative isotropic components; two largest families have similar mechanisms with opposite polarities and depths near 11.2–11.9 km. Deep mechanisms are mostly strike-slip to oblique with some thrust/normal; a polarity flip between top (~35.2 ± 0.3 km) and deeper (~36.5 ± 0.8 km) families is observed. Stress inversions indicate rotated stress fields with depth and time: near-NS compressive stress at shallow levels consistent with dike paths, and NW–SE stresses at depth slightly rotated from regional −N25°W compression. - Dynamic stages: Four stages during co-eruptive period: (1) Eruption onset with decreased seismicity and high tremor; (2) From 27 September, tremor drop and appearance of shallow cluster at ~11.8 km, with subsidence and lateral GNSS transients towards the shallow cluster implying drainage; (3) From 1 October, activation of deep cluster with largest events (Mw up to 4.1 on 3 November) and alternating seismicity bursts in shallow and deep clusters; (4) From 13 December (eruption end) to 21 December, tremor ceased, deep seismicity ended first, shallow seismicity waned later. - Aseismic transfer and reservoir system: The absence of co-eruptive earthquakes at ~18–32 km suggests an open, high-temperature aseismic conduit linking reservoirs. The system comprises three storage levels: shallow (~9–13 km), intermediate mushy (~18–32 km), and deep (~33–38 km). A shear-wave velocity anomaly near ~13 km supports a shallow reservoir. - Stress heterogeneity and focal mechanism flips: Reversed focal mechanism polarity between upper and lower parts of both shallow and deep clusters indicates strong local stress perturbations due to reservoir pressurization/drainage. Analytical models reproduce curved pre-eruptive dike paths, higher seismicity asymmetry, and mechanism reversals during depletion. - Doublets and pulsatory flow: Earthquake doublets (pairs of M ≥3.5 within seconds to a minute) occurred in the deep cluster, consistent with pulse-like magma ascent. Seismicity bursts were shortly followed by increased lava emission and new vents. - Deformation and isotropic MT: GNSS shows subsidence and horizontal motion towards the shallow cluster during drainage. Negative isotropic MT components are consistent with shallow reservoir depletion and high volatile contents causing transient pore volume reduction during co-seismic pressure increases. - Quantitative highlights: 8,488 relocations; shallow cluster 6145 events at ~10–13 km; deep cluster 1054 events at ~33–38 km; 156 MT solutions; peak magnitude Mw 4.1 (3 November 2021); eruption duration ~85 days.

The findings establish a detailed, dynamic picture of the magmatic plumbing beneath Cumbre Vieja. Pre-eruptive swarms and gas anomalies indicate progressive accumulation and pressurization within an intermediate mushy reservoir. The eruption was initiated by dike propagation from a shallow reservoir (~12 km), guided by stress heterogeneities that produced a curved path before aligning with regional NS compression away from reservoir influence. During the eruption, most seismicity localized above and below the intermediate swarms’ depths, delineating shallow and deep reservoirs separated by an aseismic conduit that enabled sustained magma transfer without wall-rock fracturing. The temporal evolution shows a top-down control early in the eruption: drainage and pressure loss at the shallow reservoir triggered delayed activation and drainage of the deep reservoir. From early November, bottom-up control emerged as deep seismic bursts preceded enhanced lava output at the surface. The cessation of deep seismicity before shallow activity after 13 December further supports vertically coupled reservoir dynamics. Reversed focal mechanism polarities within both clusters and stress inversion results highlight strong local stress rotations due to reservoir pressurization and depletion. These stress perturbations explain the observed dike trajectories, seismogenic volumes, and focal mechanism distributions. The integrated seismological, geodetic, and petrological evidence supports a three-level reservoir system, rather than a single large reservoir, with dynamic interplay governing magma supply, eruption tremor, deformation, and effusion rates. This refined plumbing model enhances interpretability of monitoring signals and supports improved early warning on oceanic island volcanoes.

This study resolves the fine-scale structure and temporal dynamics of the 2021 La Palma magmatic system using high-resolution relocations, a new moment tensor catalogue, stress analysis, and geodetic observations. Key contributions include identification of two active co-eruptive seismogenic clusters (shallow 9–13 km and deep 33–38 km) separated by an aseismic intermediate mushy reservoir (~18–32 km); documentation of stress heterogeneity and focal mechanism flips linked to reservoir pressurization and drainage; and a four-stage eruption chronology showing a shift from top-down to bottom-up control. The results argue for a three-reservoir system with an open, high-temperature conduit during eruption, explaining aseismic magma transfer and the timing of seismic and effusive signals. The conceptual model provides a framework to better interpret monitoring data and to improve early warning and response for future eruptions at La Palma and similar ocean island volcanoes.

- Stress inversion assumptions and heterogeneity: Estimated principal stress tensors are approximations because stress near reservoirs and conduits is highly perturbed and heterogeneous; true stress fields are likely more complex and time-variable. - Velocity model and depth estimates: Differences from previous depth estimates may arise from the chosen velocity model; event depths and structure are model-dependent. - MT inversion constraints: The inversion used only bodywaves, with fixed epicentral coordinates and limited sensitivity to lateral centroid shifts; lack of clear surface waves restricted data types. Manual alignment and data curation may introduce subjectivity. - Network evolution and station geometry: Temporal changes in the seismic network and differing station sets for shallow vs deep events can impact resolution and location precision, especially for early swarms and deeper events.