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

HAL
open science

# Modeling the Shape and Velocity of Magmatic Intrusions, a New Numerical Approach

S. Furst, F. Maccaferri, Virginie Pinel

► To cite this version:

S. Furst, F. Maccaferri, Virginie Pinel. Modeling the Shape and Velocity of Magmatic Intrusions, a New Numerical Approach. Journal of Geophysical Research: Solid Earth, 2023, 128 (3), 10.1029/2022jb025697. ird-04995434

HAL Id: ird-04995434

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

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

Distributed under a Creative Commons Attribution 4.0 International License
JGR Solid Earth

RESEARCH ARTICLE

10.1029/2022JB025697

**Key Points:**

* We present a new modeling scheme to compute the shape and velocity of a growing fluid-filled crack

* Our magmatic dykes show a tear drop head and open tail, on a wide range of propagation velocities

* We reproduce the velocity and fit important parameters for the 1998 Piton de la Fournaise intrusion

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

F. Maccaferri,
francesco.maccaferri@ingv.it

**Citation:**

Furst, S., Maccaferri, F., & Pinel, V. (2023). Modeling the shape and velocity of magmatic intrusions, a new numerical approach. *Journal of Geophysical Research: Solid Earth*, 128, e2022JB025697. https://doi.org/10.1029/2022JB025697

Received 28 SEP 2022
Accepted 27 FEB 2023

**Author Contributions:**

**Conceptualization:** F. Maccaferri, V. Pinel

**Funding acquisition:** F. Maccaferri, V. Pinel

**Investigation:** S. Furst

**Methodology:** S. Furst, F. Maccaferri

**Project Administration:** V. Pinel

**Supervision:** F. Maccaferri, V. Pinel

**Validation:** S. Furst, F. Maccaferri

**Visualization:** S. Furst

**Writing – original draft:** S. Furst, F. Maccaferri, V. Pinel

**Writing – review & editing:** S. Furst, F. Maccaferri

© 2023. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

Modeling the Shape and Velocity of Magmatic Intrusions, a
New Numerical Approach

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

S. Furst¹,² ID, F. Maccaferri³ ID, and V. Pinel¹ ID

¹University Grenoble Alpes, University Savoie Mont Blanc, CNRS, IRD, University Gustave Eiffel, Grenoble, France,

²GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany, ³Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Napoli - Osservatorio Vesuviano, Napoli, Italy

**Abstract** Dykes are magma-filled fractures propagating through the brittle crust. Understanding the physics of dyking process is essential to mitigate the volcanic hazard associated with the opening of new eruptive fissures at the surface. Often, physics-based models view either fracturing of the host rock or viscous flow of the magma as the dominating energy sink during dyke propagation. Here, we provide a numerical model that captures the coupling of fracturing at the crack tip and the transport of a viscous fluid. Built with the boundary element technique, our model allows for computation of the shape and velocity of a growing fluid-filled crack accounting for the viscosity of the fluid: the fluid flow induces a viscous pressure drop acting at the crack walls, and modifies the shape of the crack. The energy conservation equation provides the constraints to solve for the crack growth velocity, assuming that brittle fracturing and viscous flow are the main processes that dissipate energy. Using a parameter range that represents typical magmatic intrusions, we obtain crack shapes displaying some typical characteristics, including a tear-drop head and an open tail that depend on rock rigidity, magma viscosity, and buoyancy. We show that viscous forces significantly contribute to the energy dissipated during the propagation of magmatic dykes. Applied to the 1998 intrusion at Piton de la Fournaise (La Réunion Island), we provide ranges of dyke lengths and openings by adjusting the numerical velocity to the one deduced from the migration of volcano-tectonic events.

**Plain Language Summary** Magma is a viscous fluid that can propagate through the crust by fracturing rocks and flowing through them. These magma-filled fractures are called dykes. Magma pressure is the force ensuring the opening of the fracture and maintaining the flow of magma. Although physics-based models provide simplified but reliable representations of dykes, few consider both fracture creation and viscous flow simultaneously. Here, we present a new model for magmatic dyke propagation, implementing viscous flow equations into an existing model, such that the dyke shape, trajectory, and velocity are determined together. We show that the viscous dissipation within the magma is a major contribution to the energy balance for a range of viscosity and velocity values which are typical for magmatic intrusions. The dyke shapes obtained with our model compare well with the ones obtained with previous models. The velocity derived from the model has been compared with that derived from the spatio-temporal evolution of seismic events recorded before the 1998 eruption at Piton de la Fournaise. Our new modeling scheme may represent a step forward in predicting the timing and location of future eruptions at monitored volcanoes.

1. Introduction

The journey of magma through the crust may be complex, with dykes following complex paths, getting trapped before erupting (Pedersen et al., 2007; Roult et al., 2012) or getting only temporarily arrested before resuming their propagation with the supply of additional magma (Sigmundsson et al., 2015). When a dyke reaches the surface, it may open fissures on the volcano flank or at unspecified locations in a volcanic field. The geometry, velocity, and trajectory followed by a magmatic intrusion can be tracked by the migration of volcano-tectonic earthquakes (Battaglia et al., 2005; Duputel et al., 2019; Lengliné et al., 2016, 2021) and geodetic observations including InSAR and GNSS data (e.g., Beauducel et al., 2020; Davis et al., 2021; Maccaferri et al., 2016; Peltier et al., 2005; Smittarello et al., 2019). In particular, the velocity of propagation has been shown to be unsteady (Peltier et al., 2007; Sigmundsson et al., 2015), usually lasting for a few hours to a few days, before potentially reaching the surface. Estimating where and when magma may reach the surface is crucial to assess the volcanic risk.
Magmatic dykes can be seen as a sub-set of the larger family of fluid-filled fractures. The physics describing fluid-filled fracture growth may be rather simple to describe, but is extremely challenging to implement in an analytical (Lister, 1990; Roper & Lister, 2007; Rubin, 1998; Spence & Turcotte, 1990), and even in a numerical modeling scheme (Detournay, 2016; Dontsov & Peirce, 2015; Lecampion et al., 2018). The fracturing process is described by equations for a crack embedded into a brittle-elastic medium, while the internal flow is described by fluid dynamics equations. The pressure profile within the fluid-filled crack, the crack shape, and the velocity of crack growth, result from the solution of a combined set of these equations. The solution of the coupled elastic and fluid-dynamic problem, even in 2D and under the most simplifying assumptions, is far from trivial. In the past, two main schools have been established for modeling magmatic dykes: they have been named the “Weertman school” and the “lubrication theory school”, according to the fracture propagation regime that they target, fracture- or viscous-dominated, respectively (Rivalta et al., 2015). The “Weertman school”, after Weertman (1971), makes use of static crack theory and quasi-static approaches, disregarding the fluid motion. Such models may provide, at best, only indirect information about the crack propagation velocity, as recently shown by Pinel et al. (2022), and they can be used when the energy dissipated by viscous flow is negligible with respect to the total budget of the problem. The main advantage of this approach is that by using the Boundary Element (BE) technique they can describe complex crack shapes, account for the crack interaction with heterogeneous stress fields and crustal heterogeneities, and compute the direction of crack growth, which is not necessarily planar (Dahm, 2000a; Davis et al., 2020; Heimisson et al., 2015; Maccaferri et al., 2010, 2011). These models are typically 2D, while some compelling new approaches also describe the 3D non-planar fracture-process (Davis et al., 2021).

On the other hand, the “lubrication theory school”, simplifies the crack geometry and the crustal stress and structures, but can account for the interaction between elastic forces acting at the crack walls and viscous forces due to the fluid motion (Lister, 1990; Spence & Turcotte, 1990). Hence those models cannot compute complex magma trajectories, but they can quantify changes in crack propagation velocity due to dynamic changes in the magma source or variations of crust properties and stress (Mériaux & Jaupart, 1998; Pinel & Jaupart, 2000; Traversa et al., 2010). Such models may also be used to estimate the dyke velocity and shape along a prescribed trajectory, assuming that the propagation velocity is limited by the viscous flow, and that fluid-dynamic processes do not affect the propagation path of the intrusion (Pinel et al., 2017).

Merging these two approaches in a single modeling scheme would result in an understanding of the trajectory and velocity of a magmatic dyke, which is of central importance for the mitigation of volcanic hazards associated with the creation of new eruptive vents. A previous attempt of coupling the two approaches in a comprehensive (2D) modeling scheme has been done by Dahm (2000b), this scheme inspired the work we present here. Additionally, more recent 3D numerical approaches that couple fluid-flow equations and fracture mechanics, can now be used to simulate the propagation of buoyant magmatic intrusions (Detournay & Napier, 2019; Napier & Detournay, 2020; Zia & Lecampion, 2020). However, all these formulations remain limited to planar fracture propagation in a homogeneous elastic medium.

Here, we present a new modeling scheme, to compute the dynamic shape of a growing fluid-filled crack, built with a plane strain BE technique (2D, with infinite out-of-plane extent and no strain along the out-of-plane direction). The model computes the crack shape accounting for the fluid viscosity and the crack propagation velocity, and has been designed in order to account for heterogeneous crustal stress and complex crack paths as in Maccaferri et al. (2011). Depending on the modeling purpose, the crack velocity can be given as input parameter, or calculated as output, based on an energy conservation equation, and the assumption that the main sources of energy dissipation are the brittle fracturing and the laminar viscous flow.

We validate our model results against relevant analytical solutions available in literature for a vertical fluid-filled crack in a homogeneous infinite elastic medium, and show how some of the most important parameters characterizing magmatic intrusions determine the shape and propagation velocity. These analyses allow us to determine the parameter space within which our model is valid. We then show results in different settings, including the effect of an inclined - curved - dyke trajectory, the interaction with crustal layering, and the effect of the free surface. We discuss our model assumptions and approximations, and eventually show an application to the dyking event that fed the 1998 Piton de la Fournaise (PdF) eruption (La Réunion Island).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Scheme of the boundary element crack model, at the current propagation step *k* and at the previous one *k* − 1. The crack is made of *N* + 1 dislocation elements of constant length *l* and variable opening *h*<sub>*i*</sub>.</figcaption></figure>

## 2. Governing Equations

We consider an ascending finite volume of fluid flowing through a thin crack that grows in an isotropic, plane strain, brittle elastic medium, such that the governing equations use linear elasticity and the lubrication theory for the rock deformation and the fluid flow respectively. In the following *s* will be the curvilinear coordinate (positive upwards) along the fracture cross-section (i.e., along the intersection between a vertical plane and the crack surface, which in general can be inclined or bending), and *r* will indicate the coordinate perpendicular to the fracture walls (Figure 1). Please note that all the variables used through the manuscript are listed, along with their definitions, in the “List of Notation” at the end of the article.

### 2.1. Elasticity

A dislocation surface $\Sigma$, surrounded by a dislocation line, defines a discontinuity in the displacement field that can be represented by the Burgers vector $\mathbf{b}$, which indicates the magnitude and direction of the displacement discontinuity over $\Sigma$. Combining the linear elastic equations for (a) the balance of linear momentum (Newton's second law), (b) the strain-displacement relations and (c) the linear constitutive equations (Hooke's law), it provides the following equation governing the boundary value problem (Chapter 5 from Jaeger et al., 2007):

$$ \frac{1}{1 - 2\nu} \nabla(\nabla \cdot \mathbf{b}) + \nabla^2 \mathbf{b} + \mathbf{F} = 0 \qquad (1) $$

with $\nu$ the Poisson's ratio and $\mathbf{F}$ the body force per unit volume.

The BE approach consists in combining dislocation elements to solve the so called “crack problem”: that is computing the distributed fracture-opening and -shear, given the stress boundary conditions at the crack walls. Here, we describe the (2D) cross-section of a fracture by combining linear dislocation elements in plane strain approximation, with constant opening and shear on each of them. We used the vertical dislocation solutions for a layered elastic medium computed by Bonafede and Rivalta (1999) and by Rivalta et al. (2002), which has been generalized for inclined dislocations in Maccaferri et al. (2010).
## 2.2. Fluid Flow

The motion of the fluid within the fracture can be described using the Navier-Stokes equation for the momentum conservation of the system. According to the lubrication theory, a fluid with viscosity $\eta$, moving through a thin fracture with opening $h(s)$, in absence of turbulence, flows with velocity parallel to the fracture walls $u = [0; u_s(r)]$, if the velocity changes are small enough for the inertial force to be neglected. Neglecting fluid body forces, the Navier-Stokes equation of motion simplifies to give the Reynolds equation for the viscous pressure gradient $\frac{\partial P_{visc}}{\partial s}$:

$$ -\frac{1}{\eta} \frac{\partial P_{visc}}{\partial s} = \frac{d^2 u_s}{dr^2} \qquad (2) $$

Solving the Reynolds equation for a Poiseuille flow between two parallel walls provides the velocity $u_s(r)$ of the fluid:

$$ u_s(r) = -\frac{\partial P_{visc}}{\partial s}(s, t) \frac{r(h-r)}{2\eta} \qquad (3) $$

Note that here we introduced the fluid velocity between parallel walls because in Section 3.2 we will discretize the fracture cross-section with dislocation elements with constant opening, and the fluid velocity $u_s(r)$ at the time $t$, within each dislocation element at $s$, will be represented by Equation 3.

Using Equation 3, in the definition for the fluid flow per unit length through $h(s)$, at the time $t$: $f(s, t) = \int_0^h u_s(r)dr$, provides the equation for the flux $f(s, t)$ in a laminar regime:

$$ f(s, t) = -\frac{1}{12\eta} \frac{\partial P_{visc}}{\partial s}(s, t) h^3 \qquad (4) $$

The flow of magma through the cross-section of the crack is related to its opening by the continuity equation that ensures the conservation of mass. In 2D, for a compressible fluid, it can be written as (Lecampion et al., 2018) follows:

$$ \frac{\partial(\rho_m(s, t)h(s, t))}{\partial t} = -\frac{\partial(\rho_m(s, t)f(s, t))}{\partial s} \qquad (5) $$

with $\rho_m(s, t)$ the fluid density. Note that an assumption of the Poiseuille flow is that the density of the fluid is constant along the crack, that is, it does not vary spatially, whereas we allow the fluid density to vary in time as the continuity equation accounts for the compressibility.

## 2.3. Energy Conservation

The total energy budget for the propagation of a buoyant volume of compressible magma includes several terms:

1. The potential strain- and gravitational-energy variations, that provide a “buoyancy” term;

2. The friction dissipation due to the magma viscous flow;

3. The fracture energy needed to extend the fracture surface;

4. The kinetic energy variations due to any change in the velocity of the intrusion.

In a closed system, the total energy is conserved over time, so that the sum of all energy terms, differentiated through time, must be zero. As we will discuss later in Section 5, within the range of velocity and viscosity values that are significant for magmatic dyke propagation, the variation of kinetic energy is negligible with respect to the other terms. The rate of the potential strain energy variations $\frac{\partial W}{\partial t}$ represents the change in the work done against elastic forces during propagation (Section S1.1 in Supporting Information S1), and the rate of gravitational energy variations $\frac{\partial G}{\partial t}$ represents the work done to move the magma body within the gravitational field (Section S1.2 in Supporting Information S1). The surface fracture energy $E_f$ (the energy per unit surface needed to fracture the rocks), is related to the fracture toughness $K_c$, the Poisson's ratio $\nu$, and the shear modulus $\mu$:

$$ E_f = K_c^2 \frac{1 - \nu}{2\mu} \qquad (6) $$
The rate at which mechanical energy is dissipated by viscosity per unit mass of fluid $\Phi$, can be written as (Batchelor, 1967) follows:

$$ \Phi = -2\eta e_{mn}e_{mn} \qquad (7) $$

where $e_{mn}$ is the strain rate tensor (in plane strain approximation), and $e_{mn}e_{mn}$ is the Einstein summation notation, implying that $\Phi$ is a scalar quantity. In 2D, the total viscous dissipation rate $\frac{\partial E_v}{\partial t}$ is the integral of $-\Phi$ over the area of the magma body A:

$$ \frac{\partial E_v}{\partial t} = \int_A \Phi dsdr \qquad (8) $$

Finally, we can express the conservation of energy as follows:

$$ \frac{\partial W}{\partial t} + \frac{\partial G}{\partial t} - \frac{\partial}{\partial t}(E_f \cdot ds + E_v) = 0 \qquad (9) $$

Note that in plane strain configurations, the fracture extends infinitely in the out-of-plane direction, therefore the energy contributions that will be computed within the context of the BE scheme, will have the units of energies per unit length (in the out-of-plane direction). A 3D estimate for these quantities can be obtained by multiplying them for the out-of-plane extension of a 3D fracture (the dyke “width”).

# 3. Numerical Model Formulation

## 3.1. The Quasi-Static Approach for a BE Fluid-Filled Fracture Propagation Model

Fluid-filled fractures, from the small - centimeter - scale of laboratory experiments, up to kilometer scale magmatic dykes, can be modeled numerically by using a set of N contiguous and interacting dislocation elements (BE technique). Here, we consider the model implemented by Maccaferri et al. (2011), where the dislocation elements are plane strain (2D), and the model refers to a vertical cross-section of the fracture. The fracture is filled with a compressible fluid, trapped between the crack walls. If the fluid is less dense than the host solid, the fracture would tend to grow upwards, driven by the fluid buoyancy. In Maccaferri et al. (2011), the fracture growth and propagation are modeled with a quasi-static approach, by adding a dislocation element at the shallower tip of the crack, and solving for the static shape of the fluid-filled crack. This approach neglects the effect of the fluid viscosity on the fluid pressure profile, and cannot account for the crack propagation velocity. The fluid-filled fracture may or may not grow based on an energy budget condition: if the fracture growth corresponds to a positive energy budget (i.e., a drop of the potential energy larger than the energy needed to create the new fracture), the propagation is allowed. An important feature of this model is to compute the trajectory followed by a propagating crack, by testing different directions for the crack growth and choosing the one that maximizes the energy budget of the system.

From a numerical prospective, this modeling scheme requires solving for the elastic shape of the fluid-filled fracture (Equation 1) at each propagation step, and computing the energy release of the system along different test directions (cf. Figure 1 and Figure S1 in Supporting Information S1). In order to compute the fracture shape it is necessary to solve the so called “crack problem”, that is, computing the opening $h$ and slip $q$ of each dislocation element, given the overpressure $\Delta P$ (i.e., the fluid excess pressure with respect to the confining stress) and the shear tractions $\tau$ acting at the center of each dislocation element (stress boundary conditions). The problem is defined by a set of linear equations linking the openings of all dislocations to their stress boundary conditions, through the influence coefficients that define the mutual effect of the opening ($T^i$, $S^i$) - and slip ($T^s$, $S^s$) - of each dislocation element on the others. These influence coefficients provide at each dislocation element for the tensile (and shear) stress due to the opening (and slip) of the other elements, given their spatial position, their dip and length, and accounting for the effect of an elastic interface or free surface (Bonafede & Rivalta, 1999; Rivalta et al., 2002). At the equilibrium, the stress produced by the crack opening and slip has to balance the boundary conditions:

$$ \left\{ \begin{array}{l} \sum_{j=1}^{N} [T_{i,j}^t h_j + T_{i,j}^s q_j] = -\Delta P(i) \\ \sum_{j=1}^{N} [S_{i,j}^t h_j + S_{i,j}^s q_j] = -\tau(i) \end{array} \right. \qquad (10) $$
Actually, this formulation is still not the final one to be solved as the overpressure term (ΔP) also depends on the crack opening due to the fluid compressibility, which links any change in the fluid volume to a proportional pressure change. However, due to its linearity, the overpressure term that depends on the crack opening can be moved to the left side of Equation 10, and grouped with the T' · h terms (cf. Dahm, 2000b; Maccaferri et al., 2011). Such rearrangement of the linear system will be simply indicated here by using the symbols ΔP̃ and T̃'. Therefore we can represent Equation 10 in matrix notation:

$$ \begin{bmatrix} \tilde{T}^t & T^s \\ S^t & S^s \end{bmatrix} \times \begin{bmatrix} h \\ q \end{bmatrix} = - \begin{bmatrix} \Delta \tilde{P} \\ \tau \end{bmatrix} \qquad (11) $$

where the matrix is 2N × 2N and can be inverted to solve for the unknowns h and q (with the additional constrain of positiveness for the openings h).

In this paper, we extend the formulation of Maccaferri et al. (2011) to a “dynamic” crack solution by introducing in the elastic problem defined by Equation 11 the effect of the fluid flow on the fluid pressure profile, such that ΔP<sub>dyn</sub> represents the dynamic overpressure:

$$ \begin{bmatrix} T^t & T^s \\ S^t & S^s \end{bmatrix} \times \begin{bmatrix} h \\ q \end{bmatrix} = - \begin{bmatrix} \Delta P_{dyn} \\ \tau \end{bmatrix} \qquad (12) $$

The system of equations represented by Equation 12 is not ready to be solved, in this form. ΔP<sub>dyn</sub> depends on the openings h, and in general Equation 12 does not even represent a linear system. The solution is not trivial: Dahm (2000b) presents a numerical solution for the vertical (straight) propagation of a fluid-filled crack growing at a constant velocity. We partially follow Dahm (2000b) approach, dropping the assumption of a constant crack growth velocity, and generalizing it to complex fracture trajectories. We derive an approximate formulation for ΔP<sub>dyn</sub> that allows us to write Equation 12 as a linear system, obtaining an approximate - fast - solution for the dynamic shape of a BE fluid-filled crack propagating at a certain velocity.

In the following section (Section 3.2), we will illustrate our procedure to compute ΔP<sub>dyn</sub>, along with the definitions of all the pressure terms and their mathematical formulation, that are also summarized in the “List of notations”.

## 3.2. Computing the Dynamic Shape of a Growing BE Fluid-Filled Crack

We define a coordinate system where x is the horizontal axis and z is the vertical axis (positive downwards). Within a growing fluid-filled crack, the fluid overpressure ΔP<sub>dyn</sub>, is defined as the difference between the fluid pressure P<sub>fluid</sub>(z), and the confining pressure P<sub>conf</sub>(z). P<sub>conf</sub>(z) is the normal stress acting at the crack walls at depth z, and includes the lithostatic pressure and any other possible source of crustal stress (i.e., tectonic or topographic). P<sub>fluid</sub>(z) can be split into three components:

1. a linear, depth dependent, hydrostatic pressure profile: ρ<sub>m</sub> · g · z, where g is the gravitational acceleration;

2. a pressure term due to the fluid compressibility (K): ΔP<sub>K</sub> = K ΔV/V₀, proportional to a fluid volume change (ΔV) with respect to a reference volume at vanishing lithostatic pressure (V₀);

3. a viscous pressure change: ΔP<sub>visc</sub>, due to the flow of a viscous fluid through the crack walls (note that ΔP<sub>visc</sub> is negative, leading to a reduction of the pressure gradient along the crack).

The fluid overpressure is written as follows:

$$ \Delta P_{dyn} = \rho \cdot g \cdot z + \Delta P_{K} - P_{conf}(z) + \Delta P_{visc} \qquad (13) $$

Note that when considering an incompressible fluid, the term ΔP<sub>K</sub> is zero.

We describe the fluid flow within the crack with similar assumptions as Dahm (2000b), that solved the steady problem of a buoyant fluid-filled crack (in plane strain condition) rising vertically at a constant velocity, within a homogeneous elastic crust (hence having a fixed shape during propagation). Similarly, we assume a constant crack growth velocity during a single propagation step of our BE scheme, however we allow for changes in the
velocity and direction of crack growth between consecutive propagation steps. At each propagation step, our model selects the direction for the fracture growth by testing a number $n_α$ of different directions of propagation, defined by an angle $δα$, such that the tested directions are $α(j) = ±0.5(j - 1)δα$, with $j = 1, ..., n_α$ (Figure 1 and Figure S1 in Supporting Information S1). $n_α$ and $δα$ are prescribed as input parameters, and they should not pose any contrived limit to the direction of fracture growth. Eventually, the selected direction is the one with the maximum energy release, that will also be used to constrain the crack propagation velocity. In doing so, we are able to account for complex propagation paths, changes in the crack shape, and the interaction with heterogeneous stress fields, and crustal heterogeneities.

Following the scheme for the BE crack growth in Figure 1, we compute the fluid flux through each element $i$. For the sake of simplicity we initially write the flux crossing the upper boundary of each dislocation element $f_i^{top}$, that is the volume of fluid that flows out of the $i$th dislocation element during the time $Δt$ elapsing between the propagation step $k-1$ and $k$. During the time $Δt$, the crack extends by one dislocation element of length $l$ at its tip. Therefore, we define the propagation velocity of the crack as $v = \frac{l}{Δt}$, that is the velocity at which the crack tip is advancing along its propagation path. The volume variations ($ΔV_i$) associated with each dislocation element of length $l$ during a propagation step can be written as follows:

$$ \Delta V_i = (h_i^k - h_{i-1}^{k-1}) \cdot l \qquad (14) $$

with $h_i^k$ being the $i$th dislocation's opening at the propagation step $k$ (i.e., the normal component of the Burger vector, $h$ in Equation 12). Note that in our 2D model, the “volume” variations, are actually “cross-sectional area” variations, under the assumption that the crack is infinite in the out-of-plane direction. For an incompressible fluid, integrating the continuity equation (Equation 5) leads to:

$$ f_i^{top} = -\frac{1}{\Delta t} \sum_{j=1}^{i} \Delta V_j \quad \text{with} \quad i = 1, \dots, N \qquad (15) $$

(cf. Section S2 in Supporting Information S1 for a detailed description of the discretization of Equation 5 that led to Equation 15).

Our fluid-flow assumptions imply that the fluid density at a given time is homogeneous in space. However, the fluid density can evolve through time as the fluid-filled crack rises to shallower depths, and it is subjected to a lower confining stress, causing a fluid expansion proportional to the fluid compressibility, and a consequent density decrease (as mass is conserved). For a compressible fluid, we need to modify Equation 15 accounting for the volume variation between the propagation step $k-1$ and $k$ due to any changes in the confining stress. Assuming that such volume change will be distributed along the crack in proportion to the opening of each dislocation element, we can write:

$$ f_i^{top} = -\frac{1}{\Delta t} \sum_{j=1}^{i} \left( \Delta V_j - \frac{V^k - V^{k-1}}{V^{k-1}} \cdot h_j^{k-1} \cdot l \right) \qquad (16) $$

where $V^k = \sum_{i=1}^{N} h_i^k \cdot l$ is the (2D) crack volume at the propagation step $k$. Substituting Equation 14 in Equation 16 and introducing $v$, we obtain:

$$ f_i^{top} = -v \sum_{j=1}^{i} \left( h_j^k - \frac{V^k}{V^{k-1}} h_{j-1}^{k-1} \right) \qquad (17) $$

In order to make use of the flux condition in our BE scheme, we compute this at the middle point of each dislocation element, since we will have to prescribe the stress boundary conditions at those points. Given Equation 17, we can write the flux at the middle point of the element $i$ as follows:

$$ \begin{aligned} f_i &= \frac{1}{2} (f_{i-1}^{top} + f_i^{top}) = f_i^{top} - \frac{1}{2} (f_i^{top} - f_{i-1}^{top}) \\ &= -v \sum_{j=1}^{i} \left( h_j^k - \frac{V^k}{V^{k-1}} h_{j-1}^{k-1} \right) + \frac{1}{2} v \left( h_i^k - \frac{V^k}{V^{k-1}} h_{i-1}^{k-1} \right) \qquad (18) \\ &= -v \cdot \Psi_i \end{aligned} $$
with $\Psi_i = \sum_{j=1}^i (h_j^k - \frac{V^k}{V^{k-1}} h_j^{k-1}) - \frac{1}{2} (h_i^k - \frac{V^k}{V^{k-1}} h_i^{k-1})$. Note that combining Equation 18 with the equations for the fluid flow (Equations 3 and 4, Section 2.2) provides the relation between the crack propagation velocity $v$ and the fluid flow velocity $u$. The next step will be to write the flux as a function of the viscous pressure gradient acting within each dislocation element, and finally solve for $\Delta P_{\text{visc}}(z)$.

We assume a Hagen-Poiseuille flow within piece-wise parallel fracture walls, which results in a constant viscous pressure gradient $\Gamma_i = \frac{\partial P_{\text{visc}}(s,f)}{\partial s}|_i$, and a parabolic fluid velocity profile $u_i$ within each dislocation element $i$ (Figure 1). Equation 4 leads to:

$$f_i = -\Gamma_i \cdot \frac{(h_i^{k-1})^3}{12\eta}, \qquad (19)$$

Note that here we use the dislocation opening $h_i^{k-1}$ to compute the flux between $k-1$ and $k$, this approximation is fundamental to simplify our set of equations (as it will be clear in the following), and we expect it to work if the crack shape does not change significantly during a single propagation step (hence the time resolution is high enough).

The viscous pressure change due to the $i$th element of length $l$ is $\Delta P_{\text{visc}}^i = l \cdot \Gamma_i$. The viscous pressure change at the bottom of the $i$th element is the sum of the pressure drops of all the elements underneath, and therefore the viscous pressure change at the middle point of the $i$th dislocation element will be:

$$\Delta P_{\text{visc}}^i = \sum_{m=1}^{i} \Delta P_{\text{visc}}^m - \frac{1}{2} \Delta P_{\text{visc}}^i \\ = l \cdot \left( \sum_{m=1}^{i} \Gamma_m - \frac{1}{2} \Gamma_i \right) \qquad (20)$$

We can now use Equation 19 to substitute $\Gamma$ in Equation 20, and then use Equation 18 to substitute $f_i$. In this way we obtain an expression for $\Delta P_{\text{visc}}^i$ that depends on our model parameters and unknowns:

$$\Delta P_{\text{visc}}^i = \sum_{m=1}^{i} \frac{12\nu\eta l}{(h_m^{k-1})^3} \sum_{j=1}^{m} \left(h_j^k - \frac{V^k}{V^{k-1}} h_j^{k-1}\right) \\ -\frac{1}{2} \sum_{m=1}^{i} \frac{12\nu\eta l}{(h_m^{k-1})^3} \left(h_m^k - \frac{V^k}{V^{k-1}} h_m^{k-1}\right) \\ -\frac{1}{2} \frac{12\nu\eta l}{(h_i^{k-1})^3} \sum_{m=1}^{i} \left(h_m^k - \frac{V^k}{V^{k-1}} h_m^{k-1}\right) \\ +\frac{1}{4} \frac{12\nu\eta l}{(h_i^{k-1})^3} \left(h_i^k - \frac{V^k}{V^{k-1}} h_i^{k-1}\right) \qquad (21)$$

It is important to note that Equation 21 is valid for $i = (1, ..., N-1)$, when $i = N$ Equation 19 cannot be used to constrain $\Gamma_i$, because at the propagation step $k-1$ the dislocation element $N$ has not been opened yet. We will consider the case $i = N$ separately.

The most important feature of Equation 21 is that it is linear with respect to $h^k$, which are the unknowns of the problem at the propagation step $k$, since $h^{k-1}$ is known from the previous propagation step. Using Equation 21 for $\Delta P_{\text{dyn}}$ in Equation 12, we will still have a linear system of equations. Rearranging Equation 21 allows us to bring all the terms that are linearly dependent on $h^k$ to the left side of Equation 12. Finally we group them, and incorporate their coefficients with the influence coefficients of the matrix in Equation 12. The dynamic shape of the crack at the propagation step $k$, given the fluid viscosity and the propagation velocity, will then be simply obtained by solving a linear system. Note that when the velocity is zero, $\Delta P_{\text{visc}}$ is null all along the crack and we have the static case.

We now need to find an expression for $\Delta P_{\text{visc}}^N$. Equation 20 and Equation 18, for $i = N$ provide conditions for the viscous pressure change and the flux in $N$:

$$\Delta P_{\text{visc}}^N = \Delta P_{N-1}^{\text{visc}} + \frac{l}{2} \Gamma_N \quad \text{and} \quad f_N = \frac{1}{2} \nu h_N^k \qquad (22)$$
and we need an expression similar to Equation 19 that links $\Gamma_N$ and $f_N$, $h_N^{k-1}$ is not defined, and we cannot use $h_N^k$ otherwise Equation 12 would no longer be linear. Therefore, we provide an estimate of the flux $\Gamma_N$ approximating the opening of the tip element at the propagation step $k$ ($h_N^k$) with the opening of the tip element at the previous propagation step ($h_{N-1}^{k-1}$), and complement Equation 22 with:

$$ \Gamma_N = - \frac{12\eta}{(h_{N-1}^{k-1})^3} f_N \qquad (23) $$

In this way, we preserve the linearity of Equation 12 and we do not disregard the contribution of the Nth element of the BE crack on the viscous pressure profile. The approximations that we introduced have the advantage that the problem is simple to solve from a mathematical and computational prospective, while this resolves the viscous pressure profile accurately within the crack tail and head (as we will show in the result section). Nevertheless, our approximations underestimate the viscous pressure drop at the upper tip of the crack, and this will be critical when trying to simulate the propagation for vanishing values of rock fracture toughness.

In this section, we have shown that Equation 12 can be written as a linear system if we write $\Delta P_{visc}$ according to Equation 21, provided that we know the fluid-filled crack shape at the previous propagation step and its propagation velocity. In the next two sections, we will explain how we initiate the fluid-filled crack propagation, for example, how we compute the shape of the crack at the first propagation step, and how we can solve for the crack propagation velocity, assuming that the velocity is limited by the viscous fluid flow.

### 3.3. Initialization of the Model

In order to be able to compute the crack shape at a given propagation step, our modeling scheme requires the crack shape and velocity at the previous one. To initiate a simulation we use the method described in Dahm (2000b), which computed the shape of a fluid-filled BE crack growing at a constant velocity. Following Dahm (2000b) approach, we assume that the initial shape of the crack is stationary, that is the crack is already propagating with constant shape and velocity at the time the dynamic simulation starts. This stationary shape is computed following the iterative procedure in Dahm (2000b), starting with a static shape, estimating a viscous pressure drop associated with it, and use it to update the crack shape, and so forth. This scheme requires a few numerical iterations to reach a stable stationary solution. After computing the stationary shape, the procedure for the dynamic simulation is initiated, the crack starts growing in length, and develops a thin open tail region.

### 3.4. Computing the Fluid-Filled Crack Propagation Velocity Based on an Energy Budget Equation

We have shown how to solve Equation 12 given the crack propagation velocity $v$. In order to solve for $v$, we need to introduce an additional equation to the system of Equation 12, without introducing new unknowns to the problem. We set this additional condition using the energy budget associated with the fluid-filled crack growth.

As we already mentioned in Section 3.1, in the quasi-static model the crack propagation is allowed when the energy release associated with the crack elongation is larger than the energy needed to create the additional fracture surface. Between two consecutive propagation steps, we define the energy release ($\Delta E$) as the sum of the strain ($\Delta W$) and gravitational ($\Delta G$) energy contributions (Maccaferri et al., 2010, 2011), so that $\Delta E = \Delta W + \Delta G$. The energy needed to extend the fracture of an element of length $l$ is $E_f \cdot l$ (which is an energy per unit length, cf. Equation 6 and Section 2.3). In the quasi-static model, the crack propagation is allowed when $\Delta E \ge E_f \cdot l$, which represents the energetic criterion for fluid-filled fracture propagation (cf. Dahm, 2000b), and is equivalent to the condition of stress intensity factor larger or equal to the fracture toughness (cf. Rubin, 1998).

Such propagation condition implies that the excess energy release with respect to the fracture energy ($\Delta E - E_f \cdot l \ge 0$) is somehow spent, or dissipated, by the dynamic processes that are neglected in the quasi-static modeling scheme, such as kinetic energy variations, not-elastic processes at the crack tip, and the internal friction due to the viscous fluid flow. In the assumption that the viscous flow is the dominating force limiting the crack propagation velocity, we can write a new energy balance equation by integrating Equation 9 over time (from $t$ to $t + \Delta t$):

$$ \Delta E = E_f \cdot l + \Delta E_v \qquad (24) $$
where $\Delta E_v$ is the energy dissipated by viscous flow during $\Delta t$ and $E_f \cdot l$ results from the time-integral of $E_f \cdot ds$ in Equation 9 ($E_f$ does not depend on time, and the integral of $ds$ over $\Delta t$ is the amount of fracture growth, $l$). Equation 24 will depend on the crack propagation velocity $v$, since the fluid flow can be expressed as a function of $v$ (cf. Equation 18), and therefore represents a condition which allows us to constrain the crack propagation velocity.

In the following, we will write a formula for $\Delta E_v$, while for the calculation of $\Delta E$ we refer to Maccaferri et al. (2011).

For the Hagen-Poiseuille flow, the fluid velocity field, within each dislocation element $i$, is in the form $u = [0; u_s(r)]$ and Equation 7 simplifies to:

$$ \Phi = -\eta \left( \frac{\partial u_s}{\partial r} \right)^2 \quad \text{with} \quad i = 1, \dots, N \qquad (25) $$

For our BE crack, we can write the total viscous dissipation $\Delta E_v = \frac{\partial E_v}{\partial t} \cdot \Delta t$ introducing Equations 3 and 25 in Equation 8. It follows that:

$$ \Delta E_v = \frac{l}{12\eta v} \sum_{i=1}^{N} (\bar{h}_i)^3 (\Gamma_i)^2 \qquad (26) $$

Note that the fracture opening at each dislocation element ($h_i$) actually varies during a propagation step. Here, we choose to use the average values $\bar{h}_i = (h_i^k + h_i^{k-1})/2$ between the opening before and at the end of the $k$th propagation step ($h_i^{k-1}$ and $h_i^k$, respectively). By grouping $v$ in the equation for $\Gamma_i$ (Equations 18 and 19), the viscous dissipation $\Delta E_v$ at propagation step $k$ can be written as follows:

$$ \Delta E_v = v \cdot D \qquad (27) $$

with $D = 12\eta l \sum_{i=1}^{N} (h_i)^3 (\frac{\psi_i}{h_i^{k-1}})^2$. Using Equation 24, we can make explicit the velocity condition at each propagation step:

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

However, both $\Delta E$ and $D$ depend on $v$ (through the crack openings $h_i$ and the viscous pressure gradient $\Gamma_i$, which are velocity-dependent quantities), therefore Equation 28 must be evaluated iteratively at each propagation step by updating the right side of the equation, until convergence is eventually reached.

# 4. Results

The model allows us to simulate the propagation of a given mass of compressible buoyant fluid filling a fracture of a given initial length. In this section, we first validate our numerical scheme against existing solutions for the shape, velocity, and overpressure profile of a vertical fluid-filled fracture propagating with constant volume in a homogeneous infinite elastic medium. We then present results from numerical simulations with different model settings, showing the effect of varying several parameters characterizing magma and rocks, including the effect of the free surface, crustal layers, and an inclined - bending - trajectory followed by the intrusion when crossing the interface between two crustal layers. Finally, we perform a parametric analysis to constrain the domain of applicability of the numerical scheme.

## 4.1. Validation Against Previous Models

In order to validate our numerical results, we use analytical, numerical, and asymptotic solutions for the velocity and the shape of a vertical fluid-filled fracture, containing a finite volume of buoyant, viscous fluid. These solutions have been obtained by previous authors in the context of the so called “lubrication theory school” (as defined in Rivalta et al., 2015).

The velocity of crack propagation (and its time evolution) can be obtained from Spence and Turcotte (1990) by differentiating with respect to time their Equation 22 that links the time for the magma to reach the surface
to the intrusion's parameters. According to Spence and Turcotte (1990) formulation, the velocity $v_{th}$ is valid for non-dimensional times $t^* \gg 1$ (with $t^* = \frac{A_0^2 g \Delta \rho t}{\pi^2 (L/2)^3 \eta}$, $\Delta \rho$ being the density contrast between the fluid and the surrounding medium, $L$ the initial crack length and $A_0$ the initial cross-sectional area of the magma body) and decays as $t^{-2/3}$ (cf. Davis et al., 2023):

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

Note that the velocity is independent from the rock fracture toughness $K_c$, because this solution applies to the limit of low fracture toughness.

Here, we use typical magma and rock properties for basaltic intrusions at mid-crustal depth, considering a finite batch of magma with a viscosity $\eta = 100$ Pa·s, a buoyancy $\Delta \rho = 300$ kg/m³, and a cross-sectional area $A_0 = 0.009$ km², propagating within a homogeneous crust with Poisson's ratio $\nu = 0.25$, and shear modulus $\mu = 20$ GPa. Note that in our simulations we use a rather low compressibility for magma (bulk modulus $K = 20$ GPa), in order to compare our results with the theoretical velocity from Equation 29, which is for an incompressible fluid.

We run several simulations progressively lowering the rock fracture toughness (i.e., decreasing the fracture energy $E_f$ from 20 to 8, 7, 6, and 5 MPa·m, pink - green - blue - purple - and black lines on Figure 2a, respectively), to compare our results in the limit of low fracture toughness. Figure 2a shows that the numerical solutions (solid lines) tend to the analytical formula (dashed line) as $t^*$ grows. We calculated that the relative difference between the numerical solution for $E_f = 5$ MPa and the analytical one is 29.6% at $t^* = 10$, but goes down to 17.4% at $t^* = 20$, and 7.1% at $t^* = 40$. For the same numerical velocity profile, we computed the exponent of the best fitting power low, which can be directly compared with the analytical exponent, -2/3: we obtain -0.50 for $t^* > 20$, but for $t^* > 40$ it reaches -0.55. The fact that the numerical solution differs from the analytical one for the first part of propagation, may partially be due to the different initial conditions (Spence & Turcotte, 1990, start with an elliptical crack, our initial shape accounts for the magma buoyancy), but may also be due to the fact that our simulations cannot be computed for vanishing fracture toughness. The latter implies that at the beginning of our simulations, when the fluid-filled crack is still developing a tail region, the resistance to fracture at the crack tip still affects the crack velocity. As the crack tail grows, viscous forces progressively become the limiting factor for the crack propagation velocity, and the numerical model converges toward the analytical solution for $K_c = 0$.

Another important observation from our numerical simulations is that they become less and less stable as we lower the rock fracture toughness. In fact when lowering the fracture toughness below the values used in Figure 2a, during the first part of propagation, the velocity and crack opening start to display some oscillations, that can even prevent the simulation to continue—if they grow large enough (Figure S2 in Supporting Information S1). This limitation is very likely due to the approximation done for the viscous pressure drop at the crack tip (Equation 23), and we will discuss it further in the next section. However, it is important to notice that even though we cannot arbitrarily lower the fracture toughness, our simulations can reach a fracture propagation regime with the velocity dominated by viscous forces, as our velocity profiles are independent from the fracture toughness—for the lower values we used—and very close to the analytical velocity (cf. Figure 2a).

In general, the shape of a propagating fluid-filled fracture is characterized by a tear-drop head region (where the buoyancy forces dominate), and a thinner tail region (where the viscous forces dominate). Equation 25 in Spence and Turcotte (1990) provides an analytical formula for the shape of the fluid-filled crack below the crack head. This formula can be used for a comparison with the thickness of the tail region from our numerical simulation. In Figure 2b, we show the dyke opening from our numerical simulations (plain lines) and the analytical formula (dashed lines) for dimensionless time $t^* = 14$ and 43 (this formula from Spence & Turcotte, 1990, also works in the limit of $t^* \gg 1$). The fit between the numerical and analytical shape of the crack tail is remarkably good, especially for large $t^*$, when the numerical and analytical viscous dominated crack velocities begin to converge.

Rubin (1998) extended the work done by Spence and Turcotte (1990), providing the shape of the entire fluid-filled crack within the same set of assumptions from the lubrication theory, and introducing the elasticity and the resistance to fracturing at the tip. Rubin (1998) provided a set of equations that he solved numerically to display the full crack shape and the fluid overpressure profile for a given set of parameters (Figure 3 in Rubin, 1998). We
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Evolution of the propagation velocity for simulations with various *E*<sub>f</sub> (plain lines), as function of dimensionless time *t*\* given by Equation 4 in Spence and Turcotte (1990). The black dashed line shows the evolution of theoretical velocity from Equation 29 inferred from Spence and Turcotte (1990). The velocities can be compared for *t*\* >10, that is, after the gray hashed area. Note that for dimensionless time *t*\* >32 (black arrow) the difference between numerical velocities and the theoretical value from Spence and Turcotte (1990) is less than 10%. (b) Crack opening represented at two different dimensionless times (pink arrows on (a)) according to Equation 25 from Spence and Turcotte (1990) (dashed lines) and results from the numerical simulations (plain pink lines) for the crack propagating in a viscous dominated regime (*E*<sub>f</sub> = 5 MPa·m, pink curve in (a). In panels (c) and (d), we modified Figure 3 from Rubin (1998) and plot our results (red dotted curves) along the dimensionless dyke length *z*/L<sub>w</sub> (normalized by the Weertman length scale) of dimensionless overpressure ΔP/P (normalized by the pressure scale) (c) and dyke half thickness *w*/W (normalized by the thickness scale) (d) from the numerical simulations using the boundary scheme set with the parameters described by Rubin (1998). Plain black curves and dashed curves represent the data from Rubin (1998), respectively the overpressure and dyke half thickness at selected times during ascent of a buoyant dyke of fixed volume, and the pressure and half thickness of static Weertman crack. (e) Asymptotic solutions for crack tip shape estimated using Dontsov and Peirce (2015) for a viscous-dominated crack (pink dashed line) and for a fracture-dominated crack (black dashed line) of 7 and 7.75 km length respectively. The pink circles are the results of the numerical simulation for a crack propagating within low fracture toughness medium (*E*<sub>f</sub> = 5 MPa·m, *K*<sub>c</sub> = 516 MPa·m<sup>1/2</sup>) and the black circles are results for high fracture toughness regime (*E*<sub>f</sub> = 20 MPa·m, *K*<sub>c</sub> = 1,033 MPa·m<sup>1/2</sup>).</figcaption></figure>

use Rubin (1998) results for a direct comparison with our numerical simulations: the shape and overpressure profiles that we obtained using the same set of parameters as Rubin (1998), are plotted as red-dotted lines on top of Rubin (1998) results in Figures 2c and 2d. The solutions for the crack shapes are remarkably similar, and the overpressure profiles along the crack display a very good fit, except for the overpressure drop at the tip of the crack, which is underestimated by our numerical simulation. This effect has been anticipated in Section 3.2, when introducing the approximation on the viscous pressure drop at the crack tip (Equation 23). We also notice that, as for the comparison with Spence and Turcotte (1990), our initial conditions are not the same as Rubin (1998), which—similarly to Spence and Turcotte (1990)—starts with an elliptical crack. This is reflected in the different overpressure profiles at the crack bottom at the beginning of our simulation (cf. Figure 2c, first red-dotted curve).

The previous comparisons were performed in the viscous dominated regime, next we compare our simulations with the asymptotic solutions for the crack tip shape in both fracture- and viscous-dominated regimes. We used the Equations 2.6a and 2.6c in Dontsov and Peirce (2015), which are: $h = \beta_k \chi^{1/2}$ for the fracture dominated, and $h = \beta_v \chi^{2/3}$ for the viscous dominated regime (where $\chi$ is the distance from the crack tip and $\beta_k$ and $\beta_v$ are coefficients that depend on the intrusion parameters, as defined in Dontsov & Peirce, 2015). Using the same parameters as in Figures 2a and 2b, we plotted on Figure 2e the asymptotic solutions as black dashed line ($\beta_k = 0.0618$), and the pink dashed line ($\beta_v = 0.0161$). The dyke tips from our simulations are plotted as black and pink dots for the
fracture- and viscous-dominated regimes: $E_f$ = 20 and 5 MPa, respectively. We enhanced the tip shapes by zooming in the uppermost 0.25 km of these two cracks, that are 7.0 and 7.75 km long, respectively. Figure 2e shows a good agreement for the dyke tip shape in the fracture dominated regime (which is not surprising since for large fracture toughness the dyke shape gets closer to the static Weertman shape, which is well modeled by the static version of our BE code), and is reasonably close to the viscous dominated asymptotic, considering also that we used a positive fracture toughness, corresponding to $E_f$ = 5 MPa.

## 4.2. Numerical Model Results

Here, we present results from our numerical simulations considering some configurations that could not be addressed with previous models. We set magma and rock parameters according to the values reported in Figure 3 and Table S1 in Supporting Information S1, first row. All simulations start with the initial stationary shape described in Section 3.3, dark blue in Figure 3. For comparison we also show the corresponding static crack shapes as dashed orange lines. We represent the temporal evolution of the crack shape—sampled at regular intervals—and the velocity profiles, for six different simulations.

We first look at the effect of the elastic free-surface, Figure 3a. For comparison, here we also plot the velocity profile corresponding to the infinite elastic medium, red solid line in Figure 3a. The free surface causes an acceleration of the intrusion that starts to be appreciable when the dyke reaches a depth of ~0.5L (where L is the crack initial length). The dyke accelerates, reaching velocities of about an order of magnitude larger than the corresponding values without any free surface, but only in the very last part of propagation (≤0.1L). An acceleration of fluid-filled fractures close to the surface has also been observed in previous analog experiments (i.e., Rivalta et al., 2005; Rivalta & Dahm, 2006), and explained in terms of strain energy release rate by Pinel et al. (2022). In Figure 3b, we use the same set of parameters as in Figure 3a, except for the magma viscosity which is reduced by a factor two. We obtain a very similar solution, but with a velocity that is twice as much as the one in Figure 3a. Note that the maximum opening of the dyke and the thickness of its tail at the same depth are the same for these two simulations. A very relevant parameter for magmatic intrusions is the magma compressibility, which directly affects the buoyancy of magma, as the lithostatic pressure decreases when an intrusion rises to the surface. In Figure 3c, we decrease the magma bulk modulus by a factor two, with respect to Figure 3a: we notice that the velocity is initially slower, due to a lower magma buoyancy, with respect to the low compressibility case. However, as the magma rises, the buoyancy increases, and the velocity starts to increase. It is remarkable, how the velocity profile obtained for a compressible magma is fundamentally different with respect to the one for an incompressible fluid, that decays as $t^{-2/3}$, rather than increasing. In Figure 3d, we reduce the rigidity of the crust by a factor two with respect to Figure 3a. Here, we also reduce the rock fracture toughness such that the value of $E_f$ does not change (Equation 6). A lower rigidity produces a larger propagation velocity, which is due to a higher rate of energy release ($\Delta E$ in Equation 28), and a thicker tail (a more compliant rock deforms more easily), allowing for a faster magma flow.

In the last two simulations, we consider a layered crust. We set the interface between two layers with different rigidities and fracture toughness values at a depth of 15 km. Here, we do not consider the effect of the free surface, and stop the simulation at 5 km depth, after the dyke propagated through 10 km of crust. This is because the model allows only one analytical discontinuity in the elastic parameters, which is either a free surface or an interface between two layers. In Figures 3e and 3f, the rigidity and fracture energy decrease of a factor two across the layers (so that also the fracture toughness is reduced by the same factor, cf. Equation 6). In Figure 3e, the dyke starts vertically (dipping 90°), while in Figure 3f the initial dip angle is 60°. In both simulations, the dyke accelerates while approaching the interface with the softer layer (similarly to the free surface effects), and reaches a peak in the propagation velocity when the tip crosses the layers' interface. A peak in the energy release across such interfaces has been also described in Maccaferri et al. (2010). Once the dyke tip enters the upper layer the velocity starts decreasing but remains larger than the minimum velocity in the lower layer. In Figure 3f, the initial lower dip angle results in a lower velocity (an inclined dyke has a reduced vertical extension, and therefore buoyancy). When crossing the layers' interface, the dip angle of the dyke tip changes (red dotted line in Figure 3f), and the final dip—after the dyke head crossed the interface - is 65°. The deflection of a dyke crossing two layers' interface was already shown in Maccaferri et al. (2010), using a quasi-static BE dyke-code (i.e., neglecting the effect of magma flow). We also found that with this set of parameters, the deflection obtained with the quasi-static BE dyke-code is the same as the one obtained with the current modeling scheme, indicating that
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Temporal evolution of crack shapes (left panels) and velocity *v* (right panels) of simulations for (a) the propagation of a finite volume at a variable velocity within a rigid medium under free surface condition (reference configuration) and without the free surface effect (red line on the right panel), (b) for a less viscous magma, (c) for a more compressible magma, (d) for a less rigid medium, (e) for crustal layering with the interface represented by the dashed line at 10 km depth and (f) for an inclined dyke where the evolution of dip is shown in red dots. Initial stationary crack shape (dark blue shape, *t* = 0 s) compared with static crack shape (dashed orange shape). Representations of dynamic crack shapes are sampled at constant space intervals. The color code indicates the time. Time scale is normalized by the total duration of propagation for a more compressible magma (c). The maximum crack opening and the tail opening at 13 km depth are indicated for subfigure (a–d).</figcaption></figure>
for this magma viscosity (and generally this set of parameters), the dyke trajectory is not significantly affected by
the magma flow. However, more tests and simulations are needed, as exploring a wider range of parameters may
reveal some different effects.

4.3. Parametric Analysis and the Energy Budget

The energy balance equation (Equation 24) constrains the crack propagation velocity, and provides a measure of the relative contribution to the total energy release ($\Delta E$) of the viscous energy dissipation ($\Delta E_v$), and the fracture energy ($E_f \cdot l$). Here, we show the relative contribution of $\Delta E_v$ to $\Delta E$, and how it changes as a function of magma and rock parameters.

As we anticipated, our algorithm fails for low values of $E_f$, which is when $\Delta E_v$ tends to $\Delta E$. This parametric analysis, aims at quantifying the limit of $\Delta E_v$ to $\Delta E$ for which our model provides stable solutions, and shows how some model parameters affect the relative contribution of $\Delta E_v$ to $\Delta E$, possibly bringing our model close to its limit of applicability.

In Figure 4a, we use the set of parameters in the first row of Table S1 in Supporting Information S1 varying $E_f$ and we plot the relative contribution of $\Delta E_v$ to $\Delta E$. For large values of $E_f$, our model provides stable solutions, and the crack propagation velocity $v$ progressively approaches zero, when the contribution of $\Delta E_v$ to $\Delta E$ tends to zero ($E_f$ tends to the critical value that allows propagation with this set of parameters). We find that the lowest value of $E_f$ for which the model provides stable solutions, corresponds to a contribution of $\Delta E_v$ to $\Delta E$ of about 75%. For lower values of $E_f$ the contribution of $\Delta E_v$ to $\Delta E$ should further increase (reaching 100% for $E_f = 0$), but our model starts to display numerical instabilities (Figure S2 in Supporting Information S1) causing oscillations in the propagation velocity and eventually resulting in a null opening at the propagating tip (ending the crack propagation procedure).

In Figure 4b, we show that, for a given value of fracture energy ($E_f = 12 \text{ MPa} \cdot \text{m}$), varying the magma viscosity over several order of magnitude only affects the dyke propagation velocity without affecting the relative contribution of $\Delta E_v$ to $\Delta E$.

In Figure 4c, we show the effect of varying the rock rigidity with constant fracture energy, which implies a variation of the fracture toughness according to Equation 6. We observe that lowering the rigidity produces an increased contribution of $\Delta E_v$ to $\Delta E$, also increasing the crack propagation velocity. As a consequence, given a set of parameters, our model will progressively approaches its limit of applicability when lowering the rock shear modulus.

Similarly, in Figure 4d, we show that increasing the cross-sectional area of the dyke $A_0$ (equivalent to the magma volume for our 2D model) increases the contribution of $\Delta E_v$ to $\Delta E$ and the corresponding propagation velocity, moving toward the boundary of our model applicability.

5. Discussion About the Domain of Applicability and Model Simplifications

The “viscous dominated” and the “fracture dominated” regimes capture two end-members of magmatic dyke propagation processes. The numerical scheme that we proposed here, can access the configurations between these end-members: we have shown that our model is able to handle configurations where the contribution of viscous energy dissipation with respect to the total energy balance remains below 75%. Our numerical simulations become unstable as $K_c$ becomes small ($E_f \lesssim 0.25 \cdot \Delta E$), and this seems to be related to the underestimation of the viscous pressure drop at the crack tip (cf. Figure 2c), due to the approximation introduced in Equation 23. This directly leads to an underestimate of the viscous energy dissipation, which impacts our velocity calculation through the energy balance equation, introducing a critical source of numerical instability (Figure S2 in Supporting Information S1). These numerical instabilities rise in the initial part of the dyke propagation path, when the underestimation of the viscous pressure drop at the crack tip is larger. We noticed that a finer discretization of the BE crack actually helps delaying the rise of the instabilities (allowing for lower fracture toughness values), but does not eliminate the problem.

Dropping the assumption in Equation 23 is not trivial and would require major modifications in the modeling scheme, but may represent an important future development in order to overcome the current limitation and
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Evolution of viscous dissipation contribution with respect to the total energy produced by the system, for magmatic intrusions ascending below the surface. z* is the depth of the crack tip normalized by the initial depth for numerical simulations. Influence of (a) the fracture energy (i.e., fracture toughness) for a given set of parameters (first row in Table S1 in Supporting Information S1). Fixing the fracture energy E<sub>f</sub> = 12 MPa·m, we also vary (b) the viscosity of the magma (second row in Table S1 in Supporting Information S1), (c) the rigidity of the medium (fourth row in Table S1 in Supporting Information S1) and (d) the initial volume of dyke (fifth row in Table S1 in Supporting Information S1). Some values of K<sub>c</sub> are indicated at the bottom of the curves, along with the mean velocity of propagation estimated along the “constant part” of the propagation (defined by the gray domain).</figcaption></figure>
achieve stable solutions also for vanishing fracture toughness. A possible strategy may be to use asymptotic solutions for the crack tip, such as the ones provided in Dontsov and Peirce (2015). Such approach may possibly be coupled with non-uniform crack discretization (with higher resolution at the propagating tip), and/or using “linear” fracture elements (which are characterized by a linear opening, rather than constant), in order to improve the BE crack shape at the propagating tip.

Another numerical parameter that affects the stability of the numerical calculation is the “input velocity”, which represents an initial velocity guess (given as input), and is used to initiate the propagation and find a stable solution for the first propagation step. However, not all values of input velocity allow our algorithm to successfully compute the initial crack shape and velocity: we found that, when $E_f > 0.25 \cdot \Delta E$, the initial velocity guess should be within a range of about one order of magnitude with the actual—but yet unknown—velocity. In this case, it is relatively easy to find a good velocity guess with a simple trial-and-error approach. However, when $E_f$ approaches $0.25 \cdot \Delta E$, the acceptable interval of input velocities become progressively narrower. Finally, we verified that stable solutions obtained with different input velocities are consistent with each other, and fully overlapping velocity profiles are found after a very few initial propagation steps (Figure S4 in Supporting Information S1).

The energy balance equation is also affected by other model approximations: the 2D (plane-strain assumption) consists in an infinite dyke width (out-of-plane extension). Obviously, a 3D crack has a finite width and the effect of its lateral boundaries on the viscous dissipation are neglected here. In addition, the 3D shape of a hydro-fracture is typically rounded at the upper tip (Dahm, 2000a; Watanabe et al., 2002), such that the actual width decreases toward the upper tip of the crack. Therefore, the 2D assumption may underestimate the fluid flux at the upper tip of the crack, and consequently the viscous energy dissipation. Even though we cannot provide an estimate of how the 2D assumption affects the viscous energy dissipation, we expect that the range of viscous dissipation energies that our model is able to solve may actually extend beyond the threshold of 75% of the total energy, if 3D effects would be taken into account. We plan to investigate this aspect in the next future with the help of analog laboratory experiments, which are 3D.

Our energy budget equation (Equation 24) assumes a purely brittle-elastic fracturing process, disregarding non-elastic effects at the crack tip. Those effects may be particularly important as the size of the dyke (and its propagation velocity) increases, as suggested by previous studies that showed scale-dependent fracture toughness (cf. Rivalta et al., 2015, and reference therein). Also this point may be addressed in the future by comparing numerical simulations with laboratory experiments.

One of the main advantage of our modeling scheme is that the calculation of the fluid-filled crack shape results from the inversion of a linear system. The linearity of the system is obtained assuming that we know the initial crack shape and propagation velocity, and that the fluid flow is laminar and stationary during a single propagation step. However, at the next propagation step, we allow the crack (and fluid) velocity to “jump” to a different “stationary” value. Clearly, this assumption may become more critical if the velocity (and shape) variations between consecutive propagation steps grow large.

Other important model simplifications that we have made consist in considering magma as a homogeneous, Newtonian, single phase fluid and whose propagation within the dyke is laminar and unidirectional with no inertia. This may not be the case for complex dyke trajectories implying changes in the flow direction. Magma is rarely homogeneous and may be composed of different phases including crystals and gases. The process of gas exsolution can indeed lead to the formation of a gas-filled region at the crack tip, that would complicate the propagation dynamic, and would have implications on the energy release, trajectory, and velocity. Such process may be particularly effective at relatively shallow depths (2–3 km from the surface). Dealing with a two-phase fluid implies that the density, viscosity, compressibility, and volume of the fluid vary within the crack and will therefore depend on the depth, that is, on the confining stress. In addition, we do not account for temperature related processes, including heat conduction and convection, which may also affect the energy balance and therefore the velocity and trajectory of the dyke.

Our energy budget equation disregards the contribution of kinetic energy variations associated with acceleration or deceleration in the magma motion and rock displacements (Sections S1.3 and S1.4 in Supporting Information S1). Kinetic energy variations due to the magma motions can be estimated by computing the velocity variations of the center of mass of the intrusion. The kinetic energy variations associated with the rock displacements can be estimated considering the velocity at which the crack walls move during a propagation
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Depth of volcano-tectonic events as a function of time during the 3 days of the deep seismic crisis that started on 6 March 1998 (after Battaglia et al., 2005, with permission). The color code indicates the temporal distribution in days. Orange lines represent the velocity trend estimated over three time periods from Battaglia et al. (2005), brown line shows the average velocity used in the present paper to compare with numerical results and the red vertical line shows the onset of the eruption.</figcaption></figure>

step. We computed kinetic energy variations associated with the simulation of a magmatic intrusion (using the parameters in the first row of Table S1 in Supporting Information S1), and display our results in Figure S5 in Supporting Information S1. The larger contribution to the kinetic energy variations is provided by the magma motion rather than the fracture walls (Figure S5a in Supporting Information S1), and focuses when the dyke approaches the free surface, and the propagation velocity increases. Overall, the kinetic energy remains largely lower than the total energy release of the system $\Delta E$ (Figure S5b in Supporting Information S1), and is effectively negligible.

Despite all these simplifications, our approach allows us to describe the complex interaction between rock fracturing and magma flow within a dyke, over a wide range of model parameters which are relevant for magmatic intrusions in different volcanic contexts. Importantly, our approach—based on the BE technique—allows to simultaneously solve for the crack trajectory and its propagation velocity (Figure 3f).

# 6. Velocity of Dyke Propagation Preceding the 1998 PdF Eruption

## 6.1. Context and Seismic Data From the Magmatic Intrusion

PdF volcano has been built on the flank of two older volcanoes, Piton des Neiges and Les Alizées on La Réunion Island. It is located in the western Indian Ocean, at the youngest tip of the Réunion hotspot. This basaltic shield volcano is one of the most active on Earth with a mean time between consecutive eruptions of 9 months over the last decades (Peltier et al., 2009). Geophysical and geochemical measurements, performed since 1980 by the Volcanological Observatory of PdF, have identified three magma reservoirs beneath the volcano, located at different depths, in coincidence with structural discontinuities (e.g., Battaglia et al., 2005; Beauducel et al., 2020; Peltier et al., 2009). In particular, at sea level, magma density reaches that of the surrounding rock (level of neutral buoyancy), favoring magma stagnation in a shallow reservoir. This shallow reservoir is being occasionally recharged by magma from a deeper reservoir, located at the boundary between the volcanic edifice and the oceanic crust, at about 5–6 km below sea level (bsl). In addition, a seismic swarm from 1996 has potentially revealed another reservoir at 12.5 km bsl, at the crust-mantle boundary (Battaglia et al., 2005).

As an application of the modeling scheme we illustrated in this article, we run some simulations for the propagation of the dyke leading to the 9 March 1998 eruption of PdF volcano. Similarly to Rivalta and Dahm (2006), we assume that the migration of the hypocenters of a seismic swarm preceding the eruption indicates the approximate position of the tip of a magmatic dyke. Out of more than 3,100 seismic events detection, 583 earthquakes were located by Battaglia et al. (2005), showing a vertical migration of approximately 7 km during the three days preceding the eruption. Battaglia et al. (2005) suggested three trends of hypocenter migration: at the early stage of the seismic swarm, between 6 and 7 March (corresponding to Julian days 65–67), the seismicity propagates at a velocity of 0.4 m/min (0.024 km/hr, Figure 5). The frequency of volcano-tectonic events started to increase from the beginning of 8 March (Figure 5), with a pronounced upward migration at a rate of 1.59 m/min (0.095 km/hr). On 9 March, seismicity indicated a propagation velocity of 1.25 m/min (0.075 km/hr) up to 10.00 a.m. and then the propagation velocity drastically increased, reaching approximately 2.5 km/hr (Battaglia et al., 2005), and leading to the eruption at 11.05 a.m.

With our simulations, we aim to fit the dyke propagation velocity inferred by the hypocenter migration. Because of the scarce number of event locations before day 66.5, we focus our simulations on a time period starting at day 66.5, and we estimated an average propagation velocity of 1.82 m/min (0.11 km/hr) from day 66.5 to day 68.4, when the dyke reaches the upper part of the edifice (brown trend in Figure 5). We disregard the seismicity between day 65.5 and 66.5 as it may represent a different regime of dyke propagation (as proposed by Battaglia et al., 2005), or may be related to the preparatory process preceding the dyke nucleation, such as for instance the pressurization of the reservoir at 5–6 km depth.
<caption>**Table 1**</caption>

Parameters for All the Simulations Performed for the Piton de la Fournaise Case: Initial Tip Depth Below Sea Level ($z_0$), Initial Dyke Length (L), Rigidity ($\mu$), Rock Density ($\rho_r$), Magma Density ($\rho_m$), Viscosity ($\eta$), Bulk Modulus (K), Input Velocity ($v_0$), Cross-Sectional Area ($A_0$), Dislocation Length (l), Fracture Energy ($E_f$), Fracture Toughness ($K_c$) and Mean Velocity ($V_{mean}$)

<table><thead><tr><th>z<sub>0</sub> (km)</th><th>L (km)</th><th>μ (GPa)</th><th>ρ<sub>r</sub> (kg/m<sup>3</sup>)</th><th>ρ<sub>m</sub> (kg/m<sup>3</sup>)</th><th>η (Pa·s)</th><th>K (MPa)</th><th>v<sub>i</sub> (m/min)</th><th>A<sub>0</sub> (km<sup>2</sup>)</th><th>l (km)</th><th>E<sub>f</sub> (MPa·m)</th><th>K<sub>c</sub> (MPa·m<sup>1/2</sup>)</th><th>v<sub>mean</sub> (m/min)</th></tr></thead><tbody><tr><td rowspan="2">9.5</td><td rowspan="2">2.5</td><td rowspan="2">10</td><td rowspan="2">3,000</td><td rowspan="2">2,700</td><td rowspan="2">100</td><td rowspan="2">2.00E + 04</td><td rowspan="2">0.042</td><td rowspan="2">0.00075</td><td rowspan="2">0.05</td><td>1.5</td><td>200</td><td>0.78</td></tr><tr><td>2.2</td><td>242</td><td>0.09</td></tr><tr><td rowspan="3">9.5</td><td rowspan="3">3.5</td><td rowspan="3">10</td><td rowspan="3">3,000</td><td rowspan="3">2,700</td><td rowspan="3">100</td><td rowspan="3">2.00E + 04</td><td rowspan="3">0.042</td><td rowspan="3">0.0021</td><td rowspan="3">0.05</td><td>4.1</td><td>331</td><td>2.23</td></tr><tr><td>4.45</td><td>344</td><td>1.81</td></tr><tr><td>6.3</td><td>410</td><td>0.14</td></tr><tr><td rowspan="3">9.5</td><td rowspan="3">5</td><td rowspan="3">10</td><td rowspan="3">3,000</td><td rowspan="3">2,700</td><td rowspan="3">100</td><td rowspan="3">2.00E + 04</td><td rowspan="3">0.042</td><td rowspan="3">0.0062</td><td rowspan="3">0.05</td><td>12.9</td><td>587</td><td>12.71</td></tr><tr><td>17.1</td><td>675</td><td>1.83</td></tr><tr><td>18.7</td><td>706</td><td>0.23</td></tr><tr><td rowspan="2">9.5</td><td rowspan="2">6</td><td rowspan="2">10</td><td rowspan="2">3,000</td><td rowspan="2">2,700</td><td rowspan="2">100</td><td rowspan="2">2.00E + 04</td><td rowspan="2">0.042</td><td rowspan="2">0.0108</td><td rowspan="2">0.05</td><td>31.2</td><td>912</td><td>2.36</td></tr><tr><td>34.8</td><td>963</td><td>101.17</td></tr><tr><td rowspan="2">9.5</td><td rowspan="2">7</td><td rowspan="2">10</td><td rowspan="2">3,000</td><td rowspan="2">2,700</td><td rowspan="2">100</td><td rowspan="2">2.00E + 04</td><td rowspan="2">0.0006</td><td rowspan="2">0.017</td><td rowspan="2">0.05</td><td>50.7</td><td>1,163</td><td>3.46</td></tr><tr><td>51.8</td><td>1,175</td><td>0.80</td></tr><tr><td rowspan="2">6.8</td><td rowspan="2">2.5</td><td rowspan="2">10</td><td rowspan="2">3,000</td><td rowspan="2">2,700</td><td rowspan="2">100</td><td rowspan="2">2.00E + 04</td><td rowspan="2">0.042</td><td rowspan="2">0.00075</td><td rowspan="2">0.05</td><td>1.25</td><td>183</td><td>0.78</td></tr><tr><td>2.2</td><td>242</td><td>0.08</td></tr><tr><td rowspan="3">6.8</td><td rowspan="3">3.5</td><td rowspan="3">10</td><td rowspan="3">3,000</td><td rowspan="3">2,700</td><td rowspan="3">100</td><td rowspan="3">2.00E + 04</td><td rowspan="3">0.042</td><td rowspan="3">0.0021</td><td rowspan="3">0.05</td><td>4.1</td><td>331</td><td>2.91</td></tr><tr><td>5.4</td><td>379</td><td>1.39</td></tr><tr><td>6.3</td><td>410</td><td>0.07</td></tr><tr><td rowspan="3">6.8</td><td rowspan="3">5</td><td rowspan="3">10</td><td rowspan="3">3,000</td><td rowspan="3">2,700</td><td rowspan="3">100</td><td rowspan="3">2.00E + 04</td><td rowspan="3">0.042</td><td rowspan="3">0.0061</td><td rowspan="3">0.05</td><td>12.7</td><td>582</td><td>19.76</td></tr><tr><td>17.6</td><td>685</td><td>1.55</td></tr><tr><td>18.4</td><td>700</td><td>0.28</td></tr><tr><td rowspan="2">6.8</td><td rowspan="2">6</td><td rowspan="2">10</td><td rowspan="2">3,000</td><td rowspan="2">2,700</td><td rowspan="2">100</td><td rowspan="2">2.00E + 04</td><td rowspan="2">0.042</td><td rowspan="2">0.0082</td><td rowspan="2">0.05</td><td>31.9</td><td>922</td><td>2.76</td></tr><tr><td>34.5</td><td>959</td><td>151.02</td></tr><tr><td rowspan="2">6.8</td><td rowspan="2">7</td><td rowspan="2">10</td><td rowspan="2">3,000</td><td rowspan="2">2,700</td><td rowspan="2">100</td><td rowspan="2">2.00E + 04</td><td rowspan="2">0.0006</td><td rowspan="2">0.017</td><td rowspan="2">0.05</td><td>51.5</td><td>1,172</td><td>2.58</td></tr><tr><td>52.1</td><td>1,179</td><td>0.23</td></tr></tbody></table>

Note. $K_c$ is deduced from Equation 6 and $v_{mean}$ is the average dyke velocity excluding the initial phase of crack growth and the final acceleration due to the free surface.

## 6.2. Numerical Simulations of 1998 Magmatic Intrusion at PdF

We run several simulations considering a basaltic intrusion with a magma density of 2,700 kg/m³, viscosity of 100 Pa·s, crust rigidity of 20 GPa, Poisson's ratio of 0.25 and a rock density of 3,000 kg/m³ (Gailler et al., 2009), and we vary the initial dyke length between 2.5 and 7 km (Table 1). In addition, we tested the sensitivity of our simulations to the starting depth of the tip of the intrusion considering two different initial positions, at 9.5 and 6.8 km bsl (Figures 6a and 6b respectively), both consistent with a dyke fed by a deep reservoir, as the lower tip of the dyke would always be deeper than 10 km bsl. We set the starting time of our simulations, so that the dyke tip will be located at the depth of the deeper seismicity (~5 km bsl) when the high-rate seismic swarm began (day 66.5).

In Figure 6, we show the depth of the dyke tip as a function of time for our simulations. We display results for two different starting depths (Figures 6a and 6b), and for five different dyke lengths (solid colored lines), superposed with the hypocenter-depth migration (black spots).

The rock fracture toughness is one of the most weakly constrained parameters. Therefore, for each starting depth and dyke length, we tested different values of fracture toughness in order to produce the widest possible range of dyke propagation velocities (Table 1).

If we consider the hypocenter migration between -5 and 0 km bsl (before the final acceleration), we can identify a rather wide range of dyke lengths and rock fracture toughness that can reproduce the average dyke propagation
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Simulations of the 1998 dyke intrusion at Piton de la Fournaise volcano for two scenarios: (a) a deep initial crack at 9.5 km bsl and (b) a shallow initial crack at 6.8 km bsl. The volcano-tectonic events are represented by the black dots (after Battaglia et al., 2005, with permissions). The color code symbolizes the different crack lengths used in the numerical simulations, from 2.5 to 7 km. Dyke shapes at the initial stage of the propagation are plotted with an opening-to-length ratio of 0.0027 (the horizontal bar shows the opening scale). Colored areas show the domains of fracture toughness $K_c$ for which the associated dyke is propagating. Outside the plot frame, values of $K_c$ are displayed for each boundary of all propagation domains. The thick colored lines are the best fitting models.</figcaption></figure>

velocity: for the deeper starting depth, dyke lengths range between 3.5 and 6 km, with fracture toughness between 344 and 912 MPa·m$^{1/2}$, respectively (thick curves on Figure 6a), and for the shallower starting depth, lengths range between 2.5 and 6 km, with fracture toughness from 183 to 922 MPa·m$^{1/2}$ (Figure 6b). Note that these values of $K_c$ are all compatible with general estimates of $K_c^{eff}$ from field studies (Rivalta et al., 2015).

Our model also provides an acceleration when the intrusion approaches the surface, which depends on the dyke length. Therefore, considering the final acceleration of hypocenter migration over the last 2 km of dyke propagation, provides a further constraint on dyke length, and indicates that the best fit is obtained for dyke lengths of 6 km, and fracture toughness values of 918 and 922 MPa·m$^{1/2}$ for the deep and shallow starting depths, respectively. These dyke lengths are representative of the dyke-head length during propagation, as shown in Figure 7.

In Figures 7a–7d, we show the final crack shapes for the simulations relative to the thick colored lines in Figure 6a. For these sets of parameters, we obtain average tail thicknesses which are in the order of $10^{-1}$–$10^{-2}$ times the maximum dyke opening. Particularly, for our best fit solution (Figure 7c) the average tail over maximum dyke opening is 0.007. In Figures 7a–7d, we also show the static dyke shapes (dashed lines): in comparison with the dynamic dyke shapes, they display thicker heads in their uppermost part, and maximum openings which are about 20%–40% larger. Eventually, in Figure 7e, we display the crack velocities (and the velocity intervals that have been tested) for the simulations with different dyke lengths.

## 7. Conclusions

Our model of magmatic dyke propagation produces a dynamic shape consisting of an open tail region, which had been so far an exclusive characteristic of the “lubrication theory” school of models, overhung by a teardrop-shaped head, characteristic of the “Weertman crack” school of models. Within the range of viscosities and velocities typical of magmatic dykes, we have shown that they propagate in a regime where viscous forces are never fully negligible.

Although dyking events are of great interest and are generally analyzed by using both geodetic and seismic data, their propagation at mid crustal depths (larger than 4–5 km), may be hard to capture, and therefore are rarely documented, such that estimates of their velocity are often limited to the last part of propagation. The 1998 eruption at PdF is a great example of a vertical magma transfer from a relatively deep reservoir
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Final dyke shapes for the best fitting simulations performed with a deep initial dyke tip (thick color lines in Figure 6a) with respect to dyke lengths. The dashed colored lines show the static crack shape for each corresponding parameter set. The opening-over-length ratio is 0.0017. (a) 3.5 km, (b) 5 km, (c) 6 km, and (d) 7 km long dykes. (e) Mean velocity estimated for each Piton de la Fournaise (PdF) simulation, accounting only for the most constant part of the propagation. $v_{vt}$ is the mean velocity estimated between day 66.5 and 68.5 from the distribution of volcano-tectonic events (propagation below sea level in Figure 5).</figcaption></figure>

to the surface. The use of the dynamic model to characterize this dyke intrusion has provided additional information about the dyke shape, including length (3.5–6 km) and opening (0.37–1.24 m). We also determine a range for the fracture toughness of the medium, between 344 and 912 MPa·m$^{1/2}$ which are closer to field estimates than laboratory values, in agreement with Rivalta et al. (2015). Besides, this application has shown that the thickness of a dyke is not constant and that it can vary by one or two orders of magnitudes from head to tail. This implies that dyke thicknesses measured in the field can vary significantly even though the magma composition does not; conversely, models show that different sets of parameters and viscosities can lead to identical thicknesses, albeit at different depths. This means that the dyke thickness depends on its dynamics and where it is measured, which unfortunately is something that cannot usually be determined in the field.

Finally, we would like to stress that even though we have mostly focused our attention on the vertical ascent in a homogeneous crust, our modeling scheme has been implemented in order to account for the effect of heterogeneous crustal stress and layering (Figure 3f), and in particular to compute complex (non-straight) dyke trajectories (Maccaferri et al., 2011). In the future, we are going to use our new dynamic model to investigate non-vertically ascending dykes in interaction with complex stress configurations.

## Notation

### Rock parameters

$K_c$ fracture toughness

$\mu$ shear modulus

$\nu$ Poisson's ratio

$\rho_r$ density of the rock
Fluid-filled crack parameters

<table><tr><td>h, q</td><td>fracture opening/slip</td></tr><tr><td>K</td><td>fluid compressibility</td></tr><tr><td>u = [0, u<sub>s</sub>]</td><td>fluid flow velocity</td></tr><tr><td>v</td><td>propagation velocity of the crack</td></tr><tr><td>v<sub>th</sub></td><td>theoretical velocity from Spence and Turcotte (1990)</td></tr><tr><td>η</td><td>viscosity of the fluid</td></tr><tr><td>ρ<sub>m</sub></td><td>density of the fluid</td></tr></table>

Model parameters

<table><tr><td>A</td><td>cross-sectional area of the magma body</td></tr><tr><td>A<sub>0</sub></td><td>reference cross-sectional area of the initial magma body at vanishing lithostatic stress</td></tr><tr><td>b</td><td>Burgers vector</td></tr><tr><td>D</td><td>the term linking linearly the energy release ΔE to the crack velocity v (Equation 27)</td></tr><tr><td>e<sub>mn</sub></td><td>strain rate tensor</td></tr><tr><td>ΔE</td><td>energy release</td></tr><tr><td>E<sub>f</sub></td><td>fracture energy</td></tr><tr><td>∂E<sub>f</sub>/∂t, ΔE<sub>v</sub></td><td>rate of viscous dissipation</td></tr><tr><td>∂G/∂t, ΔG</td><td>rate of gravitational energy variations</td></tr><tr><td>∂W/∂t, ΔW</td><td>rate of the potential strain energy variations</td></tr><tr><td>f</td><td>volumetric fluid flow/flux</td></tr><tr><td>f<sub>i</sub><sup>top</sup></td><td>flux crossing the upper boundary of each dislocation element during the time Δt</td></tr><tr><td>f<sub>i</sub></td><td>flux at the middle point of ith dislocation element</td></tr><tr><td>F</td><td>the body force per unit volume</td></tr><tr><td>g</td><td>gravitational acceleration</td></tr><tr><td>h<sub>i</sub></td><td>average fracture opening at each dislocation element between propagation step k − 1 and k</td></tr><tr><td>i</td><td>ith dislocation element</td></tr><tr><td>k</td><td>propagation step</td></tr><tr><td>l</td><td>dislocation length</td></tr><tr><td>L</td><td>initial crack length</td></tr><tr><td>L<sub>w</sub></td><td>Weertman length scale from Rubin (1998)</td></tr><tr><td>n<sub>α</sub></td><td>number of directions tested at each propagation step</td></tr><tr><td>N</td><td>total number of dislocation elements</td></tr><tr><td>[T<sup>i</sup> T<sup>s</sup>]</td><td>influence coefficients matrix (Green functions) of the crack problem</td></tr><tr><td>P&#773;</td><td>pressure scale from Rubin (1998)</td></tr><tr><td>∂P<sub>visc</sub>/∂s</td><td>viscous pressure gradient (Equation 2)</td></tr><tr><td>ΔP, ΔP&#773;, ΔP<sub>dyn</sub></td><td>static, linearized static and dynamic overpressures</td></tr><tr><td>P<sub>fluid</sub></td><td>fluid pressure</td></tr><tr><td>P<sub>conf</sub></td><td>confining pressure</td></tr><tr><td>ΔP<sub>K</sub></td><td>pressure change due to the fluid compressibility</td></tr><tr><td>ΔP<sub>visc</sub></td><td>viscous pressure change</td></tr><tr><td>ΔP<sub>i</sub><sup>visc</sup></td><td>viscous pressure change due to ith dislocation element</td></tr><tr><td>ΔP<sub>i</sub><sup>visc</sup></td><td>viscous pressure change at the middle point the ith dislocation element</td></tr><tr><td>r, s</td><td>coordinate perpendicular to the fracture wall, coordinate along the crack cross-section</td></tr><tr><td>t</td><td>time</td></tr><tr><td>Δt</td><td>time elapsing between consecutive propagation steps</td></tr><tr><td>t*</td><td>non-dimension time from Spence and Turcotte (1990)</td></tr><tr><td>u<sub>i</sub></td><td>parabolic fluid velocity profile within each dislocation element i</td></tr><tr><td>v<sub>0</sub></td><td>input crack velocity</td></tr><tr><td>v<sub>mean</sub></td><td>average velocity excluding the initial phase and the final acceleration (see Table 1)</td></tr></table>
<table><tr><td>V</td><td>fluid-filled crack volume</td></tr><tr><td>w</td><td>dyke half thickness from Rubin (1998)</td></tr><tr><td>W</td><td>thickness scale from Rubin (1998)</td></tr><tr><td>ΔV</td><td>fluid-filled crack volume variations between consecutive propagation steps</td></tr><tr><td>ΔV<sub>0</sub></td><td>reference volume at vanishing lithostatic pressure</td></tr><tr><td>x, z</td><td>horizontal/vertical axis</td></tr><tr><td>z*</td><td>depth of the crack tip normalized by the initial depth</td></tr><tr><td>z<sub>0</sub></td><td>initial depth</td></tr><tr><td>α</td><td>angle between the vertical and the orientation of the dislocation element</td></tr><tr><td>β<sub>k</sub>, β<sub>v</sub></td><td>fracture and viscous coefficient for asymptotic tip solutions (Dontsov & Peirce, 2015)</td></tr><tr><td>Γ</td><td>viscous pressure gradient</td></tr><tr><td>δα</td><td>angle between each n<sub>α</sub> tested directions</td></tr><tr><td>Δρ</td><td>density contrast between fluid and medium (equivalent to the buoyancy)</td></tr><tr><td>Σ</td><td>dislocation surface</td></tr><tr><td>τ</td><td>shear traction</td></tr><tr><td>Φ</td><td>volumetric contribution of the rate of energy dissipated by the fluid flow</td></tr><tr><td>X</td><td>distance from a point inside the fracture to the tip from Dontsov and Peirce (2015)</td></tr><tr><td>Ψ</td><td>term linking linearly the velocity v to the fluid flux f (Equation 18)</td></tr></table>

## Data Availability Statement

The numerical code presented in this study is available on a repository: https://doi.org/10.5281/zenodo.7118734.

## Acknowledgments

This work received financial support from the ANR-DFG NLE 2018 MagmaPropagator projet (ANR-18-CE92-0037) and the project LOVE-CF (Progetti Dipartimentali, DIP8, INGV). The authors thank Jean Battaglia for providing the data from 1998 intrusion at Piton de la Fournaise and Aline Peltier for discussions on parameters for basaltic intrusions used in this work. The authors also thank the Editor Douglas Schmitt, the Associate Editor, Thierry Menand, Weronika Goczyk, Caitlin Chalk and two anonymous reviewers for their valuable feedback and insightful comments. Eventually, the authors would like to thank Timothy Davis for his informal review, and for the fruitful discussions and ideas that we have been regularly exchanging about this topic. Open Access Funding provided by Istituto Nazionale di Geofisica e Vulcanologia within the CRUI-CARE Agreement.

## References

Batchelor, G. K. (1967). *An introduction to fluid dynamics*. Cambridge University Press.

Battaglia, J., Ferrazzini, V., Staudacher, T., Aki, K., & Cheminée, J.-L. (2005). Pre-eruptive migration of earthquakes at the Piton de la Fournaise volcano (Réunion Island). *Geophysical Journal International*, 161(2), 549–558. https://doi.org/10.1111/j.1365-246X.2005.02606.x

Beauducel, F., Peltier, A., Villié, A., & Suryanto, W. (2020). Mechanical imaging of a volcano plumbing system from GNSS unsupervised modeling. *Geophysical Research Letters*, 47(17). https://doi.org/10.1029/2020GL089419

Bonafede, M., & Rivalta, E. (1999). The tensile dislocation problem in a layered elastic medium. *Geophysical Journal International*, 136(2), 341–356. https://doi.org/10.1046/j.1365-246X.1999.00645.x

Dahm, T. (2000a). Numerical simulations of the propagation path and the arrest of fluid-filled fractures in the Earth. *Geophysical Journal International*, 141(3), 623–638. https://doi.org/10.1046/j.1365-246x.2000.00102.x

Dahm, T. (2000b). On the shape and velocity of fluid-filled fractures in the Earth. *Geophysical Journal International*, 142(1), 181–192. https://doi.org/10.1046/j.1365-246x.2000.00148.x

Davis, T., Bagnardi, M., Lundgren, P., & Rivalta, E. (2021). Extreme curvature of shallow magma pathways controlled by competing stresses: Insights from the 2018 Sierra Negra eruption. *Geophysical Research Letters*, 48(13). https://doi.org/10.1029/2021GL093038

Davis, T., Rivalta, E., & Dahm, T. (2020). Critical fluid injection volumes for uncontrolled fracture ascent. *Geophysical Research Letters*, 47(14). https://doi.org/10.1029/2020GL087774

Davis, T., Rivalta, E., Smittarello, D., & Katz, R. F. (2023). Ascent rates of 3-D fractures driven by a finite batch of buoyant fluid. *Journal of Fluid Mechanics*, 954, A12. https://doi.org/10.1017/jfm.2022.986

Detournay, E. (2016). Mechanics of hydraulic fractures. *Annual Review of Fluid Mechanics*, 48(1), 311–339. https://doi.org/10.1146/annurev-fluid-010814-014736

Detournay, E., & Napier, J. A. (2019). Simulation of buoyancy-driven fracture propagation using the displacement discontinuity boundary element method. In A. Zingoni (Ed.) *Advances in engineering materials, structures and systems* (pp. 525–530). CRC Press/Balkema.

Dontsov, E. V., & Peirce, A. P. (2015). A non-singular integral equation formulation to analyze multiscale behavior in semi-infinite hydraulic fractures. *Journal of Fluid Mechanics*, 781, R1. https://doi.org/10.1017/jfm.2015.451

Duputel, Z., Lengliné, O., & Ferrazzini, V. (2019). Constraining spatiotemporal characteristics of magma migration at Piton de la Fournaise volcano from pre-eruptive seismicity. *Geophysical Research Letters*, 46(1), 119–127. https://doi.org/10.1029/2018GL080895

Gailler, L.-S., Lénat, J.-F., Lambert, M., Levieux, G., Villeneuve, N., & Froger, J.-L. (2009). Gravity structure of Piton de la Fournaise volcano and inferred mass transfer during the 2007 crisis. *Journal of Volcanology and Geothermal Research*, 184(1–2), 31–48. https://doi.org/10.1016/j.jvolgeores.2009.01.024

Heimisson, E. R., Hooper, A., & Sigmundsson, F. (2015). Forecasting the path of a laterally propagating dike. *Journal of Geophysical Research: Solid Earth*, 120(12), 8774–8792. https://doi.org/10.1002/2015JB012402

Jaeger, J. C., Cook, N. G. W., & Zimmerman, R. W. (2007). *Fundamentals of rock mechanics* (4th ed.). Blackwell Pub.

Lecampion, B., Bunger, A., & Zhang, X. (2018). Numerical methods for hydraulic fracture propagation: A review of recent trends. *Journal of Natural Gas Science and Engineering*, 49, 66–83. https://doi.org/10.1016/j.jngse.2017.10.012

Lengliné, O., Duputel, Z., & Ferrazzini, V. (2016). Uncovering the hidden signature of a magmatic recharge at Piton de la Fournaise volcano using small earthquakes. *Geophysical Research Letters*, 43(9), 4255–4262. https://doi.org/10.1002/2016GL068383

Lengliné, O., Duputel, Z., & Okubo, P. (2021). Tracking dike propagation leading to the 2018 Kîlauea eruption. *Earth and Planetary Science Letters*, 553, 116653. https://doi.org/10.1016/j.epsl.2020.116653

Lister, J. R. (1990). Buoyancy-driven fluid fracture: The effects of material toughness and of low-viscosity precursors. *Journal of Fluid Mechanics*, 210, 263–280. https://doi.org/10.1017/S0022112090001288
Maccaferri, F., Bonafede, M., & Rivalta, E. (2010). A numerical model of dyke propagation in layered elastic media. *Geophysical Journal International*, **180**(3), 1107–1123. https://doi.org/10.1111/j.1365-246X.2009.04495.x

Maccaferri, F., Bonafede, M., & Rivalta, E. (2011). A quantitative study of the mechanisms governing dike propagation, dike arrest, and sill formation. *Journal of Volcanology and Geothermal Research*, **208**(1–2), 39–50. https://doi.org/10.1016/j.jvolgeores.2011.09.001

Maccaferri, F., Rivalta, E., Passarelli, L., & Aoki, Y. (2016). On the mechanisms governing dike arrest: Insight from the 2000 Miyakejima dike injection. *Earth and Planetary Science Letters*, **434**, 64–74. https://doi.org/10.1016/j.epsl.2015.11.024

Mériaux, C., & Jaupart, C. (1998). Dike propagation through an elastic plate. *Journal of Geophysical Research: Solid Earth*, **103**(B8), 18295–18314. https://doi.org/10.1029/98JB00905

Napier, J., & Detournay, E. (2020). An unstructured mesh algorithm for simulation of hydraulic fracture. *Journal of Computational Physics*, **419**, 109691. https://doi.org/10.1016/j.jcp.2020.109691

Pedersen, R., Sigmundsson, F., & Einarsson, P. (2007). Controlling factors on earthquake swarms associated with magmatic intrusions; constraints from Iceland. *Journal of Volcanology and Geothermal Research*, **162**(1), 73–80. https://doi.org/10.1016/j.jvolgeores.2006.12.010

Peltier, A., Ferrazzini, V., Staudacher, T., & Bachèlery, P. (2005). Imaging the dynamics of dyke propagation prior to the 2000–2003 flank eruptions at Piton de la Fournaise, Reunion Island. *Geophysical Research Letters*, **32**(22). https://doi.org/10.1029/2005GL023720

Peltier, A., Staudacher, T., & Bachèlery, P. (2007). Constraints on magma transfers and structures involved in the 2003 activity at Piton de la Fournaise from displacement data. *Journal of Geophysical Research: Solid Earth*, **112**(B3), B03207. https://doi.org/10.1029/2006JB004379

Peltier, A., Staudacher, T., Bachèlery, P., & Cayol, V. (2009). Formation of the April 2007 caldera collapse at Piton de la Fournaise volcano: Insights from GPS data. *Journal of Volcanology and Geothermal Research*, **184**(1–2), 152–163. https://doi.org/10.1016/j.jvolgeores.2008.09.009

Pinel, V., Carrara, A., Maccaferri, F., Rivalta, E., & Corbi, F. (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), 1107–1125. https://doi.org/10.1002/2016JB013630

Pinel, V., Furst, S., Maccaferri, F., & Smittarello, D. (2022). Buoyancy versus local stress field control on the velocity of magma propagation: Insight from analog and numerical modeling. *Frontiers of Earth Science*, **32**. https://doi.org/10.3389/feart.2022.838318

Pinel, V., & Jaupart, C. (2000). The effect of edifice load on magma ascent beneath a volcano. *Philosophical Transactions of the Royal Society of London, Series A: Mathematical, Physical and Engineering Sciences*, **358**(1770), 1515–1532. https://doi.org/10.1098/rsta.2000.0601

Rivalta, E., Böttinger, M., & Dahm, T. (2005). Buoyancy-driven fracture ascent: Experiments in layered gelatine. *Journal of Volcanology and Geothermal Research*, **144**(1), 273–285. https://doi.org/10.1016/j.jvolgeores.2004.11.030

Rivalta, E., & Dahm, T. (2006). Acceleration of buoyancy-driven fractures and magmatic dikes beneath the free surface. *Geophysical Journal International*, **166**(3), 1424–1439. https://doi.org/10.1111/j.1365-246X.2006.02962.x

Rivalta, E., Mangiavillano, W., & Bonafede, M. (2002). The edge dislocation problem in a layered elastic medium. *Geophysical Journal International*, **149**(2), 508–523. https://doi.org/10.1046/j.1365-246X.2002.01649.x

Rivalta, E., Taisne, B., Bunger, A., & Katz, R. (2015). A review of mechanical models of dike propagation: Schools of thought, results and future directions. *Tectonophysics*, **638**, 1–42. https://doi.org/10.1016/j.tecto.2014.10.003

Roper, S. M., & Lister, J. R. (2007). Buoyancy-driven crack propagation: The limit of large fracture toughness. *Journal of Fluid Mechanics*, **580**, 359–380. https://doi.org/10.1017/S0022112007005472

Rout, G., Peltier, A., Taisne, B., Staudacher, T., Ferrazzini, V., & Di Muro, A. (2012). A new comprehensive classification of the Piton de la Fournaise activity spanning the 1985–2010 period. Search and analysis of short-term precursors from a broad-band seismological station. *Journal of Volcanology and Geothermal Research*, **241**, 78–104. https://doi.org/10.1016/j.jvolgeores.2012.06.012

Rubin, A. M. (1998). Dike ascent in partially molten rock. *Journal of Geophysical Research: Solid Earth*, **103**(B9), 20901–20919. https://doi.org/10.1029/98JB01349

Sigmundsson, F., Hooper, A., Hreinsdóttir, S., Vögfjörð, K. S., Ófeigsson, B. G., Heimisson, E. R., et al. (2015). Segmented lateral dyke growth in a rifting event at Bárarbunga volcanic system, Iceland. *Nature*, **517**(7533), 191–195. https://doi.org/10.1038/nature14111

Smittarello, D., Cayol, V., Pinel, V., Froger, J.-L., Peltier, A., & Dumont, Q. (2019). Combining InSAR and GNSS to track magma transport at basaltic volcanoes. *Remote Sensing*, **11**(19), 2236. https://doi.org/10.3390/rs11192236

Spence, D. A., & Turcotte, D. L. (1990). Buoyancy-driven magma fracture: A mechanism for ascent through the lithosphere and the emplacement of diamonds. *Journal of Geophysical Research*, **95**(B4), 5133. https://doi.org/10.1029/JB095iB04p05133

Traversa, P., Pinel, V., & Grasso, J. R. (2010). A constant influx model for dike propagation: Implications for magma reservoir dynamics. *Journal of Geophysical Research*, **115**(B1), B01201. https://doi.org/10.1029/2009JB006559

Watanabe, T., Masuyama, T., Nagaoka, K., & Tahara, T. (2002). Analog experiments on magma-filled cracks: Competition between external stresses and internal pressure. *Earth Planets and Space*, **54**(12), 1247–1261. https://doi.org/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), 1171–1183. https://doi.org/10.1029/jb076i005p01171

Zia, H., & Lecampion, B. (2020). PyFrac: A planar 3D hydraulic fracture simulator. *Computer Physics Communications*, **255**, 107368. https://doi.org/10.1016/j.cpc.2020.107368