In this study, we merged datasets derived by the two governmental agencies operating on La Palma prior to and during the eruption: the Instituto Geográfico Nacional (IGN; 11 stations) and the Instituto Volcanológico de Canarias (INVOLCAN; 12 stations). Both datasets cover the time span from 7 October 2017 to 13 December 2021. When merging, we assumed that both agencies had recorded the same event if the difference between the origin times in the provided solutions was < 1 s. The merged dataset contains 13,681 events recorded on 23 stations, with 140,078 P wave and 155,231 S wave picks. Event magnitudes were based on solutions provided by INVOLCAN using the Ml magnitude scale.
Tomographic inversion was based on the local earthquake tomography code LOTOS18, which has been used to investigate dozens of different volcanoes9,19–22. First, event locations in the reference 1D velocity model are determined using a grid-search method and linear approximation of ray paths9. We used the topography surface to limit the depths of events; therefore, events can be located above the sea surface. During this step, we selected events with eight or more picks and removed all data with absolute residual values of > 0.5 s for P waves and 0.7 for S waves. After removing outliers, the dataset used for the tomographic procedure included 11,349 events with 121,572 P wave and 127,766 S wave arrival times (mean of 22 picks per event).
Next, we relocated event sources using the gradient method and 3D bending algorithm for ray tracing9. In the first iteration, the relocation was conducted in the starting 1D model; in subsequent iterations, calculations were performed in the updated 3D model.
The 3D distributions of P and S wave velocity anomalies were parameterised using grid nodes irregularly distributed in the study area according to the ray density. Between the nodes, velocities were approximated continuously using bi-linear interpolation. The minimum grid spacing in both the horizontal and vertical directions was 0.7 km, which is considerably smaller than the size of resolved anomalies in our model. To reduce the influence of parameterisation on the results, we performed inversions in four grids with different basic orientations to the azimuthal direction (0, 22, 45, and 67 degrees) and then averaged them into one regular grid.
The inversion was performed simultaneously for the P and S velocity anomalies, source parameters (dx, dy, dz, and dt), and station corrections. To stabilise the solutions, we used two types of regularisation—amplitude damping and flattening—which were performed by adding the corresponding equations to the general system. The values of the damping coefficients were determined from a series of synthetic tests with realistic anomaly sizes and noise levels. The derived sparse matrix was inverted using the LSQR method23,24.
After inversion in the four grids, the resulting P and S wave velocity anomalies were recalculated in a regular grid and then used for the next iteration, which included source relocation, matrix calculation, and inversion. In total, for the experimental and synthetic data, we performed five iterations, which was a compromise between solution stability and calculation velocity. The calculations enabled considerable variance reduction. In the L1 norm, the average P wave residuals reduced from 0.2397 to 0.0860 s (64.12%) and those for S waves reduced from 0.3058 to 0.1401 s (54.17%).
Extended Data Figure 1 and 2 show the inversion results of experimental data, including horizontal and vertical sections of Vp/Vs, P wave velocity, and S wave velocity anomalies. In the context of magma-related structures, it is important to present the Vp/Vs ratio, which was calculated by division of the derived P and S absolute velocities. The adequacy of this method was determined by the similar volumes of P and S wave data, and was confirmed using synthetic tests.
The optimal 1D reference velocity model was derived after several runs of the complete tomographic procedure. After each trial, we calculated the average P and S velocities at certain depths and used them for the next iteration. As a result, we obtained a fair balance between high- and low-velocity anomalies at all depth layers.
To access the spatial resolution of the resulting models and to derive optimal values of controlling parameters for the tomographic procedure, we performed a series of synthetic tests. In all cases, the synthetic model was defined by a set of positive and negative anomalies with respect to the a priori 1D reference model. Synthetic travel times were calculated for the same source-receiver pairs as derived in the main experimental model after five iterations. The derived data were perturbed by random noise with an average deviation of 0.03 s for both P and S data, which enabled the same variance reduction as in the experimental data inversion. Before starting the synthetic model recovery, we “forgot” any information about the sources. Then, calculations were performed based on the same workflow as in the experimental data processing, including source location in the starting 1D model using the grid-search method. During synthetic modelling, we tuned the values of the controlling parameters to derive an optimal quality of the initial model recovery; then, the same controlling parameters were used for the experimental data inversion.
We separately investigated the resolution in the horizontal and vertical directions. In the first series of tests, we defined several checkerboard models with different anomaly sizes in map view. In Extended Data Figure 1, we present three tests with anomalies of 2, 3, and 4 km separated by 1 km spacing with zero anomaly values. In all cases, the amplitudes of anomalies were ±8%. We defined the opposite signs of the dVp and dVs anomalies to enable contrasted variations of the Vp/Vs ratio. The recovery results are presented in two depth sections. Anomalies of 2 km in size were only resolved in the central part of the study area, where most earthquakes were located. For the models with larger anomalies, fair resolution was observed in most parts of the study area. Based on the results of these tests, we defined a contour of the resolved area, in which the results of experimental data were plotted (Extended Data Figure 1-3). From these tests, we confirmed that the distribution of Vp/Vs was correctly recovered, demonstrating the adequacy of the method for this parameter calculation.
In another series of tests, we explored the vertical resolution. Regarding the trade-off between velocity and source coordinates in passive source tomography, as well as dominantly vertical orientations of ray paths, we could expect poorer vertical resolution compared with horizontal resolution. In particular, there was a concern about the capacity of the existing data to provide abrupt changes in velocity at ~10 km depth, as observed in the experimental data inversion. To address this problem, we defined models with alternated anomalies defined in each of the vertical sections in which the main results are presented. Along these sections, the anomalies had a size of 4 km and a spacing of 2 km. In the vertical direction, they formed two rows with an interval of zero values at depths between 6 and 10 km. This test confirmed that major anomalies in the central part of the study area were correctly recovered; however, some diagonal smearing is observed in marginal regions. This effect was considered during interpretation.
As true event coordinates were presumed to be unknown in the recovery procedure, the synthetic tests allowed us to assess the accuracy of source locations. Extended Data Figure 3. Resolution test results based on the ability to solve realistic anomalies using a previous synthetic model (left column) and observing the recovery results (right column). Extended Data Figure 4 shows the mislocations of events in the model with a vertical checkerboard in section 2 with respect to the true locations. When using the starting 1D velocity model, the average error of source coordinate determination in the L1 norm was 0.69 km. In the final model, the average error reduced to 0.38 km. As expected, the maximum errors were observed for events on the periphery of the network and at the greatest depths. These errors do not significantly affect the interpretation of the results presented in this manuscript.