Banding in the Margins of Basaltic Dykes Indicates Pulsatory Propagation During Emplacement

Basaltic fissure eruptions, which are the most common type of eruption on Earth, are fed by dykes which mediate magma transport through the crust. Dyke propagation processes are important because they determine the geometry of the transport pathway and the nature of any geophysical signals associated with magma ascent. Here, we investigate small‐scale (mm–cm wide) banding features at the margins of dykes in the Teno Massif (Tenerife, Spain) and the Columbia River Basalt Province (CRBP) (USA). Similar marginal bands have been reported for dykes in numerous localities around the world. Dyke margins record valuable information about propagation because they are the first material to solidify against the host rock at the propagating dyke tip. We find that the marginal bands are defined by cyclic variations in phenocryst concentration and vesicularity, and we infer that these cyclic variations in texture are a product of cyclic variations in magma flow rates and pressures within the dyke tip. This indicates that dyke emplacement occurs in pulses, with propagation repeatedly hindered by the rapid cooling and solidification of magma in the narrow dyke tip. Using a 1D conduction model, we estimate the time taken for each band to cool and solidify, which provides a timescale of several minutes to tens of minutes for the pulses. The occurrence of similar bands in various volcanic settings suggests that pulsatory propagation is a common, if not ubiquitous, process associated with dyke emplacement.


Introduction
Basaltic feeder dykes are planar conduits that feed fissure eruptions.Flow processes within dykes influence eruptive behavior at the surface, such as the localisation of fissure eruptions (Bruce & Huppert, 1989;Delaney & Pollard, 1982;Jones & Llewellin, 2021).In order to forecast potential propagation pathways and to determine which areas might be most at risk from eruptions, we need an understanding of dyke emplacement, from initial propagation to late-stage flow processes.Subsurface flow cannot be observed directly, but the rates and overall direction of dyke propagation have been inferred using the migration of seismic activity, with typical time-averaged propagation speeds inferred to be tens of centimeters per second (e.g., Einarsson & Brandsdóttir, 1980;White et al., 2011;Woods et al., 2019).Rock textures from ancient, solidified dykes can complement the geophysical data from active systems as an indirect means to investigate subsurface processes (e.g., Philpotts & Asher, 1994;Philpotts & Philpotts, 2007;Wada, 1992).For example, many studies have used rock textures or magnetic fabrics at dyke margins to infer initial propagation directions (e.g., Baer & Reches, 1987;Staudigel et al., 1999).
The margins of a dyke are the first material to solidify against the host rock, so they must originate at the leading edge of a dyke as it propagates through the crust.Marginal textures are therefore intrinsically related to the propagation process, and provide a record from initial crack propagation toward sustained magma flow within the dyke.Numerical models of dyke propagation, whether 2D (e.g., Maccaferri et al., 2010;Taisne & Jaupart, 2009) or 3D (e.g., Davis et al., 2021;Zia & Lecampion, 2020), operate within the framework of linear elastic fracture mechanics (LEFM), where propagation is controlled by the geometry of the dyke tip, and either the fracture toughness of the medium, or viscous losses from the fluid (Gonnermann & Taisne, 2015).However, despite LEFM models requiring a tapered dyke tip geometry to concentrate stress and instigate fracture, most overlook the role of solidification, even though the narrow tip must be prone to rapid cooling.Marginal samples provide information about the role of cooling and solidification at the tip of a propagating dyke, and about the geometry of the tip itself.This provides a framework for interpreting geophysical evidence of dyke propagation, and for informing numerical propagation models.
The aim of this work is to interpret variations in magma flow and dyke propagation from the analysis of dyke margin features, in the context of this dyke propagation framework.We focus on two samples which show banded textures at their margins: one from the Teno Massif (Tenerife, Spain) and one from the Columbia River Basalt Province (CRBP) (USA).We view these textures as a time series, with the oldest material at the dyke wall, and progressively younger material toward the dyke center.The samples therefore provide a record from initial propagation to sustained magma flow, and reveal the dyke tip flow conditions associated with propagation.
Marginal bands are typically distinguishable by eye through variations in phenocryst concentration, groundmass grain size, and/or vesicularity.Where bands are defined by phenocryst concentration, they can appear as porphyritic and non-porphyritic zones, sometimes associated with a gradation in phenocryst size (Drever & Johnston, 1959).Alternatively, where bands are defined by groundmass grain size, they appear as a series of discontinuities, involving a drop in grain size moving toward the dyke center (Holness & Humphreys, 2003).Other bands are defined by narrow, highly vesicular or amygdaloidal regions running parallel to the dyke margins (Delcamp et al., 2012;Galindo & Gudmundsson, 2012;Kavanagh et al., 2018;Kile, 1993;Platten, 2000;Thiele et al., 2021;Thordarson & Self, 1996;Walker & Eyre, 1995).Finally, some authors describe platy textures and color variations at dyke margins, the cause of which is uncertain (Coward, 1980;Healy et al., 2018;Roberts & Sanderson, 1971).
Banding is usually restricted to within approximately 15 cm of the dyke margin, especially for bands defined by phenocrysts and groundmass variations (e.g., Drever & Johnston, 1959;Holness & Humphreys, 2003).However, vesicle bands have been observed further toward the centers of intrusions; for example, in the CRBP, highly vesicular bands are present all the way from the margins to the dyke center, increasing in width from 1 to 10 cm (Thordarson & Self, 1996), and vesicular bands have also been observed across the entire widths of dykes in American Samoa (Walker & Eyre, 1995).

Formation of Bands and Layers Within Dykes
Marginal bands are distinct from the large-scale layering commonly found within dykes, which is often tens of centimeters to several meters wide.The large-scale layers are defined by distinct compositions (Eriksson et al., 2011) or jointing patterns (Gudmundsson, 1984;Ray et al., 2007), or by a textural change between layers, evidenced by chilled margins or melt-back margins.A chilled margin forms when magma cools rapidly against its surroundings, forming a glassy or fine-grained layer (e.g., Brouxel, 1991;Platten, 2000), whereas a melt-back margin forms when magma thermally erodes into older dyke material, then cools slowly, forming a coarser texture (Huppert & Sparks, 1989;Svenningsen, 1994).In either case, a substantial time gap is required between layers.As such, large-scale layering is inferred to form via multiple, discrete injections of magma; for example, where felsic magma has been injected into the center of a pre-existing basaltic dyke (e.g., Eriksson et al., 2011;Tomé et al., 2020).Large-scale layers can therefore be viewed as products of distinct intrusion episodes, with a unique emplacement history for each dyke.This is very different from marginal bands, which always show a general trend of widening and becoming less distinct toward the dyke center.Marginal bands are not, therefore, small-scale versions of a composite dyke.
The most common interpretation of marginal bands is that they represent variations in flow rates within the dyke, attributed to the repeated injection of material (Drever & Johnston, 1959;Platten, 2000).Other explanations of marginal banding include each layer representing a plane of movement at the edge of the dyke (Roberts & Sanderson, 1971).Bands defined by vesicles alone have been interpreted to result from pressure fluctuations and temporal variations in flow (Kavanagh et al., 2018;Kile, 1993;Thiele et al., 2021;Walker & Eyre, 1995), or preferential degassing pathways (Delcamp et al., 2012).However, few reports of marginal banding include observations of internal chilled margins; in fact, their absence has been noted (Platten, 2000).Therefore, although some bands may be defined by drops in grain size short of an internal chilled margin (Holness & Humphreys, 2003;Platten, 2000), it seems that band formation is too rapid to allow long cooling intervals between them.

Field Sites and Samples
Two samples displaying prominent marginal banding are investigated to determine the nature of the bands: one from the Teno Massif (Tenerife, Spain); and one from the CRBP (USA).We refer to these samples as the "Teno sample" and the "CRB sample" respectively.

Teno Sample
The Teno sample comes from a 0.5-m-wide basaltic dyke in the Teno Massif (Figure 2), and comprises the outermost 10 cm of the dyke, with six bands recognizable by eye.Banded margins were seen on several dykes in the immediate area, usually defined by a dark-to-light color variation, and often associated with variations in vesicularity (Figure 1).The Teno Massif is the eroded remains of a basaltic shield volcano active between 5.0 and 6.3 Ma (Leonhardt & Soffel, 2006;Longpré et al., 2009).Most dykes in the region are single intrusions, but multiple and composite intrusions are common, showing large-scale layering with chilled margins and distinct compositions.Dyke tips are commonly exposed, typically tapering into a crack-like geometry (Marinoni & Gudmundsson, 2000).

CRB Sample
The CRB sample comes from the margins of a 0.95 m wide basaltic dyke in SE Washington State.The dyke is one of many that fed the 14.7 Ma Roza lava flow, which was emplaced over at least 14 years and covers an area of 40,300 km 2 (Thordarson & Self, 1998).The sample comprises the outermost 12 cm of the dyke and contains nine bands distinguishable by eye (Figure 3).Many of the Roza dykes display bands of vesicles at their margins, symmetrical on either side of the dyke, and some contain up to 18 apparent bands.The inferred depth of emplacement is 500 m.

Methods
Thin sections were taken from both the Teno sample and the CRB sample, so that their textures could be analyzed from the margins inwards.For the CRB sample, which had a known orientation, thin sections were taken on the vertical plane, perpendicular to the margin.Images were collected at ×5 magnification under a petrological microscope, then stitched together using Inkscape to form image swaths extending across the width of the thin section (Supporting Information S1).Due to the different color and brightness contrasts in the thin section images, plane polarized light images were used for the Teno sample, and cross-polarized light images were used for the CRB sample.Image swaths from separate thin sections from the same sample were aligned using the bands as a guide (Figures 2 and 3), allowing measurements of phenocrysts and vesicles to be given with positions relative to the dyke margin.
Phenocrysts and vesicles were manually outlined in ImageJ (Schindelin et al., 2012).The area, maximum diameter and best-fit ellipse of the outlined features were generated using the internal functions within ImageJ.The best-fit ellipse provided measurements of an object's major and minor axes, aspect ratio and orientation.Two sets of binary images were created using the manual outlines: one set containing phenocrysts, and one set containing vesicles (Figure 4).
To quantify the textural data with distance from the dyke margin, the binary images and raw data were analyzed in strips 1,080 pixels (1 mm) wide, perpendicular to the dyke margin.This bin width was chosen as it was wider than the largest vesicles, and captured enough phenocrysts and vesicles to calculate reliable means, while remaining narrow enough to show changes within individual marginal bands.Medians and interquartile ranges were calculated from the raw data within each strip, for phenocryst and vesicle area, maximum diameter, aspect ratio, long-axis orientation, and circularity.Aspect ratio was defined as major axis/minor axis, while circularity was defined as object perimeter over the circumference of a circle with equal area.For orientation data, circular statistical methods were used to calculate means and standard deviations (Jammalamadaka & SenGupta, 2001).
The binary images were analyzed in MATLAB, where values for phenocryst area fraction and vesicle area fraction were generated based on the proportion of white pixels within a given region of the binary image, and number densities were generated from the number of independent white areas.Vesicle area fraction was calculated as vesicle area/total area, while phenocryst area fraction was calculated as the phenocryst area/(total area-vesicle area), to exclude the vesicles.Phenocryst number densities were calculated as the phenocryst number/(total area-vesicle area), and vesicle number densities were calculated as vesicle number/(total areavesicle area).
For the CRB sample, groundmass and microlites were investigated using a Zeiss Evo scanning electron microscope with 15.00 kV beam at the GJ Russell Microscopy Facility, Durham University.Three rows of isolated images were captured at 400× magnification, with each row containing images at intervals of 0.5 mm, from the margin into the sixth apparent band.At such high magnification, it was not feasible to create continuous image swaths across the sample.Within each image, plagioclase microlites were outlined manually, and oxides were captured automatically via thresholding in ImageJ.
Within the Teno sample, phenocrysts are 0.02-0.45mm long (median 0.10 mm) within a dark groundmass, which is glassy within 5 mm of the margin, but becomes progressively more crystalline toward the dyke center.Within the CRB sample, there is a greater range in crystal size; phenocrysts are arbitrarily defined as crystals longer than 0.1 mm (up to a maximum of 4.28 mm), and anything shorter than this is considered a microlite.Additionally, there are two vesicle populations within the CRB sample: angular diktytaxitic voids, all <30 μm, which were only seen using the SEM; and large, rounded, primary magmatic vesicles.The diktytaxitic voids are discussed in more detail later, in Section 5.2.

Teno Sample
Six bands are present within the Teno thin sections (Figure 5).By eye, the bands appear to be defined by dark, glassy material on their marginal side, transitioning into lighter-colored material toward the dyke center, before a relatively sharp boundary with the next band (Figure 2).The bands and their boundaries generally become wider and less distinct toward the dyke center.Here, we focus on the first four bands, which have the most distinct edges.

Teno Phenocrysts
Binary image analysis reveals that the bands are defined by variations in the concentration of phenocrysts, which are dominantly plagioclase.Nearer the dyke margin, the bands contain fewer phenocrysts than the side nearer the center, as shown by the variation in number density, which ranges from 30 to 70 mm 2 .The onset of each band is defined by a sharp drop in phenocryst area fraction (Figure 6a) and number density (Figure 6b).
The size of phenocrysts increases gradually toward the dyke center, with mean lengths increasing from around 0.1-0.15mm across the 9 cm region studied, but there is no relation to the bands (Supporting Information S1).As such, the increase in phenocryst crystallinity within each band is predominantly controlled by the phenocryst number density, not size.There is also a general increase in the crystallinity of the groundmass, which is glassy at the margin, but becomes more crystalline toward the dyke center.
Phenocryst aspect ratio shows no significant variation across the sample and no relation to the bands (Supporting Information S1).The mean orientation of phenocrysts is constant across the bands at around 10°, measured from the dyke margin.The standard deviation shows a gradual increase from around 30°near the dyke margin to values exceeding 50°nearer the dyke center, indicating a reduction in the strength of the shape preferred orientation (SPO) (Supporting Information S1).

Teno Vesicles
Binary image analysis reveals that vesicularity is generally highest near the centers of bands compared with their edges (Figure 6a).This is most evident in the first band, where the peak value of 0.22 in the band interior contrasts with values below 0.05 at the band edges, but this pattern is less distinct in other bands.In all cases, the peak vesicularity does not fall exactly in the band center, but is slightly skewed toward the dyke center.Although vesicle size increases away from the dyke margin, with mean maximum diameters increasing from 0.05 to 0.32 mm across the studied section, vesicle size shows no trend within the bands themselves.The circularity of vesicles decreases toward the dyke center, with no relation to the bands (Supporting Information S1).The mean orientation of vesicles is 10°, the same as the phenocrysts, but the strength of their SPO is lower, with a circular standard deviation of 60°(Supporting Information S1).Vesicle orientation shows no relation to the bands.

CRB Sample
Eight bands were identified in the CRB thin sections (Figure 5).By eye, the bands appear to be defined by a variation in vesicularity, but unlike the Teno sample, there are no sharp, color-defined boundaries between bands (Figure 3).By analogy with the Teno sample, we assume that vesicularity peaks in the centers of bands, and define band edges as the regions with the lowest vesicularity.The bands become wider and less distinct toward the dyke center.

CRB Phenocrysts
Binary image analysis reveals no relationship between the bands and any of the measured phenocryst characteristics.Area fraction is shown in Figure 7, and other characteristics (number density, length, aspect ratio and orientation) are presented in the Supporting Information S1.The mean orientation of phenocrysts is 13 ± 41°from the dyke margin, in the vertical plane.

CRB Vesicles
Binary image analysis reveals a clear peak in vesicularity near the center of each band (Figure 7).Peak vesicularity is typically between 0.15 and 0.25, while the lows in vesicularity, taken as the boundaries between bands, have values of 0.01-0.10.Vesicle size increases toward the dyke center, with mean maximum diameters increasing from 0.02 to 0.3 mm, and vesicle size fluctuates across each band, with peaks in maximum size coinciding with peaks in vesicularity (Supporting Information S1).Minimum vesicle circularity coincides with maximum vesicle size, but mean circularity decreases toward the dyke center.The mean orientation of vesicles is constant across the sample, with an overall mean of 24 ± 55°from the dyke margin within the vertical plane (Supporting Information S1).Vesicle orientation shows no relation to the bands.
We note that there are two populations of vesicles present in the CRB sample.In the SEM images, it is apparent that the smallest vesicles (<30 μm) form a diktytaxitic texture (Figure 8).These small voids have angular shapes with straight edges, often completely filling the space between plagioclase microlites.In some cases, the diktytaxitic voids have curved edges on their contact with the groundmass, but they always have at least one angular edge.The size distribution from small, diktytaxitic voids to large, rounded vesicles is not continuous, so they are separate populations.The diktytaxitic textures are only seen in SEM images, and there is no evidence of them connecting with the large, rounded vesicles that were observed under the petrological microscope at lower magnification.The vesicle measurements derived from binary image analysis, presented above, only apply to the larger, primary magmatic vesicle population.

CRB Microlites
SEM analysis of the CRB sample reveals the presence of plagioclase and oxide microlites.Plagioclase microlites have a median length of around 20 μm and a modal aspect ratio of 5, which are both constant across the studied region (Figure 9).No aspect of the plagioclase microlites showed any relation to the bands.Oxide microlites show a gradual increase in maximum diameter from ∼6 to 10 μm across the studied region, and a relatively sharp increase in circularity from 0.3 to 0.5 between 20 and 30 mm from the dyke margin, associated with a change from skeletal, dendritic shapes toward more equant forms (Figure 8).However, no change in oxide size or shape showed any relation to the bands (Figure 9).

Conceptual Framework for a Band-Formation Model
We first use the primary field and textural observations to create a conceptual framework for marginal band formation, based on the overall similarities between the two samples, and other examples of marginal bands from This repetitive process is significant because dyke margins comprise the earliest material to solidify within the fracture at the dyke tip, meaning that their textures are intrinsically related to propagation.The cyclical textures of the marginal bands indicate a repetitive dyke-tip process, and their prevalence within dykes across diverse volcanic settings (Brouxel, 1991;Delcamp et al., 2012;Drever & Johnston, 1959;Galindo & Gudmundsson, 2012;Holness & Humphreys, 2003;Kavanagh et al., 2018;Platten, 2000;Thiele et al., 2021;Thordarson & Self, 1996;Walker, 1987;Walker & Eyre, 1995) suggests that this process is common.
We propose that the most likely cause of marginal band formation is repeated cooling, solidification and refracturing of magma within the dyke tip.Therefore, we interpret the bands to represent sequential pulses of magma through the narrow dyke tip region, where cooling and solidification is rapid.We propose that each pulse occurs as the dyke propagates a step further, causing repeated widening of the dyke further back from the tip (Figure 10).The similarity in the banding pattern wherever these textures are found along the dyke suggests that the timescales associated with each magma pulse are likely to be roughly constant.
The marginal bands accumulate as follows.Consider a fixed point of observation, initially positioned just ahead of the dyke tip.A new crack forms in the country rock, and a very thin sheet of magma is injected at our point of observation.This is the first pulse, and it cools rapidly to form band one.Propagation stalls, causing pressure to build behind the solidified tip until it rises high enough to fracture the tip and the host rock beyond, allowing the dyke to propagate another step.This is pulse two, and at our point of observation, this forms the second band.It differs from the first in two key ways: firstly, the dyke tip has advanced further away, and so the magma-filled crack is wider at our point of observation; secondly, the magma has been injected against the previous band, which is warmer than the host rock.As such, band two is wider than band one and lacks a glassy margin.The consistent time interval between pulses is sufficiently long for band two to solidify fully, despite being wider.This process repeats, and the bands build up.Pulse three cracks through band two, and pulse four cracks through band three.Each new pulse intrudes the center of the last, because this is where material is warmest and therefore weakest (Gudmundsson, 1984).As the tapered dyke tip propagates further, the magma-filled crack gets wider, but the interval between pulses is long enough for each band to fully solidify, despite them getting wider.However, as the tip moves further from the point of observation, pressure drops and fluctuations in flow rate become less pronounced, causing bands to become less texturally distinct.Eventually, magma pulses will be wide enough that they will not solidify to their centers before the next pulse arrives.From then on, no more bands will be produced at the point of observation.Magma flow may still be pulsatory, reflecting the rhythmic propagation of the nowdistant crack tip, but this will not be recorded.
Finally, it is interesting to note that, at each step of propagation, the magma pulse responsible will form a first band at the dyke tip, a second band immediately behind the dyke tip, a third band behind that, and so on.This means that each band could, in principle, eventually be traced to the outermost dyke margin, although the extremely high aspect ratio of bands combined with the often poor and patchy exposure of most dyke margins means that this would be very challenging to test in practice.
In the remainder of this section, we explore the feasibility of this model in the context of the data presented in Section 5.

Band Widths and Dyke Geometry
If we read each band as a pulse of magma that propagated the dyke a step further, then the increase in band width toward the dyke center (Figure 5) provides us with information about the geometry of the dyke tip.The first band represents the half-width of the first pulse, from the very tip of the dyke, while each subsequent band formed as the tip migrated a step further beyond the point of observation, demonstrating that the molten width of the dyke increased with distance from the tip.We can therefore infer that the dyke tip was tapered, which is concordant with field observations (Figure 11; Clemente et al., 2007;Geshi et al., 2010;Healy et al., 2018;Marinoni & Gudmundsson, 2000;Stephens et al., 2017).
The band widths demonstrate that the tip of the propagating dyke was very narrow.Dyke exposures are sometimes found to be less than a few centimeters wide (Jolly & Sanderson, 1995;Marinoni & Gudmundsson, 2000;Stephens et al., 2017), so it is feasible for a narrow tip region to exist.Indeed, the propagation of fluid-filled, pressurized cracks in the framework of LEFM requires a crack to have a tapering point in order to concentrate stress and instigate fracturing (Rivalta et al., 2015).Our conceptual model does not require dyke tips to have a specific shape (i.e., concave or convex), only that they narrow to a point.
In our framework, the width of each band reflects the tapering dyke tip geometry, not the duration of each magma pulse.The increase in band width toward the dyke center does not indicate longer intervals between magma pulses, as this would manifest as an increase in width of successive first bands in the direction of dyke propagation.Instead, the intervals between magma pulses are likely to have been broadly consistent.Each band represents a pulse of magma through the dyke tip, causing the dyke to propagate a step further, and to inflate further upstream.Note that once a sufficient number of bands have accumulated, a pulse will not fully solidify, causing it to mingle with the next pulse.At this point, the band formation mechanism is lost.

Solidification and Emplacement Timescales
In our conceptual model, each band represents the half-width of an active section of the dyke at the time of emplacement.We can estimate the time required to solidify each band by considering the rate of heat conduction into the adjacent host rock.We construct a simple numerical model in which magma is emplaced instantaneously into a narrow dyke, then cools conductively until the dyke center reaches a temperature that we choose.A fresh pulse of magma is then emplaced instantaneously along the center of the dyke.The temperature profile is advected laterally to accommodate the new pulse, and the system continues to cool conductively.The widths of each magma pulse are defined by the measured widths of our bands.Full details of the model, including relevant equations and numerical scheme are given in Appendix A.
We apply the model to the Teno sample, using the band width data shown in Figure 5.The CRB bands are of very similar widths, resulting in almost identical timescales, and these are provided in the Supporting Information S1.We estimate the time required for each band to solidify and cool to a defined final temperature, which we take to be the band's emplacement time.We generate a lower bound for emplacement times by ignoring the latent heat released during crystallisation, and an upper bound by including latent heat, which delays solidification.The initial magma temperature, T M , is taken to be 1,200°C, the same as the eruption temperature measured in 2021 La Palma eruption (Carracedo et al., 2022;Castro & Feisel, 2022), with a solidus temperature T S = 1,000°C, and an initial host rock temperature T 0 = 0°C.We assume a specific heat capacity of c p = 1200 J kg 1 K 1 (Mollo et al., 2011;Turcotte & Schubert, 2002), a constant density of ρ = 2,700 kg m 3 for the liquid and solidified magma (Lesher & Spera, 2015), a thermal conductivity of k = 1.3 W m 1 K 1 (Sehlke et al., 2020), and where latent heat is included, a latent heat of L = 4 × 10 5 J kg 1 (Bruce & Huppert, 1990;Turcotte & Schubert, 2002).Subsequent magma pulses are injected once the center of the last reaches a defined temperature T f , and we run the model for six different values of T f : 1,000, 900, 800, 700, 600 and 500°C.
We proposed in Section 6.2 that the increase in band widths toward the dyke center reflects the tapered geometry of the dyke tip, and that the duration of each magma pulse is likely to be consistent, albeit with variability derived from the stochastic nature of tip propagation.Clearly, greater widths of material take longer to solidify, so the cooling timescales of the widest, innermost bands will be the longest, and will reflect the duration of magma pulses most closely.However, estimating cooling timescales for the inner bands requires us to model the emplacement of the earlier, outer bands, as these act as insulation for all subsequent magma pulses.We do not intend the following solidification timescales for sequential bands to imply that magma pulse duration increased over time.Instead, we view the timescales of the innermost (later-forming) bands as an indicator of potential pulse durations.
We use central temperatures T f as a constraint on our model because we infer that each band fully solidified before the next pulse arrived.Our samples show sharp, cyclic variations in phenocryst concentration, which would not exist if magma pulses had mingled, or if there had been thermal erosion on the arrival of a new pulse.As such, we know that T f ≤ T S for each band.Using this temperature constraint, we explore the range of possible magma pulse durations.
We compare the solidification timescales for each band by tracking the position of the solidification front through time, defined as the isotherm within the dyke where T = T S (Figure 12).Bands solidify faster when the center of the previous band is allowed to cool to a lower temperature before the injection of the next pulse.However, if the next pulse is injected when the center of the last has cooled to only 900°C or 1,000°C, the previous band is melted back.This causes the solidification front to migrate back toward the margin before starting to move toward the dyke center.If pulses are injected with no interval between them (i.e., as soon as the dyke center reaches the solidus, where T f = T S = 1,000°C), up to 80% of the previous band can be thermally eroded (see Figure 12, L = 0 J kg 1 K 1 , band 5).
We define the "cooling interval" for a band as the time taken for the dyke center to reach T f .In all our model runs, T f ≤ T S , so that the cooling interval is always at least as long as the time taken for the band to fully solidify.The cooling intervals for the innermost bands give us a likely timescale for band emplacement.latent heat release is neglected; (right) latent heat release is included.T f is the final temperature at the dyke center before the next magma pulse is emplaced.When the solidification front moves back toward the margin (which is most obvious in the yellow curves), it is melting back previously solidified material.
To put a lower bound on cooling intervals, we must ascertain the maximum temperature T f to which the dyke center could have cooled prior to the emplacement of the next band.Our modeling suggests that T f > 800°C leads to thermal erosion of the previous band when the next pulse arrives.However, we see no textural evidence for thermal erosion between our bands.Therefore, we conclude that each band cooled to at least 800°C before the next band arrived.
Estimating the upper bound for the cooling interval is more challenging.One line of evidence is the absence of chilled margins between bands.Internal chilled margins are characterized by a drop in groundmass grain size, created where hot magma cools rapidly against cold surroundings (Brouxel, 1991;Huppert & Sparks, 1989;Platten, 2000).Within our samples, the outermost margin is glassy, where the initial magma pulse cooled rapidly against cold host rock, but none of the subsequent bands have particularly glassy outer edges.This implies that magma pulses never cooled close to the initial host rock temperature.To produce a glass, the magma must drop below its glass transition temperature, usually around 739°C for a basaltic melt, depending on cooling rate and concentration of H 2 O (Giordano et al., 2005).The first response of material in contact with the wall is to cool to a temperature halfway between the host rock and the magma; therefore, if the initial magma temperature is 1,200°C , the wall must be ≤278°C to induce a chilled margin.As such, we would only expect the bands to have individual glassy margins if the center of the previous band had cooled to 278°C.
We estimate the cooling interval required for the second and third marginal bands to develop a glassy margin by setting T f = 278°C.These bands are the most likely to develop an internal glassy margin, as the first band will have cooled rapidly, and not heated the host rock significantly.We find that for the Teno sample, the first band must cool for 27-45 min to produce a glassy margin on the second band, and the second band must cool for 34-57 min to produce a glassy margin on the third band.For the CRB sample, where the bands are slightly narrower, the first band must cool for 15-26 min, and the second band must cool for 30-50 min.To produce glassy margins on any bands further in, the cooling intervals would be even longer.The fact that we don't see any glassy margins, even on the second and third bands, shows that these bands were being emplaced faster than these estimates, and so we take this value as the upper bound on our emplacement timescales.These timescales are similar to the solidification times of the inner bands, and so we conclude that the intervals between magma pulses were on the scale of tens of minutes.

Vesicles and Pressure Drops
Within both samples, bands are associated with peaks in vesicularity.Since dykes solidify progressively inwards from their margins, each band can be interpreted as a time series, with the oldest material at the margin, getting progressively younger toward the center.We interpret the peaks and troughs in primary magmatic vesicularity as showing a variation through time, rather than representing the spatial organization of vesicles in the active system.Since vesicularity must depend strongly on pressure, we therefore interpret peaks in vesicularity as periods of low pressure, in which bubbles grew, and troughs in vesicularity as periods of high pressure, in which bubbles shrank.
We propose that the fluctuations in pressure recorded by the vesicles within the dyke margins result from the pressure drops that occur when a new pulse fractures the host rock and causes the dyke to propagate a step further (Figure 10).We draw on extensive previous work on bubble growth in magma (Prousevitch et al., 1993) and on the LEFM of dykes (Menand & Tait, 2002;Rivalta et al., 2015) to assess quantitatively the feasibility of our model.Specifically, we test whether the pressure drop required to produce the cyclic variation in vesicularity is consistent with the pressure drop associated with fracture propagation, and whether the time required for bubble growth is consistent with the time available between magma pulses.A physical framework for this analysis is presented in Appendix B.
According to LEFM, a crack can propagate only when its geometry and internal pressure combine to produce a stress intensity at its tip which exceeds a critical threshold known as the fracture toughness, K C .We can estimate the internal pressure for fracture propagation P F via Equation B4, reproduced below: Journal of Geophysical Research: Solid Earth

10.1029/2023JB028007
where ∆ρ is the density difference between magma and host rock, and g is gravitational acceleration.Values for fracture toughness are highly uncertain; the literature reports values from 0 to 10 4 MPa m 1/2 , with laboratory measurements (e.g., Atkinson, 1984;Atkinson & Meredith, 1987;Balme et al., 2004) several orders of magnitude smaller than estimates based on dyke geometry (e.g., Daniels et al., 2012;Delaney & Pollard, 1981;Olson, 2003).However, assuming that fracture toughness is in the region of 0-500 MPa m 1/2 , the same as most lab and field values, and that the density difference between magma and its surroundings is between 100 and 500 kg m 3 (assuming that magma and host have the same density of 2,700 kg m 3 , but magma has a vesicularity up to 0.2), all estimated overpressures are less than 10 MPa.This value is taken to represent the maximum pressure drop that could occur on the instigation of a new fracture in the tip of a temporarily stalled dyke, and this pressure drop would cause the growth of bubbles and the peaks in vesicularity.
Using the estimated value of P F allows us to calculate a characteristic growth timescale t R for bubbles within the magma, following the sudden decompression when the fracture propagates a step further.We calculate t R via Equation B8, reproduced below: This equation is valid when the pressure drop is much larger than the Laplace stress arising from surface tension (i.e., when P F ≫ 2σ/R).Vesicles in our samples typically have R ≳ 10 5 m, and surface tension of a basaltic melt is typically σ = 0.2 N m 1 (Colucci et al., 2016), hence the Laplace stress is of order 0.1 MPa, which is much smaller than the fracture pressure.Taking a typical viscosity value for basalt of η = 100 Pa s (Giordano et al., 2008), using the measured values for vesicle fraction ϕ (Section 5), and assuming P F is between 1 and 10 MPa, all growth timescales are less than 4 × 10 4 s.
We also calculate the diffusional growth timescale t D for the vesicles in our samples via Equation B9, reproduced below: where L W is the characteristic thickness of the walls between bubbles, and D H2O is the diffusivity of H 2 O in the melt.We find that the diffusional timescale for the vesicles in both samples is around 20 s (for detailed calculations, see Appendix B).This is much shorter than the emplacement timescales of the bands (Section 6.3), and so we conclude that bubble size is predominantly controlled by pressure.This makes the vesicularity within the bands a reliable indicator of pressure fluctuations within the dyke tip.
We relate vesicularity to pressure using the equation of Shishkina et al. (2010), which relates the solubility of H 2 O in basalt to pressure (Equation B10, reproduced below): where S H2O is the solubility in wt.%, and P is pressure in MPa.The saturation pressure can be found when S H2O equals the initial H 2 O concentration.For any pressures lower than this, it is possible to estimate the total volume of H 2 O exsolved per kilogram of magma, which can be expressed as a gas volume fraction.This provides a means to relate pressure to vesicularity (see Appendix B for details).Working the other way around, to get pressure from vesicularity, requires a numerical solution, for which we used the solvers in MATLAB.
For both the Teno and CRB samples, areas of high and low vesicularity were selected within the binary images to provide minimum and maximum vesicularities.For the Teno sample, the maximum vesicularity is 0.22, and the mean value for a vesicularity peak is 0.16, while the minimum vesicularity is 0.04.In the CRB sample, the maximum is 0.24, the mean peak value is 0.20, and the minimum is 0.02.For the Teno basalts, initial water concentrations are likely to be between 0.8 and 1.5 wt.%, based on basalts from La Palma and El Hierro, which are neighboring volcanoes in the early stages of shield evolution, comparable to the setting of the Teno samples (Taracsák et al., 2019;Weis et al., 2015).Saturation pressures are therefore likely to have been between 8.4 and Self, 1996), with saturation pressures of 25-41 MPa.
Assuming a basaltic melt density ρ m of 2,700 kg m 3 and a temperature of 1,200°C, we calculate the equilibrium vesicularity that would result from a given pressure P, below the saturation pressure P sat .For the Teno sample, the maximum pressure drop (from maximum to minimum vesicularity) is 1.5-5.9MPa, and the mean pressure drop (from mean to minimum vesicularity) is 1.0-1.5 MPa, with the ranges dependent on H 2 O concentration.For the CRB sample, the maximum pressure drop is 7.2-13.4MPa, and the mean pressure drop is 6.1-11.4MPa.These values are consistent with the ballpark figure for fracture pressure of 10 MPa.
A final point worth mentioning is that within the CRB sample, the smallest vesicles are diktytaxitic voids (Figure 8).This texture has previously been found within dykes and lava flows and is commonly interpreted to result from late-stage degassing driven by crystallisation (Fuller, 1939;Martin et al., 2006;Rogan et al., 1996;Walker, 1989).These voids are a separate population from the larger, rounded vesicles measured under the optical microscope, as the distribution in vesicle size is not continuous.There is no evidence of diktytaxitic voids larger than 30 μm, or of any coalescing into larger, rounded vesicles.As such, it seems likely that they are a late-stage feature related to crystallisation, while the main population of larger, rounded vesicles formed earlier, before the groundmass crystallized.

Phenocrysts and Flow Differentiation
The bands in the Teno sample are defined by a cyclic variation in the concentration of plagioclase phenocrysts.Bands start with a low concentration of phenocrysts on their marginal side, then increase in concentration toward the dyke center, before a relatively sharp drop in concentration at the beginning of the next band (Figure 6).There is an accompanying variation in vesicularity, which peaks toward the side of the band closer to the dyke center.In the literature, most occurrences of marginal bands fall into one of three categories: phenocryst-defined (Drever & Johnston, 1959); groundmass-defined (Brouxel, 1991;Holness & Humphreys, 2003); or vesicle-defined (Delcamp et al., 2012;Galindo & Gudmundsson, 2012;Kavanagh et al., 2018;Kile, 1993;Platten, 2000;Thiele et al., 2021;Thordarson & Self, 1996;Walker, 1987;Walker & Eyre, 1995).The Teno sample is therefore rare in that it displays cyclic variations in both phenocryst concentration and vesicularity, and this allows us to infer that they form via the same process.
Within the Teno sample, the length of plagioclase crystals shows only a slight increase toward the dyke center (0.1-0.15 mm), with no relation to the bands.Close to the margin, plagioclase is the dominant crystal species within a glassy groundmass.As such, we have classified the plagioclase crystals as phenocrysts, and they likely grew elsewhere before being emplaced at the dyke margin.The gradients in phenocryst concentration are therefore unrelated to cooling rates, and we require another process to explain their presence.
One explanation for the gradients in phenocryst concentration, which we discount, is that there was a cyclic variation in the phenocryst load carried by the magma.The dyke solidifies inwards from its margins, and so a change in the concentration of phenocrysts within the magma flowing through the dyke would lead to a change in the concentration being captured by the solidification front.However, this would require a systematic variation in the phenocryst load over time, for which we see no rationale.An alternative explanation is that each band arrived with the phenocryst gradients already in place, and that these gradients were preserved when the magma solidified.This is in better agreement with the sharp drops in phenocryst concentration at the onset of each band.
We propose that the variation in phenocryst concentration is best explained by flow differentiation, the process by which suspended particles migrate away from steep velocity gradients at the edges of a flow.This phenomenon has previously been observed in analogue laboratory experiments and has been used to explain increased phenocryst concentrations in dyke centers (Barrière, 1976;Bhattacharji, 1967;Bhattacharji & Smith, 1964;Gibb, 1968;Komar, 1972).
Several mechanisms are known to cause lateral particle migration in flows, including wall effects (Maude & Whitmore, 1956) and the Magnus effect (Komar, 1972), but the most widely accepted mechanism for flow differentiation within dykes is caused by particle interactions (Bagnold, 1954), sometimes referred to as the Bagnold effect.When suspended particles collide, they temporarily gain momentum perpendicular to the flow direction.This momentum is higher for particles colliding with a greater velocity difference, prompting a net migration away from steep velocity gradients at the edges of flows, and causing particles to concentrate toward the flow center (Bagnold, 1954;Petford & Koenders, 1998;Phillips et al., 1992).Laboratory experiments have shown that the rate of lateral migration increases for higher shear rates and higher particle concentrations, both of which increase the frequency of collisions (Abbott et al., 1991;Karnis et al., 1963;Leighton & Acrivos, 1987).
To determine whether flow differentiation could have occurred within the bands of the Teno sample, we calculate a development length L dev using the method of Lecampion and Garagash (2014).This is the characteristic lengthscale traveled by flowing magma over which particles migrate into a steady-state concentration distribution across the channel (the actively flowing part of the dyke).The model is presented in Appendix C.
The model for L dev requires three input parameters: the initial particle fraction ϕ 0 ; the channel half-width H; and the particle radius a.Initial particle fraction is estimated by taking the total area of phenocrysts within a band, and dividing it by the area of the band excluding vesicles.We determine a mean value of ϕ 0 = 0.07, with a standard deviation of 0.01, showing that the initial particle fraction is consistent between bands.The channel half-width is taken to be the width of each band (so H ranges from 4.42 to 11.70 mm) and the particle radius is approximated as a = 0.05 mm, which is half the median phenocryst length.
Taking a mean value for the width and initial particle concentration of the first four Teno bands, we find that L dev = 1.4 km.This value is less than the length of most dykes, demonstrating that flow differentiation can be expected to occur within the narrow tips of dykes.It should also be noted that flow differentiation is an asymptotic process, and most differentiation occurs before L dev is reached.
The development length depends strongly on channel width (L dev ∝ H 3 , Equation C1), such that flows in narrower channels differentiate over shorter distances.This could partly explain why the earliest bands closest to the margin, which formed when the actively flowing width of the dyke was narrowest, are more distinct, with fewer phenocrysts remaining close to their marginal edges.
We infer that flow differentiation is a viable mechanism to produce the systematic increase in phenocryst concentration within each band in the Teno sample.The flow differentiation model presented in Appendix C does not account for the presence of vesicles, which likely affect particle migration.For the Teno sample, there is evidence that vesicles formed after flow differentiation had occurred (Section 6.6).Nonetheless, we caution that the presence of vesicles may be a confounding factor in the analysis above.For the CRB sample, the vesicle fraction greatly exceeds the phenocryst fraction (0.2 and 0.05 respectively), and the order of formation of the two phases is not as clear.Consequently, we do not attempt to analyze flow differentiation for the CRB sample, but we note the lack of textural evidence for it having occurred (Figure 7).
The SPO of phenocrysts within the Teno sample shows no variation toward the dyke center, or within the bands themselves, suggesting that the flow direction for each pulse remained constant (Philpotts & Asher, 1994;Wada, 1992).The degree of SPO is highest closer to the margins, echoing the results of previous studies (Coward, 1980;Das & Mallik, 2020), and the reduction in the degree of preferred orientation toward the dyke center could be due to an increasing groundmass crystallinity within each subsequent pulse, increasing the prevalence of crystal collisions.However, the overall consistency in SPO between the bands suggests that pulses were not substantially different from one another in terms of their fluid dynamics.In the CRB sample, plagioclase phenocrysts also display an SPO which remains constant toward the dyke center, implying the same direction of magma flow throughout the emplacement of the bands.

Relationship Between Vesicles and Phenocrysts
In the Teno sample, phenocrysts and vesicles have almost identical preferred orientations.Many vesicles are flattened against neighboring phenocrysts, or wrap around the end of a phenocryst (Figure 13).These relationships indicate that the vesicles grew within a pre-existing network of aligned phenocrysts, causing them to elongate in such a way that they attained the same direction of preferred orientation.This implies that, within the Teno sample, the vesicles experienced growth after the phenocrysts had been emplaced, and after the phenocryst textures had been established through flow differentiation.In the CRB sample, phenocrysts and vesicles also have similar preferred orientations, but we do not see textural evidence that allows us to infer whether the vesicles grew constrained by a phenocryst network.
The relationship between phenocryst and vesicle textures provides evidence for the order of events associated with our hypothesis of pulsatory magma flow.We infer that the gradients in phenocryst concentration were produced by flow differentiation as the magma pulse was emplaced, and that these concentrations were fixed in place when the pulse cooled and stalled.The vesicles then grew within this existing phenocryst network at the moment that a new fracture formed, when the dyke tip experienced a temporary pressure drop before continuing to advance.

Band Formation Model
Following the analysis presented in Sections 6.2-6.6,we can refine the model for marginal band formation presented in Section 6.1 as follows: 1.A magma pulse with a tapering geometry propagates through the host rock 2. The narrow tip solidifies rapidly, stalling the progress of the propagating crack 3. Pressure builds behind the blockage until it ruptures 4. A fresh pulse of magma causes the dyke to propagate a step further This process is cyclic, leading to the progressive build-up of bands along the dyke (Figure 14).The pulsatory flow of magma can create gradients in phenocryst concentration through flow differentiation, while the fluctuation in pressures associated with pulsatory flow and fracturing leads to the variation in vesicularity.The lack of evidence for melt-back or internal chilled margins between bands suggests that the interval between magma pulses is tens of minutes.Bands become wider toward the dyke center because they formed progressively further from the dyke tip, meaning there was a greater width of material available for solidification.The increasing distance to the dyke tip also causes bands to become less distinct toward the dyke center, because the magma becomes progressively insulated by previous bands, allowing melt-back and mixing between pulses, and causing flow and pressure fluctuations to become less severe as the tip propagates further away.Eventually, the dyke tip will be far enough away, and the dyke will be so well-insulated by previous bands, that flow past that point will become steady, and the mechanism for band formation will be lost.As such, we only see these bands at dyke margins.

Numerical Propagation Models
Dyke propagation models generally fall into two categories: those in which propagation is limited by fracture toughness, and those in which propagation is limited by viscous resistance to flow (Rivalta et al., 2015).In either case, propagation is initiated when the stress intensity at the dyke tip exceeds the fracture toughness of the host material.These models operate in the framework of LEFM, requiring dykes to have a tapering tip to instigate fracturing (e.g., Lister & Kerr, 1991).Once propagation commences, it is continuous: the dyke is arrested only on entering a material with a higher fracture toughness, or by viscous resistance, depending on the model type.
Solidification is commonly neglected (e.g., Kavanagh et al., 2006;Maccaferri et al., 2010), and so these models overlook the impact of solidification on dyke propagation.
The omission of solidification is usually justified by comparing propagation rates, indicated by seismicity to be around 1-2 km per hour (e.g., Einarsson & Brandsdóttir, 1980;White et al., 2011), with cooling rates, calculated to be a few days for a molten slab a few meters in width (Rivalta et al., 2015).However, our analysis suggests that using the width of solidified dykes to estimate cooling times is ill-posed; if a dyke has grown by inflation from several pulses, the molten width would always have been less than the final width.If the molten region is narrow, then solidification timescales become comparable to propagation timescales.Furthermore, most propagation models require the dyke tip to taper to a narrow point, which would make it prone to solidification.Even if the main dyke body remains molten, once the tip solidifies, it cannot transmit the fluid pressures required for fracture.This is likely to have a significant impact on propagation, with the potential to stall the dyke's advance at least temporarily.As such, models might consider whether resistance to propagation depends not only on fracture toughness or viscous losses, but also on the efficiency with which the cooling and solidifying dyke tip can exert pressure on its surroundings.

Experimental and Field Analogs
Analogue models used to explore dyke propagation typically use isothermal fluids (e.g., Kavanagh et al., 2017;Menand & Tait, 2002).However, pulsatory propagation has been observed by injecting paraffin wax into gelatin, as the wax cools and solidifies, leading to stalling and refracturing (Taisne & Tait, 2011).Stalling a interpreted to be caused by solidification reducing the stress intensity at the edge of the wax, meaning it could no longer fracture the gelatin.The paraffin wax intruded through a cycle of stalling, inflating, then rupturing and propagating, the same process we have inferred from the marginal bands in dykes.
Similar pulsatory behavior has also been seen in syrup flowing through chilled tubes, where the balance between cooling, viscosity and local pressures produces unstable flow even for a fixed driving pressure (Whitehead & Helfrich, 1991).This shows that pulsatory flow could occur within the dyke tip even in the absence of complete solidification, due to the effect of cooling on the magma's viscosity.
Finally, pulsatory propagation within dykes could be considered analogous to the emplacement of pahoehoe lava lobes, which propagate by cyclical rupturing, cooling and inflation.Although lava flows lose heat rapidly via radiation, and dykes lose heat slowly via conduction, the same pulsatory processes are likely to be ongoing, albeit on a different timescale.

Dyke Propagation Inferred From Seismicity
The concept of episodic or pulsatory dyke propagation is not a new idea.One of the few ways to monitor active dyke propagation is through the migration of seismicity, commonly interpreted to originate around the leading edge of the dyke (Einarsson & Brandsdóttir, 1980;White et al., 2011;Woods et al., 2019).However, dyke propagation mapped through seismicity is often seen to happen in bursts at timescales of hours to days, showing that magma transport is not a continuous process, and can occur in pulses.Common interpretations of such behavior include changes in host rock properties, interaction with faults, changing topography on the surface, and variation in pressure at the magma source (Rivalta et al., 2015).For example, two temporary halts in the vertical migration of earthquakes below Piton de la Fournaise were attributed to lithological discontinuities (Battaglia et al., 2005), while bursts in seismicity during the lateral intrusion of a dyke in Iceland were related to caldera subsidence at the magma source (Sigmundsson et al., 2015).However, these events are on a much larger scale than evidenced in the marginal bands, where the timescale for magma pulses and pressure fluctuations is measured in tens of minutes.Furthermore, these larger-scale processes are unique to each dyke, while marginal bands exhibit the same pattern in every setting and must therefore share a ubiquitous formation process.
Another example of large-scale pulsatory propagation was recorded by the migration of earthquake foci during the Bárðarbunga-Holuhraun rifting event, where rapid lateral growth was interspersed by periods of stalling, lasting hours to days.During stalled periods, earthquake foci migrated backwards from the dyke tip before propagation resumed, interpreted as a result of inflation (Woods et al., 2019).Such behavior was proposed to have arisen due to changes in topography, host rock properties, or source pressure, then exacerbated by solidification around the dyke margins (Woods et al., 2019).Bursts of seismicity related to dyke propagation have previously been interpreted to be a result of solidification and re-fracturing (Hayashi & Morita, 2003;White et al., 2011); however, this occurs on longer timescales than the marginal bands, and must be at a larger scale in order to produce detectable waves of earthquakes.
Our model for pulsatory dyke tip propagation is consistent with the seismic signals from dyke emplacement, in that it indicates the discrete nature of fracturing and dyke advance.Even "smooth" dyke propagation is characterized by a series of discrete earthquakes, each one of which must originate in a distinct fracturing event.The cyclic nature of our propagation model is likely disguised by the fact that the leading edge of the dyke will be segmented, split into regions that propagate in competition with different advance rates.As such, small, systematic, and repetitive bursts of propagation will be hidden by the signals from neighboring sectors, and by largerscale pulses transiting the system.

Relation to Dyke Arrest or Eruption Onset
The marginal bands demonstrate that solidification can influence magma flow in the narrow dyke tip.It is therefore feasible that, under the correct conditions (e.g., low magma flux into the tip, highly resistant host rock), solidification could stall propagation for longer, or even entirely.The relation between magma temperature, viscosity and advection rates can lead to positive feedback, meaning that small-scale perturbations in the dyke tip have the potential to cascade into larger events, influencing the behavior of the main dyke body.It is therefore worth considering whether large-scale layering within dykes may not arise from discrete magma pulses from the source, but from pulses generated within the dyke itself due to the influence of solidification.It is interesting to imagine the largest scale on which this behavior might occur; for example, whether centimetre-to meter-scale layering marked by internal chilled margins (e.g., Platten, 2000) or differences in columnar jointing structure (e.g., Das & Mallik, 2020;Gudmundsson, 1984;Ray et al., 2007) could result from large scale, internally generated, solidification-driven magma pulses.
Finally, it is worth noting that although the physical processes discussed here are ubiquitous to propagating batches of solidifying fluids, pulsatory propagation and band formation should not be expected in all dykes.Bands defined by varying vesicularity are more likely in magma with a high volatile content, or at shallower depth, making it more sensitive to pressure drops (Section 6.4).Similarly, variations in phenocryst concentration will only form when flow differentiation is significant, which requires a high volume fraction of phenocrysts with an appropriate size relative to the width of the dyke (Section 6.5).We also predict that, in the shallow subsurface, bands may be removed by thermal erosion if a dyke becomes an established flow pathway feeding eruptions at the surface.

Conclusion
Marginal bands are a common feature of basaltic dykes in volcanic settings around the world.They can be defined by a variation in phenocrysts and/or vesicles, and show a general trend of getting wider and less distinct toward the intrusion center.Variations in phenocryst concentration are likely to be caused by flow differentiation, while variations in vesicularity can be caused by pressure fluctuations as the dyke progressively solidifies.Both processes can be explained by pulsatory magma flow.
The margins are the earliest material to solidify, originating in the tip of the dyke, and they therefore hold information on the propagation process.The presence of banding at the margins suggests that pulsatory magma flow is associated with propagation, which is likely to arise from repeated cooling, solidification and re-fracturing of magma in the narrow dyke tip.The millimeter to centimetre widths of the bands implies rapid cooling in tens of minutes, suggesting that the pulses are local and short-lived.Material closer to the dyke center solidified at a greater distance from the tip, with flow becoming steadier over time, causing bands to become progressively wider and less distinct.
Our analysis does not put any particular constraints on the conditions under which pulsatory propagation due to solidification within the dyke tip can occur, and so we conclude that it is generally applicable, which is consistent with the fact that it can be found in diverse volcanic settings across the world.By contrast, large-scale layering, such as in composite dykes, can be attributed to other processes which will be unique to each dyke, such as fluctuations in magma flow rate from the source.However, it is possible that cooling, blocking and unblocking within the dyke can lead to fluctuations in flow rate, causing larger-scale textural layering.More importantly, the repeated solidifying and re-fracturing of the narrow tip region may have implications for dyke propagation models, especially in the way that magmatic pressures are transmitted to the host rock.

Appendix A: Solidification Model
We estimate solidification timescales based on the widths of the marginal bands, assuming them to be the halfwidths of magma pulses through the dyke tip.Following a magma pulse event, heat transport within the dyke is assumed to be purely due to heat conduction normal to the plane of the dyke.We assume that the magma and host rock have the same thermal conductivity, which is reasonable for a basalt dyke intruding layered basaltic lava flows.We then have the option of whether to include or omit the release of latent heat on solidification.Omitting latent heat is justifiable for the first band, where the groundmass is glassy; however, the groundmass becomes more crystalline toward the dyke center, suggesting that latent heat was released as the later bands solidified.We can create a cooling model that accounts for latent heat by considering the enthalpy of the system.
The model boundary condition is that the host rock is held at a fixed temperature T 0 at a distance equal to 10 times the final dyke half-width away from the dyke center, far from the warming effect induced by the dyke on the host rock (we find that the induced warming only extends ∼0.5 m into the host rock).The geometry of the model is shown in Figure A1.
Assuming the magma to be incompressible, the following energy conservation equation applies: subjected to the initial condition and boundary conditions: where ρ is mass density (kg m 3 ), h is enthalpy per unit mass (J kg 1 ), t is time after a magma pulse event (s), k is the thermal conductivity (W m 1 K 1 ), T is temperature within the magma (K), x is distance from the center of the dyke normal to the dyke plane (m), T M is the temperature of the magma pulse (K), W is the half-width of the dyke (m) and T 0 is the temperature of the host rock (K).
Assuming magma and host rock to have the same specific heat capacity, the relationship between temperature and enthalpy takes the form In response to the pressure drop instigated by fracturing, bubbles will also grow by diffusive degassing of H 2 O from the melt.The timescale for diffusive bubble growth can be approximated as where L W is the characteristic thickness of the walls between bubbles, and D H2O is the diffusivity of H 2 O in the melt (Note that Equation B9 is presented as Equation 3 in the main text).For the Teno basalts, initial water concentrations are likely to have been between 0.8 and 1.5 wt.%, based on basalts from La Palma and El Hierro, which are neighboring volcanoes in the early stages of shield evolution, comparable to the setting of the Teno samples (Taracsák et al., 2019;Weis et al., 2015).For the CRB samples, water contents are likely to have been higher, at 1.5-2.0wt.% (Thordarson & Self, 1996).Using these water contents, we estimate D H2O to be 3.6-5.3 10 m 2 s 1 for the Teno sample, and 5.3-6.7 10 m 2 s 1 for the CRB sample, based on the equations of Zhang et al. (2007).Based on average vesicle area, and the average groundmass area per vesicle, we find L W = 0.09 mm in the Teno sample, and L W = 0.11 mm in the CRB sample.The timescale for diffusive bubble growth is therefore around 20 s for both samples.This is much shorter than the emplacement times of the bands (Section 6.3), and so vesicle growth is predominantly controlled by the pressure drops, making vesicularity within the bands a reliable indicator of pressures changing through time.
The variation in vesicularity in the bands can be used to estimate the magnitude of the pressure drops occurring within the dyke tip.We relate vesicularity to pressure using the equation of Shishkina et al. (2010), which relates the solubility of H 2 O in basalt to pressure, so that S H 2 O = 0.2351 P 0.5758 .(B10) where S H2O is the solubility in wt.%, and P is pressure in MPa (Note that Equation B10 is presented as Equation 4in the main text).The saturation pressure can be found when S H 2 O equals the initial water concentration.For any pressures lower than this, it is possible to estimate the total volume of water exsolved per kilogram of magma, thereby providing a means to relate pressure to vesicularity.We can therefore estimate the pressures associated with the peaks and the troughs in vesicularity within our samples.The difference between these pressures gives us the pressure drop associated with each step of pulsatory propagation.
We can calculate the equilibrium vesicularity that would result from a given pressure P, below the saturation pressure P sat .The difference in wt.% between the initial H 2 O content S 0 and the solubility S P at a given pressure can be used to calculate the mass of H 2 O exsolved per cubic meter of magma m H2O , where m H2O = ρ m (S 0 S P )/100.(B11) This volume of this exsolved mass of H 2 O can be calculated from the ideal gas law, PV = nRT, where P is pressure in Pa, V is volume in m 3 , n is the number of moles of H 2 O, R is the universal gas constant 8.3145 J mol 1 K 1 , and T is 1473 K. Using the volume of H 2 O exsolved per cubic meter of magma, we can calculate the vesicularity ϕ as To get pressure from a given vesicularity, the steps outlined in Equation B10-B12 must be followed in reverse, which requires the use of a numerical solver (e.g., in Excel or MATLAB).
across the channel (i.e., the actively flowing part of the dyke).The development length depends on the initial, uniform particle fraction within the channel ϕ 0 , the particle radius a, and the channel half-width H, so where I is the dimensionless viscous number, S is the dimensionless inelastic storage coefficient, and κ is the dimensionless permeability.The viscous number depends on the maximum particle volume fraction, ϕ m , so that The inelastic storage coefficient S is controlled by the dimensionless friction coefficient μ, which can be calculated via where μ 1 is the value of μ at ϕ m , and β is a dimensionless compressibility coefficient.The value of S can then be calculated using By differentiating Equation C3 with respect to ϕ, we find dμ/dϕ, the reciprocal of which is which we can then substitute into Equation C4 to find S.The dimensionless permeability is calculated as where α is a constant.The constants used by Lecampion and Garagash (2014) have been experimentally verified for suspensions of hard spheres, but reliable empirical relationships do not exist for tabular shapes.Therefore, we must proceed by treating the phenocrysts as spherical; the constants for spheres are: ϕ m = 0.585, μ 1 = 0.3, β = 0.158 and α = 5.1.
The development length is then Taking a mean value for the width and initial particle concentration of the first four Teno bands, we find that L dev = 1.4 km.However, the development length is highly sensitive to the channel width, and so band 2, at only 4.44 mm width, has a development length of only 70 m, while band 3, at 18.06 mm, has a development length of 3.9 km.These are still feasible values for flow differentiation to occur, in that they are shorter than the lengths of most dykes.It should also be noted that flow differentiation is an asymptotic process, and most differentiation occurs before L dev is reached.

Figure 1 .
Figure 1.Examples of marginal banding from dykes in the Teno Massif, NW Tenerife.Bands can be found as weathered ridges (a, b), or defined by variation in color and vesicularity (c, d).

Figure 2 .
Figure 2. (left) Origin of the Teno sample in the margins of a 0.5-m-wide dyke; (right) three aligned thin sections taken from the sample, showing the first six bands.Heavily altered streaks through the sample are highlighted.The margin indicated is the true margin of the dyke, with no host rock attached.

Figure 3 .
Figure 3.The CRB sample, from the margins of a feeder dyke to the Roza lava flow, with four aligned thin sections taken from the sample in the vertical plane, perpendicular to the margin.Dashed yellow lines on the thin sections show the positions of minimum vesicularity, used to define the bands.The margin indicated is the true margin of the dyke, with no host rock attached.

Figure 4 .
Figure 4. Examples of binary images created from Teno (top left, plane-polarized light) and CRB (top right, cross-polarized light) microscope images.Binary phenocryst images (middle row) and binary vesicle images (bottom row) are created from manual outlines.Angle measurements from the vertical are demonstrated on the CRB phenocrysts.Margin is to the left in both images.

Figure 5 .
Figure 5. Widths of bands in the Teno and CRB samples.Teno bands were measured based on color variations by eye, later found to be caused by changes in phenocryst concentration.Teno band 6 overran the end of the thin section, making the value plotted a minimum estimate.CRB bands are based on distances between troughs in vesicularity.Absolute measurement error is ±1 mm for each sample, but the edges of bands are subjective.

Figure 6 .
Figure 6.(a) Variation in vesicularity (vesicle area fraction) and phenocryst area fraction across the Teno sample; (b) number densities of vesicles and phenocrysts across the Teno sample.Dashed lines show the positions of the bands determined by eye, based on color variations.Band edges coincide with drops in the phenocryst area fraction (highlighted with red arrows in panel (a)) and number density.

Figure 7 .
Figure 7. Variation in vesicularity (vesicle area fraction) and phenocryst area fraction across the CRB sample.Gray dashed lines have been placed at troughs in vesicularity.Bands are defined purely by variation in vesicularity, and show no relation to the phenocrysts.

Figure 8 .
Figure 8. SEM image of the CRB sample.Examples of plagioclase microlites are outlined in yellow, diktytaxitic voids in red, and oxide microlites have been circled in cyan.The large plagioclase crystal on the left side of the image, outlined in green, is 0.15 mm in length, so is classed as a phenocryst.

Figure 9 .
Figure 9. Variation in plagioclase microlite length and aspect ratio, and oxide circularity and maximum diameter across the first six bands of the CRB sample.The solid line shows the median value, and the dashed lines show the interquartile range.

Figure 10 .
Figure10.Schematic representation of our conceptual framework: (a) the cooling-driven cycle that leads to pulsatory propagation of the dyke tip; (b) gradual accumulation of bands.Each band represents a pulse of magma through the dyke tip, causing the dyke to propagate a step further, and to inflate further upstream.Note that once a sufficient number of bands have accumulated, a pulse will not fully solidify, causing it to mingle with the next pulse.At this point, the band formation mechanism is lost.

Figure 11 .
Figure11.The tapering tip of a dyke in the Teno Massif, Tenerife, which narrows to a width of <2 cm.Such narrow dykes were relatively common in the area, and are assumed to come from the leading edge of a magma pulse, stalled through rapid solidification.

Figure 12 .
Figure12.The position of the solidification front (the T = T S isotherm), moving inwards from the margin, based on the widths of the bands in the Teno sample: (left) latent heat release is neglected; (right) latent heat release is included.T f is the final temperature at the dyke center before the next magma pulse is emplaced.When the solidification front moves back toward the margin (which is most obvious in the yellow curves), it is melting back previously solidified material.

Figure 13 .
Figure 13.Examples of vesicles within the Teno and CRB samples.Teno: (a) vesicle showing evidence of coalescence; (b, c) vesicle shape and orientation determined by the orientation of surrounding phenocrysts; (d) vesicle enveloping ends of phenocrysts.CRB: (e) highly irregular vesicle shapes in the first band of the CRB sample caused by coalescence; (f) vesicleshapes influenced by neighboring phenocrysts (note that these are rounded, magmatic vesicles, as opposed to the small, angular, diktytaxitic voids that were found using the SEM, which are a separate population (Section 4).

Figure 14 .
Figure 14.Model for band formation.Each band is a pulse of magma which caused the dyke to propagate another step.Each pulse experienced flow differentiation before cooling and solidifying, capturing the phenocryst gradient.Pressure builds until the dyke tip ruptures, leading to a pressure drop and vesiculation.Propagation resumes with a fresh pulse of magma through the dyke tip.

(
S ≤ h ≤ c p T S + L h L c p h > c p T S + L Note that Equation B8is presented as Equation 2 in the main text).