Satellite Radar and Camera Time Series Reveal Transition from Aligned to Distributed Crater Arrangement during the 2021 Eruption of Cumbre Vieja, La Palma (Spain)

Abstract

Magma-filled dikes may feed erupting fissures that lead to alignments of craters developing at the surface, yet the details of activity and migrating eruptions at the crater row are difficult to monitor and are hardly understood. The 2021 Tajogaite eruption at the Cumbre Vieja, La Palma (Spain), lasted 85 days and developed a pronounced alignment of craters that may be related to changes within the volcano edifice. Here, we use COSMO-SkyMed satellite radar data and ground-based time-lapse photographs, offering a high-resolution dataset to explore the locations and characteristics of evolving craters. Our results show that the craters evolve both gradually and suddenly and can be divided into three main phases. Phase 1, lasting the first 6 weeks of the eruption, was characterized by a NW–SE linear evolution of up to seven craters emerging on the growing cone. Following two partial collapses of the cone to the northwest and a seismicity increase at depth, Phase 2 started and caused a propagation of the main activity toward the southeastern side, together with the presence of up to 11 craters along this main NW–SE trend. Associated with strong deep and shallow earthquakes, Phase 3 was initiated and continued for the final 2 weeks of the eruption, expressed by the development of up to 18 craters, which became dominant and clustered in the southeastern sector in early December 2021. In Phase 3, a second and oblique alignment and surface fracture was identified. Our findings that crater and eruption changes coincide together with an increase in seismic activity at depth point to a deep driver leading to crater and morphology changes at the surface. These also suggest that crater distributions might allow for improved monitoring of changes occurring at depth, and vice versa, such that strong seismicity changes at depth may herald the migration and new formation of craters, which have major implications for the assessment of tephra and lava flow hazards on volcanoes.

1. Introduction

1.1. From Dikes to Crater Alignments

Many ocean islands are affected by pronounced volcano–tectonic interactions such as dike intrusions, faulting and rifting, gravitational spreading, flank instability, and sector collapses [1,2,3,4]. Rifting and flank motion may episodically accelerate [5,6] and create the necessary space for new intrusions along dikes and rift zones [7,8]. Normal faults, central grabens, magmatic fissures, and aligned cones are common morphological features of rift volcanoes [9]. In addition, the direction of rift zones can be subject to changes and may be adjusted if the surrounding stress field conditions change [10,11].

Dike intrusions have a tabular shape, where their width usually ranges from a few tens of centimeters to several meters. In contrast, the propagation length ranges from hundreds of meters to kilometers at ocean islands [12]. The direction of dike intrusion is determined by pre-existing pathways (weak zones, faults, or lithological barriers) and/or by the principal compressive stress axes, which cause dikes to intrude perpendicular to the minimum compressive stress (σ3) and toward the maximum compressive stress (σ1) [12]. Likewise, a volcanic edifice load or changes associated with pressurized reservoirs at depth can deflect the ascent direction of dike propagation laterally [12]. Dike propagations cause faulting and earthquakes [13]. However, earthquakes produced by dike propagation reflect the differential stress and not the extension of the dike [12].

Dikes are the feeders of fissure eruptions and aligned eruptive vents. At Kilauea, a prominent crater alignment followed the seismically inferred dike intrusion in 2011 [14]. A similar example can be found in the eruption of Holuhraun in Iceland in 2014. This eruption was characterized by the movement of magma from the Bardarbunga subglacial caldera to the Holuhraun eruption site, approximately 45 kilometers away. The eruption occurred at fissures and focused at vents aligned above a feeder dike [15]. The feeder dike was inferred to be oblique to the rifting direction, leading to complex faults, step-overs, and en echelon segments of fissures and eruption sites [14,16]. Likewise, dikes at the Mount Etna NE rift zone are expressed at the surface by elongated, breached, and aligned scoria cones along the fault. This fault developed in relation to an unstable eastern volcano flank, which slides eastward toward the sea [17]. On the Canary Islands, dike intrusions and rift zones are common, associated with the governing stress field, and provide favorable conditions for magma plumbing systems [7,18]. Rift zones are characterized by crater alignments and parallel groups of dikes beneath, which are particularly common on Tenerife Island [10], El Hierro Island [11], and La Palma Island [19]. In addition, the 2021 eruption on La Palma was associated with crater alignments into a crater row, which may imply the presence of an underlying feeder dike, as our study suggests.

1.2. The Challenge of Monitoring the Development of Crater Rows

The accessibility and visibility conditions at an active volcano are often challenging and obscured due to the emission of ash, steam, and volcanic gasses, which, combined with rapidly changing weather conditions, make direct and close observations hazardous. This is also why successful monitoring of active volcanoes with passive optical remote-sensing satellites such as Landsat, Pleiades, SPOT, or Sentinel-2 rarely provide a temporal or spatial resolution sufficient for change detection. On the other hand, active sensors such as synthetic aperture radar (SAR) transmit their own radio signal and measure the amount of energy backscattered by the object, making them independent of sunlight and atmospheric conditions. Furthermore, SAR collects information in the longer wavelengths, specifically in the microwave part of the electromagnetic spectrum, thus is allowing visibility through clouds, volcanic gasses, ash, and steam [20,21]. In this study, we have used X-band radar amplitude signal images acquired from the COSMO-SkyMed (Constellation of small Satellites for Mediterranean basin Observation) SAR satellites to track surface changes and analyze the development of new eruption vents. The analysis of X-band amplitude signal has already provided results at other volcanoes, e.g., at Merapi in Indonesia [21], Soufrière Hills on Montserrat Island [22,23], or Cotopaxi in Ecuador [24].

Volcanic structures and edifices can display rapid changes during volcanic activity, and interpretation of many of these variations requires visual confirmation. Therefore, video and time-lapse camera records helped to monitor volcanic surface changes [24]. Time-lapse cameras have been used to elaborate on the transition from fissure eruptions to formations of crater rows and associated topographic changes [14,15]. It was shown that dike intrusions feed a long and continuous fissure eruption in the first hours to days, which then evolve into several aligned eruption vents with the formation of craters [14]. It is likely that the distinct craters at the surface are connected to multiple deeper dikes along finger-like conduits also during a later stage of the eruption. However, the geometric details of this connection are debated [25]. While the erupting craters grow in depth and dimension, their location often remains stable throughout the eruption [14]. This differs from observations made at La Palma, where the feeder dike opened multiple new eruptive vents during the 85 days of the 2021 eruption. At first, vents were aligned in a clear profile facilitating the monitoring of the lava fountains; however, they disaggregated during the later stages of the eruption. The details of such a deviation of crater rows from an idealized chain alignment have not been studied in previous works. In this study, we combine SAR and optical camera time-lapse data to trace craters, to locate crater positions and lava fountain activity, further comparing the results to a seismic catalogue.

1.3. Study Area

The 2021 Tajogaite eruption occurred at the upper western flank of the Cumbre Vieja (Figure 1), located on La Palma Island, which is the most active volcano of the seven Canary Islands. Repeated phases of island construction during eruptions and destruction by lateral collapses characterize the evolution of La Palma [26]. Two polygenetic volcanic complexes shape the Island: on the northern part, a large circular volcano and uplifted ocean crust complex is dominant, whereas on the southern part, the Cumbre Vieja is an elevated volcanic rift most active in written history [26]. The northern near circular volcano evolved during the Pliocene by superpositioning of four distinguished volcanic structures, from which one is a submarine and the other three are subaerial volcanoes. The three subaerial volcanoes are named Garafía (formed 1.77 to 1.20 million years ago), Taburiente-Cumbre Nueva (formed 0.89 Ma) and Bejenado (formed 0.4 Ma) [27]. The northern edifice has been affected by several mass-wasting events, possibly catastrophic, which have encarved the deeply incised morphology of the Taburiente-Cumbre Nueva (560 ka) and later affected the Cumbre Nueva rift zone. The southern volcano edifice is the youngest (123 ka) and, thus, the area where the last six historical eruptions were recorded on the island in the last 500 years [27]. Cumbre Vieja is a rift oriented north to south (N–S) with an altitude of ~1950 m above sea level. Some of the most recent events were the San Juan eruption in 1949, which created three eruptive vents (Duraznero, Llano del Banco and Hoyo Negro), and the Teneguía Eruption in 1971 [28].

Dikes and sills can be found all over the island and usually range from 0.2 to 1.2 m wide [26,29]. Shallow intrusions at Cumbre Vieja show topographical irregularities and produce prominent structural features, forming the volcanic ridge associated with faults, fissures, and vents, many of which are aligned along the N–S directed rift zone [19,27,30]. Furthermore, some less-expressed structural features, or weak rift zones, in the northern sector of the Cumbre Vieja diverge NE–SW and NW–SE [19,31]. The 2021 erupting craters follow a NW–SE direction, as will be detailed further below.

Following 50 years of quiescence, a new unrest period with intense seismic activity began on 11 September 2021 and continued through the 85 days of the eruption [32,33]. During the 85-day eruption period, several vents opened and nested on the elongated cone on the western flank of Cumbre Vieja, near to the faulting system of 1949 [32,33,34]. It was estimated that the Tajogaite cone was built up to 187 m from its base due to the accumulation of tephra (ash, lapilli and larger bombs) and lava flows [35]. The eruption also produced constant emission of ash, steam, volcanic gasses, and approximately 200 million m³ of lava [32]. However, more detailed estimates come from drone surveys that suggest a smaller volume of ~180 Mm³ [35]. The eruption had a Volcanic Explosivity Index of VEI3 [32]. The 2021 Cumbre Vieja eruption presents a unique opportunity to study crater arrangement as an indicator of dike intrusion and its possible link to volcano–tectonic processes [33].

2. Data and Methods

2.1. SAR Data and Methods

We analyzed a catalog of SAR amplitude images acquired by the Italian Space Agency (ASI) COSMO-SkyMed satellite constellation (CSK). The two values that characterize a SAR image are the azimuth and the range [36]. The SAR amplitude signal represents the intensity of the energy backscattered by all ground scatterers within a resolution cell (pixel) [37]. The signal backscatter is determined by the local incidence angle of the electromagnetic radiation, the surface roughness, the surface chemical characteristics, as well as the water content of the soil [20]. This difference in intensities makes recognizing the surface features and volcano eruption deposits possible [20]. As a result, rough surfaces facing the sensor (i.e., whose slope has a value close to the incidence angle) have a strong backscatter, which results in brighter pixels in the SAR image; conversely, surfaces with an opposite slope (backslope) are characterized by lower scattering values, resulting in darker pixels [20]. Data are displayed in greyscale maps, and white represents the highest amplitude and black the lowest. The images used in this study were acquired in both ascending and descending strip map (HIMAGE) mode with a spatial resolution of about 3 m and in HH polarization.

Geometric distortions such as layover [38], foreshortening, and shadowing effects [22] (Figure 2) characterize the SAR images. With the help of an available digital elevation model (DEM), the geometric distortions from the layover and foreshortening areas can be corrected by means of geolocation procedures [22,37]; the images can be converted into a geographic (latitude, longitude) or cartographic (easting, northing) map projection for integration into GIS systems, thus allowing for easy comparison with other data [38]. In the case of rapidly and significantly changing topography, as in the La Palma example, the DEM would need to be updated before each acquisition for the terrain correction procedures to be effective. As theses DEMs were not available, some residual geometric distortions remain visible in our data. This however does not hamper their interpretation.

In total, we studied 48 CSK synthetic aperture radar (SAR) amplitude images. The descending images (Table S1) were acquired between 25 September and 14 December 2021 (track H4-08) with a look angle of ~38°. In the ascending orbit (track H4-06, Table S2), we have 11 images acquired between 25 September and 6 December 2021. These images have a look angle of ~35.5°. Lastly, the remaining 20 images of the ascending orbit (track H4-12, Table S3) were acquired between 30 December and 9 December 2021, with a look angle of ~44°. The resulting database allows us to monitor the 85-day-long crater evolution at an unprecedented level of temporal detail.

The SAR amplitude images were co-registered to a corresponding reference image for each of these three tracks to create three series of aligned and stacked images [39]. Next, each SAR amplitude image stack was terrain corrected and projected onto a common geographic reference system [37]. Both registration and terrain correction were performed with a 5 m LiDAR DEM generated in 2009 before the eruption by the Centro Nacional de Informacion Geográfica of Spain. Lastly, the images were speckle-filtered and normalized to reduce scatters and radiometric noise. Toward the end of the eruption, a new DEM created from a UAS survey became available and was used for the terrain correction of the descending SAR data acquired on 6 December 2021 (See Section 2.3).

To trace the evolution and locations of the craters, we studied the SAR data stack in different ways. First, we generated composite maps (Figure S1) to identify general changes in the amplitude from one image to the next. Secondly, we extracted NW–SE amplitude profiles (see Figure 2c) and traced the crater centers as points with latitude and longitude through time. Finally, a post-eruption UAS DEM and orthophoto were used to validate the results of the SAR by comparing them with the geocoded SAR image on 6 December 2021. The crater center points are hence manually identified from a characteristic brightness pattern: The crater rims and slopes facing the satellite line of sight (LOS) appear brighter. The deep crater floors and shadowed crater walls are shown as darker (Figure 2a,b). Considering the amplitude profiles (Figure 2c), we see that the peaks and troughs of the profile allow for the identification of a crater, where its center point is inferred near the transition from dark to bright regions. Due to the geometric distortion, the crater walls facing the satellite are foreshortened, while the opposite walls and bottoms of the craters appear shadowed and wider (Figure 2b), which is why the brightness pattern appears non-circular and crescentic. Using this approach, we quantified and analyzed all the crater locations through time.

2.2. Camera Data and Methods

Time-lapse and video camera employment has shown to be of great use to monitor volcanic processes, as it provides visual support that helps to interpret data associated with plumes [40], erosion processes at Volcán de Colima, Mexico [41], gas pistoning at Kilauea, Hawaii [42], lava flow surface velocities [43], and surface deformation at Mount St. Helens, USA [44]. Cameras can provide data on a volcano edifice from different viewing angles, allowing one to study specific details such as the migration of eruption vents. In total, we installed seven digital time-lapse cameras (Reconyx Ultrafire XR6) located around the eruption site to record images with a resolution of 3776 × 2124 pixels on a regular interval. Here, we used three of the cameras that have a good view onto the craters (Figure 1), for which we set the time-lapse interval to 1 image every 15 min, day and night. The three cameras had two lenses each to capture images in RGB during the day, and at night, they shifted to residual light view mode. Camera 1 was installed at ~1500 m distance from the erupting vents at 28.607°N, 17.849°W at 1605 m above sea level and recorded from 2 October 2021, until the end of the eruption. The field of view (FOV) covered the SE part of the active craters, facing toward the NW direction. Camera 2 was installed at ~500 m from the erupting vents at 28.618°N, 17.860°W at 1087 m height, and operated from 11 November 2021, onward. The FOV covered the NE craters and was facing SW. Camera 3 was installed on 28.612°N, 17.863°W at 1055 m height on 4 December 2021. It looked toward the SE flank facing NW and was located at a distance of ~140 m from the closest vent.

For the time-lapse camera data processing, we used the image metadata (exif file), providing time and date (in UTC) and camera characteristics. We resampled the data by time to select only one picture every 15 min for each camera. The images were scaled from pixel to meter. For this, we measured the distance in meters using the DEM between two prominent marker points visible in each view of the cameras. Then, the images were processed in the Fiji ImageJ open source software package with the multi-kymograph tool. Kymographs are space vs. time plots that serve to analyze dynamic volcanic processes [14,25,40,45]. This allowed us to extract the pixel values (RGB or intensity) along a predefined profile and assemble the profiles from all images together (Figure 3). The graph’s x-axis corresponds to the number of pixels along the profile, and the y-axis corresponds to the time of each pixel row. To overcome differences between day (RGB) and nighttime (grayscale) images, we concentrated on the nighttime images only, acquired from 20h00 to 08h00 (UTC). To improve the contrast for better visual analysis we plotted the kymographs in a jet color theme, showing the bright nighttime lava fountains by warmer colors (red, orange, and yellow).

For the kymographs analysis, we manually picked a trend line following the highest intensity values representing the lava fountains’ activity. With this, it was possible to identify periods of stability and instability in the vent activity. Similarly, it was possible to identify intermittent activity between nearby lava vents in the NW cone and more drastic changes such as shifts in activity from the NW to SE vents. This helped identify the vent activity’s general trend during the eruption. Changes in the position of lava fountains were attributed to migrating eruptive activity, the opening of new vents, or the synchronous activity at vents. This corroborated the SAR results.

2.3. Comparative Data (SFM, Seismic Data)

Unmanned aircraft systems (UAS) are used for acquisition of near-field photogrammetric data that allow for affordable DEM generation [15,16,41,46,47]. In this study, we used a DJI Phantom 4 RTK drone (Figure S6) to collect image data on 15 January 2022, covering the entire post-eruptive main crater area. We used the Structure-from-Motion (SfM) approach implemented in the Agisoft Metashape software (vs. 1.8.0) to generate a high-quality DEM and orthophoto with a resolution of 0.10 and 0.05 m, respectively. Due to the RTK module, the geolocation accuracy was significantly improved, here achieving ~1.3 cm accuracy, although without additional ground control points, a slightly higher vertical error was expected [48]. Our data are very similar to a previous UAS survey in [35], but provide us with an even higher resolution.

We used the UAS data to validate our SAR amplitude image analysis and to explore satellite artifacts, such as geometric distortions of the satellite radar data, and camera visibility limitations that can lead to uncertainties in identifying surface structures, especially those associated with steep slopes and changing topography.

The Cumbre Vieja eruption was preceded by an increase in seismic activity recorded two years before, which motivated maintaining a dense seismic network that was still operating when the eruption started in 2021. Two sets of new seismic swarms were interpreted as indicators of sill and dike intrusions that eventually led to the opening of the first eruptive vents on 19 September 2021 [49]. According to De Luca et al. (2022) [49] and Del Fresno et al. (2022) [50], the syn-eruptive seismic activity clustered at two depths located at a shallow level of 10–14 km depth and at a deeper level of 33–39 km depth. The seismic catalog, which we use here for comparison, was produced and provided by Del Fresno and colleagues of the National Geographic Institute (IGN) of Spain, and it contains some stations operated by us [51]. The IGN catalogue is used to count the number of daily events at different depths from September to December 2021. To be consistent with our kymograph analysis, we filtered the data and selected only seismic events that occurred during the night from 20h00 to 08h00 UTC. For our analysis, the data were plotted according to the number of events per day and classified into the following depth ranges: greater than 20 km, between 9 and 20 km, and less than 9 km. We compared the seismic data, main crater changes in SAR data, and the lava fountain activity extracted from the kymographs as a track of intensity peaks.

3. Results

3.1. SAR Observations

The results obtained by analyzing the descending orbit SAR images show the details of the emerging craters and migrating crater chains. The sequence of the georeferenced descending SAR amplitude images is presented in Figure 4.

Prior to the 2021 eruption, the data allowed to locate an old (Holocene) eruption crater open to the NW. Day 1 of the eruption was 19 September 2021, where a localized eruption vent appeared at the northern flank of this crater. On Day 7 of the eruption (25 September 2021), the presence of new craters could still not be identified in our SAR amplitude images, although the presence of an extensive cone was evident by a moderate and smooth appearance implying predominantly fine-grained tephra deposits (Figure 4a). SAR amplitude images imply that the initial eruption was fed from a single point-like vent and was not associated with a buildup of large craters or fissures. On Day 11 of the eruption (29 September 2021), the SAR image showed five new craters (Figure 4b). Four were in the NW part of a newly evolving cone and one on the SE flank at the southern rim of the pre-existing Holocene crater. In this image, the craters of the main cone are very well-defined by the bright (high reflectivity) upper rim of the crater walls, while the deepest parts and one of the inner crater walls are shown in shadow (low reflectivity). On the other side, the southeastern crater has only one visible brighter rim, but its location can be confirmed by the intensity profile and in the optical images. Furthermore, it can be seen in the SAR amplitude images that the location of the new craters is almost perfectly linear, and the shape of the craters is round (Figure 4b). In the following SAR images, linearity is still observed, but the number of craters varies. First, fewer craters are visible, which may be attributed to intense coverage and obscuration by tephra fall. On Day 23 (11 October 2021), the SAR image shows that a new open crater evolved on the NW continuation of the crater row associated with a crater wall collapse. In addition, other craters formed both to the NW and SE of the main cone (Figure 4f). Likewise, the volcanic edifice appeared more circular (in map view) after 11 October (day 23 of the eruption). Until Day 27 (15 October 2021), the craters maintained a circular outline and were NW–SE aligned; thus, they evolved following a stable linear trend for approximately 20 days (Figure 4g).

On Day 39 of the eruption (27 October 2021), we see that new craters began to develop in the SE, where two new parallel craters emerge (Figure 4h). This generates two concentrations of craters: a linear one to the NW and a more disaggregated one in the SE sector of this structure. On Day 43 of the eruption (31 October 2021), we observe the appearance of a new crater between the concentration of craters in the NW and the SE (Figure 4i). In the subsequent two weeks (captured by SAR images on 3, 4, and 12 November 2021), we see a small variation in the number of visible craters while they maintained their linear distribution to the NW and were more disaggregated to the SE.

On Day 71 (28 November 2021), we observe an outbreak of craters on the northern flank and in the southern part of the volcanic edifice (Figure 4m). This date is associated with major structural changes, as we observed the appearance of other craters, some only half visible in a crescent form, which are located mostly on the NW and SE flanks of the cone. This change from a linear crater evolution to a more distributed cluster continues afterward, especially in the last phase of the eruption, which was declared over on Day 85 (13 December 2021).

Looking at the general distribution of the craters seen in the SAR amplitude images, we can visualize three phases: The first one (Phase 1), lasting from the beginning of the eruption until 15 October 2021 (day 27 of the eruption), is expressed by a clear alignment of all craters (Figure 4a–g). Furthermore, the regular spacing between craters is notable on 29 September 2021 (Figure 4a) and varies slightly on the following dates. The second phase (Phase 2), from 27 October 2021 (day 39 from the eruption), is characterized by the migration of craters to the SE, maintaining its linearity in the NW and disaggregated to the SE (Figure 4h–l). Again, we can see the equal spacing of craters, especially in the NW-aligned craters. Finally, in the third phase (Phase 3) from 28 November (day 71 of the eruption), we see an outbreak and migration of the craters to the NE and SE of the main cone (Figure 4m). In this phase, the new crater positions suggest the formation of a broader zone or even development of a second row of craters oblique to the main row, striking roughly N–S (Figure 4m–p).

While the descending data provide very good insight into the evolving crater row, the ascending radar data (Figure S2) add temporal details and allow for independent confirmation of the three-phase crater evolution. Due to the higher geometric distortion of the ascending data and the combination of two orbits, the linearity of the craters is not as clear as in the descending data. This is because the ascending data LOS direction is orthogonal with respect to the crater row, which forms across the evolving cone, resulting in a rather curved crater row arrangement. In addition, other geometric distortions prevail. For instance, in the ascending data (Figure S2), the SW flank looks brighter and tilted, especially in the more central craters. Furthermore, the SE craters show stronger distortion and are barely recognizable. The craters are close to each other but do not show an alignment; instead, they can be traced along a curve. Despite this, these data complement the results of the first dataset. We see that Phase 1 did not end on 27 October (day 38 of the eruption) but earlier, around 16 October (day 27 of the eruption) (Figure S2j). Here, we can also observe the beginning of the migration and disaggregation of the craters toward the southeastern part of the cone. Moreover, it also supports the results of the descending data about Phase 2 ending around November 25 (day 68 of the eruption) (Figure S2z). Finally, at the beginning of Phase 3 (Figure S2aa), the appearance of the second row of craters in the northeastern flank of the cone, oblique to the first, is detected more clearly. Here, we see how the NE craters that emerge on this flank on 28 November (Figure S2aa) align with the craters that emerge on 31 October and connect the linear to the disaggregated group (Figure S2p). We note that the number of craters identified in the ascending data is lower (14 craters) than the descending data (18 to 20 craters) by 6 December. However, it shows similar disaggregation characteristics concentrated in the cone’s NW and SE parts.

Likewise, we validated the SAR data analysis results by comparing and identifying the craters in an orthophoto and the DEM generated from the drone survey realized on 15 January 2022. Here, we observed the same structure of the volcanic edifice and found a similar number of craters, which maintain the same structure and distribution (Figure S6).

3.2. Camera Observations

Figure 5 shows the aligned kymographs from the three cameras. Due to the location of Camera 1, in the SE part of the cone, the row of craters is seen from the front, which limits the exact location of the vents. Even so, it is possible to identify the phase shift’s most drastic and decisive changes. In the kymograph of Camera 1 (Figure 5a), we observed that from 3 to 18 October 2021 (day 14 to 29 of the eruption), the activity tended to move slightly to the right in the NW direction, approximately 150 m, and remained relatively stable afterward. On the 18th (day 29 of the eruption), we observed the first drastic shift to the left (SE), approximately 440 m, which lasts at most one day. After that, the activity migrated again to the right (NW) approximately by 290 m, where it remained stable for about five days. On 24 October 2021 (day 36 of the eruption), we observed that the magmatic activity fluctuated between the NW and SE craters, which was maintained for about eight days. After this, the activity stabilized again in the central vents. The next change occurred at the end of November and the beginning of December, when we observed that the activity migrated to the SE and remained active there until the end of the eruption.

Because Camera 2 had a lateral view of the craters’ row, it facilitated the differentiation of the activity across the craters. The data from 11 November (day 54 of the eruption) showed small activity jumps between the NW central craters (Figure 5b). On 17 November (day 60 of the eruption), we see a drastic and brief jump of around 250 m toward the left (SE) to then stabilize in the central vents. This is also somewhat visible in the kymograph of Camera 1. Then again, we can see how from 28 November (day 71 of the eruption), the activity migrated toward the left (SE) vents. Because of the camera acquisition angle, the activity of the second crater row, identified in the SAR data, might have been overlooked. Finally, the kymograph of Camera 3 (Figure 5c), operating only for a shorter period, validates the results since the migration of activity to the SE craters is also visible.

3.3. Comparison of SAR, Photogrammetric and Seismic Observations

The camera kymographs support the results of the SAR data since they help to establish the phase changes and the appearance of new craters that coincide with the dates observed in the SAR data more accurately. In some examples (Figure 5d–h), we can explore how craters identified in the SAR data are also visible in the photographs and are established when they first emerged. This is the case at the beginning of the second phase with the emergence of craters in the SE flank (Figure 5d). Here, we see two craters, which we can identify in the photographs as well (Figure 5f,g). Similar to on 31 October (Figure 5e), we note the emergence of a new crater that connects the NW–SE linear crater row to the disaggregated ones in the SE. The emergence of this crater was visible on the photographs (Figure 5h). Finally, at the beginning of Phase 3, the development of the craters in the NE flank (Figure 5i) can be identified in the photographs (Figure 5j), as can be seen for the craters on 5 December 2021 (Figure 5k) and the corresponding photographs (Figure 5l).

By comparing the peaks of intensity extracted from kymographs of Camera 1 and Camera 2 with the number of seismic events, we noted that around 24 October 2021 (day 36 of the eruption), an increase in seismicity occurred. This is located at a shallow depth of 9 to 12 km, immediately followed by a shift of lava fountain activity toward the SE (Figure 6a,b). A similar seismic increase followed by a shift in eruptions from one to the other craters occurred again on 17 November (day 60 of the eruption). This date concurs with a rapid increase in the number of earthquakes from 34 on 16 November to 158 on 17 November. We note that this is one of the very few episodes where seismic events were shallow between 9 and 12 km depth only. In addition, in the lava fountain activity recorded by Camera 2, we see an immediate shift of activity toward a SE vent. This phenomenon occurred again on 30 November (day 73 of the eruption) and 1 December 2021 (day 74 of the eruption), with a rapid increase in the number of earthquakes from 67 on 29 November to 184 and 143 on 30 November and 1 December, respectively. On these dates, the seismic events were of medium (9 to 12 km) and very shallow (<9 km) depth. We noted a strong directivity, as the following eruption activity uniformly and rapidly shifted toward the SE vents, before gradually moving back again. However, this time, it took longer (three to four days) until the changes were visible as activity in lava fountains. In addition, it was evident that from 10 to 27 November 2021 (day 53 to 70 of the eruption), when the activity of the lava fountains looked stable in both cameras, the number of seismic events of greater depth (>20 km) was higher than the ones of medium depths (9–12 km).

4. Discussion

The 2021 Tajogaite eruption was the longest eruption ever recorded at La Palma and on the Canary Islands, lasting 85 days, and it caused major infrastructure damage in the pathways of the lava flows and during ash fall on the western and southwestern flanks of the Cumbre Vieja. A new cone rapidly emerged within the first days, continued growing and was incrementally dissected by craters that became nested and aligned, forming a central depression in the direction NW–SE. During later stages, the eruption further shifted to the SE and produced a number of small and clustered craters feeding new lava flows. Our study uses synthetic aperture radar (SAR) satellite data to describe the chronologic evolution and distribution of craters developing during the eruption, with important implications for lava flow hazards and possible drivers of sudden shifts of eruption sites discussed below.

4.1. Data and Method Limitations

Over the past 20 years, SAR data have been increasingly available for studying active volcanoes and their hazards. The SAR interferometry approach (InSAR and time series) was successfully applied to constrain surface deformations at volcanoes [52,53]. During the 2021 Tajogaite eruption, however, near-field InSAR were not successful due to the loss of coherence during the ash-rich eruption, which is why deformation data remained limited. Therefore we followed earlier studies that exploited the amplitude signal of SAR data, which helped to identify surface changes due to their ability to operate throughout the day and night, independent of visibility and weather conditions [54]. However, as many volcanic processes, such as eruption changes and formation of craters, occur rapidly during hours or days only [22,36], we could not possibly capture all changes occurring at La Palma. Although we used a large number of satellite images, on average one SAR image every two days, monitoring faster changes occurring on Tajogaite remained challenging.

Another important limitation of the SAR data we used is related to its geometric distortion. During the Tajogaite eruption, the topography was rapidly evolving, and the elevation of the cone was growing and partly collapsing again. The eruption cone was built by the accumulation of scoria material, increasing its relative height to 187 m [35]. We can identify a geometric distortion increasing during the course of the eruption and new cone buildup, which is largely pronounced in the ascending data, and to a smaller degree also seen in our descending data. Furthermore, the geometric distortions varied depending on the look angle, which was evident when analyzing our ascending data for two tracks: foreshortening and shadowing were found to be stronger in track 12 than in track 06. We find that this makes the crater flanks in track 12 look approximately half the size of those in track 06. A similar difference was found for the shadowing effect of track 12, where craters look more elongated, and consequently, the crater row appears wider in scale (see supplementary material Figure S2). Therefore, for correct interpretation of the location of craters, comparing SAR images acquired from different tracks (looks) helps to avoid false interpretations.

Due to these geometric distortions, we also combined and validated findings with independent optical imagery. By comparing our findings with high-resolution drone orthophotos and a post-eruptive digital elevation model, we found the number of craters to be the same as the ones found in the SAR data. Similarly, the general orientation of the aligned structure is comparable, although the strike locally changes. We also performed a correction on the later SAR data (e.g., the 6 December 2021 image) with the post-eruptive digital elevation model and found some distortion effect compared to the data terrain corrected with a pre-eruptive digital elevation model (Supplementary Material Figure S7). This implies that a regularly updated digital elevation model would be required to achieve improved results.

Wadge et al. (2006) [23] and Arnold et al. (2018) [20] pointed out that for X-band, a height variation of a few millimeters only can cause surfaces to look brighter in the SAR data. The overprinting of materials of different grain sizes might cause variations in the amplitude signal. This could explain the expression of visible craters, especially in Phase 1, where we see five craters on 29 September 2021 and only three on 1 October 2021. This may suggest that finer material covered the new cone and craters, resulting in lower contrasts of the SAR amplitude and, therefore, less defined or partially indistinguishable crater rims.

The camera and kymography data were found to be highly instructive, but were limited by the single angle perspective, which prevented us from tracking the activity of vents perpendicular to the field of view. Camera data hence face the limitations similar as optical satellite images, such as reduced visibility due to bad weather conditions and the presence of gasses and ash from the volcano [40,43,44]. These limitations also hindered the interpretation of the data in our kymographs. This is the case for the daytime images, which had more cloud cover and haze; thus, the activity of the lava fountains was not visible; only the time-lapse photography from 20h00 to 08h00 UTC was used. However, the combination of SAR images and camera data complement each other well, resulting in a whole dataset with good spatial coverage and nearly daily temporal resolution.

4.2. Spatial Distribution of Craters and Implications for a Feeding Dike

The SAR images from Phase 1 suggest that the new Tajogaite cone was already visible on day 7 of the eruption, which is allegedly linked to a dike at depth as inferred by geodetic data [49]. On day 11 of the eruption, the Tajogaite cone was found dissected by five pronounced craters. This allowed us to speculate that the eruption initially was fed by a rather localized finger or conduit connected to the deeper dike, which during the following days of Phase 1, propagated and developed a shallow feeder dike leading to the observed alignment of craters dissecting the cone.

During Phase 2 of the eruption, this trend was further established, fed by a NW–SE oriented dike. We interpreted the appearance of new craters in the SE at the end of October by further propagation/opening of the underlying dike intrusion. This is also supported by the information from the kymographs, as brief activity of the lava fountains was spotted in the SE flank. A reason for this migration could be attributed to the volcanic edifices load, which can alter the local stress field and thus affect the course of nearby intruding magma [12,22,55]. Likewise, migration of the eruption site can also indicate a change in the supply rate [56]. The orientation of the dike feeding the crater row of Phase 1 and 2 is in general agreement with a previously inferred rift zone [27], and with the orientation of the fissures that were active during the 1949 eruption [19].

In Phase 3, the appearance of new craters in the NE and in the SE coincide with the formation of a pronounced set of fractures striking approximately N–S [57], which is a new azimuthal direction that cuts the prevailing dike trend and hence suggests a change in volcano–tectonic dynamics starting from 28 November 2021. We were also able to identify a new arrangement of craters, which suggests the fractures to be related to the intrusion of a new dike that is oblique with respect to the trends identified in Phase 1 and Phase 2. A number of separated craters developed, together with two eruption vents, one in the north and one in the south of this dike. We conjectured that this last dike is driven by late changes in the stress field, just 2 weeks before the eruption came to an end.

Observations made at La Palma significantly differed from crater and vent developments analyzed elsewhere. During the 2014 Holuhraun eruption in Bardabunga (Iceland), a gradual development from fissure eruption to venting at isolated craters was observed [15]. Above the feeding dike, the number of active vents decreased from 57 to 10 in the first 5 days of the eruption, explained by the solidification of magma channels with less heat transport due to less activity [15]. In contrast, during the 2021 Tajogaite eruption on La Palma, the number of craters and active vents increased from 1–8 in Phase 1 to 6–11 in Phase 2, and eventually to 18–20 in Phase 3. This might suggest that at Cumbre Vieja, the feeding dike developed multiple paths to the surface, either due to reactivation of existing weakness zones, or due to a complex segmentation of the dike in the presence of a complex stress field, or by other local effects such as material and structural heterogeneities in the newly built cone over the preexisting topography.

4.3. Temporal Coupling between Crater/Vent Activity and Deep Seismicity

Over the years, seismic data have been largely used to predict eruptive cycles, eruption sizes, and other changes during eruptions [12,56,58], as well as the coupling between deep, shallow or distal processes. Similar coupling of the deep and shallow storage regions was also inferred during the 2021 La Palma eruption [59], propagating from 10 km depth to the Tajogaite eruption site in less than one week. We compared our remote sensing data and the emergence of venting activity and crater formation to seismic activity. We noted a period of approximately 14 days with a higher number of seismic events, during which the highest seismicity peak (24 to 26 October) concurred with the transition from Phase 1 to Phase 2. This happened again toward 28 November, which marks a pronounced seismicity peak, and the transition from Phase 2 to Phase 3.

In addition, other short and distinct peaks of seismic activity occurred almost simultaneously with newly emerging craters toward the SE flank as well as high lava fountain activity, which we see in the kymographs, as well as newly emerging craters toward the SE flank. For a number of these changes, we can infer that seismicity increased prior to changes observed at the surface, although we noted that the SAR satellite data barely have the temporal resolution to further substantiate such a causal relationship. By closer examining the seismicity, it may be possible to distinguish more than three phases, albeit these did not have a clear impact on the eruption processes on the surface, which we observe with SAR and camera data.

Of particular interest is the timing of significant changes other than seismicity. During Phase 1, we could identify a strong fluctuation of lava fountain activity occurring toward the NW of the crater alignments. Within Phase 1, a number of spectacular gravitational collapses occurred, involving parts of the Tajogaite cone and the NW crater walls on 25 September [60], as well as on 4, 9, 11 and 26 October. These collapses were excavating a deep open crater on the NW flank also visible in our SAR data from 11 October (Figure 4f). A closer look at both our SAR intensity peak profiles and the kymographs reveals that after every collapse, the lava fountain activity changed. Specifically, we found fountaining activity that apparently alternated between closer vents on the NW craters. This suggests that the crater wall collapses influenced the dynamics and location of lava fountain activity at those vents located close by. However, looking at the seismic data of these dates, it is evident that the latest crater wall collapse is coincident with a peak in seismic activity and with the transition from Phase 1 to Phase 2 that we identified from the crater analysis. This may lead to speculations about possible joint-driving mechanisms. We conjecture that seismicity changes may be related to pressure changes at depth. As there is only a small delay between these observations, we believe that pressure surges in the underground may be responsible for the observed conjoint changes at the eruption site. Similar concepts were applied for explaining the communication between eruption sites and the caldera at Kilauea [14,61], between the Holuhraun eruption and the Bardabunga earthquakes [62], and between neighboring volcanoes in Kamchatka [63,64]. Such a pressure surge may have initiated the collapse of a crater wall, which then changed the venting activity in the craters, which marked the transition from Phase 1 to Phase 2. Similarly, the oblique dike intrusion that dissects the previously aligned crater row and the development of new and clustered craters may have been initiated by a pressure surge, marking the transition from Phase 2 to Phase 3.

On some occasions, we identified simultaneous activity at more than one crater. On 17 November 2021, a rapid increase in the number of shallow seismic events occurred, escorted by a rapid jump of activity to the SE crater visible in the kymographs (Figure 5b). This is also seen in the SAR ascending data, where we identified two craters only visible in this image (see Supplementary Material Figure S2). From 18 to 28 November 2021, we see that seismic activity decreased in a number of earthquakes and depth. Likewise, we see that the activity of the lava fountains was stable around the central vents. This coincides with the SAR image results, as only one or two new craters were identified within this period. This trend changed on 28 November (Figure 6a), when the seismic activity started to increase again, reaching a peak of 184 events on 30 November 2021. Here, seismic activity reached shallow depths of less than 9 km until 2 December, when activity drastically decreased from 143 to 29 earthquakes. This indicates that our pressure hypothesis may also explain sudden waning activity, when a new dike intrusion took place.

Our SAR results support this; the appearance of the new craters on the NE flank was identified on 28 November (Figure 4m). The activity of these lava fountains is not visible in the kymographs due to the angle at which the images were taken. However, visual inspection of the series of images confirms the observations. In early December, new craters in the SE flank of the cone were identified, and activity of the lava fountains was visible, occurring only days after the peak of seismic activity.

Therefore, the shifts between phases during the 2021 Cumbre Vieja eruption were presumably related to pressure surges and were preceded by increases in seismic activity at depth.

4.4. Conceptual Model

Based on our measurements from SAR amplitude images, time lapse cameras and comparison to seismic data, we developed a conceptual model illustrating the alignment of craters. The eruption was characterized by the opening of vents that produced craters that increased in number and aligned over a length of approximately 650 m in the direction NW–SE. The eruption went through three phases, according to the evolution and opening of the vents illustrated in Figure 7. In Phase 1 (Figure 7a), opening and evolution of the aligned craters were concentrated on the NW side of the Tajogaite cone, initially also hosting a crater in the SE. During this phase, further craters developed to form an alignment of vents dissecting the cone. Just prior to Phase 2 (Figure 7b), collapse events occurred in the NW, followed by the opening of new craters toward the SE during Phase 2. Phase 3 began with the formation of less aligned and more clustered craters, associated with the intrusion of an oblique dike oriented in a north–southernly direction. This generated the connection of the two groups of craters thus far separated and the appearance of another crater row oblique to the main one on the SE side.

5. Conclusions

Our study of the 2021 Tajogaite crater alignment suggests that profound changes of the location and arrangement of eruption craters and lava fountains occur, which is important to understand, as they also strongly affect the direction of lava flows. The 2021 Tajogaite eruption was distinguished by the opening of vents, which developed in an aligned direction NW–SE. We studied a large 48-image dataset of synthetic aperture radar (SAR) amplitude data acquired by the CSK satellite. We found an increase in the number of craters and a decrease in the level of linear arrangement thereof. We compared profound changes occurring at the crater distribution and identified three main phases of activity consisting of (1) an initial linear of craters above the erupting dike in NW–SE orientation, (2) renewed intrusion forming separate craters along the same lineament, and (3) disaggregation of craters due to a new dike intrusion in N–S orientation. Notably, the transition of these phases is marked or preceded by pronounced seismic activity increases. We speculated about the occurrence of possible pressure surges explaining these observations and developed a conceptual model summarizing the surface expression of the craters. Specifically, we conjectured that changes in the crater row are linked to changes in the underlying feeder dike, which may propagate, open or close, or newly orient to form an oblique dike during the late stages of the eruption.

Our study has relevance for other volcanoes as well, as through the analysis of seismic data in relation to the activity of lava fountains, the opening of erupting fissures and emergence of lava flows at the surface can be foreseen in seismic data days to hours prior to occurrence. We did not attempt to forecast locations and expressions of these changes, which appeared rather complex and rapidly evolving during the transition between the phases identified. Importantly, understanding of such a deep to shallow communication may in the future lead to the development of tools where seismic data can be used as an indicator of changes at the surface, and where SAR data rapidly depict crater arrangements and lava flow direction in near real time.