Dyke apertures record stress accumulation during sustained volcanism

The feedback between dyke and sill intrusions and the evolution of stresses within volcanic systems is poorly understood, despite its importance for magma transport and volcano instability. Long-lived ocean island volcanoes are crosscut by thousands of dykes, which must be accommodated through a combination of flank slip and visco-elastic deformation. Flank slip is dominant in some volcanoes (e.g., Kilauea), but how intrusions are accommodated in other volcanic systems remains unknown. Here we apply digital mapping techniques to collect > 400,000 orientation and aperture measurements from 519 sheet intrusions within Volcán Taburiente (La Palma, Canary Islands, Spain) and investigate their emplacement and accommodation. We show that vertically ascending dykes were deflected to propagate laterally as they approached the surface of the volcano, forming a radial dyke swarm, and propose a visco-elastic model for their accommodation. Our model reproduces the measured dyke-aperture distribution and predicts that stress accumulates within densely intruded regions of the volcano, blocking subsequent dykes and causing eruptive activity to migrate. These results have significant implications for the organisation of magma transport within volcanic edifices, and the evolution and stability of long-lived volcanic systems.

. Stress accumulation assuming a Maxwell viscoelastic rheology for a volcanic edifice with Young's modulus of 2 GPa, Poisson's ratio of 0.25 and viscosity (µ) of 1-5×10 22 Pas. The constant rate of dyke-injection ( F d ), dyke height h and initial excess pressure P 0 are determined from observed maximum and mode dyke apertures (which relate to the initial and equilibrium excess pressure) and estimates of the bulk tangential strain using a least-squares solver. The results for viscosities of 3-4×10 22 Pas match observed dyke aperture distributions and reach an equilibrium state after 400 ka, which corresponds to the amount of time it took for activity to localise onto the southern flank of Volcan Taburiente.

Photogrammetric survey methodology
Seventeen UAV surveys were conducted in May 2017 with approval and assistance from Parque National Caldera de Taburiente, using a DJI Phantom 4 Pro and its integrated 20-megapixel camera. Survey methodology was tailored to the topography being captured: cliff faces were flown manually using horizontal flight lines and horizontal + 30 degree downward oriented viewing angles, while MapPilot was used to survey flatter areas, using a nadir viewing angle and overlaps of 80% along flight lines and 70% between them. The target ground sampling distance (GSD) was a compromise between desired survey area and the features of interest, and so varies between surveys (Table 1). Generally, large cliff-exposures (e.g. Risco Liso) were flown at a distance of >50 m to enable a large survey area (at lower GSD), while smaller areas were flown at distances of 10-50 m to resolve smaller features.
Poor GPS reception and limited or impossible access to survey areas precluded the employment of accurately surveyed ground control points in most cases (Table 1). Instead, the surveys were approximately georeferenced using on-board GPS data and then aligned to the publically available 2014 Spanish LiDAR survey series (which has a resolution of 0.5 to 0.75 points/m2) after photogrammetric reconstruction. This alignment was done using the iterative closest point (ICP) algorithm in Cloud-Compare. All of the surveys were co-registered to this LiDAR data to within 1-2 m, suggesting no systematic errors or distortions in the final digital outcrop models. Structure-from-motion multi-view-stereo (SfM-MVS) photogrammetric reconstruction was performed using Agisoft Photoscan Professional (now called Agisoft Metashape) version 1.4.3. Image contrast and colour balance was enhanced using a batch-script in Adobe Photoshop (such that the same operation was applied to all images in a survey) prior to reconstruction. If radial dykes are emplaced into a volcanic edifice that deforms visco-elastically, the tangential (hoop) stress that accumulates during repeated dyke emplacement can be investigated using a 1D Maxwell visco-elastic model. Here we detail the assumptions and derivation of the model beginning with a generalized three-dimensional viscoelastic model.

Assumptions and boundary conditions
First, we assume that: (1) Ascending dykes are deflected towards areas of low stress, such that long-term averaged circumferential stress σ θ is spatially uniform (at least to some critical distance). This assumption means that the intrusions induce a change in the circumferential normal strain θ that is spatially constant.
(2) Dyke intrusion does not change vertical stress (due to the free surface). Adopting a notation where σ and are taken as changes in stress and strain, respectively, due to the dyke intrusion, this assumption is stated as σ z = 0.
(3) The edifice is free to move in the r direction (due to its conical shape), meaning σ r = 0.
(4) The r, θ, z directions define the principal axes for the induced stresses and strains, that is, there are no induced shear stresses or strains in this coordinate system.
From these assumption, we can express the change in the stress tensor as: Using Hooke's law, a purely elastic response to this stress would be: where E and ν are the Young's modulus and Poisson's ratio, respectively. It is often convenient in constitutive modeling to separate stress and strain into their volumetric and deviatoric parts. So, by decomposition we can write that: where δ ij is the Kronecker delta and s ij and e ij are the deviatoric components of the stress and strain tensors.

3D visco-elasticity:
A general form of a 3D viscoelastic behaviour is then given by (e.g. [1]): where the { } notation indicates differential operators of the form Expanding this for first-order (Maxwell) viscoelasticity (m = 1), and accounting for the uniaxial stress and zero shear stresses/strains implied by the assumptions outlined above gives, and considering the circumferential normal stress strain relationship (i, j = 2, 2) from the deviatoric part of the constitutive law, leads to Rearranging Eq. 8. givesė Next we will prescribe the parameters of the constitutive law in order to recover Hooke's law under rapid deformation and viscous relaxation under slow deformation. The material does not relax due to instantaneous strain, which causes a purely elastic response, so from Hooke's law (Eq. 4): Here G is the shear modulus. Similarly, if the strain rate is very slow then strain will be entirely viscous and it will not induce a change in stress, soė When viewed in this way it is evident that e θ is an initial value before the strain starts slowyly changing. Hence, as we are interestied in solving for stress-changes relative to an initial value before the intrusions form, we can use e θ = 0. Similarly, we can use µ to describe the viscosity of the rock (under uniaxial strain), such that: Substituting into the generalized constitutive law leads to a constitutive equation for the material under uniaxial circumferential stress given byė Hence, if∆ is the strain rate induced by the emplacement of the dykes, theṅ Thus we can relate the stress induced by the intrusions to the rate at which they induce strain using the 1D Maxwell-type viscoelastic equation∆ Relationship to overpressure In our instance, the induced strain rate∆ will depend on two factors: the rate of dyke injections F d and the aperture a of those dykes. Assuming the dykes are propagating laterally (as suggested by field observations) and have longer strike-length than height h, their aperture can be linked to the stress when they form: If we assume a constant rate of dyke injections F d along a circumference with length L = 2πr, strain rate becomes: Hence, as per the constitutive relationship derived above: where, for convenience: This expands to give a linear differential equation in standard form: which can be solved analytically to give stress as a function of time assuming σ t=0 = 0:

ANALYSIS AND IMPLEMENTATION
All computations for this paper, including analysis of the structural data and computation of the Maxwell models, were implemented in Python notebooks. These have been included below (1) for reference and (2) to highlight how the pycompass package can be used to analyse huge structural datasets extracted from digital outcrop models.

Structural analyses
Dyke orientation and thickness measurements extracted from the UAV surveys using the Compass plugin in CloudCompare have been exported to a .xml format that can be loaded and analysed using pycompass. These .xml files are avaliable here.
The following notebook provides a complete record of this workflow and the logic behind our structural analyses. Specifically, it contains the code required to create Figures 2, 3, 4, and S1.

Load data
Intrusion traces, orientation estimates, thicknesses measurements etc. have been exported into an .xml format for analysis. For convenience, we load all these datasets and group them into classes that make access convenient. During this step we also classify the intrusions into a steeply dipping (dyke) and shallow-dipping (sheet-intrusion) set.

Estimate swarm centre
Many of the dykes in the above dataset are radial to an unknown swarm centre. We can constrain how well candidate swarm centres explain the observed SNEs by grid-searching the model space and evaluating a likelihood function for each candidate centre location. First, we define our model searchspace. As the radial model has only two parameters (the x and y coordinate of the swarm center) then we can simply grid-search this whole space!

Estimate vertical strain
We repeat all of the above but using vertical rather than horizontal scan lines, to estimate vertical strain.

Plot strain rose (Fig. 6)
Plot all the strain estimates on a "rose diagram" to show the distribution of strain through Caldera Taburiente.

Thickness and aperture
Finally, we analyse the distribution of intrusion thickness (all measurements) and aperture (measurements of the thickest portions of individual intrusions). Although somewhat arbitrary, we define aperture as the 75-90th percentile of all the thickness measurements for each dyke.

Estimate and plot overpressure (Fig. 5)
Some of the UAV surveys contain dykes that are exposed from tip to tip, meaning the aperture-height scaling relationship can be measured and used to infer magma overpressure. This assumes that: 1. The dykes propagate horizontally 2. They are significantly longer than they are high (i.e. have a blade-shape) 3. Vertical overpressure variations are negligible Intrusions that are completely exposed are stored in a separate xml node, so we can extract them from the datasets as follows: #loop through SNEs belonging to this GeoObject and gather thickness estimates SNEs = d. Data from the steeply dipping intrusions (dykes) are then used to estimate the magma overpressures that could explain the observed aspect ratios. A variety of Young's modulii (from 1 -5 GPa) have been used for these calculations, and the top 25% of pressure estimates treated as outliers (as they result from unusually short/thick dykes).

Parmeter optimization
Assuming we can estimate reasonable values for the material properties of our volcano (from other sources), there are three main parameters in the above model that are unknown: h, P 0 and f d . These parameters can be directly related to observations of: 1. the maximum dyke aperture prior to a buildup of accomodation stress, 2. the final (equillibrium) dyke aperture and, 3. the final bulk-strain.

Stochastic model
The aperture distributions produced by the previous model have a shape that is generally similar to the observed aperture distribution, however the model cannot produce dykes with apertures less than the asymptote aperture of 1 m (resulting in the "spike" in aperture frequency at 1 m).
To better fit the observed data, we assume that the overpressure in ascending dykes is not uniform (it is very unlikely that it would be!), and add a random element to the model by sampling dyke overpressures from a normal distribution: P = N(µ = P 0 , σ) − σ accom Due to the added random element, this can no-longer be solved analytically, so we instead use a numerical solution: In [32]: """ Model the accomodation of 1d radial dykes by maxwell visco-elasticity using a time-stepping numerical approach. This allows the introduction of randomly varying magma pressure, determined the the parameter sigma.
**Arguments**: -E = Young's modulus (Pa) -v = poissons ratio -mu = viscosity (Pa.sec) -Pmax = initial fluid overpressure (Pa) -r = the distance from the radial swarm centre (m) -h = the height of the dykes (m) -Fd = the rate at which dyking events occur (dykes/year) -sigma = the standard deviation of the normal distribution (with mean Pmax) to sample overpressure from.

Parameter Optimization
As with the analytical solution, we optimise the model parameters against observations. As the aim of this model is to reproduce the observed dyke aperture distribution, we simply optimise the unknown parameters (the mean and standard deviation of P 0 , dyking frequency f d and height h) by maximizing the overlap between the observed and predicted aperture distributions.