Late complex tensile fracturing interacts with topography at Cumbre Vieja, La Palma

Volcanic eruptions are often preceded by episodes of inﬂation and emplacement of magma along tensile fractures. Here we study the 2021 Tajogaite-Cumbre Vieja eruption on La Palma, Canary Islands, and present evidence for tensile fractures dissecting the new cone during the terminal stage of the eruption. We use synthetic aperture radar (SAR) observations, together with drone images and time-lapse camera data, to determine the timing, scale and complexities associated with a fracturing event, which is diverging at a topographic ridge. By comparing the ﬁeld dataset with analogue models, we further explore the details of lens-shaped fractures that are characteristic for faults diverging at topographic highs and converging at topographic lows. The observations made at Cumbre Vieja and in our models are transferrable to other volcanoes and add further evidence that topography is substantially affecting the geometry and complexity of fractures and magma pathways, and the locations of eruptions.


INTRODUCTION
The propagation of basaltic magma at shallow depths commonly occurs along tensile fractures, forming dikes and sills.These sheet intrusions are of high relevance for shaping a volcano, control the eruption location and build up of topography, as well as the endogenous growth of the edifice [Pollard et al. 1983;Rubin and Pollard 1987;Annen et al. 2001;Segall et al. 2001;Ito and Martel 2002;Kühn and Dahm 2008].Conversely, dike intrusions are directed, oriented or arrested by stress field, anisotropies and topographic loading [Muller et al. 2001;Gudmundsson and Brenner 2004].Therefore, complex interdependencies exist between propagating dikes and the evolving volcano edifice.Dike intrusions occur on large scales (>100 km) associated with plate boundaries and swarm intrusions [Ernst et al. 1995;Gudmundsson and Brenner 2004], on mid-scales (1-10 km) affecting entire volcanoes and volcanic islands [Fiske and Jackson 1972;Pollard et al. 1983;Walker 1987;Carracedo 1999;Walter and Troll 2003], and on very local scales (<1 km) in direct vicinity of an eruption vent [Bonaccorso et al. 2002;Geshi et al. 2002].
Dike intrusions are widespread on the Canary Islands, an archipelago with seven main islands located in the Atlantic Ocean.Dikes in the Canaries have mainly been studied in deeply eroded valleys and excavated road cuts and tunnels, indicative of ancient principal stress fields, rift zones and flank movements [Carracedo 1994;Day et al. 1999;Gudmundsson et al. 1999;Ablay and Martı 2000;Walter and Schmincke 2002;Ancochea et al. 2003;Fernández et al. 2006;Ancochea et al. 2008;Galindo and Gudmundsson 2012;Thiele et al. 2020].Mapping of hundreds to thousands of dikes on La Palma showed that they have a median thickness of ∼1 m at depth [Thiele et al. 2020], which is in agreement with measurements done at other Canary Islands [Marinoni and Gudmundsson 2000].More recently, geodetic and geophysical monitoring of unrest episodes helped to identify dike propagation and widening processes under the most western Canary Islands of El Hierro [González et al. 2013;Longpre and Felpeto 2021] and La Palma [Fernández et al. 2021;De Luca et al. 2022].
Dikes that propagate close to the surface may develop aligned eruption vents and associated subparallel fractures [Witt et al. 2018], and a graben subsidence above the opening dike [Rubin 1992;Müller et al. 2017].The details of these fractures are not fully understood, with models and data indicating that the fracture propagation may be upward [Grant and Kattenhorn 2004] or downward [Mastin and Pollard 1988;Trippanera et al. 2015b].The graben is defined by the inward dipping normal faults, the trace of the erupting vents and the underlying feeder dike [Müller et al. 2017;Silva et al. 2018].At topographic highs, the graben and normal faults may diverge, as prominent examples on Iceland (e.g., Laki Mountain) and experimental studies suggest [Trippanera et al. 2015b].Displacement on the normal faults is enhanced at topographic highs [Trippanera et al. 2015a].Similarly, other fault types, such as strike-slip faults, dip-slip faults and even caldera faults, may also evolve complexities in strike and dip, with sigmoidal features and divergence at topographic highs [Lagmay et al. 2003;Belousov et al. 2005;Wooller et al. 2009;Mathieu and de Vries 2011;Grosse et al. 2020], suggesting a topographic effect for different modes of failure.Faults and especially those associated with tensile fractures and dike intrusions may reactivate pre-existing structures, leading to trend changes and slip reversals [Gudmundsson et al. 2008;Trippanera et al. 2015a;Müller et al. 2017].Some normal faults are also associated with sinkholes, which may form due to unconsolidated materials falling into voids that open during the faulting process [Ferrill et al. 2004].
La Palma is a volcanic island rich in dike intrusions [Thiele et al. 2020], although the first active dike emplacement was instrumentally measured just prior to and during the early stage of the 2021 eruption at the Cumbre Vieja [De Luca et al. 2022].La Palma hosts the most active volcanoes and number of historical eruptions of the Canary Islands [Longpre and Felpeto 2021].Its geological history includes stages of island growth and destruction, so that geologists could identify the main formation of the Garafía and Taburiente volcanic edifices in the north (Figure 1A), succeeded by a southward migration of volcanic activity and, eventually, the emergence of the Cumbre Vieja edifice [Carracedo et al. 1999;Klügel et al. 2005a].The Cumbre Vieja hosts a N-S directed rift zone, the last eruptions occurred in 1971 further south, and in 1949 at the central Cumbre Vieja ridge [Klügel et al. 1999].The 1949 eruption was also associated with formation of a local graben-like fracture (Figure 1B) that has led to speculations about flank instability and a possible collapse [Carracedo et al. 1999;Day et al. 1999].
A new eruption started on 19 September 2021 at the Tajogaite volcano, which is at the western flank of the Cumbre Vieja, just 1.5 km to the north of the vents of the 1949 eruption, and terminated after 85 days on 13 December 2021 [González 2022].The location of the 2021 eruption at the Cumbre Vieja was not foreseen, although a diffuse unrest was identified years before already [Torres- González et al. 2020;Fernández et al. 2021;Longpre and Felpeto 2021].At the location of the eruption, the slope of the edifice was gentle, but a number of older vents, mostly open to the west, were evident [González et al. 2010;González 2022].An important structural trend directed NW-SE, diverging away from the main N-S axis of the Cumbre Vieja, was hypothesised based on structural mapping [Carracedo et al. 1999;Day et al. 1999] and geophysical imaging methods [Camacho et al. 2009;García and Jones 2010;González et al. 2010;Di Paolo et al. 2020].Satellite radar interferometry (InSAR) showed an intrusion just prior the 8-16 September 2021 eruption, followed by shallow dikes near to the eruption location [De Luca et al. 2022].Significant ground displacement was identified mainly in the very early stages of the eruption in September 2021 by geodetic networks and InSAR [De Luca et al. 2022].As no surface ruptures or graben systems evolved, these September 2021 dikes might have been too deep [0.95 to 1.4 km depth b.s.l., De Luca et al. 2022], or the surface expressions might have been obscured by volcanic deposits at the surface.Here, we use satellite and field observations to show that also during the late stage of the eruption, in November-December 2021, a further fracturing event occurred.This initially remained undetected due to the local dimension, but was associated with major changes of the dynamics of the terminating eruption.As we show, an associated graben structure is strongly affected by the volcano topography, with fracture divergence and structural deviations that highlight the complex interaction of topography and magma pathways, and controlled locations of eruption and incandescence.

DATA AND METHODS
In this work, we consider high resolution satellite, drone and time-lapse camera observations to identify changes associated with a local fracturing event at La Palma in early December 2021, just a few days prior to the termination of the eruption on 13 December 2021.Furthermore, we realised sandbox analogue models to simulate tensile fracturing beneath a complex topography, which helps to interpret the surface observations.

Satellite observations
As a reference, a pre-eruptive 5-m-resolution LiDAR digital elevation model (DEM) is considered, which is freely available from the National Center for Geographic Information Systems (CNIG * ), which is part of the National Geographic Institute (IGN) (Figure S1 in Supplementary Material 1).A post-eruptive DEM was derived from Pleiades tri-stereo images acquired on 31 December 2021.The Pleiades images were photogrammetrically processed by us in the Agisoft Metashape (vs.1.8) software package.The rational polynomial coefficient data were employed for georeferencing.The relative image orientation was performed with 5,184 tie points with the further generation of a 36-million-point dense cloud with a total error of 15 cm.The final DEM we used on large scale has 1.5 m resolution, and is shown in Figure 1C.
To access the local structural changes at the Tajogaite eruption site, we used high-resolution Synthetic Aperture Radar (SAR) observations acquired by the COSMO-SkyMed (CSK) mission and provided by the Italian Space Agency (ASI).The SAR instruments onboard CSK satellites are operating in Xband (3.1 cm wavelength).The three most populated CSK datasets for Cumbre Vieja were acquired in HIMAGE mode (stripmap), one on descending orbit, beam H4-08, the other two on ascending orbits, beam H4-06 and H4-12.Those data allowed us to achieve a very high temporal sampling (sometimes daily).In particular, we here use the CSK data before and after an early December 2021 fracturing event.Using the same CSK data catalogue allowed us to explore formation of the crater row [Muñoz et al. 2022].With a spatial resolution of about 3 m (in range and azimuth) of the radar-amplitude data, we are able to look even through eruption ash plumes and clouds.CSK data were coregistered and geocoded using the CNIG high resolution digital elevation model, which, as recalled earlier, is pre-eruptive.However, in the area of interest, the topography changed substantially (over several hundred metres) during the observed period because of the growth of the new eruptive cone; therefore, the terrain corrected data still show some distortion effect (such as layover).Using the post-eruptive DEM would give even worse results; although this would better correct the latest images, it would introduce strong terrain artifacts in the earlier images, thus making their identification in the image sequence more difficult.In spite of such residual image geometrical distortions, it was still possible to visually analyse the relevant stack of data, and to perform a change detection achieved by bandcoloured image stacks, where two channels (red and blue) are represented by two consecutive SAR images.These composite maps allow a graphical representation of the decrease (red channel) and increase (blue channel) in radar amplitude.We also compare these composite maps to difference maps calculated by subtracting the second from the first image amplitudes.Due to widespread ash coverage, radar interferograms showed very high decorrelation effects in the area of interest, thus being of little use and are hence not further discussed in this work.Conversely, the amplitude analysis turned out to be very helpful, to contextualise the timing of any fracture formation.

Time-lapse cameras
Time-lapse cameras were used to record occurrence and location of eruptions and morphological features.The cameras we used were Reconyx UltraFire trail cameras, which we programmed to record images every 5-15 minutes at a resolution of 2304×1296 pixels, 72 dpi.In total, we have set up seven cameras, but due to malpositioning of four of the cameras and theft of one of the cameras, we can only use two cameras, Cam-4 (Figure 1B) and Cam-5 (Figure 1C) for this work.The 5 mm focal length (the 35 mm equivalent focal length is 37 mm) was fixed.The cameras have a night vision in black-and-white (BW) and a day vision in full red-greenblue (RGB) colour.For night eruptions generating strong incandascence, the night vision was set to automatically change to the day colour vision, in order to avoid any overexposure.The cameras were powered by internal lithium batteries, and time-controlled by synchronization to a GNSS (Global Navigation Satellite System) station before the installation, which was controlled again at the end of the experiment for possible time drift correction.By picking up the cameras on 15 Jan 2022 we found a time drift of 31 s and 16 s for Cam-4 and Cam-5, respectively, which is well below the sampling interval and therefore considered insignificant.We strapped the cameras to trees (Cam-4 at 28.606330°N, 17.851424°W, height of 1507 m, heading NW to N60°W; and Cam-5 is at 28.617742°N, 17.860086°W, height of 1087 m, heading SW to N140°W).Cam-4 was positioned at 1500 m distance and high on the often cloudy Cumbre Vieja, which is why only few images were of clear view.A near continuous view is provided by Cam-5 at <700 m distance.For the Cam-5 data we did the analysis by using kymograph recording [Witt et al. 2018], which is found to be an effective way of analysing long photo sequences.To this aim, we generated an image data stack and chose a vertical line through the flank crater in all images (section x-y in Figure 4).We then plot the RGB pixel values of all these sections together along a time axis.The kymograph is therefore a time-space-plot that allows us to quickly identify occurrence, location, and expression of eruptions, as well as morphological changes as seen in the camera images.As kymographs are calculated for day and night images together, the result also shows strong diurnal effects of brightness and colourisation.

Drone observations
Drone surveys were realised when weather (rain, wind strength and direction) and ash fall permitted the safe operation during daylight.Here we use data collected during drone surveys made during the eruption on 12 November 2021 and then again on 5 December 2021, and after the eruption on 15 January 2022.While the first two were done by a Mavic-2 drone, for the latter we used a Phantom 4 RTK.The Phantom 4 RTK survey was also the reference dataset, which due to the high precision of the dual frequency base and radio module yields absolute horizontal position errors on the order of a few centimetres, even in the absence of ground control points [Nota et al. 2022].
The Mavic 2 drone is equipped with a 1" CMOS and 20 MPixel optical camera manufactured by Hasselblad (5472×3648 pixels) and a 35 mm focal equivalent of 28 mm, set to record images at 0.5 Hz interval while flying manually over the fracture.The flight height was ∼80 m above the ground surface.The Phantom 4 RTK drone is equipped with a 1" CMOS and 20 MPixel optical camera as well, a single frequency GPS+GLONASS+Galileo receiver, plus a dual frequency L1/L2 GNSS receiver connecting to our base station fixed near the drone launch site.The drone's positioning accuracy is up to 1.5 cm vertically and horizontally 1 cm.We preprogrammed the flight path of the RTK drone to guarantee constant flight height of 80 m, image overlap of 60-90% and full coverage of the summit craters and the fracture.
In total, we recorded over 4,000 RGB images using the optical drone cameras.The images were manually preselected, excluding repetitive and lower-quality shots.We then processed the data in a Structure-from-Motion (SfM) computer vision approach [Westoby et al. 2012] by using the Agisoft Metashape (vs.1.8) software package.We performed a point cloud reconstruction for each overflight on separate data chunks to reconstruct the three-dimensional morphology and structure.Geoinformatics analysis is then performed in ArcGIS (vs.10.2.1), allowing us to investigate topography, difference, hill shade, and slope, and measure locations, distributions and offsets of fractures and sinkholes (see also Figure S2 in Supplementary Material 1).
Due to the difficult and hazardous terrain while the eruption was still ongoing, we did not collect ground control points to create a reliable georeferencing.Instead, we rely on the GNSS RTK drone for georeferencing, and associate all models with respect to this best and largest dataset acquired after the eruption.This reference model was based on 661 high quality nadir photographs, yielding 175,000 tie points and a dense cloud consisting of 114 million points used to generate a 10 cm digital elevation model and a 5 cm orthomosaic and a total error of 3 cm.This dataset from 15 January 2022 is comparable and acquired with similar drones (Phantom 4 RTK) as the surveys realised by other colleagues between 24 and 28 January 2022 [Civico et al. 2022].
The analysis of our Mavic-2 drone flights was based on 330 (13 November 2021) and 135 (5 December 2021) images acquired at the new fracture zone, yielding 260,000 and 95,000 tie points, 23 and 55 million dense points, DEMs of 16 and 4 cm and orthomosaics of 8 and 2 cm, respectively.Using artificial ground control points picked from the RTK-flight images, these models could be co-located to the RTK model with an error of 45 and 20 cm, respectively.

Analogue Modelling
To better understand the formation of tensile fractures affecting complex cratered topography and the divergent fault trends at the southern slopes of the La Palma eruption site, we performed a suite of analogue experiments.Analogue experiments allowed simulation of the geometry and path of dike intrusions into an elastic material [McLeod and Tait 1999;Takada 1999;Muller et al. 2001;Walter and Troll 2003;Rivalta et al. 2015] and to study faults related to intrusions into a brittle (often sandbox) material [Klügel et al. 2005b;Mathieu et al. 2008;Kervyn et al. 2009;Le Corvec and Walter 2009;Abdelmalak et al. 2012;Trippanera et al. 2015b;Poppe et al. 2019].The sandbox models [Mastin and Pollard 1988;Trippanera et al. 2015b], in particular, were associated with a "process zone" developing brittle faults at the upper dike tips, capable of forming lateral widening and a graben-like subsidence above the tensile fracture.Dike tips can be simulated by a wedge or by a rectangular extension [Trippanera et al. 2015b], both of which develop similar graben faults at the surface.At Cumbre Vieja, as we do not have constraints on simulating a dike of uncertain geometry, here we keep the number of unknown variables low by employing a simplified experimental approach by simulating tensile failure only.This aims to represent the extension near to the upper dike tip, that were quantified in the field and in models to be horizontally dominantly [Baer 1991;Rubin 1993;Gudmundsson and Loetveit 2005].Our deformation driving mechanism was thus horizontal extension of two basal plates, mimicking tensile faulting.Models were performed by simulating both constant height (flat) topography and ridge-like topographies that are pulled apart.Only the ridge-like topographies yielded divergent fault pattern similar to the field observations.Our experimental table had dimensions of 50×50 cm plane, and allowed topographic height simulations up to 12 cm, at mean slopes of 25°(Figure S3 in Supplementary Material 1).We also tested steeper slopes close to the repose angle, but neglected those due to frequent sliding processes.Since we use a Mohr-Coulomb material, the extension rate does not significantly affect the faulting pattern.Experiment interruption and restarting also had no visible effect on the fault growth due to the (largely) strain rate independence and our simplified design with a lack of viscous materials.Two DSLR cameras monitored horizontal and vertical displacement every 1 mm extension.The images with 4928×3264 pixels from a distance of ∼40 cm yield a resolution of 20 px mm −1 .
As analogue materials, we used dry rounded quartz sand with ∼400 µm grain size diameter [Rosenau et al. 2018].In order to approach higher cohesion we mixed this sand with 10 wt% plaster powder <250 µm diameter, creating a commonly used analogue for crustal volcanic rocks [Brothelande and Merle 2015;Brothelande et al. 2016], which are also found ideal to visualise brittle faulting and structural features under compression [Zorn et al. 2020] and extension [Poppe et al. 2021].Sand-plaster mixtures allow us to control material mechanics, larger plaster amounts decrease density, but increase tensile strength and cohesion [Poppe et al. 2021].This effect comes with more pronounced vertical fracture heights, allowing the experimental study of fault-related offsets in greater detail.For mixtures with <20 wt% of plaster, the friction coefficient typically ranges from 0.71 to 0.81 [Poppe et al. 2021], equivalent to an angle of internal friction of 35-39°.
The experiments were scaled to allow geometric and kinematic similarity between modelled and natural processes.Our granular material scales rate-independent and similarity therefore mainly depends on the friction coefficient and the cohesion [Merle 2015].If the cohesion is very low, such as for our pure sand material, scaling can be reduced to a geometric ratio [Brothelande and Merle 2015].For mixed cohesive materials, by employing a scale factor using the stress ratio σ * [Merle 2015] one may scale the cohesion of the analogue model to be c analogue = σ * c nature , with σ * = ρ *  * ℎ * , where ρ * ,  * , and ℎ * are the density, gravity, and length ratios, respectively.By unchanging gravity, the length ratio and the density ratio become σ * = ρ * ℎ * .In this respect, for the E-W elongate topographic ridge with a dimension of 140 m along N-S and a height of 40 m from its base, and a fracture 50 m wide and 170 m long, the scaling factor is 1×10 3 , so that 1 m in nature scale to 1 mm in experiments yielding an analogue ridge of 140 mm in N-S and 40 mm height.The stress scaling factor accordingly suggests that the analogue ridge is approximately 1000 times weaker than in nature, which for the chosen materials is in the scaled range as also shown in Zorn et al. [2020].Further details on the scaling and mechanical tests of the used materials are provided in Zorn et al. [2020] and in Poppe et al. [2021].
The experiments recorded by cameras were further analysed using the digital image correlation (DIC) approach [Pan et al. 2009].By this, we are able to track the deformation in an arbitrary reference (scaled) frame, and monitor displacement and strain rate.In DIC the images are searched for changes in a moving sub-window, where a correlation peak is found if corresponding sub-windows match.Changing positions and dimensions of these allow determination of displacement and strain (Figure S3 in Supplementary Material 1).Similar techniques have been previously applied in analogue models and in volcano monitoring [Le Corvec and Walter 2009;Burchardt and Walter 2010;Walter 2011;Trippanera et al. 2015b;Zorn et al. 2020].We are using LaVision DaVis (vs 10.0.5), and map the 2D displacement field in x (Ux) and in y (Uy) direction, and calculate the horizontal axial strain Exx to identify tensile fractures.

Timing of graben subsidence
The emergence of a set of new fractures were noted and first reported by one us (PJG) to the volcano emergency team (Plan de Emergencias Volcánicas de Canarias -PEVOLCA) around 14:00 local time on 03 December 2021.Identifying the exact timing of the graben formation was challenging in the field, however, because the location was not visible from common viewpoints and often inaccessible due to frequent air fall and gas hazards.
The chronology of events can be well accessed by considering both our CSK and time-lapse camera data.The regular recording of the CSK satellite radar images shows the evolution of the near-linear NW-SE aligned crater row until 12 November 2021 [Muñoz et al. 2022].We can identify the first ure 3A).The first visible change on that slope is identified in an image on 02 December 2021 at 15:15 (UTC), showing a zone of broad subsidence forming a gentle depression affecting the ridge (Figure 3B).The next image is taken only 10 minutes later on 02 December 2021 at 15:25, and reveals the presence of a sharp fault trend (Figure 3C).Later images basically confirm the sharp fault trend, indicative of a rapid fracture event.Therefore, Cam-4 reveals that the fractures evolved gradually, firstly causing broad subsidence and (10-minutes) later showing trends of fractures.The fractures converge on the left (south) and diverge uphill.A number of nighttime images of Cam-4 shows an incandescent spot on the southern part of the fracture (Figure 3D), indicative of high temperature or even lava close to the surface (no lava flow emitted at this location though).
The other camera located close to the northern open vent is Cam-5, and reveals strong changes at this northern site (Figure 4A).This, we can best study using kymographs along a vertical line at the location of the later flank open vent eruption (Figure 4B).The timeline starts in November 2021 and shows a generally low level of activity during day and night images.Occasionally eruptions were strong enough to leave bombs on the slopes (e.g., 20 November 2021), but activity along this line seems low.This abruptly changes in the night of 28 November 2021, when the first flank eruption starts.The eruption then becomes incandescent enough to cause the camera to switch from nighttime to day view settings, seen by the red glow in Figure 4.The flank eruption lasted six days, starting on 28 November 2021 at 03:15 and ending on 03 December 2021 at 08:00.The kymograph also shows the silhouette of the growing cone (dashed line in Figure 4B), indicating a growth of ∼40 m, although the main growth episode is over on 02 December 2021.Image by image inspection reveals that at the beginning of this episode, the eruption occurs at several distributed short-lived vents (Figure S5 in Supplementary Material 1), implying fracturing of the northern edifice slopes four days prior to the formation of the lens-shaped fracture seen 400-500 m further to the south (Figure 4C).

Fracture geometry from drones
The fracture geometry developing on the southern slopes was hardly distinguishable in every detail by time-lapse cameras.
Drone data allows a much more detailed morphological and structural characterization.Our drone data acquired before the fracture formation were also before the northern open vent formation, showing the aligned craters and slopes of mainly smooth surface due to ash and lapilli air fall in this proximal region (<400 m from the main eruption vent).On the southern slopes, we identify a west-east trending morphological ridge (Figure 5), which is already seen in pre-eruption digital elevation data (Figure S1 in Supplementary Material 1).Closer comparison of both sets of drone data, 12 November 2021 and 05 December 2021, reveals that prominent trees and boulders can hardly be traced suggesting that the morphology of the southern ridge, although being largely unchanged, has been covered by loose clastic materials due to air fall (which makes co-alignment of the data challenging).
The 12 November 2021 drone photomosaic data show the westerly directed plume and red glow in the central crater (Figure 5).All active craters are aligned NW-SE.The northern and southern outer slopes appear smooth in hill shade images due to air fall deposits, although pre-existing topography of a buried scoria cone is still seen in the south (Figure S1 in Supplementary Material 1).The 05 December 2021 drone photomosaic shows a much stronger and wider cratered region.Most prominent also in the hill shade is a flank crater open to the NE in the north, and the fracture in the south.Both are aligned NNE-SSW, but a direct link can not be determined from the morphologic drone data.
A closer view of the fracture for the 05 December 2021 data is provided in Figure 6.Shaded relief, DEM difference and slope maps show the details of the lens-shaped fracture, exposing sinkholes, vertical subsidence and sharp associated fault traces (Figure 6A-C).The area of the lens-shaped fracture and sinkholes is 50 m wide and 170 m long.In the northernmost part, the fracture was notably the site of a small eruption vent, which was long enough to feed a lava flow on 4-5 December 2021 (LF), but too short to build up a new crater or cone (Figure 6B).In fact, the fracture is defined as a divergent fracture system, being narrow at topographic lows, and wider at the elevated region of the topographic ridge.DEM difference maps (Figure 6B,E,H) reveal large vertical throw of the fracture system.We find maximum vertical fault offset to be 4 m, then decreasing to <1 m at the northern and southern edges (Figure S2 in Supplementary Material 1).Thus, higher values are found at higher parts of the topographic ridge and where co-located sinkholes can be identified (Figure 6A-B).
Slope maps better reveal the arrangement of the fracture system (Figure 6C,F,I), showing that the main fractures are escorted by sub-parallel small fractures that are on the order of 10 m long.Some of these arrange as mini-grabens (Figure 6D), where the central part 2-4 m wide has subsided.Numerous sinkholes established in the very southern part of the fracture system (Figure 6D-F).They are 3-4 m in diameter and even located outside the graben, becoming larger and nested within the graben, there exceeding 10 m diameter.The largest sinkholes are found on the northern part of the divergent graben, with diameters of 25 and 32 m, respectively (Figure 6G-I), and depths of up to 4.5 m (Figure S2 in Supplementary Material 1).A view of the main faults (Figure 6H) reveals other complexities, such as left or right-stepping en-echelon fracture sets (Figure 6I).These en-echelon fracture sets are aligning and present only at the lower proportions of the ridge, and may [G, H, I] Zoom to the northeastern edge reveals that the fault is crossed by en-echelon-like fracture sets, here with minor vertical displacement only.
imply involvement of strike slip faulting during subsidence.Later drone surveys show that the relief is eroding quickly [Civico et al. 2022].

Graben evolution from modelling
The morphological analysis of the drone data reveals the presence of a pronounced graben, diverging towards the highs of a pre-existing topographic ridge.While associated local deformations seem complex, no major vertical offsets of the graben shoulders were identified, suggesting a horizontal widening source mechanism.To this aim we designed analogue models, simulating horizontal widening beneath a topographic ridge, which may represent a widening mechanism associated with tectonic faulting, dike intrusion or gravity driven extension.Our analogue models allow us to investigate chronological and structural changes in homogeneous materials (Figure 7) and also in a layered material to increase sinkhole formation (Figure S4 in Supplementary Material 1).The first tensile fractures visible in the analogue experiments were near the base, whereas fractures appeared at the higher ridge at a later stage only.The deformation maps show that already at this early stage, although no fractures are visible, a divergent lens-shaped region is mechanically decoupled from the western and eastern block.This is implying that faults are pronounced at depth, but not yet visible at the surface.With continuation of extensional faulting at the base, we find that the lens-shaped region becomes narrower, and better defined by fractures.These fractures have clear vertical throw, and connect to the linear N-S fault experimentally simulated at the base.Therefore, the lens-shaped faults are inward dipping normal faults.An N-S displacement calculation (Uy) reveals also a movement orthogonal to the main driving process.The wedges enclosed by the converging fractures in the north and in the south are laterally displacing.The E-W strain calculation reveals the exact position of the faults reaching the surface, for which we find en-echelon side stepping faults in the northern and southern wedge, in agreement with a strike-slip component.We repeated the experiments with cohesive materials and overlain by low-cohesive granular materials, and find a stronger expression of fault step-overs and local direction changes, and an increase of the dimension and number of sinkholes (Figure S4 in Supplementary Material 1).

DISCUSSION
We identified a fracture event in the late stage of the 2021 Cumbre Vieja eruption on La Palma and investigated details on timing and expression, proposing a link to the eccentric 6day eruption in the north side.This period was also showing strong changes as recorded in other data, such as at seismic data a peak of the number of located events and shallow seismicity occurred from 29 November 2021 * , and also the opening of new short-lived vents downslope away from the main cone and towards the west [González 2022].

Limitations
Continuous observations of surface structures and morphology were challenging due to rapidly changing eruption dy- * https://www.ign.es/web/vlc-serie-palmanamics and weather effects hindering clear views.We could produce repeat drone flights, but due to high winds and ash dispersion, and unflexible drone permits often limiting the flights to a single group at a certain time and area, acquisition of complete data and well planned nadir photographs remained challenging.Therefore, the flights realised by Mavic drones may also contain larger drifts and processed data contain gaps and variable quality.For instance, the drone flight realised on 05 December 2021 was during a day with high winds, near ground wind speeds exceeding 30 km h −1 , and higher winds at drone flight elevations.Therefore, even a hovering drone was drifted by the wind, so that proper flight planning was elusive during these days.
Crowd-sourcing observations contributed to better constrain the evolution of the eruption [Wadsworth et al. 2022], but as the fracture opening we describe occurred on the backside hidden to most observers, instrumental monitoring and remote sensing was most valuable for our study.Satellite radar has the benefit of looking through clouds, even eruption clouds, providing regular observation of Cumbre Vieja during bad weather and strong ash producing episodes.However, spatial resolutions are low and repeat intervals of the satellites are relatively poor compared to the fast transitions occurring at volcanoes, so that details of the occurrence of the fractures remain vague.Here, we showed that time-lapse camera records provide added value, especially if views from multiple directions can be realised.We installed cameras, and found them highly useful for identifying details of the evolving fracture and for determining the timing of the event.
We run analogue experiments to simulate fracture divergence caused by topographic effects.This idea and approach is not new, and was used to explain large scale structural changes at cone-shaped topographic mounts [Acocella and Neri 2009;Wooller et al. 2009].Earlier studies showed that dikes change their path towards topographic highs [Maccaferri et al. 2017] and they curve around topography scars [Acocella and Tibaldi 2005].Fault divergence can cause volcano cones to become unstable and produce collapses [Belousov et al. 2005].Moreover, graben fractures evolve above dikes [Rubin and Pollard 1988] and diverge at topographic cones [Müller et al. 2017], where also asymmetric fault patterns may evolve [Trippanera et al. 2019].Recent analogue models could reveal that the evolving graben system may be rather asymmetric and complex [Abdelmalak et al. 2012] and strongly diverge at topographic cones to form lens-shaped fracture sets [Trippanera et al. 2015b].In this work we build on these previous studies, and considered realistic topographies, to show fracture divergence and convergence depending on topographic highs and lows, respectively.Our models are held simple, simulate not the dike but the extension above, and hence do not consider possible complexities at the dike itself.Our models are constrained at depth, implying a constant depth of the top of the dike, which may be too simplistic.However, our model approach is producing an extension and structures that are comparable to the flat topped dikes described in Trippanera et al. [2015b].We simulate extension faulting only, i.e., any larger dip slip and strike slip component of the driving fault was not further elaborated.Furthermore, our models ig- nore rheological complexities and time dependence.Despite this simplistic approach, the geometric similarities found to the real case on La Palma are striking, giving confidence that the basic principle is correct.

Importance of graben for the 2021 Cumbre Vieja eruption
The structural features seen in satellite radar, drone and timelapse camera data during the late La Palma fracture are (i) a gradual deepening and faulting (Figure 3), (ii) association to eruption and incandescence (Figures 3 and 4), (iii) divergent geometry (Figure 5), and (iv) sinkholes and strike slip mechanisms (Figure 6).Even though our analogue models are simplified, we could show that normal faults develop at depth and migrate to the top, while tensile fractures form at the surface and migrate downwards, which is seen by a firstly diffuse lens-shaped deformation before inward dipping faults cut through the models.Strike slip indicators develop in response to the topographic effect, and do not necessarily imply the source (dike) to be transtensional.Our models show frac-ture divergence only if a topographic high is present, which we simulated in the form of a topographic ridge trending E-W.At topographic lows, the fracture converges again to form a single tensile fault.These are the locations of incandescence and eruption, suggesting magma presence and propagation is at the level of the tensile fracture (where the divergent graben faults merge) rather than at the graben faults above.All in all, these findings strongly suggest that the fractures formed due to a dike intrusion, mainly with tensile widening, while the complex topography has led to complications and fault divergences as described by us.However, we cannot rule out the involvement of alternative mechanisms, such as a pure tectonic or gravitational failure that is simply being used as a magma pathway.
Conjecturing that the interpretation of a shallow dike intrusion occurrence on 02 December 2021 is correct, the temporal and structural association to the eccentric eruption at the northern flank becomes particularly interesting.At the surface, the sequence of events started with an eruption change in the north, and opening of multiple eruption vents on the former flanks of the cone (Figure S5 in Supplementary Material 1).The eruption change to the north was followed by propagation of the dike to the south, as the lens shaped graben and localised eruption sites (incandescence) occurred.Interestingly, the general trend of the eccentric vent open to the north is almost the same as the dike induced fracture in the south.Deviations might have been caused by the strong topographic expression not only in the south, but also in the north, but also by trend changes of this dike.
Hypothesising that this interpretation is correct, that the eccentric eruption and the southern lens shaped fractures were caused and fed by the same dike intrusion, a north to south progression or propagation may be indicated.It took four days from the eruption in the north on 28 November 2021 03:15 to the fracture formation, incandescence and small eruption in the south on 02 December 2021 15:15, implying that the dike length in total might have been 600-700 m.
Comparison to other observations underlines that the timing of the fracturing event outlined in this work is also temporally associated with seismic activity changes.Seismicity occurring at two main depth levels, at 34 km and at 12 km depth, might be associated with magma storage regions.This behaviour is sporadically interrupted by bursts of larger magnitude earthquakes, occurring mainly in the shallow region at 12 km depth [del Fresno et al. 2022].During the period we ob-served the occurrence of the fracture set, also shallower earthquakes were identified.In fact, a comparison to available seismic catalogues [del Fresno et al. 2023] * reveals a chronologic succession, with large earthquakes at the deeper reservoirs at 34 and 12 km depths, followed in the end of November 2021 and early December 2021 by more shallow earthquakes at <12 km depth, which is followed by the new eccentric crater on the northern flank of the main eruptive crater and then by the divergent fracture on the southern side as described in our work.In this view, the vertical time lag between earthquake and eruption changes is less than one day, whereas the horizontal time lag between the northern vent and the fracture emergence in the south is up to four days.
Possible drivers may be (i) increased magma flux causing earthquakes and then reaching the eruption site, there, causing a late dike intrusion, or (ii) tectonic faulting, possibly associated with the growing load of the edifice, which is then serving as a pathway for propagating magma.Understanding the particularities of the drivers and interrelations of the deep and lateral changes observed are of fundamental importance for understanding eroded outcrops of older volcanoes and also for improving rapid assessment of changes during volcano eruptions elsewhere.This is essential also for assessing new lava flow hazards, that had completely different pathways for the eruption sites on the northern flank and at the southern flank.

Complex topography
Volcanoes evolve a highly complex topography that is continuously changing during the course of an eruption.Also, the Cumbre Vieja 2021 eruption led to a topography much more complex than simulated by our lens-shaped models.Therefore, in an attempt to consider the complex topography, we repeated the same experiments described above, but incorporated a topography that included the main morphological features, such as the topography ridge in the south, the central eruption cone, and the central eruption crater.The experiment material is homogeneous, the morphology is manually constructed in our sandbox, which means that the same exact geometry can hardly be reproduced.Results are recorded again by camera and analysed using DIC to test the strain normal to the tensile fracture (En) and the displacement along the fracture (U).Results are shown in Figure 8, and demonstrate the presence of the lens-shaped fracture across the ridge in the south.This is slightly asymmetric due to the oblique trend of the tensile fault with respect to the ridge, with length and strength, causing a parallelogram-shaped fault geometry.To the north, the faults converge at the 4-5 December eruption site, and diverge again towards the main eruption cone.At the location of the central crater, a tensile fault geometry curves around the crater walls, and bends on the northern outer rim to converge again.We need to emphasise, this complex fault geometry is driven by a single linear tensile fracture at depth.The very northern margin converges, there locally mismatching the open vent, although we note that the northern continuation of this was the site of numerous short term flank eruption sites and fractures (seen in the time-lapse camera on 28 November 2021 at 07:45 UTC; Figure S6 in Supplementary Material 1).Although we believe that the models can help understanding and locating the feeders and new vents, eruptive deposition and scoria cone or crater formation was not attempted to be simulated.

CONCLUSIONS
We identified a fracturing event during the late stages of the La Palma eruption.By combining satellite radar images and time-lapse camera records in the field, we could determine the timing of the fracture formation occurring on a topographic ridge on 02 December 2021 between 15:15 and 15:25.In fact, we find that initial broad subsidence occurs, which then develops into main fractures.During the same period, a flank eruption occurred on the other side of the edifice, which may also be genetically related to the same fracture set that lead to incandescence and a small lava flow on the south side.This leads us to speculate the fracture was driven by a tensile dike opening.Using high resolution drone imagery we were able to determine subsidence and strike slip kinematic indicators as well as sinkholes at the lens-shaped fracture in the south.We designed analogue models to show that a single and linear tensile fracture at depth diverges upwardly to form broad subsidence and then well-defined lens-shaped inward dipping normal faults at the top of the ridge.We further speculate if this fracture to the south and the eccentric eruption to the north are related.In this view, we show that the complex topography in a cratered volcanic summit region is leading to very complex fracture geometries, which is of high relevance for understanding volcanoes and their evolution, stability and propagating or reorganizing vents elsewhere.

AUTHOR CONTRIBUTIONS
T.R.W. initiated the study, contributed to writing the manuscript, performed the drone and analogue modelling analysis.E.Z. contributed to writing the manuscript, processed drone data, and contributed to analogue modelling.P.J.G. contributed to writing the manuscript and coordinated the field activities.E.S. contributed to writing the manuscript and together with D.R. processed the SAR data.V.M. contributed to the field work.A.V.S. and S.P. analysed and processed the Pleiades data.T.R.W., E.Z., P.J.G., N.R., V.M., A.V.S. contributed to the field work, drone flights, and camera setups.

Figure 1 :
Figure 1: Study area at western Cumbre Vieja, La Palma (Canary Islands).[A] La Palma is located at the western Canary Islands, the 2021 eruption occurred on the southern edifice Cumbre Vieja.[B] Pre-eruption lidar shaded relief map, final lava flow outline (red shaded area, Copernicus data) and 1949 fault trace (black line).[C] Zoom-in to eruption location, post-eruptive (31 December 2021) Pleiades tri-stereo shaded relief map.Note the location of the crater row, the December 2022 fracture and open vent.

Figure 2 :
Figure 2: Satellite radar change analysis from image dataset acquired by COSMO-SkyMed (orbit H4 08D) just before and after the December 2021 fracturing event.[A] Composite maps of two consecutive radar images, the primary image is shown in red channel, the secondary image shown in blue channel, dates are provided above the images.[B] Radar image difference in gray scale, for the same dates.Analysis suggests fracture occurrence in early December 2021 (see also Supplementary Material 1) was associated with a shift in eruptive activity and new lava flows to the south and north side of the main crater row.Note the orientation of the fracture aligns with the open vent to the north.R stands for the satellite range direction, LF -new lava flow.

Figure 3 :
Figure 3: Camera records, realised by fixed installed time-lapse trail Cam-4 (for location, see Cam-4 in Figure 1B).The selected images allow us to determine approximate timing and evolution of the fracture.[A] View onto the eruption crater, also visible is the steaming of the new open vent to the right.Close up shows gentle slope.[B] Limited visibility, but close up (right column) onto the gentle southern slope illustrates the emergence of a broad subsidence region.No fractures are visible yet on 02 December 2021 at 15:15 (UTC).[C] Strongly degassing volcano, close up (right column) shows the presence of knife-sharp fractures on 02 December 2021 at 15:25.[D] Night time image shows incandescent spots on southern part of the fracture.Due to limited visibility, other images for the period were not found useful.

Figure 4 :
Figure 4: Camera records, realised by fixed installed time-lapse trail Cam-5.[A] Example images in November and December 2021.Note the outline shows the edifice growth.The section x-y indicates the line for the kymograph.[B] Kymograph shows night (greyscale) and day images (showing the eruption in red color).The eruption abruptly changed on 28 November 2021 to a flank eruption that abruptly stopped again on 03 December 2021.Associated with this short 6-day episode, a 40-m topography built up.The stop of this episode correlates with emergence of the lens-shaped fracture in the south.[C] Example images for the flank eruption highlights a larger number of vents at the beginning, which eventually develops into a major fountaining eruption site.

Figure 5 :
Figure 5: Drone-based orthomosaic (left) and shaded relief map (right).[A] Data acquired before the fracture event reveals the crater row and presence of the southern ridge that is a pre-eruption morphology.[B] Data acquired after the eruption shows the lens-shaped fracture set on the south side, and the flank eruption site to the north, located eccentric to the crater row.

Figure 6 :
Figure 6: Drone-based results close-up to the lens-shaped fracture as of 05 December 2021, showing the shaded relief maps (left panels), the DEM difference (centre) and the slope maps (right panels).[A, B, C] The fracture is a divergent feature, the largest widening at high topography.Sinkholes appear.[D, E, F] Zoom to the southern edge (for location, see black box), revealing the presence of small sinkholes and subparallel mini-grabens.Vertical difference (middle column) shows that largest subsidence occurs at this place.Slope maps reveal 4-5 subparallel fractures following the trace of the main fault, but with only minor offset.[G,H, I] Zoom to the northeastern edge reveals that the fault is crossed by en-echelon-like fracture sets, here with minor vertical displacement only.

Figure 7 :
Figure 7: Modelling results on evolution of lens-shaped fractures.[A] Experimental setup, showing that a topographic ridge 40 mm high and wide, 140 mm wide (in y-direction) and 300 mm long (elongated in x-direction of the experiments) was subject to tensile faulting at the base.[B] After 2 mm of extension, horizontal displacement in x-direction (Ux) shows diffuse graben.Displacement in y-direction (Uy) shows no relevant deformation.Strain in x-direction (Exx) reveals fault locations.[C] After 4 mm of extension, faults become more clear, the graben narrows, Uy shows strike slip movement, and Exx reveals the presence of en-echelon faults and inactive faults in the outer region.

Figure 8 :
Figure 8: Conceptual model of the late La Palma fracture and complex topography effect.[A] Drone derived high resolution shaded relief map and the tensile fracture (dike) hypothesised in this work, connecting the northern eccentric eruption and the lens shaped fracture in the south.[B] Analogue model with a complex morphology, considering a ridge in the south, and a high cone/crater in the centre, may explain the structures observed.Northern continuation (question mark in [A]) showed multiple short term vents.Color contours show normal strain (En) in the direction perpendicular to the tensile fracture.See discussion section for details.