In this study, we use the 15 experiments described in Maccaferri et al. (2019). Whereas this previous study focused on the description of the parameters influencing the crack path, we here focus on the description and interpretation of the velocity of the propagating cracks. A detailed description of the experiments is provided in Maccaferri et al. (2019), here we summarise the experimental protocol and provide new information on the velocity measurements. We performed experiments of air-filled crack propagation into a transparent brittle-elastic gelatin block, whose stress field is perturbed by a load applied at the surface (Figure 1 for the experimental setup description). A rectangular Plexiglas tank of 40 (length) × 20 (depth) × 20 (height) cm was filled with 16 L of liquid gelatin (concentration of 2% by weight) and was put in a fridge for 20 h at a temperature of 5°C, where it solidified while remaining in a fully hydrostatic state of stress. We assume that the Poisson’s ratio (ν) for gelatin is 0.5 (Kavanagh et al., 2013). The gelatin rigidity (Young’s modulus, E) was quantified for each tank by measuring the maximum vertical displacement at the surface due to the loading. This displacement is interpreted using a 3D Finite Element Model of the strain and stress induced inside the gelatin block by the rigid load, as detailed in Smittarello et al. (2021). We estimated E to be 2,150 ± 230 Pa with no significant variations between various experiments. At room temperature, we injected a controlled volume of air with a syringe from a given hole at the base of the tank (Figure 1). Several injections (up to 4) were performed successively in the same tank, using different holes and ensuring that the new paths were not affected by the previous ones. A sheet-like, air-filled crack formed and started propagating upward due to buoyancy. When the intrusion reached a chosen depth (Z _s ), the loading mass was put on the surface of the gelatin block, at the chosen horizontal distance from the air-filled crack (X _s ). The load applied at the surface was a metal plate with rectangular base of 6 (length) × 14 (depth) cm and mass ranging between 25.4 and 262.9 g (the load applied for each experiment is given in Table 1). In order to match the experimental setting used by Watanabe et al. (2002), we tried to set the starting depth and horizontal distance at, respectively, 2.7 times and 2.8 times the half-length of the load (i.e., 8.1 and 8.4 cm). We applied the load manually based on a visual estimation of the distance. In addition, the application of the load sometimes induces vibrations inside the gelatin for a few seconds and we start the analysis after the gelatin has reached a state of equilibrium. It follows that in practice, both Z _s and X _s were slightly different from the prescribed values and experimental values that are reported in Table 1. This load applied at the surface induces a heterogeneous stress field within the underlying gelatin. A compressive stress is created under the load, the magnitude of which varies with the magnitude of the applied pressure and is slightly influenced by the position of the load relative to the tank walls.

Three perpendicular cameras recorded the fluid-filled crack shape and path. Two spotlights illuminated the tank from the back and right sides (Figure 1). At the end of each experiment, we took several pictures of a ruler at different locations inside the gelatin block, in order to calculate the calibration factor F needed to scale the videos. The Software TRACKER is used to measure the dimensions of the air-filled cracks and to extract their velocities and trajectories. The whole propagation is monitored, from the initiation of the crack up to the surface, but we hereafter consider only the propagation after adding the load on the gelatin. The positions of the crack tip as a function of time obtained with TRACKER are analysed with MATLAB to derive both the horizontal and vertical component of the velocity along the propagation path. Outliers are removed and velocity is filtered using a one dimensional median filter. In fact, first velocity estimations may show some outliers due to perturbations induced by applying the load; however the velocity does not vary significantly in the first centimetre following this initial perturbation. We use this part of the velocity evolution (grey rectangle in Figure 2) to estimate a reference value of the velocity characterising the behaviour of the vertical crack (Figure 2).

The three dimensional stress field inside the gelatin is calculated for each experiment using the exact position of the load with COMSOL Finite Element model in order to take into account the effect of the rigid walls of the tank (zero displacement condition applied to the lateral and bottom boundaries of the gelatin to reproduce the adherence of the gelatin to the tank walls). We use a mesh made of about 330,000 tetrahedral units, refined in a vertical plane centred below the load as well as on the upper surface around the load (maximum size of the mesh is set to 2 mm on these surfaces). The upper surface is considered as a free surface except where the load is applied. We apply a condition of zero horizontal displacement to the lateral edge of a rigid plate (rigidity of 10⁹ Pa and thickness of 4 mm) to simulate the load [see Smittarello et al. (2021) for detailed explanations.]. In practice, performing a 2D calculation in plane strain approximation for the vertical plane centred below the load, results in a similar stress field estimation with increased resolution, such that in the following we use the 2D approximation for stress field estimations along the crack path. Comparison between 3D and 2D stress field is presented in Supplementary Figure S1 of Supplementary Material.

To complement our dataset of an air-filled crack initially rising at some lateral distance from the applied load, we also use the data published by Watanabe et al. (2002), which provide the velocity of oil-filled cracks rising just below the center of a surface load. These experiments were performed in a rectangular tank (58.5 cm in length, 26 cm in depth and 35 cm in height), the applied load had a rectangular base of 9 (length) × 20 (depth) cm and the rigidity of the gelatin is reported to be 270 Pa.

The reference velocity (Figure 2) corresponding to the vertical ascent in absence of stress field perturbation is reported in Table 1 for all the experiments considered. This velocity depends on the density contrast between the injected fluid and the gelatin, the elastic properties of the gelatin and the crack length (Takada, 1990; Heimpel and Olson, 1994). When considering the 15 experiments of air injection, characterized by the same gelatin physical properties, this velocity is proportional to the volume of fluid injected (Supplementary Figure S2 of supplementary material), which is consistent with the finding of Takada (1990), the volume being proportional to the crack length to the power of four (Smittarello et al., 2021).

Figure 3 shows the relative velocity variations with respect to the reference velocity along the path followed by cracks propagating inside the heterogeneous stress field induced by the load applied at the surface. When a load is applied with a lateral offset over a vertically ascending crack, this crack tends to be deflected towards the load. This deflection increases with the ratio between the amplitude of the applied load and fluid pressure (Watanabe et al., 2002; Maccaferri et al., 2019) and is more pronounced for shorter cracks (Maccaferri et al., 2019). On Figure 3 we can observe a significant velocity decrease along the two most deflected trajectories (experiments EXP61 and EXP34). The same velocity decrease is also seen for the crack rising vertically below the load as previously described by Watanabe et al. (2002). For the crack ascending below the load, the observed velocity decrease is explained by the crack propagation inside the compressive stress field induced by the load as modelled by Pinel et al. (2017), whereas for the deflected ones this velocity decrease occurs at a normalized depth around 1, that is to say well before the crack reaches the compressive stress field (represented in dark brown on Figure 3). A significant velocity increase is also observed close to the surface for several of the trajectories reaching the surface at some distance from the load: experiments EXP60, EXP45 and EXP31, for instance, show a clear acceleration, consistently with the free-surface effect described by Rivalta et al. (2005) and Rivalta and Dahm (2006). However, not all the experiments show this effect (cf. EXP49 in Figure 3). This may be due to the stiffening of the gelatin surface which may have occurred because of water evaporation (Kavanagh et al., 2013). This process is expected to affect only the surface of the gelatin, it does not induce any change in the Young’s modulus of the gelatin at depth (below a few mm from the surface) so that no significant change is recorded in the Young’s modulus measurements by surface loading.

In order to investigate the cause of the velocity decrease observed for cracks rising vertically below the load and for cracks deflected towards the load, we simulated the crack ascent for the vertically ascending crack of Watanabe et al. (2002) and for our experiment EXP61 (cf. Tab. 1, WAT and EXP61, respectively). Figure 4 shows the evolution of the buoyancy (Figure 4B) and of the normal component of the stress field induced by the load (Figure 4C) along the experimental paths (trajectories 1 and 4 in Figure 4A). Figure 4B shows that buoyancy is constant along the vertical crack trajectory, and it always decreases along the trajectory followed by experiment EXP61. This decrease is due to the dip of the trajectory, and reaches its maximum gradient at a normalized depth around 0.5. Just beneath the surface, the dip of the trajectory being close to 45°, the buoyancy is reduced by 50%. Regarding the normal component of the stress field (Figure 4C), it decreases continuously along the trajectory of experiment EXP61 due to the fact that the crack remains outside the most compressive area. This decrease is expected to favor the crack opening and its propagation. On the contrary, the normal compressive stress increases along the trajectory for the crack ascending vertically below the load. We can thus conclude that whereas the velocity decrease observed in experiment WAT is clearly due to the compressive stress field for the vertical crack, which corresponds to the effect modeled in Pinel et al. (2017), for the deflected experimental path (experiment EXP61) the velocity decrease recorded only results from the reduction of buoyancy induced by the horizontalization of the magma trajectory. In this case, the effect of the stress field is only indirect, inducing deflection of the trajectory, which leads to a reduction in buoyancy. Note that for a crack rising vertically away from the surface load, no change in buoyancy or significant change in the stress field occurs along the path, so the velocity does not change significantly (see EXP29 in Figure 3). We performed the same analysis along two additional trajectories relative to experiment EXP61: one was obtained using the BE model to compute the energetically-preferred path, and the other was obtained following the direction perpendicular to σ ₃ and aligned with σ ₁ (which is widely used as proxy for the propagation path of hydrofractures). The results are shown in Figure 4 and labeled as curves three and two, respectively. We always start with an initially vertical crack. When the BE crack follows the direction of σ ₁ (trajectory number two in Figure 4A), the upper tip of the crack will be subject to a very sharp deflection at the beginning of the propagation. This causes a faster decrease of the normalized buoyancy parameter with respect to the other simulations (Figure 4B). This fast decrease in σ _buoy (curve number two in Figure 4B) continues until the BE crack has fully entered the curved trajectory, or—in other worlds—when the lower tip of the crack reaches the starting depth of the upper tip, Z _s . After this point the buoyancy parameter starts to increase again, due to the overall verticalization of the crack orientation. Along the trajectory derived from the BE model (curve three in Figure 4B), the buoyancy decreases faster than in the experimental trajectory EXP61, and it starts to increase again before reaching the surface load (consistently with the trajectory from the BE model being deflected at deeper depth with respect to the experimental trajectory EXP61). The normal stress experienced by the crack approaching the loading (curve 3 Figure 4C) is also remarkably different from the one in experimental trajectory EXP61 (curve 4, Figure 4C): in fact, even if the BE trajectory ends approximately in the same location as the experiment EXP61, the dip angle of these two cracks are very different, hence the normal component of the stress field they experience are different. These differences are even larger when comparing the σ ₁ trajectory with the experimental trajectory EXP61, curves two in Figures 4A–C. Given these results, it appears to be important to use the parameters extracted along the experimental trajectory (rather that the ones obtained with the BE model or, worst, with the σ ₁ trajectory) in order to compare them with the velocity variations observed for the analog experiments. We will thus follow this strategy, and in the following we will always refer to the parameters extracted from the BE model results along the trajectory that was actually followed by the experimental crack.

The BE model is based on a quasi-static approach which cannot provide any insight into the velocity of propagation. However, it provides the energy release associated with the crack propagation. The energy release represents the excess between the energy gained at a propagation-step of the BE model (crack elongation), and the energy spent to fracture the host medium at the crack tip. However, a physical process should conserve the total energy of the system, so what is the meaning of such positive energy release? What is not considered in this model are all the possible forms of energy dissipation related with the dynamic of the fluid-filled fracture propagation process (fluid viscosity, plastic effects at the crack tip, elastic wave emission, and possibly more): the energy dissipation should equal the energy release to balance the total energy budget. This implies that the larger is the energy release associated with the crack propagation, the larger should be the energy dissipated by dynamic processes (which are velocity dependent). Therefore the energy release may correlate with the velocity of crack propagation. In order to check this hypothesis, we compare both values, energy release and velocity, for experiments showing a significant velocity variation along the trajectory: the vertical crack ascending below the load as well as the most deflected path (EXP61). We also perform this analysis on a less deflected crack showing a marked acceleration at shallow depth (EXP60) and a crack rising vertically away from the surface load with no significant velocity variation (EXP29). Results are presented in Figure 5. For experiments WAT and EXP61 the energy release (red line in Figures 5C,F) decreases when the crack rises towards the surface, which appears consistent with the observed velocity decrease (Figures 5B,E). However at shallow depth, there is an increase in the energy release, which is not observed in the velocity. This increase of the energy release is due to the free surface, which causes an effective reduction of the fracture toughness (Rivalta and Dahm, 2006). However, the free surface condition may not be realistic for a crack approaching a rigid load at surface. To correct for it, we run the BE simulations removing the free surface (circles in panels c and f of Figure 5), and in fact we obtain more consistent results for the energy release when compared with the recorded velocity. For EXP60, there is a slight velocity increase due the crack inclination but this time the most important feature in the velocity evolution is the strong velocity increase at shallow depth consistent with the energy release when considering the free surface (plain line in Figure 5I). For EXP29, no significant velocity variation is observed even if a velocity increase at shallow depth seems to be observed consistently with the free surface effect. From our numerical results, the free surface effect occurs at larger depth for longer vertical cracks (Supplementary Figure S3 of Supplementary Material. We also notice that for a given crack length, the free surface influence occurs at shallower depth for inclined cracks when the angle between the trajectory and the vertical becomes larger than 40° (Supplementary Figure S4 of Supplementary Material).

Finally, using the Pearson correlation coefficient, we were able to quantify the relation between velocity variations (Figure 6A) and energy release variations, the energy release being numerically calculated either considering the free surface or removing its effect (Figures 6B,C). Figure 6D shows correlation coefficients for all experimental trajectories, circles and inverted triangles being for energy release considering the free surface or removing it, respectively (note that the calculations performed without free surface are meant to better explain the experiments where the crack approaches the gelatin surface in correspondence of the loading plate). For most of the experiments performed, one of the Pearson coefficient is larger than 0.5, which indicates a positive linear correlation between the energy release along the trajectory and the velocity. For all cracks reaching the surface close to the load (at a distance from the load smaller than the length of the load, x* < 2), the Pearson correlation coefficient calculated is always larger when removing the free surface effect, which confirms that cracks reaching the surface below or close enough to the load experience a reduced free surface effect.

In addition to the BE model for inviscid fluids, we use a dynamic viscous fluid-filled fracture propagation model along the experimental trajectories (Pinel et al., 2017). This model provides velocity evolution along the trajectory based on the assumption that this velocity of the fracture tip is controlled by the Newtonian viscous fluid flow within the crack. We simulate the fracture propagation along the most deflected experimental path (EXP61) either taking into account the external stress field induced by the load or in absence of external stress field in order to separate the effect of the external stress and the trajectory inclination (solid and dotted lines in Figure 7). Based on the experimental record, we set Z _s to 7.88 cm. The fracture half-breadth a is set to 1.4 cm. In experiment EXP61 the driving pressure is given by σ _buoy and equal to 95 Pa while the load is 268 Pa, we accordingly set N ₂, which is the ratio between the driving fluid pressure and the load applied, to 0.35. We made simulations with N ₁ equal to 0, 1 and 3, considering no-buoyancy, buoyancy equal to the bottom overpressure, and buoyancy 3 times larger than the bottom overpressure, respectively. In Figure 7, we show the evolution of the dimensionless velocity [velocity scaled by [v] as defined in Eq. 6] along the analog path as a function of depth.

The curves in Figure 7 shift towards higher velocities for increasing values of N ₁ (from black to red in 7), consistently with larger buoyancy values producing faster fluid-filled cracks (given a constant total magma pressure). For all the simulations the velocity rapidly decreases during the first (deeper) part of the propagation path. In fact, this initial - rapid - velocity decrease is due to the bottom overpressure condition, and is not displayed in the experimental velocity measurements. Such effect is slightly smaller for larger N ₁ (i.e., when the contribution of buoyancy is larger), and indicates that the condition N ₁ = 3 is the most suitable to reproduce the velocity profiles observed in the air-filled crack experiments.

Focusing on the velocity decrease near the surface, for the simulations with N ₁ = 3, we get the most interesting features for our comparison with experiment EXP61: we notice that the velocities estimated without the loading stress field (red dotted line in Figure 7), overlap with the velocity obtained with the loading stress (red solid line in Figure 7). This proves that the shallow velocity decrease (also observed in the experiments) is fully due to the trajectory inclination, and is not directly due to the normal stress change along the path (which is negligible).

This is confirmed by the velocities obtained for fluid-filled fractures propagating vertically in a lithostatic stress field (dashed lines in Figure 7) which are larger than for the deflected path, the difference increasing with the buoyancy (larger differences for larger dimensionless numbers N ₁).