# Fracturing around magma reservoirs can explain variations in surface uplift rates even at constant volumetric flux

Juliet Biggs, Timothy Rafferty, Jonathan Macha, Edna W Dualeh, Gregor Weber, Alain Burgisser, Finbar Carroll, Lauren Hart, Alison C Rust, Mark Gilbertson, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Juliet Biggs, Timothy Rafferty, Jonathan Macha, Edna W Dualeh, Gregor Weber, et al., Fracturing around magma reservoirs can explain variations in surface uplift rates even at constant volumetric flux. Journal of Volcanology and Geothermal Research, 2024, 452, pp.108129. 10.1016/j.jvolgeores.2024.108129. hal-04764557

HAL Id: hal-04764557

https://hal.science/hal-04764557v1

Submitted on 4 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Fracturing around magma reservoirs can explain variations in surface uplift rates even at constant volumetric flux

Juliet Biggs<sup>a, *</sup>, Timothy Rafferty<sup>b</sup>, Jonathan Macha<sup>b</sup>, Edna W. Dualeh<sup>a</sup>, Gregor Weber<sup>a</sup>, Alain Burgisser<sup>c</sup>, Finbar Carroll<sup>b</sup>, Lauren Hart<sup>b</sup>, Alison C. Rust<sup>a</sup>, Mark Gilbertson<sup>b</sup>, Alexandra Morand<sup>a</sup>

## ARTICLE INFO

### Keywords:

Surface deformation

Analogue models

Caldera unrest

## ABSTRACT

Many volcanoes show continuous but variable deformation over timescales of years to decades. Variations in uplift rate are typically interpreted as changes in magma supply rate and/or a viscoelastic response of the host rock. Here we conduct analogue experiments in the laboratory to represent the inflation of a silicic magma body at a constant volumetric flux, and measure the chamber pressure and resulting surface displacement field. We observe that dyke intrusions radiating from the magma body cause a decrease in the peak uplift rate, but do not significantly affect the spatial pattern of deformation or spatially averaged uplift rate. We identify 4 distinct phases: 1) elastic inflation of the chamber, 2) a gradual decrease in the rate of uplift and pressurisation, associated with the formation of visible cracks 3) propagation of a dyke by mode 1 failure at the crack tip and 4) a pressure decrease within the chamber. Phase 2 can be explained by either a) crack damage, which reduces the elastic moduli of the surrounding rock or b) magma filling pre-existing cracks. Thus these experiments provide alternative mechanisms to explain observed variations in uplift rate, with important implications for the interpretation of deformation patterns at volcanoes around the world.

## 1. Introduction

Volcanoes exhibit deformation with a variety of patterns and rates (Ebmeier et al., 2018; Biggs and Pritchard, 2017; Pritchard et al., 2022). At basaltic caldera systems, the rate of magma flux determines whether an eruption will occur (Galetto et al., 2022) while many silicic caldera systems show continuous deformation over timescales of years to decades without culminating in eruptions (Biggs et al., 2014; Acocella et al., 2015; Sandri et al., 2017). The rate of deformation at these silicic systems varies not only between volcanoes, but also as a function of time. Deformation is most often analysed using the simplest model, that of a pressure change in a point source embedded within an elastic half-space (Mogi, 1958). However, this is largely due to computational simplicity rather than realism. State-of-the-art numerical models consider a) the physical and chemical state of the chamber (e.g. pressure, temperature, gas and crystal content), b) fluxes into and out of the system and c) interactions with the surrounding host rock (e.g.

Degruyter and Huber, 2014; Anderson and Segall, 2011; Segall, 2016). Finite Element Models allow the investigation of a wider range of crustal heterogeneities and temperature-dependent mechanics (e.g. Gregg et al., 2012; Hickey et al., 2016; Alshembari et al., 2022). While these models can produce a range of complex behaviours and are thus capable of explaining variability in rates of surface uplift, they typically assume that the source geometry and material properties remains constant until the stresses in the host rock reach a failure criterion, at which point an eruption occurs (e.g. Gregg et al., 2012; Galetto et al., 2022). Laboratory studies have shown that not all dykes propagate to the surface (Taisne et al., 2011), and here we use laboratory-based analogue models to test whether the assumptions of a constant source geometry and material properties are realistic.

Experimental volcanology involves representing real-world processes through analogue modelling in a laboratory setting and is widely used in conjunction with field and theoretical works to investigate the physics of magmatic systems and to validate analytical and numerical

<footnote><sup>a</sup> COMET, School of Earth Sciences, University of Bristol, Bristol BS8 1RJ, UK</footnote>

<footnote><sup>b</sup> School of Electrical, Electronic and Mechanical Engineering, University of Bristol, Bristol BS8 1TR, UK</footnote>

<footnote><sup>c</sup> ISTerre, Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, Grenoble 38000, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: juliet.biggs@bristol.ac.uk (J. Biggs).</footnote>
models (Roche and Carazzo, 2019). McLeod and Tait (1999) conducted analogue experiments of magma chamber pressurisation and found that the propagation of pre-existing magma-filled cracks is the dominant mechanism of dyke nucleation, and can occur before the hoop stresses exceed the tensile strength of the surrounding rocks. In this case, magma viscosity determines the flux into the crack and hence the rate at which the crack pressurizes and the timescale over which an eruption will occur. Alternatively, rock mechanics experiments have shown that crack damage caused by locally high stresses can reduce the Young's modulus of volcanic rocks inducing nonlinear variations in magma overpressure and surface displacement (Heap et al., 2010; Got et al., 2017). Neither of these mechanisms are typically considered when modelling surface deformation, but have the potential to alter the interpretation of monitoring data.

We hypothesise that fracture-related processes can explain the variations in long-term deformation rate observed at many silicic systems. In this study, we explore some possible models using both analogue and numerical approaches and discuss their implications for the interpretation of satellite datasets. First, we build on the experiments of McLeod and Tait (1999) using a constant volumetric flux and monitor the pattern of surface deformation and chamber pressure that results (Bertelsen et al., 2021). We then explore two possible numerical approaches based on a) a pseudo-elastic medium in which damage reduces the elastic modulus of the surrounding rock and b) pre-existing cracks governed by fluid flow and fracture mechanics rather than the behaviour of the surrounding medium. We then qualitatively compare our results to satellite measurements of deformation at silicic calderas, and consider the role of cracks in causing variable deformation rates.

## 2. Analogue experiments

We adapt our experimental design from the studies of McLeod and Tait (1999) and Canon-Tapia and Merle (2006) who conducted analogue experiments to understand the growth of dykes from magma chambers, but perform our experiments at a constant volumetric flux. In addition, we measure the pattern of surface deformation following Tortini et al. (2014); Galland et al. (2016); Bertelsen et al. (2021); Rincón et al. (2022) which can then be compared to satellite observations.

### 2.1. Case study

Although the purpose of our models is to investigate fundamental processes and we make only qualitative comparisons to satellite observations, we use the geometry of Laguna del Maule as a foundation for our experimental design (Le Mével et al., 2021). Laguna del Maule is a rhyolitic caldera system in the Southern Volcanic Zone of Chile and has been uplifting continuously since around 2005 (Singer et al., 2014; Le Mével et al., 2021). The rate of uplift is variable: the rate increased during 2007–2010 reaching a peak of ~20 cm/year, decreased during 2010–2017 and increased again during 2017–2020 (Feigl et al., 2014; Le Mével et al., 2015, 2021). Novoa et al. (2019) model the 2010–2017 data assuming mafic magma recharge at the base of a crystal rich reservoir that behaves viscoelastically, but their timeseries ends before the 2017–2020 increase in uplift rate. Le Mével et al. (2021) explain the variations in uplift rate between 2007 and 2020 using a numerical model of viscous magma flowing through a conduit into a fluid-filled reservoir surrounded by a heterogeneous, viscoelastic crust and require two periods of magma injection, each causing an inlet pressure increase of ~9 MPa/year. The magma flow rate varied between 0 to $37 \times 10^6 \text{ m}^3\text{yr}^{-1}$, with an average of $24.3 \times 10^6 \text{ m}^3\text{yr}^{-1}$ for 2005–2020 (Le Mével et al., 2021).

### 2.2. Dimensional analysis

Scaling is a critical consideration for experimental volcanology as it

is necessary to ensure kinematic, dynamic and geometric similarity between the small-scale experiments and the large-scale, real-world volcanic phenomena. For the setup shown in Fig. 1a, surface deformation *d* can be described using ten parameters: density of the rock and magma ($\rho_r$, $\rho_m$); magma viscosity ($\mu$); volume flow rate (*Q*); chamber width (*W*) and depth (*D*); gravity (*g*); time (*T*); Young's modulus (*E*); Poisson's ratio ($\nu$); (Table 1). These can be represented using seven dimensionless numbers. $\Pi_1 = W/D$ governs the geometric similarity and scaling of the displacement between the model and the real world. $\Pi_2 = QT/W^3$ governs how the time frame of the experiment is related to the real world. $\Pi_3 = \frac{\rho_m Q}{\mu W} = \frac{\rho_m u W}{\mu}$ is the Reynolds number, which defines the balance between inertial and viscous forces in a fluid flow. $\Pi_4 = \frac{\rho_m}{\rho_r}$ governs if the system has positive, negative or neutral buoyancy. $\Pi_5 = \frac{\rho_r D_g}{E}$ governs the relationship between the weight from the overburden on the chamber and the country rock material stiffness. $\Pi_6 = \frac{E W^3}{\mu Q} = \frac{E W}{\mu \nu}$ is the ratio of fluid flow viscous stresses to elastic stresses. $\Pi_7 = \nu$ is the Poisson's ratio which is already dimensionless. If the material is assumed to be isotropic the Poisson's ratio acts as the scalar quantity between the shear, bulk and Young's modulus.

#### 2.2.1. Geometric scaling

To produce a feasibly-scaled laboratory experiment, we use a length-scaling of $10^{-5}$ and a cylindrical tank with a diameter of 29 cm and gelatin depth of 25 cm, representing a domain of $29 \text{ km} \times 25 \text{ km}$. We base the magma source geometry on the system at Laguna del Maule but even then the source geometry depends on the method, time period and model chosen. Le Mével et al. (2021) model the deformation at Laguna del Maule with an ellipsoidal source with principal axes $a = 7.06 \text{ km}$, ($a \in [5.70, 8.05] \text{ km}$), $b = 1.90 \text{ km}$, ($b \in [1.55, 2.74] \text{ km}$), and $c = 1.35 \text{ km}$, ($c \in [0.96, 2.09] \text{ km}$) at a depth $d = 4.78 \text{ km}$, ($d \in [5.48, 3.89] \text{ km}$). In contrast, Feigl et al. (2014) used a rectangular dislocation to represent a sill intrusion (Okada, 1985) and found a best-fit source with a length of $9.0 \pm 0.3 \text{ km}$, width of $5.3 \pm 0.4 \text{ km}$ and depth of $5.2 \pm 0.3 \text{ km}$. Other geophysical techniques provide an ever wider range of possible source geometries (e.g. Pritchard et al., 2018). Given these uncertainties, we conduct experiments using two simplified sources: a spherical source with a radius of 25 mm (which scales to 2.5 km) and a horizontal prolate source with a rough dimensions of $9 \text{ cm} \times 4 \text{ cm} \times 4 \text{ cm}$ (which scales to $9 \text{ km} \times 4 \text{ km} \times 4 \text{ km}$). The total volume of these chambers are $65 \text{ km}^3$ and $75 \text{ km}^3$ respectively.

Previous analogue and numerical experiments have set magma input rates to be either constant pressure (McLeod and Tait, 1999) or constant velocity (Merle and Vendeville, 1995; Guldstrand et al., 2017). In this study, we chose instead to use a constant volumetric flow rate, because it can be estimated for the real-world case using satellite data, unlike pressure and velocity. We use a pump flow rate of $6 - 10 \times 10^{-8} \text{ m}^3\text{s}^{-1}$. Based on our dimensional analysis ($\Pi_2$) and an average rate of $24.3 \times 10^6 \text{ m}^3\text{yr}^{-1}$ from Le Mével et al. (2021), this means that 1 s in the lab represents ~2 years in the real world. Thus a 20 min-long experiment corresponds to 5000 years, which is much longer than the few decades for which satellite measurements are available. In the natural world, it is unlikely that the magma supply rate would remain constant over such long time periods. We also note also that Laguna del Maule is one of the fastest deforming systems observed globally (Feigl et al., 2014), and lower rates may be more representative.

#### 2.2.2. Analogue materials

Gelatine is an ideal analogue for modelling homogeneous, isotropic and elastic materials and is often used to study the propagation of magma-filled fractures in the Earth's brittle and elastic crust (Kavanagh et al., 2013, 2015, 2017; Bertelsen et al., 2021; Smittarello et al., 2021). Kavanagh et al. (2013, 2017) and Smittarello et al. (2021) provide detailed studies of the material properties in this context. Although gelatine is viscoelastic, at low temperatures (5–10 °C) and for short
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Experimental setup. a) General setup with the pump on the left and the tank on the right. The surface deformation *d* depends on ten key parameters: density of the rock and magma (*ρ*<sub>r</sub>, *ρ*<sub>m</sub>); magma viscosity (*μ*); volume flow rate (*Q*); chamber width (*W*) and depth (*D*); gravity (*g*); time *T*; Young's modulus (*E*); Poisson's ratio (*ν*). b) Adaptions to the input system to incorporate syphon and pressure sensor, as described in section 2.3 c-e) illustrative photographs showing c) experimental setup with pump and tank; d) the balloon mould and e) the start of the experiment with the golden syrup added and the balloon removed.</figcaption></figure>

<table><caption>**Table 1**
Key system parameters and dimensionless numbers.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Real World</th><th>Analogue Model</th></tr></thead><tbody><tr><td>Magma Density (kgm<sup>-3</sup>)</td><td>ρ<sub>m</sub></td><td>2350–2600</td><td>1240</td></tr><tr><td>Rock Density (kgm<sup>-3</sup>)</td><td>ρ<sub>r</sub></td><td>2300–2600</td><td>1000</td></tr><tr><td>Magma Flow Rate (m<sup>3</sup>s<sup>-1</sup>)</td><td>Q</td><td>1.8</td><td>6 – 10 × 10<sup>-8</sup></td></tr><tr><td>Chamber Width (m)</td><td>W</td><td>9000</td><td>0.09</td></tr><tr><td>Chamber Depth (m)</td><td>D</td><td>5200</td><td>0.052</td></tr><tr><td>Magma Viscosity (Pa s)</td><td>μ</td><td>10<sup>5</sup> – 10<sup>10</sup></td><td>0.32</td></tr><tr><td>Young's Modulus (Pa)<sup>*</sup></td><td>E</td><td>10<sup>8</sup> – 10<sup>10</sup></td><td>3200 – 6200<sup>*</sup></td></tr><tr><td>Poisson's Ratio</td><td>ν</td><td>0.3</td><td>0.4–0.5</td></tr><tr><td>Timescaling<sup>**</sup></td><td>T</td><td>~ 2 yrs</td><td>1 s</td></tr><tr><td>Π<sub>1</sub></td><td>W/D</td><td>1.8</td><td>1.8</td></tr><tr><td>Π<sub>2</sub></td><td>QT/W<sup>3</sup></td><td>8.7 – 17 × 10<sup>-5</sup></td><td>8.7 – 17 × 10<sup>-5</sup></td></tr><tr><td>Π<sub>3</sub> (Reynolds Number)</td><td>ρ<sub>m</sub>Q/μW</td><td>10<sup>-11</sup> – 10<sup>-6</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Π<sub>4</sub> (Buoyancy)</td><td>ρ<sub>m</sub>/ρ<sub>r</sub></td><td>1.00–1.13</td><td>1.24</td></tr><tr><td>Π<sub>5</sub><sup>*</sup></td><td>ρ<sub>r</sub>Dg/E</td><td>0.01–0.13</td><td>0.08–0.16<sup>*</sup></td></tr><tr><td>Π<sub>6</sub><sup>*</sup></td><td>EW<sup>3</sup>/μQ</td><td>10<sup>9</sup> – 10<sup>16</sup></td><td>10<sup>7</sup> – 10<sup>8</sup><sup>*</sup></td></tr><tr><td>Π<sub>7</sub> (Poisson's Ratio)</td><td>ν</td><td>0.3</td><td>0.4–0.5</td></tr></tbody></table>

* The values of *E* and hence Π<sub>5</sub> and Π<sub>6</sub> are determined retrospectively in section 3.2.

** Applying Π<sub>2</sub> to the real and analogue values of *Q*, means that 1 s in the lab represents ~ 2 years in the real world.

periods of time (tens of minutes) the viscous component at the crack tips is negligible, so the gelatine can be considered purely elastic (Kavanagh et al., 2013). Under these conditions, failure occurs via mode I (opening) fracture but does not follow the Mohr-Coulomb failure criterion (Bertelsen et al., 2021). Young's modulus *E* and Poisson's ratio *ν* are a function of the concentration of gelatine, *C*<sub>gel</sub>. We selected *C*<sub>gel</sub> = 4 wt% using the dimensionless numbers Π<sub>5</sub>, Π<sub>6</sub> and Π<sub>7</sub>, which is consistent

with the recommendations of Kavanagh et al. (2013).

Water and oil are the most commonly used magma analogues (Kavanagh et al., 2018), and their viscosity scales appropriately for basaltic magma (μ = 1 – 10<sup>4</sup> Pas). Here we focus on rhyolitic magma which has a higher viscosity (μ = 10<sup>5</sup> – 10<sup>10</sup> Pas), thus we select golden syrup as a more appropriate analogue material based on dimensionless parameters Π<sub>3</sub> and Π<sub>6</sub>. The viscosity of golden syrup is temperature dependent, varies between batches and can be further modified by adding water (Beckett et al., 2011). We therefore test the sample viscosity under typical lab temperature conditions (18–22 °C) and at various dilutions using a concentric cylindrical rheometer. We found a concentration by volume of 80% golden syrup to provide the best fit to our dimensional analysis, with the additional benefit of having the lowest variation in viscosity with temperature. Like oil and water, golden syrup is a Newtonian fluid, which is appropriate for a pure melt, but does not consider bubble or crystal components that would be expected in a multi-phase magma and would affect magma compressibility and hence surface deformation (Kavanagh et al., 2018; Yip et al., 2024).

These analogue materials give good agreement between real world and analogue models for all dimensionless numbers with three exceptions (Table 1). The Reynolds number (Π<sub>3</sub>) is significantly higher in the analogue model than in nature, but both are small enough to lie well within the laminar regime. The Poisson's ratio (Π<sub>7</sub>) governs the relationship between Young's Modulus (*E*) and bulk (*K*) and shear modulus (*G*). The real world values of these moduli vary by orders of magnitude, so the precise value of the Poisson's ratio is not critical to these experiments. Finally, Π<sub>5</sub> indicates the relative importance of the weight from the overburden on the chamber and the Young Modulus, *E*, of the host material. *E* is sensitive to setting conditions (Kavanagh et al., 2013) and was hence determined retrospectively by fitting the measured pressures and surface displacements (see Section 3.2).
## 2.3. Experimental setup

We adapt our experimental design from the studies of McLeod and Tait (1999) and Canon-Tapia and Merle (2006). Our initial testing identified three areas where modification was required 1) the method of creating a chamber 2) a method for measuring pressure within the chamber and 3) the connection between the conduit and the chamber (Fig. 1). Similar solutions may have been employed by McLeod and Tait (1999) and Canon-Tapia and Merle (2006) but were not described in the resulting publications, thus we include details here in the interests of reproducibility.

In order to create a magma chamber, we use a balloon as a mould and fill the tank with liquid gelatin. We then cover the tank with plastic film to reduce dehydration and then place it in the fridge to set (Fig. 1d). We cover the balloon in vaseline to stop it sticking to the gelatine and use a counterweight system to hold the balloon in place while the gelatine sets. For a spherical balloon, the counterweight is attached to the balloon knot using fine fishing wire and for ellipsoidal balloons, we use a lasso to avoid creating geometric stress intensifiers. Once the gelatine has set, we pop the balloon with a needle and remove the balloon by pulling on the counterweight. The gelatine sets while the balloon is inflated with air, but once the balloon is popped and the cavity refilled with syrup, the chamber's pressure is no longer in equilibrium with lithostatic pressure in the gelatine. For this reason, we normalise our measurements to the starting time of the experiment.

Previous experiments struggled to remove trapped air produced when initially filling the chamber and found it necessary to limit the size of the magma chamber to reduce buoyancy effects (Canon-Tapia and Merle, 2006; Menand and Tait, 2001). To overcome this, we insert a secondary closeable pipe to fill the magma chamber without pressurizing it as this acts as a syphon allowing the gas to escape (Fig. 1b). Altering the height of the syphon allows for the control of the volume of air trapped in the analogue magma chamber, and removes the constraints on magma chamber size.

We use an Adafruit MPRLS flow pressure sensor linked to an Arduino Micro logger to measure the pressure within the chamber. This has an estimated accuracy of ~10 Pa in the expected pressure range of 0.4 – 1.2 kPa (McLeod and Tait, 1999). We use a push-fit pipe adaptor to permit an easy change between pressure sensors (Fig. 1b). The pipe was also fitted with a silicon compression olive and a smoothed adaptor, which allowed for the height of pipe, and as such the chamber depth, to be easily varied (Fig. 1b). Finally, we attached a hyperbolic paraboloid disk to the end of the pipe to create a better connection between the pipe and the chamber, avoiding unrealistic dynamics.

Finite Element analysis conducted in Abaqus showed that boundary effects caused by the wall of the tank (29 cm diameter) would have little effect on the deformation caused by pressurisation of the chamber. However, the analysis did not consider the boundary effects on the propagation paths of the dykes and these are likely significant. Therefore, we do not attempt to analyse the trajectories of the dykes.

The experiments were conducted at room temperature (18–22°), using gelatin which was initially at a temperature of ~4°C and syrup at room temperature. This setup means that there are likely to be temperature gradients within the system, which in turn will affect the Young's modulus of the gelatine. Over the short timescales of the experiments (<30 mins), this is unlikely to affect the fundamental processes that we aim to observe, but, among other factors, limits the quantitative interpretation of the models.

## 2.4. Measurements of surface deformation

A variety of techniques have been used in the past to measure surface deformation in analogue volcanic experiments, each with strengths and weaknesses (Tortini et al., 2014; Rincón et al., 2022; Galland et al., 2016). Motion sensor cameras were originally designed for the gaming industry, but have been used more recently in experimental applications

in Earth Sciences (Rincón et al., 2022; Tortini et al., 2014). The Microsoft Kinect v2 has been used for previous analogue experiments in volcanology, but its low accuracy (> 1 mm) makes it unsuited to these experiments (Rincón et al., 2022; Tortini et al., 2014; Lun and Zhao, 2015). Therefore, we use the Intel RealSense Depth Camera D435, a wide-angled depth camera designed for fast-moving applications (e.g. drones), which uses photogrammetry enhanced with an IR projector to track the surface. The RealSense D455 claims twice the precision of the Kinect v2 at ±0.5 mm and an improved maximum rate of capture of 90 Hz. The range of the RealSense permits closer recording than the Kinect2, which further increases precision. Unlike other methods which measure the full 3-D deformation field (e.g. Digital Image Correlation), motion sensor cameras only measure the vertical component of deformation, but that is sufficient for the experiments reported here.

For motion sensor cameras to work effectively, the surface of interest needs to be rough, opaque and well-lit, but too much light can produce blindspots (Tortini et al., 2014). Since gelatine is smooth and transparent, we first add a layer of wholemeal flour to reduce glare, followed by a layer of mixed, dyed poppy seeds to generate a speckle pattern. Two LED illumination units provided consistent lighting; however, for some experiments, reflections meant that certain areas of the surface needed to be masked during post-processing. To reduce random noise in the data, groups of ten files were averaged together to create a data file for every second (Tortini et al., 2014). To remove coherent noise, the data was transformed into the frequency domain and a Gaussian filter was applied in both the x and y directions.

## 2.5. Experimental results

Initially, we describe the results from the experiments with an ellipsoidal chamber (Fig. 3a), and then compare our observations to those from a similar experiment with a spherical chamber (Fig. 3b). Neither experiment resulted in an eruption at the surface, but both produced stalled dykes with a width of approximately 4 mm. The first dykes formed on the base of the magma chamber and grew radially (green outline in Fig. 2) while secondary dykes formed, at an oblique angle (red and blue outlines in Fig. 2). Dykes initiating from the top of the chamber towards the top surface either stagnated or took a parabolic trajectory to avoid the free surface. In contrast, dykes propagating from the base of the magma chamber grew more rapidly, accelerating downwards away from the chamber. The effect of bouyancy is discussed further in section 2.6.

For the ellipsoidal chamber, we identify 5 distinct phases of activity based on the timeseries of surface displacement and chamber pressure (Fig. 3a):

* **Phase 0 (0–60s):** Initially, the deformation is beneath the detection limit of the instrument of ±0.5 mm (section 2.2).

* **Phase 1 (60–120 s):** The chamber pressure and surface displacement increase roughly linearly, consistent with elastic inflation of the chamber.

* **Phase 2 (200 s–300 s):** The rate of displacement and pressure increase slows, and small cracks start to appear on the chamber surface.

* **Phase 3 (~250–350 s):** There is a sudden pressure drop, consistent with the propagation of a dyke (Kavanagh et al., 2017). This occurred at ~270 s for the ellipsoidal chamber and ~350 s for the spherical chamber.

* **Phase 4 (350–500 s).** Following the dyke intrusion, the chamber pressure drops rapidly and the rate of surface deformation falls further.

We calculate the deformation rate in each of these phases (Table 2). The total volumetric flux into the system is constant, and this is reflected by the mean rate of deformation across the whole image, which remains roughly constant at 1.9 – 2.1 µms⁻¹. However, the rate of maximum displacement gradually decreases from 9.6 µms⁻¹ in phase 1, to 7.2 µms⁻¹ in phase 2, and 4.0 µms⁻¹ in phase 4. Using Π₂, the measured rates of uplift (4.0–9.6 µms⁻¹) correspond to 20–50 cm/yr in real
<figure><img src="image_5.png" /><figcaption>Fig. 2. Photographs from the experiment with the ellipsoidal chamber showing the paths taken by the dykes from two different orientations.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Time series of surface displacement and chamber pressure. a) Experiment with ellipsoidal chamber. Mean displacement is calculated over the area shown in Fig. 4. Rates are given in Table 2. b) Comparison between displacement time series for ellipsoidal (blue) and spherical chambers (red). Observed dyke propagation is indicated by an arrow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

**Table 2**

Measured rates of displacement and pressure increase during phases 1,2 and 4 of the ellipsoidal chamber experiment. The time periods used to estimate the rates are given in the first row.

<table><thead><tr><th></th><th>Units</th><th>Phase 1</th><th>Phase 2</th><th>Phase 4</th></tr></thead><tbody><tr><td>Measurement Time</td><td>s</td><td>60–120</td><td>140–200</td><td>350–500</td></tr><tr><td>Max Disp. Rate</td><td>µms<sup>-1</sup></td><td>9.6</td><td>7.2</td><td>4.0</td></tr><tr><td>Mean Disp. Rate</td><td>µms<sup>-1</sup></td><td>2.1</td><td>1.9</td><td>2.1</td></tr><tr><td>Pressure Rate</td><td>Pas<sup>-1</sup></td><td>6.5</td><td>3.6</td><td>-1.0</td></tr></tbody></table>

systems, slightly larger than but on a similar order of magnitude to those observed at Laguna del Maule (Le Mével et al., 2021). Despite the differences in rate and mechanism, phase 1 and phase 2 both behave elastically with the maximum surface deformation directly proportional to the chamber pressure.

To compare the pattern of surface displacement in each of these 3 phases, we taken time slices at 100 s during the initial elastic phase, at 250 s during the crack opening and at 450 s after the dyke propagation (Fig. 4). The pattern of displacement is very similar both qualitatively (Fig. 4a-c) and quantitatively (Fig. 4d-f). The variations (±2 mm) are significantly larger than the estimated accuracy of the measurement technique (±0.5 – 0.6 mm), but do not show any systematic patterns, thus likely reflect mechanical heterogeneities in the gelatine.

We repeat the experiment using a spherical chamber and see the same phases of behaviour (Fig. 3b). In both cases, the average rate of surface displacement is constant, confirming the constant rate of magma

supply. Initially, the maximum surface displacement also increases linearly (phase 1), but the rate gradually slows (phase 2) prior to the dyke injection (phase 3). In this case, two dykes propagated very close together in time. This is followed by a decrease in the rate of maximum surface displacement (phase 4).

## 2.6. Buoyancy

In our experiments, dykes initiating from the top of the chamber towards the top surface either stagnated or took a parabolic trajectory to avoid the free surface whereas dykes propagating from the base of the magma chamber grew more rapidly, accelerating downwards away from the chamber. In contrast, similar experiments by Canon-Tapia and Merle (2006) found that the dykes propagating from the base of the chamber stagnated and those out of the top grew. The crucial difference between the investigations is the fluid density, as the trajectory of the dykes depends on their buoyancy, $\Pi_4 = \frac{\rho_m}{\rho_s}$, where $\Pi_4 < 1$ represents positively buoyant magma and $\Pi_4 > 1$ represents negatively buoyant magma. We use golden syrup which is negatively buoyant ($\Pi_4 = 1.24$) while Canon-Tapia and Merle (2006) used oil which is positively buoyant.

We ran two additional experiments using vegetable oil cooled to 4 °C ($\Pi_4 = 0.89$) and castor oil at room temperature ($\Pi_4 = 0.95$) as the magma analogues, creating positively buoyant and roughly neutral systems respectively. Both were measured to have a viscosity of the same order of magnitude as the golden syrup solution. For the vegetable oil experiment, the dykes propagating from the bottom of the chamber stagnated and those from the top grew continuously and erupted
<figure><img src="image_7.png" /><figcaption>Fig. 4. Spatial pattern of deformation measured using a RealSense D355 camera for an ellipsoidal chamber. a) time slice at 100 s, during the elastic inflation of the chamber (phase 1), b) time slice at 250 s during opening of pre-existing cracks (phase 2) and c) time slice at 450 s, after dyke propagation (phase 4). Differences between spatial pattern in each stage, calculated by linearly scaling the pattern at the second time slice to match the magnitude of the first time slice and subtracting d) time slice 250 s minus time slice 100 s (phase 2 - phase 1); e) time slice 450 s minus time slice at 250 s (phase 4 - phase 2); and time slice 450 s minus time slice at 100 s (phase 4 - phase 1). Contours are every 1 mm and measurement uncertainties are ± 0.5–0.6 mm and tank diameter is 29 cm.</figcaption></figure>

through the surface. For the castor oil experiment, the dykes initiated from the base of the chamber but then grew only upwards around the chamber and erupted through the top surface. These results are consistent with the observations of Canon-Tapia and Merle (2006) and confirm that the trajectory of the dykes in our experiments is determined by magma buoyancy.

Previous studies have shown that numerous mechanisms can lead to the arrest of positively buoyant dykes (Davis et al., 2020; Taisne et al., 2011; Kavanagh et al., 2006). Therefore, stalled dykes would be expected to occur even in a positively buoyant system, especially if the pressure drop within the chamber was sufficient to cut off the magma supply before the critical volume was reached, or if the host rocks were mechanically heterogeneous.

## 3. Physical models

Our experiments were conducted at a constant magma flux and initially show a linear response of chamber pressure and surface uplift (phase 1), followed by decrease in the rate of pressurisation and peak uplift (phase 2), prior to dyke formation (phase 3). Our interpretation is that phase 1 corresponds to a pure elastic response of the gelatine to the reservoir inflation. The decrease in uplift rate during phase 2 could be explained with either of two mechanisms: a) crack damage, which can be described by a pseudo-elastic medium in which the elastic moduli of the surrounding rock are reduced (Got et al., 2019) or b) magma filling pre-existing cracks, which are governed by fracture mechanics rather than the behaviour of the surrounding medium (McLeod and Tait, 1999). In this section, we compare our experimental results to the models of Got et al. (2019) and McLeod and Tait (1999), respectively.

### 3.1. Magma-filled cracks

McLeod and Tait (1999) conducted analogue experiments to investigate the pressurisation and growth of liquid-filled cavities. Their theoretical analysis showed that dyke nucleation by the growth of pre-existing magma-filled cracks can occur before the hoop stresses exceed the tensile stress of the surrounding rocks. Pre-existing cracks in the wall rocks are to be expected considering the ubiquitous pervasive fracturing of the crust and thermal influence of the magma storage. Thus growth of pre-existing magma-filled cracks is expected to be the dominant mechanism of dyke nucleation in nature as well as analogue experiments.

In this model, flow of magma from the chamber into the crack gradually pressurizes and widens the crack. The evolution of the crack pressure, P, is governed by the differential equation,

$$ \frac{dP}{dt} = \frac{(P - \sigma_t)^3 (P_c - P)}{3m^2\mu} $$

where $P_c$ is the chamber pressure, $\sigma_t$ is the hoop stress, $m$ is the elastic stiffness of the host rock and $\mu$ is the magma viscosity. Thus the rate of crack pressurisation is independent of crack length and lithostatic loading, but there is a direct dependence on magma viscosity, which can vary by orders of magnitude. Crack growth will occur once the crack-tip stress intensity factor, $K = 1.12(P - \sigma_t)\sqrt{l}$, exceeds the fracture toughness of the medium, $K_c$, where $l$ is the crack length. For the case of instant pressurisation, dyke nucleation is predicted to occur on the order of hours to days for low-viscosity basaltic magmas, but could take several years for higher-viscosity rhyolites.

Although McLeod and Tait (1999) noted surface deformation in their experiments, their focus was on the timescales of dyke formation and they did not attempt to quantitatively measure or model the associated surface deformation. The analytical solutions consider either an instant step change in chamber pressure or a gradual pressurisation at a constant rate, but do not consider the case of constant magma supply rate and are therefore not directly comparable to our experiments. However, their model is qualitatively consistent with our observations. As magma from the chamber slowly widens the pre-existing cracks, the rate of chamber pressurisation and hence surface uplift gradually decreases. Initially, the low flux of magma into the dyke does not completely counteract the background pressurisation within the chamber, and the chamber pressure continues to rise, but at a reduced rate (phase 2). However, once the dyke starts to propagate, the magma flux out of the chamber dominates and there is a net depressurization of the chamber (phase 3).

### 3.2. Damage model

In an alternative model, the small cracks that appeared during phase 2 of our experiments are consistent with the development of a damaged shell around the pressurized reservoir (Got et al., 2017). The pressure and displacement of our experiments can be modelled using the work of Got et al. (2017), which assumes that damage weakens the Young modulus of the host medium according to:

$$ E(t) = E_0 \left( 1 - \frac{t}{t_c} \right)^a $$

where $a$ characterizes the rate of incremental damage (the purely elastic limit is recovered when $a = 0$), $E_0$ is the undamaged Young modulus and $t_c$ is the critical time for rupture. The relation between the injected
volume, ΔV, and the overpressure variation in the reservoir, ΔP, is (Reverso et al., 2014; Got et al., 2017):

$$ \Delta V(t) = \frac{3V\gamma}{4} \frac{\Delta P(t)}{E(t)} \qquad (1) $$

where V is the reservoir volume, γ = 1 for a sphere and γ = 8(1 - ν)/(3π) for an oblate reservoir. In our experiments, the fluid flow rate, Q, is constant, which precludes us from using the analytical solution of Got et al. (2017) because they considered constant overpressure instead. Deriving (1) with respect to time and using that ∂tΔV = Q yields:

$$ Q = \frac{3V\gamma}{4} \partial_t \left( \frac{\Delta P(t)}{E(t)} \right) $$

The rate of overpressure change is thus:

$$ \partial_t(\Delta P(t)) = \frac{4Q}{3V\gamma}E(t) - \frac{a\Delta P(t)}{t_c - t} \qquad (2) $$

Surface displacements, u, taking into account the free-surface effect can be calculated with Got et al. (2017):

$$ u(t) = f \Delta P \left( \frac{t_c}{t_c - t} \right)^a \qquad (3) $$

where f = u₀/P₀ is the ratio of initial displacement, u₀, over initial reservoir pressure, P₀.

We performed a least-square fit of the maximum displacement and pressure data of Fig. 3 by using the values of ν = 0.45, t_c = 290 s for the

ellipsoidal chamber, t_c = 350 s for the spherical chamber, and the bounds of Q = 6 - 10 × 10⁻⁸ (Table 1). Initial time was set to 0 but the data from Phase 0 (≤ 60 s) were left out of the fitting procedure. The free parameters a, E, and f were constrained by their physical limits whenever possible: 0 ≤ a ≤ 1, E₀ > 0, f > 0, 0 < γ ≤ 0.47 for the ellipsoidal chamber, and γ = 1 for the spherical chamber.

Gelatine is widely used in analogue studies of magmatic systems and its elastic properties have been well documented (Kavanagh et al., 2013). Such studies suggest that for a gelatine concentration, Cgel, of 4 wt %, the Young's modulus, E, should be ~ 14 kPa. However, this is strongly dependent on setting time and conditions, so we checked these values using the relationship between pressure and displacement during phase 1 of the experiments given by the analytical solutions for deformation in an elastic half-space (Mogi, 1958; Yang et al., 1988). Using this method, our estimates of E were 3.2 kPa and 6.2 kPa for the ellipsoidal and spherical chamber models, respectively. If the fitting procedure is done with E constrained at values > 1 kPa, the lowest bound of E is systematically the best-fit value. Such behaviour generally indicates that the bracket of tested E values should include values < 1 kPa. As such values are significantly below the first-order estimate done above, we test two sets of conditions: 1) E is unconstrained and 2) E is in the range 3–6 kPa (Fig. 5).

The pressure and displacement fits for the ellipsoidal chamber are visually comparable, regardless of whether E is constrained or not (Fig. 5a-b). This is also the case for the displacement fit for the spherical chamber, but the pressure fit with E constrained to ≥ 3 kPa is poor (Fig. 5c-d). The model output of the displacement is independently

<figure><img src="image_8.png" /><figcaption>Fig. 5. Best-fits of displacement and pressure for the damage model. Grey curves represent data, solid black curves represent unconstrained E, and dashed black curves represent E = 3 kPa. Best-fit values of the free parameters are listed for each geometry with the unconstrained E values in parenthesis. A) Displacement of the ellipsoidal chamber. B) Pressure of the ellipsoidal chamber. C) Displacement of the spherical chamber. D) Pressure of the spherical chamber.</figcaption></figure>
scaled to that of the pressure with the variable *f* (ratio of initial displacement over initial reservoir pressure). The good displacement fit of the spherical chamber with constrained *E* is thus explained by a *f* value that varies by a factor 2.5 between the cases of free and constrained *E* (in comparison, *f* varies by only 13% for the ellipsoidal chamber). The pressure fits for the ellipsoidal chamber are better at high *E* if the earlier (< 150 s) part of the pressure evolution is considered, and better at low *E* if the whole duration of phases 1 and 2 is considered. So either the spherical chamber model had an (unexplained) issue with the pressure gauge, or the Young modulus of our gelatin was lower than usually reported. For simplicity, we favor the low Young's modulus estimates given by the unconstrained fitting procedure, which we analyse below. Even using the lowest obtained values of *E* = 0.7 ± 0.1 kPa and *E* = 0.833 ± 0.066 kPa, the dimensionless numbers Π₅ = 0.08 – 0.16 and Π₆ = 10⁷ – 10⁸ are still consistent with the natural ranges.

The unconstrained values for *E* and *a* overlap (Fig. 5), which indicates that, as expected, gelatine properties did not vary between experiments. The different ranges covered by *f* are expected because the two chamber geometries yield different initial pressures. The limit between Phases 1 and 2 is ~ 250 s for the ellipsoidal chamber and ~ 180 s for the spherical chamber, which is consistent with our observations (Fig. 5A-C). The low *a* values indicate that damaging is minimal during deformation because its value is close to the elastic limit of *a* = 0. The damaged region was thus a relatively thin shell around the chamber walls, which is consistent with the small cracks observed during Phase 2.

# 4. Discussion

## 4.1. Magmatic architecture

For the last century, the dominant paradigm has been that magma is stored in a melt-dominated, elastic-walled magma chamber. However, there is now ample geophysical and petrological evidence that melt-filled pockets are only one component of a much larger magmatic system with a heterogeneous distribution of melts, crystals and exsolved volatiles (Cashman et al., 2017; Edmonds et al., 2019). We must therefore consider whether these experiments still provide an appropriate analogue to our current conceptual understanding of magmatic systems.

Transcrustal magmatic systems likely develop over periods of tens to hundreds of thousands of years as magma rises in discrete, short-lived pulses and stalls at rheological or lithological boundaries in the upper crust, where it cools and crystallises forming composite bodies (Annen et al., 2015). Numerical models of heat diffusion can therefore be used to study the temperature structure and rheology of magmatic systems. These suggest that relatively young systems (<10 kyr) are more likely to create the high temperature gradients necessary to juxtapose a high melt fraction with country rock at a low enough temperature to behave elastically (Sparks et al., 2019).

In Fig. 6 we present 2D thermal models depicting pulsed magma injection into a 20 km crustal section. These models were generated using the open-source Julia package 'MagmaThermoKinematics' (Schmitt et al., 2023). In these simulations, magma is injected into the system as dykes and sills at a temperature of 1000 °C. The model simulates the injection of sills and dikes with fixed dimensions of 3 km by 0.2 km over a random depth range between 5 and 15 km. The orientation of the dikes is randomly selected between 80 and 100°, while the

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Modelled temperature and viscosity structure of a two-dimensional crustal section that has experienced repeated injection of magmatic dykes and sills. A) Initial temperature distribution for a geothermal gradient of 25 °C. B) Temperature distribution after 400 ka of magma injection. The 700 °C contour is shown, corresponding to the modelled solidus temperature. C) Temperature field after 800 ka of magma injection. D) Viscosity structure in logarithmic units for time 0 ka. Contours are shown for a viscosity of 10¹⁸, 10²⁰, and 10²⁵ Pa s, corresponding to temperatures of 700, 500 and 250 °C, respectively. E) Viscosity structure corresponding to 400 ka, and F) for 800 ka of magma injection.</figcaption></figure>
sills are randomly oriented within –10 to 10°. Sills are injected at depths >12 km, while dikes are injected at shallower depths. Overall, this corresponds to a volumetric magma flux of 1.11 × 10⁻⁵ km³ km⁻² yr⁻¹. We show the resulting temperature (Fig. 6A-C) and viscosity (Fig. 6D-F) structure of the crust after 0, 400 and 800 thousand years (ka) of magmatic activity. Viscosity (μ) is determined from the temperature (T) field in Kelvin using an Arrhenius relation: μ = Ae^(H/RT), where A is the Dorn parameter with a value of 10⁹ Pa s (Del Negro et al., 2009), H represents an activation energy of 150 kJ mol⁻¹ (Ranalli, 1995), and R is the universal gas constant.

Initially, before the onset of magma injection, viscosity values of 10¹⁸ and 10²⁵ Pa s correspond to 20 km and 10 km depth, respectively. This suggests that the upper crust would display elastic behaviour when subjected to deformation. As the crust undergoes progressive thermal maturation after 400 and 800 ka of magma intrusion, a mush body develops, surrounded by a viscoelastic shell with viscosities ranging approximately from 10¹⁸ – 10²⁰ Pa s. Although the ongoing injection of magma over time leads to the expansion of the viscoelastic shell, crustal rocks shallower than 5 km maintain relatively high viscosities even after 800 ka of magma injection. This is due to the thermal buffering effect exerted by the Earth surface. Consequently, in young and low magma flux systems, as well as at shallow crustal conditions, the dominant tendency is for elastic deformation of magma reservoir walls.

However, in a high temperature, mush system, it is also unlikely that the source geometry will remain fixed with time. Localised bands of high strain also occur in partially-molten aggregates, especially under extensional stresses. Tension across a weak, high-porosity region leads to low pressure that, in turn, causes convergence of melt flow into that region, raising its porosity and further weakening it (Katz et al., 2006). The result is the formation of a high porosity, high strain melt band. In this case it is not pre-existing cracks that localise strain into narrow bands, but small heterogenities in melt fraction within the mush. However, the end result would likely be a similar effect on surface deformation to the stalled dykes discussed here.

## 4.2. Comparison to observations

Our experiments show that anelastic behaviour associated with

fracturing is able to explain variations in uplift rate, but it does not mean that this mechanism is responsible for all such variations. Variability in supply rate and viscoelastic behaviour are both physically-realistic mechanisms that are equally able to explain the observation. Unfortunately, the spatial and temporal patterns of deformation produced by fracturing would be very similar to those associated with a change in rate of magma supply or a viscoelastic rheology. Nonetheless, we consider two types of observations that may allow us to distinguish the mechanism responsible for variations in uplift rate: the spatially-averaged rate of displacement and the pattern of seismicity.

### 4.2.1. Spatially-averaged displacement

The average displacement rate in our experiments remains constant with time because the volumetric flux is held constant. However, for viscoelastic relaxation or variations in supply rate, we would expect variations in the spatially-averaged rate to mirror those in the peak rate. Here we test this hypothesis using deformation patterns from six silicic systems: Laguna del Maule, Chile; Campi Flegrei, Italy; Domuyo, Argentina; Corbetti, Ethiopia; Tullu Moye, Ethiopia and Suswa, Kenya (Fig. 7). The deformation time series cover 2014–2023 and were produced from Sentinel-1 InSAR data using the LICS system (Lazecký et al., 2020; Morishita et al., 2020). The time-series of peak displacement illustrate the range of observed patterns: Domuyo uplifted at a rate of 14.3 cm/yr in 2014–2019, but has been subsiding even since; Tullu Moye experienced 9.6 cm of uplift in 2015–2017, followed by slower uplift at 1.9 cm/yr and the rate of uplift at Laguna del Maule is steady over this time interval, but significant variations have previously been reported (Le Mével et al., 2021).

For each of these 6 timeseries, we also calculate time series of the spatially-averaged displacement over an area of 250 km² (Fig. 7). Uncertainties are dominated by the choice of reference pixel (Albino et al., 2022) and are estimated by calculating all timeseries for a single pixel (within area of deformation) using all suitable reference pixels, those having a root mean square error within the 85th percentile. We find that variations in the average displacement mirror those in the maximum displacement, for all five of the six time series. The only example for which this is not statistically true is Suswa, but the deformation here is small and the variations appear to be due to low signal-to-noise ratio rather than systematic variations. Thus based on the correlation

<figure><img src="image_10.png" /></figure>

Fig. 7. Deformation timeseries for six silicic systems showing variations in the mean and maximum rates of uplift. The timeseries span 2014–2023 and were produced using data from the Sentinel-1 satellite using the LICS system (Lazecký et al., 2020; Morishita et al., 2020). Uncertainties are dominated by the choice of reference pixel (Albino et al., 2022) and are estimated by calculating all timeseries for a single pixel (within area of deformation) using all suitable reference pixels. Reference pixels are determined based on having a RMS within the 85th percentile.
between time-series of peak and spatially-averaged deformation, we conclude that processes other than the dyking are likely to be responsible for the changes in rate at these examples.

### 4.2.2. Seismicity

Either brittle damage or the formation of stalled dykes might be expected to cause an increase or change in the spatial pattern of seismicity. At our case study, Laguna del Maule, frequent seismic swarms occur around the Troncoso Fault to the SW of the modelled magma source (Le Mével et al., 2021) and small local faults in the hanging wall of the Troncoso Fault facilitated repeated, small volume rhyodacitic and andesitic eruptions during the post-glacial period (Garibaldi et al., 2020), supporting the idea that these structures control magmatic pathways. The rate of seismicity was high in 2012–2016 during the deceleration of uplift consistent with a decrease in pressurisation caused by the opening of pre-existing cracks. There was also a large increase in seismicity in 2018–2020, at least a year after the uplift rate began to accelerate, which Le Mével et al. (2021) interpret as fluid motion or reactivation of local faults caused by stress redistribution associated with the pressurizing magma reservoir. The complexity of the behaviour illustrates the challenge of producing a unique interpretation for geophysical observations where multiple processes occur.

## 5. Conclusions

The analogue experiments described here demonstrate that variations in surface uplift rate can occur even when the volumetric flux is constant due to fracturing of the surrounding rock. The time-series measurements of surface deformation and chamber pressure identify 4 distinct phases: 1) elastic inflation of the chamber, 2) crack opening/formation associated with a gradual decrease in the rate of uplift and pressurisation, 3) propagation of a dyke by mode 1 failure at the crack tip and 4) a pressure decrease within the chamber. Importantly, these variations do not significantly affect the spatial pattern of uplift and fracturing would therefore be difficult to detect in satellite-based observations, even at high spatial resolution.

While this suggests an alternative mechanism for variations in surface uplift rate there is, as yet, little evidence that this actually occurs over the decadal timescales of the available observations. When appropriately scaled, our experiments show that dyke propagation would occur after several hundred years, but this would be significantly less for magma with a lower viscosity. Although brittle failure is not necessarily compatible with the current conceptual understanding of mushy magma systems, other mechanisms of strain and melt localisation exist in partially molten aggregates and would be expected to have a similar effect on observations of surface uplift. Thus, while we have demonstrated that this new mechanism is plausible, further work is required to investigate whether it actually occurs. Nonetheless, the implications for forecasting are important. The formation of stalled dykes do not significantly affect the spatial pattern of deformation and can form at lower overpressures than previously considered. In fact, they may be occurring undetected on a regular basis. Thus, the conditions for dyke propagation rather than formation determine whether an eruption will occur and models of dyke propagation and pathways rather than just reservoir failure are needed to assess volcanic hazard.

### Open data statement

Custom code to analyse the RealSense data is available in the project GitHub: https://github.com/TMecEng/MAST_analogue_modeling. Processed Sentinel-1 data is available from the COMET Volcano and Magmatic Deformation Portal: https://comet.nerc.ac.uk/comet-volcano-portal/. The 'MagmaThermoKinematics' package is available from https://github.com/boriskaus/MagmaThermoKinematics.jl).

### CRediT authorship contribution statement

**Juliet Biggs:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.
**Timothy Rafferty:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation.
**Jonathan Macha:** Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation.
**Edna W. Dualeh:** Writing – review & editing, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation.
**Gregor Weber:** Writing – review & editing, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation.
**Alain Burgisser:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.
**Finbar Carroll:** Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis. **Lauren Hart:** Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation.
**Alison C. Rust:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization.
**Mark Gilbertson:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Investigation, Conceptualization.
**Alexandra Morand:** Methodology, Writing – review & editing.

### Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

The manuscript includes an open data statement which includes links to all data and codes.

### Acknowledgements

JB, AR, AB, GW and ED received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (MAST; grant agreement No. 101003173). JB, GW and ED are part of COMET, the NERC Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics, a partnership between UK universities and the British Geological Survey TR, JM, FC, LH and MG contributed as part of an MSci Group Industrial Project.

### References

Acocella, V., Di Lorenzo, R., Newhall, C., Scandone, R., 2015. An overview of recent (1988 to 2014) caldera unrest: Knowledge and perspectives. *Rev. Geophys.* 53, 896–955.

Albino, F., Biggs, J., Lazecký, M., Maghsoudi, Y., 2022. Routine processing and automatic detection of volcanic ground deformation using sentinel-1 Insar data: insights from african volcanoes. *Remote Sens.* 14, 5703.

Alshembari, R., Hickey, J., Williamson, B.J., Cashman, K., 2022. Poroelastic mechanical behavior of crystal mush reservoirs: Insights into the spatio-temporal evolution of volcano surface deformation. *J. Geophys. Res. Solid Earth* 127 e2022JB024332.

Anderson, K., Segall, P., 2011. Physics-based models of ground deformation and extrusion rate at effusively erupting volcanoes. *J. Geophys. Res. Solid Earth* 116.

Annen, C., Blundy, J.D., Leuthold, J., Sparks, R.S.J., 2015. Construction and evolution of igneous bodies: Towards an integrated perspective of crustal magmatism. *Lithos* 230, 206–221.

Beckett, F., Mader, H., Phillips, J., Rust, A., Witham, F., 2011. An experimental study of low-Reynolds-number exchange flow of two newtonian fluids in a vertical pipe. *J. Fluid Mech.* 682, 652–670.
Bertelsen, H.S., Guldstrand, F., Sigmundsson, F., Pedersen, R., Mair, K., Galland, O., 2021. Beyond elasticity: are coulomb properties of the earth's crust important for volcano geodesy? J. Volcanol. Geotherm. Res. 410, 107153.

Biggs, J., Pritchard, M.E., 2017. Global volcano monitoring: what does it mean when volcanoes deform? Elements 13, 17–22.

Biggs, J., Ebmeier, S., Aspinall, W., Lu, Z., Pritchard, M., Sparks, R., Mather, T., 2014. Global link between deformation and volcanic eruption quantified by satellite imagery. Nat. Commun. 5, 3471.

Canon-Tapia, E., Merle, O., 2006. Dyke nucleation and early growth from pressurized magma chambers: Insights from analogue models. J. Volcanol. Geotherm. Res. 158, 207–220.

Cashman, K.V., Sparks, R.S.J., Blundy, J.D., 2017. Vertically extensive and unstable magmatic systems: a unified view of igneous processes. Science 355 eaag3055.

Davis, T., Rivalta, E., Dahm, T., 2020. Critical fluid injection volumes for uncontrolled fracture ascent. Geophys. Res. Lett. 47 e2020GL087774.

Degruyter, W., Huber, C., 2014. A model for eruption frequency of upper crustal silicic magma chambers. Earth Planet. Sci. Lett. 403, 117–130.

Del Negro, C., Currenti, G., Scandura, D., 2009. Temperature-dependent viscoelastic modeling of ground deformation: Application to Etna volcano during the 1993–1997 inflation period. Phys. Earth Planet. Inter. 172, 299–309.

Ebmeier, S., Andrews, B., Araya, M., Arnold, D., Biggs, J., Cooper, C., Cottrell, E., Furtney, M., Hickey, J., Jay, J., et al., 2018. Synthesis of global satellite observations of magmatic and volcanic deformation: implications for volcano monitoring & the lateral extent of magmatic domains. J. Appl. Volcanol. 7, 1–26.

Edmonds, M., Cashman, K.V., Holness, M., Jackson, M., 2019. Architecture and Dynamics of Magma Reservoirs.

Feigl, K.L., Le Mével, H., Tabrez Ali, S., Córdova, L., Andersen, N.L., DeMets, C., Singer, B.S., 2014. Rapid uplift in laguna del Maule volcanic field of the andean southern volcanic zone (Chile) 2007–2012. Geophys. J. Int. 196, 885–901.

Galetto, F., Acocella, V., Hooper, A., Bagnardi, M., 2022. Eruption at basaltic calderas forecast by magma flow rate. Nat. Geosci. 15, 580–584.

Galland, O., Bertelsen, H.S., Guldstrand, F., Girod, L., Joannessen, R.F., Bjugger, F., Burchardt, S., Mair, K., 2016. Application of open-source photogrammetric software micam for monitoring surface deformation in laboratory models. J. Geophys. Res. Solid Earth 121, 2852–2872.

Garibaldi, N., Tikoff, B., Peterson, D., Davis, J., Keranek, K., 2020. Statistical separation of tectonic and inflation-driven components of deformation on silicic reservoirs, laguna del Maule volcanic field, Chile. J. Volcanol. Geotherm. Res. 389, 106744.

Got, J.L., Amritano, A., Marsan, D., Brothelande, E., Peltier, A., 2017. An analysis of the nonlinear magma-edifice coupling at Grimsvötn volcano (Iceland). J. Geophys. Res. Solid Earth 122, 826–843. URL https://doi.org/10.1002/2016JB012905 eprint. https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2016JB012905.

Got, J.L., Amitrano, D., Stefanou, I., Brothelande, E., Peltier, A., 2019. Damage and strain localization around a pressurized shallow-level magma reservoir. J. Geophys. Res. Solid Earth 124, 1443–1458. URL https://doi.org/10.1029/2018JB016407. https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2018JB016407.

Gregg, P., De Silva, S., Grosfils, E., Parmigiani, J., 2012. Catastrophic caldera-forming eruptions: thermomechanics and implications for eruption triggering and maximum caldera dimensions on earth. J. Volcanol. Geotherm. Res. 241, 1–12.

Guldstrand, F., Burchardt, S., Hallot, E., Galland, O., 2017. Dynamics of surface deformation induced by dikes and cone sheets in a cohesive coulomb brittle crust. J. Geophys. Res. Solid Earth 122, 8511–8524.

Heap, M., Faulkner, D., Meredith, P., Vinciguerra, S., 2010. Elastic moduli evolution and accompanying stress changes with increasing crack damage: implications for stress changes around fault zones and volcanoes during deformation. Geophys. J. Int. 183, 225–236.

Hickey, J., Gottsmann, J., Nakamichi, H., Iguchi, M., 2016. Thermomechanical controls on magma supply and volcanic deformation: application to aira caldera, Japan. Sci. Rep. 6, 32691.

Katz, R.F., Spiegelman, M., Holtzman, B., 2006. The dynamics of melt and shear localization in partially molten aggregates. Nature 442, 676–679.

Kavanagh, J.L., Menand, T., Sparks, R.S.J., 2006. An experimental investigation of sill formation and propagation in layered elastic media. Earth Planet. Sci. Lett. 245, 799–813.

Kavanagh, J., Menand, T., Daniels, K.A., 2013. Gelatine as a crustal analogue: determining elastic properties for modelling magmatic intrusions. Tectonophysics 582, 101–111.

Kavanagh, J., Boutelier, D., Cruden, A., 2015. The mechanics of sill inception, propagation and growth: Experimental evidence for rapid reduction in magmatic overpressure. Earth Planet. Sci. Lett. 421, 117–128.

Kavanagh, J., Rogers, B., Boutelier, D., Cruden, A., 2017. Controls on sill and dyke-sill hybrid geometry and propagation in the crust: the role of fracture toughness. Tectonophysics 698, 109–120.

Kavanagh, J.L., Engwell, S.L., Martin, S.A., 2018. A review of laboratory and numerical modelling in volcanology. Solid Earth 9, 531–571.

Lazecký, M., Spaans, K., González, P.J., Maghsoudi, Y., Morishita, Y., Albino, F., Elliott, J., Greenall, N., Hatton, E., Hooper, A., et al., 2020. Licsar: an automatic Insar tool for measuring and monitoring tectonic and volcanic activity. Remote Sens. 12, 2430.

Le Mével, H., Feigl, K.L., Córdova, L., DeMets, C., Lundgren, P., 2015. Evolution of unrest at laguna del Maule volcanic field (Chile) from Insar and gps measurements, 2003 to 2014. Geophys. Res. Lett. 42, 6590–6598.

Le Mével, H., Córdova, L., Cardona, C., Feigl, K., 2021. Unrest at the laguna del Maule volcanic field 2005–2020: renewed acceleration of deformation. Bull. Volcanol. 83, 39.

Lun, R., Zhao, W., 2015. A survey of applications and human motion recognition with microsoft kinect. Int. J. Pattern Recognit. Artif. Intell. 29, 1555008.

McLeod, P., Tait, S., 1999. The growth of dykes from magma chambers. J. Volcanol. Geotherm. Res. 92, 231–245.

Menand, T., Tait, S.R., 2001. A phenomenological model for precursor volcanic eruptions. Nature 411, 678–680.

Merle, O., Vendeville, B., 1995. Experimental modelling of thin-skinned shortening around magmatic intrusions. Bull. Volcanol. 57, 33–43.

Mogi, K., 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Earthq Res Inst 36, 99–134.

Morishita, Y., Lazecký, M., Wright, T.J., Weiss, J.R., Elliott, J.R., Hooper, A., 2020. Licbas: an open-source Insar time series analysis package integrated with the licsar automated sentinel-1 Insar processor. Remote Sens. 12, 424.

Novoa, C., Rémy, D., Gerbault, M., Baez, J.C., Tassara, A., Cordova, L., Cardona, C., Granger, M., Bonvalot, S., Delgado, F., 2019. Viscoelastic relaxation: a mechanism to explain the decennial large surface displacements at the laguna del Maule silicic volcanic complex. Earth Planet. Sci. Lett. 521, 46–59.

Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 75, 1135–1154.

Pritchard, M., De Silva, S., Michelfelder, G., Zandt, G., McNutt, S.R., Gottsmann, J., West, M., Blundy, J., Christensen, D., Finnegan, N., et al., 2018. Synthesis: Plutons: investigating the relationship between pluton growth and volcanism in the Central Andes. Geosphere 14, 954–982.

Pritchard, M., Poland, M., Reath, K., Andrews, B., Bagnardi, M., Biggs, J., Carn, S., Coppola, D., Ebmeier, S., Furtney, M., et al., 2022. Optimizing Satellite Resources for the Global Assessment and Mitigation of Volcanic Hazards—Suggestions from the Usgs Powell Center Volcano Remote Sensing Working Group.

Ranalli, G., 1995. Rheology of the Earth. Springer Science & Business Media.

Reverso, T., Vandemeulebrouck, J., Jouanne, F., Pinel, V., Villemin, T., Sturkell, E., Bascou, P., 2014. A two-magma chamber model as a source of deformation at Grímsvötn Volcano, Iceland. J. Geophys. Res. Solid Earth 119, 4666–4683. URL http://onlinelibrary.wiley.com/gaelnomade-1.grenet.fr/doi/10.1002/2013JB010569 abstract. https://doi.org/10.1002/2013JB010569.

Rincón, A., Márquez, A., Herrera, R., Galland, O., Sánchez-Oro, J., Concha, D., Montemayor, A., 2022. Monitoring volcanic and tectonic sandbox analogue models using the kinect v2 sensor. Earth Space Sci. 9 e2020EA001368.

Roche, O., Carazzo, G., 2019. The contribution of experimental volcanology to the study of the physics of eruptive processes, and related scaling issues: a review. J. Volcanol. Geotherm. Res. 384, 103–150.

Sandri, L., Acocella, V., Newhall, C., 2017. Searching for patterns in caldera unrest. Geochem. Geophys. Geosyst. 18, 2748–2768.

Schmitt, A., Sliwinski, J., Caricchi, L., Bachmann, O., Riel, N., Kaus, B.J., de Léon, A.C., Cornet, J., Friedrichs, B., Lovera, O., et al., 2023. Zircon age spectra to quantify magma evolution. Geosphere 19, 1006–1031.

Segall, P., 2016. Repressurization following eruption from a magma chamber with a viscoelastic aureole. J. Geophys. Res. Solid Earth 121, 8501–8522.

Singer, B.S., Andersen, N.L., Le Mével, H., Feigl, K.L., DeMets, C., Tikoff, B., Thurber, C. H., Jicha, B.R., Cardona, C., Córdova, L., et al., 2014. Dynamics of a large, restless, rhyolitic magma system at laguna del Maule, southern Andes, Chile. GSA Today 24, 4–10.

Smittarello, D., Pinel, V., Maccaferri, F., Furst, S., Rivalta, E., Cayol, V., 2021. Characterizing the physical properties of gelatin, a classic analog for the brittle elastic crust, insight from numerical modeling. Tectonophysics 812, 228901.

Sparks, R., Annen, C., Blundy, J., Cashman, K., Rust, A., Jackson, M., 2019. Formation and dynamics of magma reservoirs. Phil. Trans. R. Soc. A 377, 20180019.

Taisne, B., Tait, S., Jaupart, C., 2011. Conditions for the arrest of a vertical propagating dyke. Bull. Volcanol. 73, 191–204.

Tortini, R., Bonali, F.L., Corazzato, C., Carn, S.A., Tibaldi, A., 2014. An innovative application of the kinect in earth sciences: quantifying deformation in analogue modelling of volcanoes. Terra Nova 26, 273–281.

Yang, X.M., Davis, P.M., Dieterich, J.H., 1988. Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing. J. Geophys. Res. Solid Earth 93, 4249–4257.

Yip, S.T.H., Biggs, J., Edmonds, M., Liggins, P., 2024. The role of pre-eruptive gas segregation on co-eruptive deformation and so2 emissions. Earth Planet. Sci. Lett. 626, 118548.