# Dynamics of magmatic intrusion: what can we learn from the comparison of analog and numerical models?

Séverine Furst, Virginie Pinel, Francesco Maccaferri

► To cite this version:

Séverine Furst, Virginie Pinel, Francesco Maccaferri. Dynamics of magmatic intrusion: what can we learn from the comparison of analog and numerical models?. Volcanica, 2024, 7 (1), pp.67-87. 10.30909/vol.07.01.6787 . ird-04995427

HAL Id: ird-04995427

https://ird.hal.science/ird-04995427v1

Submitted on 18 Mar 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamics of magmatic intrusion: what can we learn from the comparison of analog and numerical models?

Séverine Furst *$\alpha, \beta$, Virginie Pinel$\beta$, and Francesco Maccaferri$\gamma$

## ABSTRACT

This study investigates the dynamics of magmatic intrusions based on the joint analysis of analog and numerical models. By injecting different fluids from the bottom of a solidified gelatin block, we simulate the propagation of magmatic intrusions through the crust and record their shapes, trajectories, and velocity as they rise towards the surface. Additionally, we make use of a 2D fluid-filled crack propagation model constrained by our experimental observations. The numerical simulations demonstrate that our viscous fluid-filled crack experiments, conducted with silicon-oil injections, propagate in the same regime as typical basaltic intrusions. The comparison between analog and numerical results allow us to define the domain of validity of the numerical model and its limit of applicability. This study provides new insights into the processes that control the propagation of magmatic intrusions and our ability to reproduce them using analog and numerical models.

## RÉSUMÉ

Cette étude combine modélisation analogique et numérique pour étudier la dynamique des intrusions magmatiques. Des huiles de silicone visqueuses sont injectées pour simuler la forme, la trajectoire et la vitesse des intrusions magmatiques. Les paramètres et les observations des expériences alimentent le modèle 2D d'éléments frontière de propagation dynamique, ce qui permet de caractériser le régime de propagation des injections d'huiles. Nous montrons que ce régime est semblable à celui d'une intrusion basaltique et que notre modèle numérique peut reproduire avec précision le comportement de ces intrusions, lorsque les variations de vitesse sont réduites. La comparaison des résultats analogiques et numériques permet donc d'étendre le domaine de validité du modèle numérique. Notre étude apporte de nouvelles informations sur les mécanismes régissant la propagation des intrusions magmatiques, et met en évidence le potentiel de l'utilisation de modèles analogiques et numériques pour améliorer notre compréhension du processus d'intrusion magmatique.

**KEYWORDS:** Viscous-fluid flow; Analog experiments; Numerical modeling; Magmatic intrusion; Propagation regimes.

## 1 INTRODUCTION

Magma is a complex, multiphase, viscous fluid that propagates through the brittle crust by creating and opening fractures. Magmatic intrusions can lead to volcanic eruptions occurring in new, and sometimes unexpected, locations within a volcanic field. The associated hazards can have severe consequences for over 1 billion people (14.3 % of the global population) living within 100 km of a Holocene volcano in 2015 [Freire et al. 2019]. Therefore, improving our understanding of the dynamics of magmatic intrusions is a significant motivation.

Direct observations of magmatic intrusions (or dykes) are only possible when they solidify and become exposed through erosion. As they travel through the crust, they can be tracked by the induced seismicity, deformation, and gravity changes. While these techniques provide accurate information for relatively shallow intrusions, their accuracy decreases for deeper intrusions. Understanding dyke dynamics therefore relies on the use of analog and numerical models in which the intrusions are generally described as growing fluid-filled fractures governed by the competition of two main dissipative processes, namely the viscous flow and the rock fracturing [Rivalta et al. 2015]. In the past, two main modeling approaches

have developed, each considering one of the two dissipative processes as dominant, thereby addressing two end-member regimes for the propagation of fluid-filled cracks: the viscous- and fracture-dominated regimes. The viscous-dominated regime is well described by the lubrication theory, and analytical solutions were first given by Lister [1990] and Spence and Turcotte [1990]. It refers to the propagation of a viscous fluid-filled crack connected to a pressurized reservoir, whose dynamics are driven by the reservoir's excess pressure and the buoyancy of the fluid, neglecting the fracturing process at the tip. The fracture-dominated regime is described by Weertman's theory [Weertman 1971], which considers the propagation of a fluid-filled crack isolated from a reservoir and driven solely by the fluid's buoyancy. However, Weertman's crack theory does not account for the effect of fluid viscosity and, consequently, does not provide information on the crack propagation velocity. When neither of the two primary dissipative processes can be neglected, more complex approaches are necessary, which consider both fracture and viscous forcing [Dontsov and Peirce 2015; Detournay 2016; Lecampion et al. 2018; Zia and Lecampion 2020; Möri and Lecampion 2021; Furst et al. 2023].

So far, analog experiments of dyke propagation have involved various experimental set-ups, fluids with different densities and viscosities (air, water, glycerine, synthetic, and veg-

<footnote>$\alpha$ GEOMAR Helmholtz Centre for Ocean Research Kiel, 24148, Kiel, Germany.</footnote>

<footnote>$\beta$ Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTerre, 38000 Grenoble, France.</footnote>

<footnote>$\gamma$ Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Napoli - Osservatorio Vesuviano, Via Diocleziano 328, 80124, Napoli, Italy.</footnote>

<footnote>* s.furst@geomar.de</footnote>
etable oils), and host-materials that deform either elastically, as solidified gelatin, or plastically, as compacted granular materials [Galland et al. 2018; Kavanagh et al. 2018b]. The choice of the set-up depends on the parameters, conceptual models, and configurations being tested [Rivalta et al. 2015; Galland et al. 2018; Kavanagh et al. 2018b]. Studies have investigated crack geometries and their kinetics under the influence of different density contrasts between fluids and hosts [e.g. Takada 1990; Lister 1991; Taisne and Jaupart 2011], stress fields [e.g. Watanabe et al. 2002; Acocella 2005; Menand et al. 2010; Maccaferri et al. 2019; Pinel et al. 2022], free surface effects [e.g. Rivalta and Dahm 2006], and mechanical layering of host materials [e.g. Rivalta et al. 2005; Kavanagh et al. 2006]. Yet a fundamental parameter for a magmatic intrusion is its propagation velocity. A compelling asset of laboratory experiments using gelatin as host material is its transparency, which allows for the direct observation of the injected fluids and their flow patterns [Kavanagh et al. 2018a; Williams et al. 2022; Pansino et al. 2023], along with crack propagation velocities [Takada 1990; Heimpel and Olson 1994; Rivalta et al. 2005; Rivalta and Dahm 2006; Taisne and Jaupart 2011; Pinel et al. 2022]. Several studies have focused on the velocity of air-filled cracks [e.g. Takada 1990; Heimpel and Olson 1994; Rivalta et al. 2005; Rivalta and Dahm 2006; Taisne and Jaupart 2011], while Watanabe et al. [2002] described velocity variations of oil-filled crack under heterogeneous stress fields. Recently, Maccaferri et al. [2019] and Pinel et al. [2022] carried out numerical simulations of fluid-filled crack propagation experiments using a two-dimensional boundary element model for non-viscous fluids [Maccaferri et al. 2011]. These simulations yielded further understanding of the dominant physical processes observed in the experiments, characterizing the relative influence of buoyancy and a heterogeneous external stress field on the propagation path and velocity.

In most cases, experiments focusing on fluid velocity have used fluids with relatively low viscosity, leading to propagation in the fracture-dominated regime. In contrast, magma viscosities may vary by several orders of magnitude [McLeod and Tait 1999; Takeuchi 2011]. Basaltic magmas typically range from a few tens of Pa·s for 2021 Cumbre Vieja (Spain) eruption [Castro and Feisel 2022], to $10^4$ Pa·s for magma from Mt Etna (Italy) [Harris and Allen 2008], with typical values of $10^2$ Pa·s for Piton de la Fournaise (France) [Grasso and Bachèlery 1995]. Rhyolitic magmas can reach more than $10^4$ Pa·s [Spera 2000; Takeuchi 2011]. However, the propagation regime depends on the balance between the viscous dissipation, which increases with magma viscosity, and the energy required for fracturing the rocks, which increases with fracture toughness. While magma viscosity can be relatively well constrained, effective rock fracture toughness values are more uncertain, and the assessment of the propagation regime of magmatic intrusions may sometimes be very challenging [Rivalta et al. 2015]. Nevertheless, it is generally assumed that in most cases, the propagation regime of magmatic intrusions may not be fully dominated by fracturing processes [e.g. Einarsson and Brandsdottir 1980; Buck et al. 2006; Traversa et al. 2010]. To fill the lack of fluid-filled crack propagation experiments in the viscous-dominated regime, we adopted a multimethod approach

combing a recently developed numerical model [Furst et al. 2023] with new laboratory experiments, using silicon oils characterized by significantly higher viscosities (up to 9.7 Pa·s) compared with those in previous studies. Our primary objective is to analyze the influence of viscosity on the shape, velocity, and trajectory of fluid-filled cracks. Additionally, to illustrate the influence of a heterogeneous stress field, we also considered the stress perturbation induced by surface loading as a case study.

The laboratory experiments were analyzed using the 2D boundary element model of Furst et al. [2023], that simulates the behavior of fluid-filled cracks as observed in our analog experiments, encompassing their shape, trajectory, and velocity. The physical parameters we measured during the analog experiments were employed to set and constrain the numerical simulations. By quantitatively comparing the outcomes of our oil-filled crack experiments with the results from our numerical simulations, we were able to characterize the propagation regimes of our experiments based on the relative contribution of viscous- and fracture-energy dissipation rates. Finally, we discuss these results in the light of parameters associated with typical magmatic intrusions.

## 2 ANALOG EXPERIMENTS OF FLUID-FILLED CRACK PROPAGATION IN A SOLIDIFIED GELATIN BLOCK

Please note that all the variables used through the manuscript are listed, along with their definitions, in the "List of notation" (Table 1).

### 2.1 Experimental set-up

Analog experiments were conducted in the new experimental laboratory at ISTerre in Le Bourget-du-Lac, France. We performed fluid-filled crack propagation experiments, using air and silicon oils injected into a transparent brittle-elastic gelatin block (Figure 1). We used a plexiglas tank with a rectangular base and dimensions $L_t \times l_t \times H_t = 40 \times 20 \times 35$ cm (Figure 1) in which we poured 18 L of salty gelatin (corresponding to a height $H_{gel}$ of 22.5 cm, cf. Figure 1) by dissolving the gelatin powder into salty water at 65 °C (2 wt.% gelatin, 15 wt.% salt). The addition of salt was intended to increase the density of the gelatin and provide greater buoyancy for our oil-filled cracks [Acocella 2005; Kavanagh et al. 2006; Ritter et al. 2013]. The gelatin tank was cooled down in a refrigerator for 20 to 30 hours to reach a temperature of 5 °C. In order to maintain a uniform stiffness of the gelatin surface, a thin layer of oil was applied over the surface before the gelatin cooled and solidified, preventing water evaporation which could otherwise cause the gelatin to strengthen. Throughout the cooling and during the experiments, we used two thermocouples to continuously record the ambient temperature and the temperature inside the gelatin block. It was crucial to keep the gelatin within the range of 5 to 15 °C to ensure its elasticity [Kavanagh et al. 2013; van Otterloo and Cruden 2016]. For each tank, we measured the density ($\rho_{gel}$) and rigidity ($\mu_{gel}$) of the gelatin before starting the experiment, when the temperature $T_{gel}$ reached its lowest value (Table 2). To determine density, we left 100 mL of liquid gelatin to solidify alongside the tank. After solidification the gelatin volume was reduced,

<figure><img src="image_2.png" /></figure>
<table><caption>Table 1: List of notation.</caption><thead><tr><th>Symbol</th><th>Definition</th><th>Unit</th></tr></thead><tbody><tr><td colspan="3"><strong>General parameters and experimental set-up</strong></td></tr><tr><td>g</td><td>gravity acceleration</td><td>m·s<sup>-2</sup></td></tr><tr><td>t</td><td>time</td><td>s</td></tr><tr><td>x, z</td><td>coordinate system</td><td>cm</td></tr><tr><td>L<sub>t</sub> × l<sub>t</sub> × H<sub>t</sub></td><td>tank dimensions</td><td>cm</td></tr><tr><td>L<sub>l</sub> × l<sub>l</sub></td><td>load dimensions</td><td>cm</td></tr><tr><td>m<sub>l</sub></td><td>mass of the applied load</td><td>g</td></tr><tr><td>x<sub>s</sub>, z<sub>s</sub></td><td>position of the crack tip when applying a load</td><td>cm</td></tr><tr><td>d<sub>l</sub></td><td>horizontal distance at the surface between crack tip and center of the load</td><td>cm</td></tr><tr><td>M<sub>tot</sub></td><td>mass of gelatin and water</td><td>g</td></tr><tr><td>ρ<sub>w</sub></td><td>density of water</td><td>kg·m<sup>-3</sup></td></tr><tr><td colspan="3"><strong>Gelatin parameters</strong></td></tr><tr><td>ρ<sub>gel</sub></td><td>density of the gelatin</td><td>kg·m<sup>-3</sup></td></tr><tr><td>T<sub>gel</sub></td><td>temperature during the experiments</td><td>°C</td></tr><tr><td>H<sub>gel</sub></td><td>height</td><td>cm</td></tr><tr><td>M<sub>gel</sub></td><td>mass of gelatin</td><td>g</td></tr><tr><td>V<sub>gel</sub></td><td>volume of solidified gelatin</td><td>mL</td></tr><tr><td>μ<sub>gel</sub></td><td>shear modulus</td><td>Pa</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>-</td></tr><tr><td>K<sub>I</sub></td><td>stress intensity factor</td><td>Pa·m<sup>1/2</sup></td></tr><tr><td>K<sub>c</sub></td><td>fracture toughness</td><td>Pa·m<sup>1/2</sup></td></tr><tr><td>K<sub>c</sub><sup>exp</sup></td><td>2D estimation by air injections</td><td>Pa·m<sup>1/2</sup></td></tr><tr><td>K<sub>c</sub><sup>sim</sup></td><td>numerical estimation based on oil injections</td><td>Pa·m<sup>1/2</sup></td></tr><tr><td>E<sub>f</sub></td><td>fracture energy</td><td>Pa·m</td></tr></tbody></table>

<table><thead><tr><th>Symbol</th><th>Definition</th><th>Unit</th></tr></thead><tbody><tr><td colspan="3"><strong>Fluid crack characteristics</strong></td></tr><tr><td>ρ<sub>f</sub></td><td>density of the fluid</td><td>kg·m<sup>-3</sup></td></tr><tr><td>ρ<sub>air</sub></td><td>density of the air</td><td>kg·m<sup>-3</sup></td></tr><tr><td>Δρ</td><td>= ρ<sub>gel</sub> - ρ<sub>f</sub></td><td>kg·m<sup>-3</sup></td></tr><tr><td>η</td><td>viscosity of the fluid</td><td>Pa·s</td></tr><tr><td>V</td><td>volume of injected fluid</td><td>mL</td></tr><tr><td>v</td><td>crack velocity</td><td>cm·s<sup>-1</sup></td></tr><tr><td>L<sub>head</sub></td><td>crack head length</td><td>cm</td></tr><tr><td>h<sub>max</sub></td><td>maximum crack head opening</td><td>cm</td></tr><tr><td>d</td><td>depth of the transition between the head and the tail of the crack</td><td>cm</td></tr><tr><td>w</td><td>width</td><td>cm</td></tr><tr><td>u</td><td>fluid velocity</td><td>cm·s<sup>-1</sup></td></tr><tr><td>Re</td><td>Reynolds number</td><td>-</td></tr><tr><td>δ</td><td>dip of the crack</td><td>°</td></tr><tr><td>γ</td><td>strike of the crack</td><td>°</td></tr><tr><td>L<sub>c</sub></td><td>critical head length in the fracture-dominated regime</td><td>cm</td></tr><tr><td>L<sup>*</sup></td><td>critical head length in the viscous-dominated regime</td><td>cm</td></tr></tbody></table>

<table><tbody><tr><td colspan="3"><strong>Numerical model parameters and outputs</strong></td></tr><tr><td>N</td><td>number of dislocation elements</td><td>-</td></tr><tr><td>l</td><td>dislocation length</td><td>cm</td></tr><tr><td>A<sub>0</sub></td><td>area of the initial fluid body</td><td>cm<sup>2</sup></td></tr><tr><td>v<sub>i</sub></td><td>input velocity</td><td>cm·s<sup>-1</sup></td></tr><tr><td>v<sub>th</sub></td><td>theoretical velocity from Spence and Turcotte [1990]</td><td>cm·s<sup>-1</sup></td></tr><tr><td>x<sub>i</sub>, z<sub>i</sub></td><td>initial tip position</td><td>cm</td></tr><tr><td>z<sup>*</sup></td><td>normalized depth</td><td>-</td></tr><tr><td>h<sub>tail</sub></td><td>mean tail opening</td><td>mm</td></tr><tr><td>ΔE</td><td>energy release</td><td>Pa·m</td></tr><tr><td>ΔE<sub>v</sub></td><td>viscous dissipation</td><td>Pa·m</td></tr><tr><td>D</td><td>term linking linearly ΔE to v (Equation 3)</td><td>Pa·m·s</td></tr><tr><td>R</td><td>proxy for the regime of propagation</td><td>-</td></tr></tbody></table>

<figure><img src="image_3.png" /></figure>
and characterized by mass $M_{gel}$. By adding water until the gelatin returned to its original volume, we also measured the total weight ($M_{tot}$). Knowing the density of the water ($\rho_w$), we computed the volume of the added water as $(M_{tot} - M_{gel})/\rho_w$, which corresponds to the volume loss during gelatin solidification, providing $V_{gel}$. We consistently obtained very similar density values of $\rho_{gel} \sim 1109 \pm 7 \text{ kg} \cdot \text{m}^{-3}$ for all experiments. The rigidity $\mu_{gel}$ was estimated by measuring the subsidence induced by applying a load at the surface of the gelatin block, following the numerical method from Maccaferri et al. [2019] and Smittarello et al. [2021]. Note that the gelatin adheres to the tank walls, which is taken into account in the numerical model by applying a zero displacement condition to the tank walls. The rigidity depended on the temperature of the gelatin, with lower temperatures resulting in higher rigidity values (Table 2). For the solid gelatin, a Poisson's ratio (v) close to 0.5 was assumed [Kavanagh et al. 2013].

We prepared five gelatin blocks and conducted a series of injections using holes located at the bottom of the tank and spaced 2 cm apart (Figure 1). We performed 16 air injections with different volumes, 7 injections with M1000 silicon oil, and 3 with M10000 silicon oil (Table 2). In the case of air injections, we carried out several injections in the same tank to estimate the gelatin's fracture toughness, which corresponds to the critical value of the stress intensity factor ($K_{I}^{2D} = \Delta\rho g[L_{head}/2]\sqrt{\pi[L_{head}/2]}$) as the crack velocity approaches zero [Smittarello et al. 2021]. For silicon oil injections, no more than three injections were performed in the same tank, in order to increase the total number of fluid-filled fracture propagation experiments. This ensures very similar gelatin conditions for all the experiments carried out within the same tank, which helps comparing thoroughly the outcomes of those experiments. In fact, although our experiments are generally consistent in terms of gelatin's physical properties, variations in the rigidity of different tanks (Table 2) may occur based on different cooling times, subsequently affecting fracture toughness [Kavanagh et al. 2013].

Both M1000 and M10000 silicon oils have the same density, $\rho_f = 970 \text{ kg} \cdot \text{m}^{-3}$ (supplier specifications), but they differ in terms of viscosity: $\eta = 0.97 \text{ Pa} \cdot \text{s}$ and $\eta = 9.7 \text{ Pa} \cdot \text{s}$ at room temperature (~20 °C), respectively. The density of air, $\rho_{air} \sim 1.2 \text{ kg} \cdot \text{m}^{-3}$, is negligible compared to the density of gelatin. Once the fluid is injected in the gelatin, a fluid-filled crack forms and propagates upwards due to buoyancy. Air-filled crack shapes are characterized by a tear-drop profile, pinch-closed at the bottom of the crack (Figure 2A), and their cross-sectional shape can be predicted by Weertman's theory [Weertman 1971]. Oil-filled cracks display a more complex shape, with a concentrated volume of oil at the crack head and an open crack tail maintained by the viscosity of the oil (Figure 2B to D), which can be described by the lubrication theory [Lister 1990; Rubin et al. 1998; Roper and Lister 2007]. The crack exhibits a tear-drop shaped head region where buoyancy forces prevail, and a narrower tail region where viscous forces dominate.

Due to the viscosity of the oil and the relatively small density difference between the oil and gelatin ($\Delta\rho \sim 135 \text{ kg} \cdot \text{m}^{-3}$), the propagation velocity of oil-filled cracks is relatively slow

($\le 0.009 \text{ cm} \cdot \text{s}^{-1}$), which can sometimes result in oil leakage from the bottom of the gelatin tank during the injection process. Such leakage may occur when the upper tip of the developing crack moves too slowly upward, leaving insufficient space for the complete formation of the oil-filled crack within the gelatin block. To mitigate this issue, we adjusted the injection rate of the oil to a slow speed of approximately 0.05–0.2 mL·s$^{-1}$. Furthermore, to enhance the visualization of the shape of the oil-filled cracks, we colored a large volume of oil (>50 mL) by mixing a few drops of green dye to M1000 silicon oil and blue dye to the M10000 oil (Figure 2). We did not notice any significant effect of the small amount of dye added on the viscosity of the oil.

In four experiments performed using the M1000 silicon oil, we placed a loading mass on the surface of the gelatin to induce a heterogeneous stress field within the block. This allowed us to investigate the impact of this stress field on the path and velocity of crack propagation, following the approach of Watanabe et al. [1999] and Maccaferri et al. [2019]. We used a rectangular loading plate with dimensions $L_l \times l_l = 14 \times 6$ cm and a mass $m_l = 30.0 \text{ g}$ (Figure 1B). Additional weights were added to reach a total mass up to 110.4 g. We define a coordinate system with origin at the center of the gelatin surface, a horizontal axis (x) rightwards oriented, and a vertical axis (z) downwards oriented (Figure 1B). During the experiments with the loading plate, once the crack had reached a certain depth $z_s$, the load was added on the gelatin surface at a horizontal distance $d_l$ from the crack tip. In order to be able to process the results, and perform 2D numerical simulations of these experiments, we required for the crack plane to be perpendicular to the x-axis of our reference frame.

Among all the experiments, two were excluded from analysis. The first one (EXP2101-I3) was discarded due to time gaps in the recorded videos, while the second one (EXP2103-I1) was excluded because of the orientation of the crack plane, that did not allow for the application of our 2D numerical model.

<table><caption>Table 2: Physical properties of the gelatin tanks prepared for this study: tank number, concentrations of gelatin and salt, density ($\rho_{gel}$), temperature at the beginning and at the end of the injections ($T_{gel}$), number of injections and rigidity ($\mu_{gel}$).</caption><thead><tr><th>Tank #</th><th>Gelatin (%)</th><th>Salt (%)</th><th>ρ<sub>gel</sub> (kg·m<sup>-3</sup>)</th><th>T<sub>gel</sub> (°C)</th><th>Nb of Inj</th><th>μ<sub>gel</sub> (Pa)</th></tr></thead><tbody><tr><td>2101</td><td>2.0</td><td>15</td><td>1107.4</td><td>8.9–11.4</td><td>5</td><td>188</td></tr><tr><td>2102</td><td>2.0</td><td>15</td><td>1121.4</td><td>6.6–11.1</td><td>5</td><td>318</td></tr><tr><td>2103</td><td>2.0</td><td>15</td><td>1106.4</td><td>8.0–13.5</td><td>5</td><td>195</td></tr><tr><td>2104</td><td>2.0</td><td>15</td><td>1105</td><td>7.6–13.7</td><td>5</td><td>188</td></tr><tr><td>2105</td><td>2.0</td><td>15</td><td>1105</td><td>5.7–15.1</td><td>6</td><td>313</td></tr></tbody></table>

## 2.2 Processing tools and methods

We define the propagation velocity ($v$) of the fluid-filled crack as the velocity at which the upper tip of the crack moves along its path. To measure $v$, we tracked the position of the crack tip using the front-view camera footage (Figure 1B) and analyzed the data using TRACKER software.

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 1: [A] Photography of the experimental set-up with two cameras capturing the injections from the front and the side of the tank. Another camera is hanging from the support out of the field of vision, and monitors the propagation from top. Two lamps light the tank from the left and the back sides through white screens. Temperature inside the gelatin box is monitored by a thermocouple (box hooked on the support to the left above the tank). [B] Scheme of the experimental set-up, location of the three cameras and the 15 injection holes at the bottom of the tank. For some experiments, we added a load at the surface of the gelatin, represented by the blue rectangle. The section below the load shows the induced normalized stress.</figcaption></figure>

TRACKER software was also used to measure the length and dip angle of air-filled cracks during propagation. By tracking the positions of the upper and lower edges of the crack, we were able to determine these parameters, as air-filled cracks close at their bottom end. The lengths and velocities of air-filled cracks were then used to estimate the gelatin's fracture toughness, following the methodology outlined by Smittarello et al. [2021] (see Section 4.4). However, the same procedure could not be applied to estimate the length of the crack head for the oil-filled cracks, as the presence of an open tail behind the crack head made it challenging to identify the lower end of the crack head. To address this issue, we developed an automatic processing procedure using the MATLAB Image Processing toolbox, based on the methodology described by Galetto et al. [2021]. This procedure facilitated the differentiation between the crack head and its tail. The process involved comparing time-lapsed photos taken during the propagation of oil-filled cracks, with a reference "background" photo captured before the injection. The reference photo was subtracted from each photo of the propagating oil-filled crack and converted to grayscale. These grayscale images were then binarized using a threshold value, which was chosen in order to return non-zero values in the area of the crack head. Eventually, we used the function regionprops from MATLAB toolbox to obtain the position, orientation, and length of major and minor axes of the ellipses that share the same normalized second central moments with the non-zero region of the image. With this approach, we measured the crack head length $L_{\text{head}}$ and dip angle $\delta$, and estimated the average value and standard deviation of $L_{\text{head}}$ measurements, considering the part of the propagation path where it remains constant (excluding the injection phase and the interaction with the free-surface).

The same procedure was used to measure the width ($w$), opening ($h_{\max}$) and strike ($\gamma$) of the oil-filled cracks using time-lapse photos from the top view camera (Figure 2B–D bottom row). However, it is important to note that due to the refraction of light in the gelatin block, accurate estimates of these parameters were only possible during the final stages of crack propagation, when the thickness of gelatin above the crack is small enough to allow for a clear view of the crack and to minimize the distortions. For vertical propagation, we determined these parameters within the last 0.5 cm of propagation and computed their mean values along with uncertainties. In experiments involving surface loads, we could only evaluate the opening, strike, and width before applying the loading plate at the surface (at the crack tip position $z_s$), when the crack was still propagating along a straight path. To correct for distortions introduced by measurements at different depths, we computed correction factors as the ratios between the crack opening near the surface and the opening at depth $z_s$. Using data from two vertical propagation experiments (EXP2102-I3 and 2105-I4), we obtained correction factors ranging from 1.26 to 1.43 for EXP2102-I3 and from 1.14 to 1.22 for EXP2105-I4, depending on the value of $z_s$. When determining uncertainties associated with these parameters, we considered errors from direct measurements as well as errors from the correction factors. Unfortunately, we were unable to directly measure the openings of the crack tails due to their small thicknesses and to the strike angle of the crack, which were on the order of tenths of millimeters, making them beyond the resolution capabilities of our laboratory set-up.

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Figure 2: Screenshots from front views (top row) and top views (bottom row) monitored during four different injections for air and silicon oils. Crack head lengths ($L_{head}$), widths ($w$) and openings ($h_{max}$) are illustrated for each crack. For silicon oil cracks, an open tail remains while propagating, represented by the *crack tail*. Vertical cracks induced by an injection of [A] of 5 mL of air, [B] 20 mL of M1000 silicon oil (the crack separated into two cracks of different volumes: only the biggest one reached the surface of the gelatin) and [C] 30 mL of M10000 silicon oil. [D] Deflected trajectory of a silicon oil crack, resulting from the injection of 30 mL of M1000 silicon oil, under a surface load (black rectangle on top of gelatin).</figcaption></figure>

# 3 NUMERICAL MODEL FOR THE PROPAGATION OF FLUID-FILLED CRACKS

## 3.1 Numerical model description

We used the 2D boundary element model recently developed by Furst et al. [2023] to simulate the shape and velocity of fluid-filled fracture propagation. This numerical approach combines the fracturing at the crack tip with the transport of a viscous fluid, enabling the computation of the shape, trajectory, and propagation velocity of a fluid-filled crack, that can be compared to the measurements obtained from our analog experiments. The numerical model is particularly well-suited for analyzing the propagation of viscous oil-filled crack. On the other hand, air-filled crack experiments, where fluid viscosity is negligible, can be effectively described by the Weertman crack theory and modeled using a previous version of the 2D boundary element model for fluid-filled crack propagation [Maccaferri et al. 2011; Pinel et al. 2017; Pinel et al. 2022]. This model is applicable for inviscid fluids but does not provide

direct information on the crack propagation velocity. For a comprehensive description of the approach, we refer to Furst et al. [2023], while here we will briefly summarize the key principles of the model necessary for comparison with the analog experiments.

Built with the boundary elements technique in plane strain approximation, this model considers a half-space with a traction-free surface in which a finite volume of fluid can propagate driven by buoyancy. The cross-section of the fluid-filled crack is discretized using *N* dislocation elements of fixed length. By imposing prescribed boundary conditions at the center of each dislocation element, the model solves the “crack problem”. This enables the computation of the opening and slip of each dislocation, given the stress boundary conditions, which include fluid overpressure and shear tractions acting at the element’s center. The overpressure is defined as the excess fluid pressure with respect to the confining stress. The fluid pressure profile accounts for the contribution of three components: 1) the depth-dependent hydrostatic pressure profile, 2)

<figure><img src="image_9.png" /></figure>
the pressure term due to the compressibility of the fluid, and 3) the viscous pressure change induced by a Hagen-Poiseuille flow. The confining pressure is a combination of an isotropic, depth-dependent lithostatic pressure and the normal component of any other possible source of crustal stress (i.e. tectonic or topographic). To simulate the propagation of the fluid-filled crack, a new dislocation element is added at the crack tip, and the model computes the new overpressure profile and crack shape accordingly.

At each propagation step, the model computes the change in potential energy considering the contributions from strain and gravitational energies. This computation establishes the conditions for crack propagation or arrest, as well as constraints on the crack propagation velocity. If the amount of potential energy released by the system (ΔE), is sufficient to fracture the medium, crack growth is allowed [Dahm 2000], and the propagation procedure is iterated. Mathematically, this means that ΔE > E<sub>f</sub> ⋅ l, where

$$E_f = K_c^2 \frac{1 - v}{2\mu} \qquad (1)$$

is the fracture energy, l is the crack lengthening (i.e. the length of a dislocation element), μ is the shear modulus, and K<sub>c</sub> is the fracture toughness of the medium. By testing different directions of propagation for each step, the model computes the path followed by the fluid-filled crack, always selecting the direction that maximizes the potential energy release [Maccaferri et al. 2011].

The propagation condition implies an excess of potential energy release with respect to the energy spent to fracture the embedding medium (Equation 1). To maintain energy conservation, the difference between ΔE and E<sub>f</sub> ⋅ l must be dissipated through dynamic processes. Considering viscous flow through the crack walls as the dominant dissipative process [Furst et al. 2023], the sum of fracture energy and energy dissipated by viscous flow (ΔE<sub>v</sub>) during a propagation step must balance the energy release:

$$\Delta E = E_f \cdot l + \Delta E_v \qquad (2)$$

where ΔE<sub>v</sub> can be described as the product of the velocity v and a function D that depends on the crack geometry, magma inflow, and viscosity. This function can be numerically estimated at each propagation step [Furst et al. 2023]: ΔE<sub>v</sub> = v ⋅ D. Therefore, the crack propagation velocity v is constrained by the amount of energy produced by the system, such that

$$v = \frac{\Delta E - E_f \cdot l}{D} \qquad (3)$$

Note that there is a dependency between the crack propagation velocity and the fracture energy: increasing E<sub>f</sub> results in a smaller value of v.

The model has been designed to account for a heterogeneous background stress field interacting with the fluid-filled crack. The background stress can be provided as an input file, containing the stress tensor components on a grid over the model domain [Maccaferri et al. 2019]. For the purpose of this study, we have also implemented the possibility to change

the stress field during the crack propagation, allowing for the simulation of experiments with a surface load placed at surface during the propagation of the oil-filled crack. From a numerical perspective, this means that the computation of the crack opening only considers an isotropic (depth-dependent) lithostatic stress when the crack tip is deeper than z<sub>s</sub>, while for 0 < z < z<sub>s</sub>, the crack shape is solved considering the heterogeneous stress field due to the loading plate at surface.

The numerical model for fluid-filled crack propagation requires several inputs: the mechanical properties of the fluid and the embedding medium, the initial position of the crack tip (x<sub>i</sub>, z<sub>i</sub>), an initial velocity guess (v<sub>i</sub>), the cross-sectional area of the fluid-filled crack at a reference pressure (A<sub>0</sub>), and the fracture energy (E<sub>f</sub>). Note that since the model is in plane strain approximation, i.e. extending to infinity in the out of plane direction, A<sub>0</sub> may be used to provide an estimate of the 3D volume of the crack only assuming a finite crack width. As outputs, the model provides the propagation path followed by the fluid-filled crack, its propagation velocity, the crack shape, and the displacements and stress fields induced in the surroundings.

## 3.2 Numerical model set-up

To conduct numerical simulations of our oil-filled crack experiments, we need to input all the necessary parameters into the boundary element model for fluid-filled crack propagation. While we have precise measurements for most of the gelatin and silicon oil parameters (ρ<sub>gel</sub>, μ<sub>gel</sub>, ρ<sub>f</sub>, η, Tables 2 and 3), the fracture energy of the gelatin (E<sub>f</sub>) and the cross-sectional area of the crack (A<sub>0</sub>) are not well constrained [Smittarello et al. 2021]. Estimating the crack area accurately is challenging as the cracks are thin and their orientation is not always perpendicular to the camera view. Therefore, we varied A<sub>0</sub> and E<sub>f</sub> to identify values that can simultaneously provide the best fit for the observed propagation velocity and for the crack length or maximum thickness, which are the most reliable geometric parameters. Additionally, considering that the measurements of crack length and opening come with different uncertainties, we further examined the impact of these uncertainties on the estimation of fracture energy to determine which geometric measurement provides the most reliable constraint for our simulations (Section 3.3).

The boundary element model also requires the initial guess for the crack propagation velocity (v<sub>i</sub>) and the initial depth of the fluid-filled crack tip (z<sub>i</sub>). These are necessary for the model to start the numerical computation of the initial crack shape. Although the propagation velocity guess does not affect the results of the simulation, it does impact the numerical performance of the model during the initial propagation steps [Furst et al. 2023]. For our simulations, we used the average oil-filled crack propagation velocity measured during the experiments as the initial velocity guess, which yielded satisfactory results for all the simulations. Since the simulations do not consider the fluid injection phase and start with the fluid volume already in place within the crack, we set the initial depth of the crack tip (z<sub>i</sub> in Table 4) to the depth reached by the tip of the oil-filled crack after the oil injection phase was completed and the crack was fully developed.

<figure><img src="image_10.png" /></figure>
Lastly, to simulate the experiments with a loading plate applied at the surface of the gelatin, we computed the stress field induced by the loading using a finite element modeling tool (COMSOL multiphysics). The model accounts for the effect of a rigid loading plate applied to the surface of a gelatin block confined within a tank with rigid boundaries [Maccaferri et al. 2019; Smittarello et al. 2021, Figure 1B]. The stress field computed with COMSOL was then given as input to the boundary element model to simulate the propagation of the fluid-filled crack in the presence of the loading plate at the surface.

## 3.3 Sensitivity of the numerical model to crack length and opening

We conducted several numerical simulations for a fluid-filled crack ($η = 0.97$ Pa·s, $ρ_f=970$ kg·m⁻³) with a head length of 10 ± 2 cm and a maximum opening of 1.0 ± 0.1 cm, propagating within a gelatin block ($μ_{gel}=300$ Pa, $ρ_{gel}=1120$ kg·m⁻³) at a velocity of 0.0035 ± 0.003 cm·s⁻¹, within a depth range of 12 to 2 cm (red dashed line and light red area in Figure 3). The crack propagation velocity was constrained within a limited depth range due to the initial deceleration during the crack growth phase and its acceleration at shallow depths as it approaches the free surface. This effect is more pronounced at deeper depths for longer crack heads [Rivalta and Dahm 2006; Pinel et al. 2022].

We varied $A_0$ from 1 to 9 cm², resulting in crack head lengths $L_{head}$ ranging from 5.92 cm and 12.32 cm, and crack openings between 0.2 to 1.2 cm. For each simulation, we searched for the fracture energy $E_f$ that would yield the target propagation velocity within the selected depth range (red dashed line in Figure 3A). Within the range of tested crack lengths, the fracture energy varied from 0.097 to 1.30 Pa·m. In Figure 3B and C, we plot the crack lengths and maximum openings as function of the fracture energy. In this example, the range of values for $E_f$ corresponding to a crack length of 10 ± 2 cm would be between 0.35 and 1.22 Pa·m (Figure 3B). Alternatively considering the maximum opening of 1.0 ± 0.1 cm, the numerical simulations provide values of $E_f$ ranging between 0.86 and 1.13 Pa·m (Figure 3C). Hence, the fracture energy would be better constrained by the measurement of the maximum crack opening rather than its length, considering the errors associated with those measurements. We found that, given our measurement errors, the maximum crack opening consistently provides the best constraint for the fracture energy.

The procedure described in this section has been used to determine the uncertainties on $E_f$ (i.e on $K_c$, see Equation 1) that will be provided for each experiment in the following sections.

# 4 RESULTS

In this section, we present the results of the vertical propagation of air- and oil-filled cracks (Section 4.1). We then focus on analyzing the relative contributions of viscous- and fracture-energy (Section 4.2). Subsequently, we present the results from experiments and numerical simulations involving surface loading (Section 4.3). Finally, we provide estimates for the gelatin's fracture toughness obtained from the different

<figure><img src="image_11.png" /><figcaption>Figure 3: Numerical simulations of a viscous fluid-filled crack ($η = 0.97$ Pa·s, $ρ_f=970$ kg·m⁻³) propagating in a gelatin block ($μ_{gel} = 300$ Pa, $ρ_{gel} = 1120$ kg·m⁻³). [A] Evolution of crack velocity along vertical propagation simulated with the numerical model with varying fluid amount $A_0$ and energy threshold $E_f$ (color scale). The red dashed line and the light red area indicate the target velocity and the associated uncertainties respectively. [B] and [C] relation between crack head length and opening respectively, with energy threshold for the target velocity of 0.0035 cm·s⁻¹ (red dashed line in [A]).</figcaption></figure>

sets of experiments (air- and oil-filled cracks) as well as their corresponding simulations (Section 4.4). For each experiment, we gather all the parameters measured in the laboratory and extracted from the video records, which are summarized in Table 3. Our measurements also confirm that these cracks exhibit low Reynolds numbers ($Re = ρ_f uL_{head}/η$, with $u$ the fluid velocity, see Furst et al. [2023]), making inertial terms negligible, as assumed in the numerical model. Additionally, Table 4 displays the input and output parameters related to the numerical simulations.

## 4.1 Vertical crack propagation

The injection of a finite volume of air is extremely rapid, resulting in the fast formation of a crack that initially decelerates until it reaches a stable propagation velocity. The air-filled crack then ascends at a constant velocity (grey area in Figure 4A) and accelerates as it approaches the surface, similar to that described by Rivalta and Dahm [2006]. We injected finite volumes of air ranging from 2 to 30 mL, resulting in cracks propagating with mean velocities ranging from 0.07 to 4.34 cm·s⁻¹ (depending on the gelatin's rigidity). Due to their large buoyancy, air-filled cracks required smaller volumes for propagation, compared to silicon oil injections. We used volumes ranging from 15 to 40 mL for the M1000 oil, and from

<figure><img src="image_12.png" /></figure>
<table><caption>Table 3: Parameters of the different injections performed in this study. Tank number, injection number, type of injected fluid, fluid viscosity (η), volume of injected fluid (V), holes used for injections (as displayed in Figure 1), average dip at the end of the injection (δ), strike (γ), density contrast between gelatin and fluid (Δρ), coordinates of the crack tip when putting the load at the surface (x_s, z_s), mass of the load (m_l), horizontal distance between the crack tip and the center of the load (d_l), crack head length (L_head), maximum width (w), maximum opening (h_max) and average velocity of propagation (v) excluding the initial phase of crack growth and the final acceleration due to the free surface. Grey rows correspond to injections performed with a load applied at the surface, all others were performed within homogeneous stress field. Brown, green, and blue rows indicate the injections displayed in Figure 4 for Air, M1000 and M10000 respectively.</caption><thead><tr><th>Tank #</th><th>Inj #</th><th>Fluid</th><th>η Pa-s</th><th>V (mL)</th><th>Hole</th><th>δ (°)</th><th>γ (°)</th><th>Δρ (kg·m<sup>-3</sup>)</th><th>x_s (cm)</th><th>z_s (cm)</th><th>m_l (g)</th><th>d_l (cm)</th><th>L_head (cm)</th><th>w (cm)</th><th>h_max (cm)</th><th>ν (cm·s<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="5">2101</td><td>1</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>2</td><td>5</td><td>4.9</td><td>-</td><td>-</td><td>1107.4</td><td>-</td><td>-</td><td>-</td><td>3.04 ± 0.05</td><td>-</td><td>-</td><td>0.07 ± 0.03</td></tr><tr><td>2</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>4</td><td>1</td><td>4.8</td><td>-</td><td>-</td><td>1107.4</td><td>-</td><td>-</td><td>-</td><td>3.86 ± 0.12</td><td>-</td><td>-</td><td>0.08 ± 0.10</td></tr><tr><td>3a</td><td>M1000</td><td>0.97</td><td>15</td><td>4</td><td>5.2±0.1</td><td>-</td><td>-</td><td>137.4</td><td>-</td><td>-</td><td>-</td><td>5.13 ± 0.35</td><td>-</td><td>-</td><td>0.004 ± 0.001</td></tr><tr><td>3b</td><td>M1000</td><td>0.97</td><td>20</td><td>4</td><td>-</td><td>-</td><td>-</td><td>137.4</td><td>-</td><td>-</td><td>-</td><td>7.23 ± 0.23</td><td>-</td><td>-</td><td>0.004 ± 0.002</td></tr><tr><td>3c</td><td>M1000</td><td>0.97</td><td>27</td><td>4</td><td>-</td><td>-</td><td>-</td><td>137.4</td><td>-</td><td>-</td><td>-</td><td>8.1 ± 0.4</td><td>-</td><td>-</td><td>0.006 ± 0.002</td></tr><tr><td rowspan="5">2102</td><td>4</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>30</td><td>2</td><td>-2.9</td><td>-</td><td>-</td><td>1107.4</td><td>-</td><td>-</td><td>-</td><td>5.76 ± 0.09</td><td>-</td><td>-</td><td>4.353 ± 0.024</td></tr><tr><td>5</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>30</td><td>3</td><td>5.1</td><td>-</td><td>-</td><td>1107.4</td><td>-</td><td>-</td><td>-</td><td>5.76 ± 0.10</td><td>-</td><td>-</td><td>4.34 ± 0.06</td></tr><tr><td rowspan="3">2103</td><td>1</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>2</td><td>5</td><td>8.0</td><td>-</td><td>-</td><td>1121.4</td><td>-</td><td>-</td><td>-</td><td>2.97 ± 0.03</td><td>-</td><td>-</td><td>0.010 ± 0.003</td></tr><tr><td>2</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>5</td><td>1</td><td>1.7</td><td>-</td><td>-</td><td>1121.4</td><td>-</td><td>-</td><td>-</td><td>4.19 ± 0.05</td><td>-</td><td>-</td><td>0.190 ± 0.007</td></tr><tr><td>3</td><td>M1000</td><td>0.97</td><td>30</td><td>4</td><td>5.5±0.3</td><td>21.9 ± 0.1</td><td>-</td><td>151.4</td><td>-</td><td>-</td><td>-</td><td>11.1 ± 1.7</td><td>7.10 ± 0.04</td><td>1.03 ± 0.03</td><td>0.005 ± 0.002</td></tr><tr><td rowspan="5">2104</td><td>4</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>30</td><td>2</td><td>5.6</td><td>-</td><td>-</td><td>1121.4</td><td>-</td><td>-</td><td>-</td><td>6.52 ± 0.04</td><td>-</td><td>-</td><td>1.97 ± 0.02</td></tr><tr><td>5</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>30</td><td>3</td><td>6.6</td><td>-</td><td>-</td><td>1121.4</td><td>-</td><td>-</td><td>-</td><td>6.43 ± 0.09</td><td>-</td><td>-</td><td>1.97 ± 0.02</td></tr><tr><td rowspan="3">2105</td><td>1</td><td>M1000</td><td>0.97</td><td>30</td><td>11</td><td>6.5±0.1</td><td>8.1 ± 0.7</td><td>-</td><td>135</td><td>-5.08</td><td>7.80</td><td>70.2</td><td>9.36</td><td>9.72 ± 0.27</td><td>5.83 ± 0.12</td><td>1.21 ± 0.13</td><td>0.009 ± 0.001</td></tr><tr><td>2</td><td>M1000</td><td>0.97</td><td>40</td><td>10</td><td>6.5±0.4</td><td>40.7 ± 0.2</td><td>-</td><td>135</td><td>-</td><td>-</td><td>-</td><td>9.8 ± 0.7</td><td>6.48 ± 0.03</td><td>0.98 ± 0.03</td><td>0.006 ± 0.001</td></tr><tr><td>3</td><td>M1000</td><td>0.97</td><td>30</td><td>9</td><td>10.1±0.2</td><td>-6.7 ± 0.1</td><td>-</td><td>135</td><td>3.50</td><td>8.04</td><td>30.0</td><td>8.32</td><td>9.5 ± 0.6</td><td>5.23 ± 0.10</td><td>1.16 ± 0.19</td><td>0.007 ± 0.002</td></tr><tr><td rowspan="5">2105</td><td>4</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>3</td><td>1</td><td>5.3</td><td>-</td><td>-</td><td>1106.4</td><td>-</td><td>-</td><td>-</td><td>3.17 ± 0.11</td><td>-</td><td>-</td><td>0.24 ± 0.03</td></tr><tr><td>5</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>5</td><td>5</td><td>0.3</td><td>-</td><td>-</td><td>1106.4</td><td>-</td><td>-</td><td>-</td><td>3.61 ± 0.04</td><td>-</td><td>-</td><td>0.77 ± 0.04</td></tr><tr><td>6</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>5</td><td>13</td><td>10.5</td><td>-</td><td>-</td><td>1106.4</td><td>-</td><td>-</td><td>-</td><td>3.45 ± 0.12</td><td>-</td><td>-</td><td>0.57 ± 0.09</td></tr><tr><td>1</td><td>M1000</td><td>0.97</td><td>30</td><td>11</td><td>6.5±0.1</td><td>8.1 ± 0.7</td><td>-</td><td>135</td><td>-5.08</td><td>7.80</td><td>70.2</td><td>9.36</td><td>9.72 ± 0.27</td><td>5.83 ± 0.12</td><td>1.21 ± 0.13</td><td>0.009 ± 0.001</td></tr><tr><td>2</td><td>M1000</td><td>0.97</td><td>40</td><td>10</td><td>6.5±0.4</td><td>40.7 ± 0.2</td><td>-</td><td>135</td><td>-</td><td>-</td><td>-</td><td>9.8 ± 0.7</td><td>6.48 ± 0.03</td><td>0.98 ± 0.03</td><td>0.006 ± 0.001</td></tr><tr><td rowspan="5">2105</td><td>3</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>3</td><td>3</td><td>-</td><td>-</td><td>-</td><td>1105</td><td>-</td><td>-</td><td>-</td><td>3.36 ± 0.05</td><td>-</td><td>-</td><td>0.210 ± 0.009</td></tr><tr><td>4</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>5</td><td>15</td><td>-</td><td>-</td><td>-</td><td>1105</td><td>-</td><td>-</td><td>-</td><td>3.79 ± 0.03</td><td>-</td><td>-</td><td>0.487 ± 0.015</td></tr><tr><td>5</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>5</td><td>13</td><td>8.0</td><td>-</td><td>-</td><td>1105</td><td>-</td><td>-</td><td>-</td><td>4.66 ± 0.02</td><td>-</td><td>-</td><td>0.354 ± 0.010</td></tr><tr><td>6</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>5</td><td>13</td><td>8.0</td><td>-</td><td>-</td><td>1105</td><td>-</td><td>-</td><td>-</td><td>3.51 ± 0.02</td><td>-</td><td>-</td><td>0.068 ± 0.039</td></tr><tr><td>7</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>7</td><td>13</td><td>-</td><td>-</td><td>-</td><td>1105</td><td>-</td><td>-</td><td>-</td><td>4.11 ± 0.04</td><td>-</td><td>-</td><td>0.28 ± 0.04</td></tr><tr><td rowspan="5">2105</td><td>8</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>3</td><td>3</td><td>-</td><td>-</td><td>-</td><td>1105</td><td>-</td><td>-</td><td>-</td><td>11.3 ± 2.0</td><td>8.22 ± 0.01</td><td>0.99 ± 0.04</td><td>0.006 ± 0.003</td></tr><tr><td>9</td><td>Air</td><td>1.81 × 10<sup>-5</sup></td><td>5</td><td>15</td><td>-</td><td>-</td><td>-</td><td>1105</td><td>-</td><td>-</td><td>-</td><td>10.8 ± 0.9</td><td>4.94 ± 0.15</td><td>0.90 ± 0.09</td><td>0.003 ± 0.001</td></tr><tr><td>10</td><td>M1000</td><td>9.7</td><td>50</td><td>5</td><td>1.0±1.0</td><td>58.9 ± 0.1</td><td>-</td><td>135</td><td>-</td><td>-</td><td>-</td><td>7.0 ± 0.4</td><td>5.36 ± 0.16</td><td>0.75 ± 0.04</td><td>0.002 ± 0.001</td></tr><tr><td>11</td><td>M1000</td><td>9.7</td><td>30</td><td>1</td><td>11.7±0.5</td><td>23.8 ± 1.0</td><td>-</td><td>135</td><td>-</td><td>-</td><td>-</td><td>7.0 ± 0.4</td><td>5.36 ± 0.16</td><td>0.75 ± 0.04</td><td>0.002 ± 0.001</td></tr><tr><td>12</td><td>M1000</td><td>0.97</td><td>20</td><td>2</td><td>0.1±1.4</td><td>38.8 ± 0.9</td><td>-</td><td>135</td><td>-</td><td>-</td><td>-</td><td>7.0 ± 0.4</td><td>5.36 ± 0.16</td><td>0.75 ± 0.04</td><td>0.002 ± 0.001</td></tr></tbody></table>

<figure><img src="image_13.png" /></figure>
<table><caption>Table 4: Parameters used to set up the numerical simulations in order to fit with analog observations: experiment number, gelatin rigidity (μ<sub>gel</sub>), gelatin density (ρ<sub>gel</sub>), fluid density (ρ<sub>f</sub>), viscosity (η), input velocity (v<sub>i</sub>) and depth (z<sub>i</sub>), cross-sectional area (A<sub>0</sub>) and maximum crack head opening (h<sub>max</sub>). Output of the numerical simulations: crack head length (L<sub>head</sub>), fracture energy (E<sub>f</sub>) and fracture toughness (K<sub>c</sub><sup>sim</sup>) which is deduced using Equation 1. Unshaded rows correspond to injections with a homogeneous stress field, whereas shaded rows correspond to injections performed with a load applied at the surface.</caption><thead><tr><th>Exp #</th><th>μ<sub>gel</sub><br>(Pa)</th><th>ρ<sub>gel</sub><br>(kg·m<sup>-3</sup>)</th><th>ρ<sub>f</sub><br>(kg·m<sup>-3</sup>)</th><th>η<br>(Pa·s)</th><th>v<sub>i</sub><br>(cm·s<sup>-1</sup>)</th><th>z<sub>i</sub><br>(cm)</th><th>A<sub>0</sub><br>(cm<sup>2</sup>)</th><th>h<sub>max</sub><br>(cm)</th><th>L<sub>head</sub><br>(cm)</th><th>E<sub>f</sub><br>(Pa·m)</th><th>K<sub>c</sub><sup>sim</sup><br>(Pa·m<sup>1/2</sup>)</th></tr></thead><tbody><tr><td>2102-I3</td><td>317.7</td><td>1121.4</td><td>970</td><td>0.97</td><td>4.80 × 10<sup>-3</sup></td><td>14.5</td><td>7.1–7.9</td><td>1.00–1.06</td><td>11.5–11.2</td><td>1.03–1.15</td><td>36.13–38.8</td></tr><tr><td>2103-I2</td><td>194.7</td><td>1106.4</td><td>970</td><td>0.97</td><td>1.00 × 10<sup>-2</sup></td><td>8.8</td><td>6.9–9.4</td><td>1.12–1.38</td><td>10.0–11.1</td><td>0.91–1.25</td><td>26.55–31.18</td></tr><tr><td>2104-I1</td><td>188</td><td>1105</td><td>970</td><td>0.97</td><td>7.40 × 10<sup>-2</sup></td><td>8.8</td><td>6.5–8.9</td><td>1.09–1.34</td><td>9.8–10.8</td><td>0.85–1.17</td><td>25.22–29.66</td></tr><tr><td>2104-I2</td><td>188</td><td>1105</td><td>970</td><td>9.7</td><td>6.20 × 10<sup>-3</sup></td><td>15.0</td><td>6.5–7.0</td><td>0.95–1.01</td><td>9.8–10.0</td><td>0.65–0.72</td><td>22.11–23.24</td></tr><tr><td>2104-I3</td><td>188</td><td>1105</td><td>970</td><td>0.97</td><td>7.20 × 10<sup>-3</sup></td><td>9.0</td><td>5.5–8.9</td><td>0.97–1.35</td><td>9.2–10.8</td><td>0.73–1.17</td><td>23.16–29.66</td></tr><tr><td>2105-I4</td><td>313</td><td>1105</td><td>970</td><td>9.7</td><td>5.70 × 10<sup>-3</sup></td><td>15.5</td><td>8.3–9.2</td><td>0.95–1.03</td><td>12.6–13.0</td><td>0.84–0.97</td><td>32.47–34.81</td></tr><tr><td>2105-I5</td><td>313</td><td>1105</td><td>970</td><td>9.7</td><td>3.20 × 10<sup>-3</sup></td><td>15.5</td><td>6.2–8.1</td><td>0.81–0.99</td><td>11.4–12.4</td><td>0.68–0.93</td><td>29.24–34.16</td></tr><tr><td>2105-I6</td><td>313</td><td>1105</td><td>970</td><td>0.97</td><td>2.20 × 10<sup>-3</sup></td><td>15.6</td><td>4.5–5.3</td><td>0.71–0.79</td><td>10.2–10.8</td><td>0.58–0.69</td><td>26.90–29.31</td></tr></tbody></table>

30 to 50 mL for the M10000 oil to investigate the propagation velocity. While the velocity behavior of oil-filled cracks is similar to that of air-filled cracks (with an initial deceleration, a constant velocity phase, and a final acceleration, see Figure 4B and 4C), the velocity of oil-filled cracks is almost three orders of magnitude slower, and ranges from 0.002 to 0.009 cm·s<sup>-1</sup> (Table 3).

Next, we focus on comparing the results obtained from the experiments with those obtained from the numerical simulations. Figure 5 presents the temporal evolution of different parameters for M1000 (Figure 5; panels A and B) and M10000 injections (Figure 5; panels C, D, and E). The approach described in Section 3.3 allowed us to determine ranges of uncertainties for the different parameters presented in this section. The resulting range of acceptable values for each parameter is summarized in Table 4 and displayed as dashed lines in panels from columns 2 to 5 of Figure 5.

Crack shapes are represented at regular depth intervals (Figure 5; A1 to E1), starting from a stationary shape (dark shape) and propagating with a dynamic shape characterized by a thin open channel developing beneath the tear-drop region at the crack head, which is due to the fluid's viscosity. To be able to qualitatively compare crack shape parameters, we need to objectively characterize the boundary between the crack head and tail: we define it as the depth *d* at which the second derivative of the crack opening—with respect to depth—is maximum (red dashes in Figure 5; A1 to E1). Above this depth, the crack head maintains a relatively constant length during the propagation (Figure 5; A2 to E2). Although we did not use the crack length as a constraint in our model, the obtained lengths fall within the range of measured lengths (pink area in Figure 5; A2 to E2), except for EXP2105-I6. In this case, the crack initially propagated as one main fracture but then a small portion of the oil volume split from the main crack (Figure 2B). As a result, the aspect ratio of the 3D crack is biased by the splitting of the crack compared with a crack that remains in one part throughout its entire propagation.

The opening of the crack head was used as a fitting parameter in the numerical simulations to constrain the crack surface A<sub>0</sub> (cf. Section 3.2, thin pink line on Figure 5; A3 to E3). As stated in Section 2.2, the measurement of the crack opening

was only feasible near the surface (or at depth z<sub>s</sub>, for the loading experiments). Therefore, we represented it as a constant value in the third column of panels of Figure 5. For M10000 silicon oil simulations, we can observe a slight decrease in the maximum opening at the beginning of the propagation that coincides with the formation of a thicker tail (region below *d*), as shown by Figure 5; A4 to E4 [Davis et al. 2023]. Note the difference in the horizontal axis scale between M1000 and M10000 silicon oil.

The last column of Figure 5 compares the simulated velocity (colored circles) with the observed velocity (pink dots). From the analog experiments, the velocity of silicon oil cracks propagating in homogeneous stress field can be described by three stages (pink dots in Figure 5; A5 to E5): initiation, propagation, and acceleration, as previously described in case of air crack injections [e.g. Takada 1990; Heimpel and Olson 1994; Rivalta and Dahm 2006]. The acceleration observed in the final stage, when approaching the surface, is due to the elastic free-surface effect, which has been well documented for air injections [e.g. Rivalta et al. 2005; Rivalta and Dahm 2006; Pinel et al. 2022]. The evolution of the simulated velocity captures the early stages of crack propagation, with a deceleration due to tail formation, followed by a more constant propagation velocity. Eventually a rapid acceleration is simulated, which occurs deeper and with greater amplitudes than during analog experiments. For a given fluid viscosity, the discrepancy between the observation and the simulation during this final stage diminishes when considering smaller cracks (see differences in Figure 5; A5 and B5 or C5 and D5), similar to what is observed for air cracks [Pinel et al. 2022].

We can also compare our simulated velocity to the estimate provided by the lubrication theory, using the analytical formula for the velocity of a 2D batch of fluid propagating with a constant area (no compressibility effect) in an infinite medium as given by Spence and Turcotte [1990]:

$$v_{th} = \left( \frac{A_0^2 \Delta \rho g}{48 \eta t^2} \right)^{1/3} \quad (4)$$

where *t* is the time and *g* the gravity acceleration. Using the time of propagation associated with the simulated velocity, we estimate the analytical velocity for each analog experiment

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>

Figure 4: Measured fluid-filled crack velocity (colored dots) for [A] air (EXP2102-I4), [B] M1000 (EXP2102-I3), and [C] M10000 (EXP2105-I5) silicon oil injections. The mean velocity (thick black line) and standard deviation (dashed lines) are estimated for the gray area, i.e. excluding the initial phase of crack growth and the final acceleration due to the free surface. Note the different velocity scales for each experiment.

(blue curves on Figure 5; A5 to E5). The analytical velocity shows a continuous decrease throughout the propagation process, reflecting the progressive transfer of fluid from the (buoyancy-dominated) crack head, to the elongating (viscously-dominated) tail. However, it does not account for the free-surface effect, such that it is more appropriate to compare it with the velocity obtained from the simulations without the free-surface condition (red curves on panels from columns 2 to 5), which indeed show a similar trend. Note that the fitting of both experimental and analytical velocities is more accurate for the most viscous oil (M10000).

## 4.2 Fracture and viscous dissipation energy contributions

For each numerical simulation of the analog experiments, we investigated the relative contributions of fracture energy ($E_f \cdot l$) and viscous dissipation ($\Delta E_v$) to the total energy release ($\Delta E$) along each crack propagation path (cf. Equation 2 and Figure 6). For M1000 silicon oil injections (dashed curves on Figure 6), the relative contribution of $\Delta E_v$ to $\Delta E$ remains low, around 2 %, for most of the propagation. The increase in contribution is influenced by the acceleration due to the free-surface condition, which is significantly larger in the numerical simulations, with respect to the experiments. For the most viscous fluid (plain curves on Figure 6), the relative contribution of $\Delta E_v$ to $\Delta E$ increases up to 20 %. In both cases, the propagation of cracks filled with silicon oil is primarily dominated by the fracturing process. In our previous work [Furst et al. 2023], we validated our numerical model for $\Delta E_v$ contributions ranging from 35 to 75 % (gray area in Figure 6). The lower limit of the gray area corresponds to validation against the analytical model proposed by Rubin et al. [1998], while the upper limit is based on comparison with

the study by Spence and Turcotte [1990]. When the contribution of $\Delta E_v$ exceeds this range (white region to the right of the gray area in Figure 6), we enter a regime where viscosity plays a dominant role. In this regime, the contribution of $E_f \cdot l$ becomes negligible. For this regime Furst et al. [2023] have shown that numerical simulations become unstable due to an underestimation of the viscous pressure drop at the crack tip. The analog experiments performed in this study allow us to further validate the numerical velocity profiles in a regime with a smaller contribution of $\Delta E_v$ (blue area in Figure 6). As shown in the previous section, we achieved a reasonably good agreement between our numerical simulations and the analog experiments performed with the M10000 oil (Figure 5; C1 to E5), which corresponds to a range of viscous energy dissipation of ~10–20 %. When the contribution of $\Delta E_v$ decreases to as low as 2 % in the M1000 experiments, we observed a decline in the ability of our model to accurately reproduce the observed velocity profiles. Despite this limitation, our model still successfully replicates the shapes of the oil-filled fractures and their velocities far away from the free surface (Figure 5; A1 to B5).

Figure 6 also includes the relative contributions of energy for simulations of magmatic dykes. The dotted area represents the range of energy contributions for the three simulations with propagation velocity fitting the migration of volcanotectonic events accompanying the 1998 dyke-intrusion at Piton de la Fournaise [Furst et al. 2023]. Depending on the volume of the intrusion, the model produces dykes with relative contribution of $\Delta E_v$ to $\Delta E$ ranging from 2 % up to 20 %. This indicates that our oil-filled fracture experiments cover a comparable range of energy contributions to that of magmatic intrusions.

<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /><figcaption>Figure 5: (Caption next page.)</figcaption></figure>

<figure><img src="image_18.png" /></figure>
Figure 5: (Previous page.) Results of the simulations for the vertical propagation of silicon oil injections: for M1000 silicon oil (panels from rows A and B) and for M10000 silicon oil (panels from rows C to E). Panel A1 shows the temporal evolution of the crack shape simulated for EXP2102-I3. Initial stationary crack shape (dark blue shape, $t = 0$ s) compared with the static crack shape (dashed orange shape). Representations of dynamic crack shapes sampled at constant space intervals. The color code indicates the time. The crack head length ($L_{head}$), the maximum crack opening ($h_{max}$) and the mean tail opening ($h_{tail}$) are defined on A1 and their evolution along the depth of propagation are displayed on A2, A3, and A4 respectively with associated uncertainties (dashed lines). The pink lines within the pink areas indicate the average values of the crack length and opening measured from the experiments and their uncertainties. The evolution of the propagation velocity $v$ is represented on A5, along with the observed velocity from the analog experiment in pink. The theoretical value of velocity $v_{th}$ defined by Equation 4, is indicated by the blue curve on A5. The red curves on A2, A3, A4, and A5 represent results of simulations removing the effect of free surface. Panels B1 to B5, C1 to C5, D1 to D5, and E1 to E5 display results for EXP2105-I6, EXP2104-I2, EXP2105-I4, and EXP2105-I5 respectively.

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

Figure 6: Evolution of viscous dissipation contribution with respect to the total energy produced by the system. Results for simulations fitting analog experiments of vertical propagation, compared with Piton de la Fournaise simulations from Furst et al. [2023]. $z^*$ is the depth normalized by the initial depth for numerical simulations.

## 4.3 Influence of a topographic load

In the experiments involving surface loading, the oil-filled cracks deflect due to the interaction with the heterogeneous stress field induced by the surface load. It is important to note that the load is applied to the surface of the gelatin after the crack has attained a stable shape and velocity (grey horizontal lines in Figure 7). For the three experiments, the propagation velocity decelerates rapidly after the injection and reaches a constant velocity (pink dots below $z_s$ in Figure 7; A1, B1, and C1). Once the load is applied at the gelatin surface, the crack trajectory starts to bend towards the load, and no major change in the velocity is observed within the range of measurement uncertainties. We initially expected a deceleration due to the bending of the trajectory, as previously observed in air injections [Pinel et al. 2022] or when oil cracks ascend vertically beneath a load [Watanabe et al. 2002]. However, such behavior was not observed here. Instead, we observed an acceleration due to the free surface effect when the crack approaches the surface.

In Figure 7, we also present the results from numerical simulations of these experiments. We conducted both “unconstrained” simulations, where the numerical model solves for both the trajectory and velocity of the oil-filled cracks (shown as blue shapes in Figure 7; A2, B2, and C2), and “constrained” simulations, where the trajectory is imposed as an input and the numerical model solves for the crack velocity only (shown as green shapes in Figure 7; A3, B3, and C3). The parameters characterizing each experiment are listed in Table 4.

For all simulations, we present results obtained with and without accounting for the free surface, as we have shown in the previous section that M1000 simulations provide much larger accelerations due to the free surface, with respect to the experiments. Overall, our model results show a good fit with the trajectories observed in the M1000 experiments, with minimal differences between simulations with or without free surface. However, the model fails to reproduce the observed velocity variations, except for the experiment with a small deflection, and without considering the effect of free surface (Figure 7A4). For this simulation, there is little difference between the velocity and energy release obtained from the “unconstrained” simulation and those obtained by imposing the observed trajectory. Furthermore, the velocity correlates well with the energy release, as previously observed for air injec-

<figure><img src="image_20.png" /></figure>
<figure><img src="image_21.png" /><figcaption>Figure 7: Observations and results from numerical simulations from the less deflected to the most deflected case: EXP2104-I3 (row A), EXP2104-I1 (row B), and EXP2103-I2 (row C). Column 1) observed velocity of propagation: the grey area indicates the propagation within the heterogeneous stress field induced by the surface load. Column 2) initial and final crack shapes (dashed shape and blue filled shape respectively) obtained from the numerical simulation of each experiment in full space (no traction-free surface). The grey-scale indicates the stress field induced by the load. The pink dashed line shows the observed trajectory and the grey rectangle at the surface corresponds to the load. The black transparent shape is the result of the simulation in half-space (with a free surface condition). Column 3) is the same as in column 2 for simulations imposing the trajectory from the experiment. Column 4) simulated velocities over-imposed on the observations (pink dots), in blue for the unconstrained simulation and in green for the ones imposing the trajectory. Plain (dashed) curves are for simulations without (with) the free surface condition. Column 5) energy release for the different simulations.</figcaption></figure>

<figure><img src="image_22.png" /></figure>
tion experiments by Pinel et al. [2022]. For the other two ex-
periments, the simulations exhibit velocity variations that are
not observed in the experiments, and considering a free sur-
face in the simulations appears to worsen the results.

4.4 Fracture toughness comparison

We can consider magmatic dykes in terms of fracture-
dominated and viscous-dominated regimes. In the fracture-
dominated regime, dyke propagation velocity is limited by the
velocity at which the fracture front extends, given the stress in-
tensity factor at the crack tip and the rock fracture toughness.
In the viscous-dominated regime, dyke propagation velocity
is limited by viscous flow, preventing magma from flowing
through the fracture as fast as the fracture front would extend,
if it were filled with an inviscid fluid. Air injections are typ-
ical examples of cracks propagating in a fracture-dominated
regime, where the propagation is limited by the fracture tough-
ness of the material. This latter property provides a quanti-
tative measure of the resistance of brittle material to crack
propagation.

Assuming that viscous forces are fully negligible in air-filled
crack experiments, we can use Weertman's crack theory to
compute the stress intensity factor ($K_I$) as a function of the air-
filled crack length, in 2D [Secor and Pollard 1975] and in 3D
[Heimpel and Olson 1994]. The gelatin's fracture toughness
($K_c$) is the critical value of $K_I$ that allows for crack propaga-
tion. Therefore, we computed $K_I$ for each air injection per-
formed within the same gelatin block, and assuming a linear
dependency between $K_I$ and the crack propagation velocity
[Heimpel and Olson 1994; Smittarello et al. 2021], we extrap-
olated the critical value of $K_I$ when the velocity approaches
zero, which is $K_c$ for that gelatin block.

Experimental values of 2D fracture toughness ($K_c^{exp}$) are displayed for each gelatin tank in Figure 8 using dashed colored lines, with the 3D values being approximately 1.103 times larger [Smittarello et al. 2021]. In addition, we obtained estimates of fracture toughness values from the numerical simulations ($K_c^{sim}$) of the oil-filled crack propagation experiments (vertical and deflected) using Equation 1 (Table 4 and symbols on Figure 8). While $K_c^{exp}$ values for each tank are based on air injections, silicon oil injections yield a single value of $K_c^{sim}$. The values obtained from the numerical simulations are always lower than those inferred from air injections, although they are of the same order of magnitude. Note that the values derived from the air injections would not allow for the propagation of the oil-filled fractures, and that $K_c$ estimates for each gelatin block (in Figure 8 the symbols and lines with the same colors) seem to correlate with the fracture propagation velocity of the experiments used for their estimate. This suggests that effective, possibly velocity-dependent values of fracture toughness should be considered for the gelatin.

5 DISCUSSION

The characterization of oil-filled fractures is more challenging
than that of air-filled fractures, which have been extensively
studied in the past. This is especially true when determin-
ing the fracture head length, as the presence of an open tail
complicates the measurement process (Section 2.2). In fact,

measurements of crack shape may be influenced by various
biases. The relatively strong mismatch between the refrac-
tive indexes of gelatin and air, and gelatin and silicone oil,
may have a significant effect on the measurements of the frac-
ture geometry. To address this issue, following each series
of injections, a sequence of photographs was captured with a
graduated ruler positioned within the gelatin tank, at the loca-
tion of each fluid injection. Subsequently, TRACKER software
was employed to gauge the apparent distance between fixed
markers, and compute a calibration coefficient for the videos.
However, the contrast between the crack and gelatin at the
transition between the crack-head and tail, may not be suf-
ficiently sharp for accurate and consistent estimations of the
crack head lengths using the automated process. Moreover,
the sharpness of photographs or videos may impact the auto-
mated detection of crack velocity and head length. The com-
bined influence of these possible biases results in the relatively
large uncertainties associated with our length measurements.

Due to the oil viscosity, tiny bubbles trapped within the
oil inside the syringe were difficult to eliminate, and some-
times they were injected into an oil-filled crack, where they
slowly gathered at its tip. We estimated that the air volume
(when present) was always less than 2 % of the total crack
volume. Based on five M1000 oil filled crack experiments
with the same volume (EXP2102-I3, EXP2103-I1, EXP2103-I2,
EXP2104-I1, EXP2104-I3), we could make a direct comparison
between the propagation velocity of an oil-filled crack with an
air bubble trapped at its tip (EXP2102-I3), and the velocities of
the other four experiments without air bubbles (before apply-
ing the surface load). All velocities were very similar, except
for EXP2102-I3, which was slightly slower because of a stiffer
gelatin (μgel = 317.7 Pa). This indicates that the effect of the
air bubbles on the crack propagation velocity is negligible in
comparison to other factors, such as gelatin rigidity. In ad-
dition, we observed that the presence of the bubble slightly
sharpens the tip of the crack, but this effect is highly localized
and has an insignificant impact on the overall length of the
crack head and opening.

Other possible sources of uncertainties related to our mea-
surements might have occurred in the case of a crack propa-
gating at the edges of the tank: in such a case the top camera
may have captured not only the opening of the crack, but also
a portion of its body, leading to an overestimation of the open-
ing. Also, the presence of oil residue from previous injections
on the gelatin surface, as well as the dip angle of the crack,
may have affected our measurements. Despite these exper-
imental challenges, we could assess the repeatability of our
experiments through the comparisons of two pairs of exper-
iments: EXP2103-I1 and EXP2104-I1, as well as EXP2103-I2
and EXP2104-I3. These experiments shared key parameters,
including injected volumes, injection hole placements, fluid
types, and gelatin parameters. Even though these experiments
were conducted under the influence of different surface loads,
we could compare crack characteristics prior to the introduc-
tion of the loads. Crack characteristics across these experi-
ments revealed a high level of consistency in terms of veloci-
ties, $h_{max}$ values, and widths (Table 3). However, it is worth
noting that the head lengths displayed a slightly wider range

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /><figcaption>Figure 8: Different fracture toughness estimates as function of the velocity. For each tank, the value of fracture toughness in 2D is deduced from air injections following Smittarello et al. [2021] (dashed colored lines). The symbols with the vertical line represent the results from the numerical simulations with associated uncertainties, by fitting the opening and velocity of each injection. Triangles are for M1000 silicon oil injections and diamonds for M10000.</figcaption></figure>

of values, possibly due to the measurement bias that we have previously discussed.

The vertical propagation of oil-filled cracks exhibits similar velocity behavior to air-filled cracks, characterized by deceleration, stable propagation, and acceleration due to the free surface. However, this acceleration is less pronounced in experiments than in numerical simulations. The larger acceleration obtained with the simulations may be attributed to the 2D plane-strain approximation used, which assumes an infinite crack width and neglects the finite rounded shape of the crack tip. This approximation may result in an underestimation of the oil flux at the crack tip, leading to a decrease in viscous energy dissipation [Furst et al. 2023]. Since the propagation velocity is directly related to the energy dissipation according to Equation 2, any underestimation of viscous dissipation would consequently lead to an overestimation of the propagation velocity.

In our study, we acknowledge the influence of different boundary conditions between the experiments and the numerical model. As our model does not account for the potential effect of the zero-displacement boundary condition at the tank wall on crack opening, we conducted validations using a finite element numerical model. We verified that even in cases where cracks were located at the smallest distance from the tank wall (EXP2105-I4), the stress field induced by the intrusion remained largely unperturbed by the presence of the tank wall. To do this, the stress field induced by the opening of the crack head was calculated by considering a zero-displacement condition at the tank wall using the commercial COMSOL software and by artificially increasing the size of the tank by a factor of two. The maximum compressive stress field induced by the crack opening in both cases is presented in Supplementary Material 1, together with their differences. The closest boundary induces an asymmetry of the

amplitude of $\sigma_1$ induced by the crack opening with larger values on the side towards it, consistently with the plot of the differences between the two simulations showing some significant positive differences at the right of the crack tip, and negative differences at its left. Most importantly, along the crack the differences are very small, except for a small region near its upper tip, which should not significantly influence the energy release calculated over the whole crack and, consequently, the calculated velocity. In addition, preexisting cuts within the gelatin block, due to the propagation of previous fluid-filled cracks, were not considered in our stress models. Preexisting cuts may interact with both the loading stress and the crack-induced stress [Le Corvec et al. 2013]. Therefore, in order to minimize their possible influence on our fluid-filled crack experiments, we always kept a distance of at least 4 cm between oil injections in the same tank [cf. Le Corvec et al. 2013]. Following Maccaferri et al. [2019], we also made sure that no cuts were present between the loading mass at surface and the propagating fluid-filled crack, which could potentially produce a disturbance in the loading stress field at the location of the fluid-filled crack.

When subjected to a surface load, cracks filled with viscous fluid deflect similarly to air-filled cracks, but their velocity does not appear to be affected, although further experiments are needed to confirm this. On the contrary, we obtained significant velocity variations in the numerical simulations of M1000 oil-filled cracks. This may be due to the fact that M1000 silicon oil experiments (and simulations) fall within the fracture-dominated domain (Figure 6). In our numerical model, the computation of the velocity is based on the assumption that the viscous dissipation is the main process controlling the fracture propagation velocity. Therefore, we expect that loading experiments performed with higher-viscosity oil injections may be better reproduced by our numerical model.

<figure><img src="image_25.png" /></figure>
However, as the fracture approaches the fracture-dominated regime, the velocity becomes increasingly influenced by the fracturing process. In order to capture the velocity variations observed in this regime, it may be necessary to consider effective, velocity-dependent values of fracture toughness in the simulations. This hypothesis is also supported by the results obtained for the gelatin's fracture toughness, which seems to correlate with the fracture propagation velocity. An effective fracture toughness could account for possible inelastic processes occurring at the crack tip, such as plastic deformation or damage, which are expected to become more significant as the crack propagation approaches the fracture-dominated regime [Rivalta et al. 2015].

In Section 4.2, we evaluated the relative contributions of viscous dissipation and fracture energy to determine the propagation regime of our experiments and compare it to that of magmatic intrusions. Previously, Roper and Lister [2007] and Rivalta et al. [2015] introduced the ratio $R = L_c/L^*$ between the critical head lengths of a crack in the fracture- and viscous-dominated regime, respectively:

$$L_c = \left( \frac{K_c}{\Delta\rho g} \right)^{2/3} \qquad (5)$$

$$L^* = \left( \frac{\mu_{gel}h_{tail}}{2(1 - \nu)\Delta\rho g} \right)^{1/2} \qquad (6)$$

When $R \ll 1$ the zero-toughness assumption from the lubrication theory is valid, while for $R \gg 1$, the Weertman's theory applies [Equation 22 in Rivalta et al. 2015]. Values of R falling within the range $0.1 \le R \le 10$ indicate a crack propagation regime where both viscous forces and fracture resistance are not fully negligible. Here, we compute R values for our experiments and compare them with the results we obtained for the relative contributions of viscous- and fracture-energy dissipation.

For the air-filled crack experiments, the condition $R \gg 1$ is satisfied because $h_{tail} \to 0$ (Equation 6). However, this is not the case for the silicon oil experiments. We computed the values of R following Equations 5 and 6, using the values directly measured in the laboratory for $\Delta\rho$ (Table 3), and $\mu_{gel}$ (Table 4), and the output from the numerical simulations, for $K_c$ and $h_{tail}$ (Table 5), assuming $\nu = 0.5$ for the solidified gelatin. For the most viscous oil (M10000), we obtained R values ranging between 4.1 and 7.1, while for the less viscous oil (M1000), the values ranged between 9.2 and 14.6. Roper and Lister [2007] have demonstrated that for $R \ge 4$ the head length scales with $L_c$ and the maximum opening of the crack head differs significantly from that of the tail, which aligns with the results from our silicon oil injection, where R is always larger than 4. These values of R are consistent with the relative values of fracture- and viscous-energy contributions computed with the numerical simulations (Table 5 and Figure 6), and they confirm that for the oil injections, the viscous dissipation is significant, even though it is not the dominant force. Injecting higher-viscosity oils would be technically challenging but could be interesting in order to reach fully viscous-dominated regimes.

Estimating crustal fracture toughness remains challenging, yet it is a crucial parameter for modeling magma propagation within the crust. Currently, there are significant discrepancies between laboratory measurements and field estimates of fracture toughness. In laboratory studies, $K_c$ measured for centimetric rock samples is generally of the order of $1 \text{ MPa} \cdot \text{m}^{1/2}$ [Atkinson 1984] to be compared with estimates derived from field observations, which vary from 40 to $4000 \text{ MPa} \cdot \text{m}^{1/2}$ [Rivalta et al. 2015, and references therein]. This discrepancy has been attributed to the dependence of effective fracture toughness on fracture length, resulting from the increasing size of the process zone at the tip with increasing lengths [Delaney et al. 1986]. However, the interpretation of field data can be biased by the model considered [Gill et al. 2022]. From the estimate of $K_c$ obtained by our numerical model applied to Piton de la Fournaise data [Furst et al. 2023], which ranges from 344 to $912 \text{ MPa} \cdot \text{m}^{1/2}$, we found that the propagation regime of basaltic dykes was similar to that experienced by our oil-filled cracks (see Figure 6 and Section 4.2).

The choice of different fluids in these experiments corresponds to magmas with varying viscosities, making them relevant for studying a range of magma types. Although air injections have commonly been used in analog experiments, they are better suited for modeling dykes containing low-viscosity magmas. Given the substantial diversity in magma viscosities found in nature [McLeod and Tait 1999], the experimentation expanded to include alternative fluids when working with gelatin [Kavanagh et al. 2018b]. For instance, silicon oils used by Takada [1990] and Watanabe et al. [2002], provide a more suitable representation for studying higher-viscosity magmas propagating in the form of dykes. In our study, we have demonstrated that viscous silicon oils effectively represent basaltic magmas, while even higher viscosity silicon oils are required to simulate phenocryst-rich andesite or rhyolite magmas [Takeuchi 2011].

Adopting a multimethod approach for investigating magma transfer offers several key benefits, allowing us to advance our understanding and gain a more holistic perspective on magma dynamics. In the context of this work, integrating analog experiments and numerical modeling offered us the opportunity to gain from the complementarity of these approaches: experiments are intrinsically 3D, while our numerical model is 2D; analogue models are very time-consuming and can be performed in a limited number for practical reasons, numerical simulations offer the possibility to increase the number of parameter configurations that we can test; analogue experiments can be used to validate the numerical approach, and evaluate the impact of model assumptions and simplifications; numerical models allow for investigating experimental parameters (or quantities) that would otherwise be difficult to assess (for instance an estimate of the crack-tail thickness, or the viscous energy dissipation for our experimental set-up). By continuously bridging the gap between theory and empirical data through multimethod studies, we can refine and enhance numerical models, making them more accurate and predictive tools for assessing volcanic hazards and informing hazard mitigation strategies.

<figure><img src="image_26.png" /></figure>
<table><caption>Table 5: Parameters to calculate the ratio *R* for vertically propagating oil-filled cracks. Experiment number; fluid type; parameters from the numerical simulations: minimum and maximum average tail opening (*h*<sub>tail</sub>) extracted from Figure 5; optimal fracture toughness (*K*<sub>c</sub>) providing parameter in plain line in Figure 5; optimal crack head length (*L*<sub>c</sub> = *L*<sub>head</sub>); critical head length in the viscous-dominated domain (*L*<sup>*</sup>) estimated from Equation 6; ratio *R*.</caption><thead><tr><th>EXP</th><th>Fluid</th><th><i>h</i><sub>tail</sub><br>(mm)</th><th><i>K</i><sub>c</sub><br>(Pa·m<sup>1/2</sup>)</th><th><i>L</i><sub>c</sub><br>(cm)</th><th><i>L</i><sup>*</sup><br>(cm)</th><th><i>R</i> = <i>L</i><sub>c</sub>/<i>L</i><sup>*</sup></th></tr></thead><tbody><tr><td>EXP2102-I3</td><td>M1000</td><td>0.39–0.19</td><td>37.1</td><td>9.12</td><td>0.96–0.67</td><td>9.5–13.6</td></tr><tr><td>EXP2104-I2</td><td>M10000</td><td>1.8–0.88</td><td>22.6</td><td>6.63</td><td>1.60–1.12</td><td>4.1–5.9</td></tr><tr><td>EXP2105-I4</td><td>M10000</td><td>1.06–0.88</td><td>33.6</td><td>8.63</td><td>1.58–1.44</td><td>5.5–6.0</td></tr><tr><td>EXP2105-I5</td><td>M10000</td><td>1.29–0.58</td><td>31.7</td><td>8.31</td><td>1.75–1.17</td><td>4.8–7.1</td></tr><tr><td>EXP2105-I6</td><td>M1000</td><td>0.30–0.12</td><td>28.6</td><td>7.75</td><td>0.84–0.53</td><td>9.2–14.6</td></tr></tbody></table>

# 6 CONCLUSIONS

We performed a series of experiments involving the injection of air and two types of silicon oils with different viscosities into solidified gelatin blocks. These experiments and comparison with numerical simulations aimed at studying the effect of viscosity on the propagation of fluid-filled cracks under different conditions. We conducted experiments in homogeneous stress conditions, where the crack propagated vertically under lithostatic stress. Additionally, we applied surface loads to induce heterogeneous stress fields, causing the fluid-filled cracks to follow curved trajectories towards the load. We also used a numerical model [Furst et al. 2023] to simulate the oil-filled cracks and infer some parameters that are otherwise difficult to access, such as the thickness of the crack tail, and the gelatin's fracture toughness. We showed that the propagation of a viscous fluid results in the formation of a tail under the crack head, the thickness of which increases with the viscosity. This tail leads to a loss of volume of fluid transported within the crack head during ascent, leading to a deceleration in ascent, that is greater with higher viscosities. We also estimated the gelatin's fracture toughness using an independent method based on the air-filled crack experiments. By comparing the two independent measurements of fracture toughness $K_c$, we identified significant differences, suggesting the presence of effective fracture toughness values that may depend on scale and/or velocity. Furthermore, the numerical simulations provided insights into the relative contributions of viscous and fracture energy dissipation during the crack propagation, which are crucial mechanisms in determining the propagation regime of a fluid-filled fracture, whether it is dominated by viscosity or fracture processes. We showed that the less viscous oil (M1000) exhibits the characteristic shape of fluid-filled fracture described by the lubrication theory, with a thicker tear-drop shaped head followed by a thin tail. However, we also showed that the M1000 oil-filled cracks propagate close to a fracture-dominated regime ($R \sim 10$, and $\Delta E_v \sim 2 \%$), and their velocity (and velocity variations) cannot be accurately described by the lubrication theory. The numerical model of Furst et al. [2023], takes into account both: the contribution of viscous forces (similarly to the lubrication theory), and the brittle fracturing, improving the fit with observed velocities, but only at a significant distance from the free surface, and it is not able to reproduce the velocity variations caused

by loading stress. This suggests that when the propagation is fracture-dominated, more complex rheology—considering velocity-dependent, elasto-plastic effects—might be needed in order to describe these velocity variations in solidified gelatin. Additionally, we demonstrated that even our most viscous intrusions (M10000 silicon oil) do not propagate in a fully viscous dominated regime ($1 < R < 10$, and $\Delta E_v \sim 20 \%$), but rather in an intermediate regime where both viscous forces and fracture resistance are non-negligible. Consequently, the velocity profiles observed in these experiments cannot be fully described by lubrication theory equations, but they can be fairly well described by the simulations performed with the numerical model from Furst et al. [2023].

We also established that the range of viscous energy dissipation of our oil experiments aligns well with that of typical basaltic intrusions, such as those observed at Piton de la Fournaise. Therefore, to investigate higher-viscosity magmatic intrusions, it will be necessary to further increase the relative contribution of viscous energy dissipation in future analog experiments, testing even higher viscosity fluids.

# AUTHOR CONTRIBUTIONS

V.P. and F.M. conceptualized the study and acquired the funding. S.F., F.M. and V.P. developed the numerical methodology. S.F. and V.P. performed the analog experiments. S.F. processed the data from the experiments, ran the numerical models and displayed the figures. All the authors contributed to the writing of the manuscript.

# ACKNOWLEDGEMENTS

We would like to express our gratitude to Editor John Browning, Sam Poppe, and an anonymous reviewer for their valuable and constructive feedback. This work was supported financially by the ANR-DFG NLE 2018 MagmaPropagator project (ANR-18-CE92-0037) and the LOVE-CF project (Progetti Dipartimentali, DIP8, INGV). The figures in this manuscript were created using the Free and Open Source Generic Mapping Tool [Wessel et al. 2019], as well as Adobe Illustrator. We thank Samuel McGowan for his assistance with the analog experiments.

<figure><img src="image_27.png" /></figure>
DATA AVAILABILITY

The Fortran scripts for the dynamic code can be downloaded from https://zenodo.org/record/7118734. The experimental data and numerical simulation examples presented in this study are available on a repository: https://doi.org/10.5281/zenodo.8073619. Supplementary Material is available alongside the online version of this article.

COPYRIGHT NOTICE

© The Author(s) 2024. This article is distributed under the terms of the **Creative Commons Attribution 4.0 International License**, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

REFERENCES

Acocella, V. (2005). "Dike Propagation Driven by Volcano Col lapse: A General Model Tested at Stromboli, Italy". *Geo physical Research Letters* 32(8), page L08308. doi: 10 . 1029/2004GL022248.

Atkinson, B. K. (1984). "Subcritical Crack Growth in Geological Materials". *Journal of Geophysical Research: Solid Earth* 89(B6), pages 4077–4114. doi: 10.1029/JB089iB06p04077.

Buck, W. R., P. Einarsson, and B. Brandsdóttir (2006). “Tec tonic Stress and Magma Chamber Size as Controls on Dike Propagation: Constraints from the 1975-1984 Krafla Rift ing Episode: DRIVING DIKE PROPAGATION”. *Journal of Geophysical Research: Solid Earth* 111(B12), n/a–n/a. doi: 10.1029/2005JB003879.

Castro, J. M. and Y. Feisel (2022). "Eruption of Ultralou Viscosity Basanite Magma at Cumbre Vieja, La Palma, Ca nary Islands". *Nature Communications* 13(1), page 3174. doi: 10.1038/s41467-022-30905-4.

Dahm, T. (2000). "Numerical Simulations of the Propagation Path and the Arrest of Fluid-Filled Fractures in the Earth". *Geophysical Journal International* 141(3), pages 623–638. doi: 10.1046/j.1365-246x.2000.00102.x.

Davis, T., E. Rivalta, D. Smittarello, and R. F. Katz (2023). “As cent Rates of 3-D Fractures Driven by a Finite Batch of Buoyant Fluid”. *Journal of Fluid Mechanics* 954, A12. doi: 10.1017/jfm.2022.986.

Delaney, P. T., D. D. Pollard, J. I. Ziony, and E. H. McKee (1986). "Field Relations between Dikes and Joints: Emplace ment Processes and Paleostress Analysis". *Journal of Geo physical Research: Solid Earth* 91(B5), pages 4920–4938. doi: 10.1029/JB091iB05p04920.

Detournay, E. (2016). "Mechanics of Hydraulic Fractures". An nual Review of Fluid Mechanics 48(1), pages 311–339. doi: 10.1146/annurev-fluid-010814-014736.

Dontsov, E. V. and A. P. Peirce (2015). “A Non-Singular Inte gral Equation Formulation to Analyse Multiscale Behaviour in Semi-Infinite Hydraulic Fractures”. *Journal of Fluid Me chanics* 781, R1. doi: 10.1017/jfm.2015.451.

Einarsson, P. and B. Brandsdóttir (1980). *Seismological Evi dence for Lateral Magma Intrusion during the July 1978 Deflation of the Krafla Volcano in NE-Iceland*. Technical report UI-79-9-7, 890964, UI-79-9-7, 890964. doi: 10.2172/ 890964.

Freire, S., A. Florczyk, M. Pesaresi, and R. Sliuzas (2019). “An Improved Global Analysis of Population Distribution in Proximity to Active Volcanoes, 1975–2015”. *ISPRS Inter national Journal of Geo-Information* 8(8), page 341. doi: 10.3390/ijgi8080341.

Furst, S., F. Maccaferri, and V. Pinel (2023). “Modeling the Shape and Velocity of Magmatic Intrusions, a New Numer ical Approach”. *Journal of Geophysical Research: Solid Earth*. doi: 10.1029/2022JB025697.

Galetto, F., A. Bonaccorso, and V. Acocella (2021). “Relating Dike Geometry and Injection Rate in Analogue Flux-Driven Experiments”. *Frontiers in Earth Science* 9, page 665865. doi: 10.3389/feart.2021.665865.

Galland, O., E. Holohan, B. van Wyk de Vries, and S. Bur chardt (2018). “Laboratory Modelling of Volcano Plumbing Systems: A Review”. *Physical Geology of Shallow Mag matic Systems: Dykes, Sills and Laccoliths*. Edited by C. Breitkreuz and S. Rocchi. Advances in Volcanology. Cham: Springer International Publishing, pages 147–214. ISBN: 978- 3-319-14084-1. doi: 10.1007/11157-2015_9.

Gill, S. P. A., R. J. Walker, K. J. W. McCaffrey, and C. Greenfield (2022). “No Unique Scaling Law for Igneous Dikes”. *Journal of Geophysical Research: Solid Earth* 127(9), e2022JB024120. doi: 10.1029/2022JB024120.

Grasso, J. R. and P. Bachéléry (1995). “Hierarchical Organiza tion as a Diagnostic Approach to Volcano Mechanics: Val idation on Piton de La Fournaise”. *Geophysical Research Letters* 22(21), pages 2897–2900. doi: 10.1029/95GL01786.

Harris, A. J. L. and J. S. Allen (2008). “One-, Two- and Three- Phase Viscosity Treatments for Basaltic Lava Flows”. *Jour nal of Geophysical Research* 113, B09212. doi: 10.1029/ 2007JB005035.

Heimpel, M. and P. Olson (1994). “Chapter 10 Buoyancy- Driven Fracture and Magma Transport through the Litho sphere: Models and Experiments”. *International Geo physics*. Volume 57. Elsevier, pages 223–240. ISBN: 978- 0-12-605070-7. doi: 10.1016/S0074-6142(09)60098-X.

Kavanagh, J. L., T. Menand, and K. A. Daniels (2013). “Gela tine as a Crustal Analogue: Determining Elastic Properties for Modelling Magmatic Intrusions”. *Tectonophysics* 582, pages 101–111.

Kavanagh, J. L., A. J. Burns, S. Hilmi Hazim, E. P. Wood, S. A. Martin, S. Hignett, and D. J. Dennis (2018a). “Challenging Dyke Ascent Models Using Novel Laboratory Experiments: Implications for Reinterpreting Evidence of Magma Ascent and Volcanism”. *Journal of Volcanology and Geothermal Research* 354, pages 87–101. doi: 10.1016/j.jvolgeores. 2018.01.002.

Kavanagh, J. L., S. L. Engwell, and S. A. Martin (2018b). “A Re view of Laboratory and Numerical Modelling in Volcanol ogy”. *Solid Earth* 9(2), pages 531–571. doi: 10.5194/se- 9-531-2018.

<figure><img src="image_28.png" /></figure>
Kavanagh, J. L., T. Menand, and R. S. J. Sparks (2006). "An Experimental Investigation of Sill Formation and Propagation in Layered Elastic Media". *Earth and Planetary Science Letters* 245(3-4), pages 799-813. doi: 10.1016/j.epsl.2006.03.025.

Le Corvec, N., T. Menand, and J. Lindsay (2013). "Interaction of Ascending Magma with Pre-existing Crustal Fractures in Monogenetic Basaltic Volcanism: An Experimental Approach". *Journal of Geophysical Research: Solid Earth* 118(3), pages 968-984. doi: 10.1002/jgrb.50142.

Lecampion, B., A. Bunger, and X. Zhang (2018). "Numerical Methods for Hydraulic Fracture Propagation: A Review of Recent Trends". *Journal of Natural Gas Science and Engineering* 49, pages 66-83. doi: 10.1016/j.jngse.2017.10.012.

Lister, J. R. (1990). "Buoyancy-Driven Fluid Fracture: The Effects of Material Toughness and of Low-Viscosity Precursors". *Journal of Fluid Mechanics* 210, pages 263-280. doi: 10.1017/S0022112090001288.

- (1991). "Steady Solutions for Feeder Dykes in a Density-Stratified Lithosphere". *Earth and Planetary Science Letters* 107(2), pages 233-242. doi: 10.1016/0012-821X(91)90073-Q.

Maccaferri, F., M. Bonafede, and E. Rivalta (2011). "A Quantitative Study of the Mechanisms Governing Dike Propagation, Dike Arrest and Sill Formation". *Journal of Volcanology and Geothermal Research* 208(1-2), pages 39-50. doi: 10.1016/j.jvolgeores.2011.09.001.

Maccaferri, F., D. Smittarello, V. Pinel, and V. Cayol (2019). "On the Propagation Path of Magma-Filled Dikes and Hydrofractures: The Competition Between External Stress, Internal Pressure, and Crack Length". *Geochemistry, Geophysics, Geosystems* 20(4), pages 2064-2081. doi: 10.1029/2018GC007915.

McLeod, P. and S. Tait (1999). "The Growth of Dykes from Magma Chambers". *Journal of Volcanology and Geothermal Research* 92(3), pages 231-245. doi: 10.1016/S0377-0273(99)00053-0.

Menand, T., K. A. Daniels, and P. Benghiat (2010). "Dyke Propagation and Sill Formation in a Compressive Tectonic Environment". *Journal of Geophysical Research* 115(B8), B08201. doi: 10.1029/2009JB006791.

Möri, A. and B. Lecampion (2021). "Arrest of a Radial Hydraulic Fracture upon Shut-in of the Injection". *International Journal of Solids and Structures* 219-220, pages 151-165. doi: 10.1016/j.ijsolstr.2021.02.022.

Pansino, S., A. Emadzadeh, and B. Taisne (2023). "Magma Flow Patterns in Dikes: Observations From Analogue Experiments". *Journal of Geophysical Research: Solid Earth* 128(3), e2022JB025463. doi: 10.1029/2022JB025463.

Pinel, V., A. Carrara, F. Maccaferri, E. Rivalta, and F. Corbi (2017). "A Two-Step Model for Dynamical Dike Propagation in Two Dimensions: Application to the July 2001 Etna Eruption". *Journal of Geophysical Research: Solid Earth* 122(2), pages 1107-1125. doi: 10.1002/2016JB013630.

Pinel, V., S. Furst, F. Maccaferri, and D. Smittarello (2022). "Buoyancy versus Local Stress Field Control on the Velocity of Magma Propagation: Insight from Analog and Numerical Modelling". *Frontiers in Earth Science*, page 32.

Ritter, M. C., V. Acocella, J. Ruch, and S. L. Philipp (2013). "Conditions and Threshold for Magma Transfer in the Layered Upper Crust: Insights from Experimental Models". *Geophysical Research Letters* 40(23), pages 6043-6047. doi: 10.1002/2013GL058199.

Rivalta, E., M. Böttinger, and T. Dahm (2005). "Buoyancy-Driven Fracture Ascent: Experiments in Layered Gelatine". *Journal of Volcanology and Geothermal Research*. The Tectonics and Physics of Volcanoes 144(1), pages 273-285. doi: 10.1016/j.jvolgeores.2004.11.030.

Rivalta, E., B. Taisne, A. Bunger, and R. Katz (2015). "A Review of Mechanical Models of Dike Propagation: Schools of Thought, Results and Future Directions". *Tectonophysics* 638, pages 1-42. doi: 10.1016/j.tecto.2014.10.003.

Rivalta, E. and T. Dahm (2006). "Acceleration of Buoyancy-Driven Fractures and Magmatic Dikes beneath the Free Surface". *Geophysical Journal International* 166(3), pages 1424-1439. doi: 10.1111/j.1365-246X.2006.02962.x.

Roper, S. M. and J. R. Lister (2007). "Buoyancy-Driven Crack Propagation: The Limit of Large Fracture Toughness". *Journal of Fluid Mechanics* 580, pages 359-380. doi: 10.1017/S0022112007005472.

Rubin, A. M., D. Gillard, and J.-L. Got (1998). "A Reinterpretation of Seismicity Associated with the January 1983 Dike Intrusion at Kilauea Volcano, Hawaii". *Journal of Geophysical Research: Solid Earth* 103(B5), pages 10003-10015. doi: 10.1029/97JB03513.

Secor, D. T. and D. D. Pollard (1975). "On the Stability of Open Hydraulic Fractures in the Earth's Crust". *Geophysical Research Letters* 2(11), pages 510-513. doi: 10.1029/GL002i011p00510.

Smittarello, D., V. Pinel, F. Maccaferri, S. Furst, E. Rivalta, and V. Cayol (2021). "Characterizing the Physical Properties of Gelatin, a Classic Analog for the Brittle Elastic Crust, Insight from Numerical Modeling". *Tectonophysics* 812, page 228901. doi: 10.1016/j.tecto.2021.228901.

Spence, D. A. and D. L. Turcotte (1990). "Buoyancy-Driven Magma Fracture: A Mechanism for Ascent through the Lithosphere and the Emplacement of Diamonds". *Journal of Geophysical Research* 95(B4), page 5133. doi: 10.1029/JB095iB04p05133.

Spera, F. (2000). "Physical Properties of Magmas". *Encyclopaedia of Volcanoes*. San Diego : Academic Press, c2000.

Taisne, B. and C. Jaupart (2011). "Magma Expansion and Fragmentation in a Propagating Dyke". *Earth and Planetary Science Letters* 301(1-2), pages 146-152. doi: 10.1016/j.epsl.2010.10.038.

Takada, A. (1990). "Experimental Study on Propagation of Liquid-Filled Crack in Gelatin: Shape and Velocity in Hydrostatic Stress Condition". *Journal of Geophysical Research: Solid Earth* 95(B6), pages 8471-8481. doi: 10.1029/JB095iB06p08471.

Takeuchi, S. (2011). "Preeruptive Magma Viscosity: An Important Measure of Magma Eruptibility". *Journal of Geophysical Research* 116(B10), B10201. doi: 10.1029/2011JB008243.

<figure><img src="image_29.png" /></figure>
Traversa, P., V. Pinel, and J. R. Grasso (2010). "A Constant In- flux Model for Dike Propagation: Implications for Magma Reservoir Dynamics". *Journal of Geophysical Research* 115(B1), B01201. doi: 10.1029/2009JB006559.

van Otterloo, J. and A. R. Cruden (2016). "Rheology of Pig Skin Gelatine: Defining the Elastic Domain and Its Thermal and Mechanical Properties for Geological Analogue Experiment Applications". *Tectonophysics* 683, pages 86-97. doi: 10.1016/j.tecto.2016.06.019.

Watanabe, T., T. Koyaguchi, and T. Seno (1999). "Tectonic Stress Controls on Ascent and Emplacement of Magmas". *Journal of Volcanology and Geothermal Research* 91(1), pages 65-78. doi: 10.1016/S0377-0273(99)00054-2.

Watanabe, T., T. Masuyama, K. Nagaoka, and T. Tahara (2002). "Analog Experiments on Magma-Filled Cracks: Competition between External Stresses and Internal Pressure". *Earth, Planets and Space* 54(12), pages 1247-1261. doi: 10.1186/BF03352453.

Weertman, J. (1971). "Theory of Water-filled Crevasses in Glaciers Applied to Vertical Magma Transport beneath Oceanic Ridges". *Journal of Geophysical Research* 76(5), pages 1171-1183.

Wessel, P., J. F. Luis, L. Uieda, R. Scharroo, F. Wobbe, W. H. F. Smith, and D. Tian (2019). "The Generic Mapping Tools Version 6". *Geochemistry, Geophysics, Geosystems* 20(11), pages 5556-5564. doi: 10.1029/2019GC008515.

Williams, K., J. Kavanagh, and D. Dennis (2022). "Focused Flow during the Formation and Propagation of Sills: Insights from Analogue Experiments". *Earth and Planetary Science Letters* 584, page 117492. doi: 10.1016/j.epsl.2022.117492.

Zia, H. and B. Lecampion (2020). "PyFrac: A Planar 3D Hydraulic Fracture Simulator". *Computer Physics Communications* 255, page 107368. doi: 10.1016/j.cpc.2020.107368.

<figure><img src="image_30.png" /></figure>
Supplementary Material for:
# Dynamics of magmatic intrusion: what can we learn from the comparison of analog and numerical models?

Séverine Furst *$\alpha^\beta$, Virginie Pinel$^\beta$, and Francesco Maccaferri$^\gamma$

INFLUENCE OF BOUNDARY CONDITION AT THE TANK WALL
ON CRACK OPENING

We used the commercial software COMSOL to compute the stress field generated by the opening of the crack head, considering a boundary condition of zero displacement at the tank wall. Specifically, we calculated this stress field using the parameters of EXP2105-I4, where 50 mL of oil was injected into hole 5, located 12 cm from the lateral wall. The calculations were conducted with the crack located 2 cm below the surface, under two conditions: one considering the actual dimensions of the experimental tank and another simulating a tank twice as large.

We modeled the crack opening of a static Weertman crack [Secor and Pollard 1975], matching the observed length, which is always larger than the opening of a dynamic crack head [Furst et al. 2023]. The resulting maximum compressive stress around the crack is oriented perpendicularly to the opening. Figure S1a and S1b illustrate this maximum compressive stress within the vertical plane perpendicular to the crack (Y=0) for the two distinct tank sizes, while Figure S1c displays the stress difference between these cases. The impact of the gelatin's fixed displacement at the tank walls on the internal stress field around the crack is significant only in a small area localized above the crack tip. Hence, we maintain confidence in the limited influence of this effect on crack propagation.

<footnote>$^\alpha$ GEOMAR Helmholtz Centre for Ocean Research Kiel, 24148, Kiel, Germany.</footnote>

<footnote>$^\beta$ Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERRE, 38000 Grenoble, France.</footnote>

<footnote>$^\gamma$ Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Napoli - Osservatorio Vesuviano, Via Diocleziano 328, 80124, Napoli, Italy.</footnote>

<footnote>* s.furst@geomar.de</footnote>
<figure><img src="image_31.png" /><figcaption>Figure S1: Amplitude of the maximum compressive stress (positive value for compression) induced by the opening of a Weertman crack (white shape) representing the largest crack from the analog experiments. The crack tip is located at 2 cm below the gelatin's surface. COMSOL simulations were performed for a gelatin block with dimensions a) $L_t \times l_t \times H_t = 40 \times 20 \times 22.5$ cm ("Normal Tank") and b) $L_t \times l_t \times H_t = 80 \times 40 \times 50$ cm ("Large Tank"). c) Difference between the "Large tank" and the "Normal Tank".</figcaption></figure>

## COPYRIGHT NOTICE

© The Author(s) 2024. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## REFERENCES

Furst, S., F. Maccaferri, and V. Pinel (2023). "Modeling the Shape and Velocity of Magmatic Intrusions, a New Numerical Approach". In: Journal of Geophysical Research: Solid Earth. DOI: 10.1029/2022JB025697.

Secor, D. T. and D. D. Pollard (1975). "On the Stability of Open Hydraulic Fractures in the Earth's Crust". In: *Geophysical Research Letters* 2(11), pages 510-513. DOI: 10.1029/GL002i011p00510.

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