Modeling Dike Propagation in Both Vertical Length and Horizontal Breadth

We present analog experiments on dike propagation, followed by a numerical model of horizontal and vertical growth, which is partially analytical and partially based on empirical observations. Experimental results show that the growth rates are similar until buoyancy becomes significant and, afterward, vertical growth dominates. The numerical model is defined for different conditions in a homogeneous medium: (a) constant flux, fracture‐limited propagation; (b) constant flux, viscous‐limited propagation; and (c) variable flux dependent on the driving pressure and dike dimensions. These conditions distinguish between cases when the influx depends on the deeper source of magma (e.g., a conduit, independent of the dike geometry) and when it depends on the dike, so the influx can change as it grows. In all cases, the ratio of vertical to horizontal propagation is proportional to the ratio of buoyancy pressure to source pressure, in which buoyancy drives vertical propagation. We test the numerical model on dikes observed at Piton de la Fournaise, in which the dimensions were estimated using geodetic and seismic data. The results show that the final dimensions can be reproduced using magma‐crust density differences of 50–300 kg/m3, viscosities of 30–300 Pa·s, influxes of 50–750 m3/s and shear moduli of ∼10 GPa. The modeled magma and host rock parameters agree with previous studies of the volcano, while the flux is higher than what is typically observed during eruption. This implies a variable injection condition, in which the flux peaks during propagation and diminishes by the onset of eruption.

The dike's dimensions contain information of the dominant pressures, which may change as it grows, in that fracture pressure decreases and viscous pressure drop increases (Lai et al., 2016). Here we will define "length" as the primary direction of propagation (thus vertical), "breadth" as the secondary, transverse direction of propagation (thus horizontal), and "thickness" as the narrowest dimension (also horizontal). The elastic pressure always balances with one of these other pressures and is characterized via the aspect ratio (thickness to breadth; ). By neglecting the breadth dimension in a model, one is left to estimate pressure via the thickness and vertical length, potentially overestimating the aspect ratio and underestimating the pressure. By fixing the breadth to a constant value, we might improperly influence the pressure balance evolution. Furthermore, the length to breadth aspect ratio contains information of its driving pressures, so that horizontal growth diminishes as buoyancy surpasses the source pressure (Menand & Tait, 2002).
In this article, we lay out a framework for both vertical and horizontal dike growth, for a simple scenario of dike ascent in a homogeneous half-space. In this scenario, the dimensions all evolve with time, which enables assessment of the overall dike geometry, in a computationally efficient way. To do so, we review the basic theory on which numerical models are formulated. We then show analog experiments whose results indicate that this theory is upheld for certain conditions. Using these results, we develop a new analytical model for dike growth in vertical length and horizontal breadth. Finally, we test this model on a series of dikes from Piton de la Fournaise, which erupted from 2000 to 2003.

Background Theory
Dike propagation is under the influence of both the magma, which can resist flow if it is too viscous, and the surrounding solid, which can resist fracture and deformation if it is too strong or too rigid. Propagation can be driven by source pressure (ΔP) and/or buoyancy pressure (P b ). Resistance to propagation is quantified by the fracture pressure (P f ), the elastic pressure (P e ), and the viscous pressure drop (P v ) (Traversa et al., 2010). These pressure scales are largely dependent on dike geometry and defined as (Dahm, 2000;Lister, 1990;Lister & Kerr, 1991;Rivalta et al., 2005;: The symbols represent the liquid to medium density contrast (Δρ), gravitational acceleration (g), the dike vertical, horizontal and thickness characteristic dimensions (L, B, and H), fracture toughness of the host rock (K c ), shear modulus (G), Poisson's ratio (ν), liquid dynamic viscosity (μ), and influx (Q). Please see the final table (Table 1) at the end of the document for a full list of symbols. The denominators in Equations 2 and 3 are the smaller of the L and B dimensions (here it is B, for vertically propagating dikes), which characterizes the magnitude of their 3 of 23 respective stresses, following . Note the factor of 12 in Equation 4 is a constant associated with fissure flow. The K c is related to G via the Young's modulus, E, via K c = (2Eγ) 1/2 (Menand & Tait, 2002) and E = 2G(1 + ν), where γ is the surface energy, which has an approximate value of 1 J/m 2 in analog experiments with gelatin (Kavanagh et al., 2013), although this can be reduced by adding salt during preparation (Smittarello et al., 2021). We assume gelatin can be approximated as an incompressible material and therefore has a Poisson's ratio of 0.5 (Pansino & Taisne, 2020).
Subsets of these pressures can dominate the propagation behavior, due to the physical properties of the magma or host rock and lead to different analytical models (Rivalta et al., 2014). For example, the Weertman model assumes fracture dominated propagation for a dike of constant volume, which ascends only due to buoyancy and must overcome the host rock rigidity and fracture toughness (Weertman, 1971), a bit like a flat teardrop ascending through the crust. Once a dike grows to a critical length scale, L bf , here referred to as the buoyancy length, it can self-propagate without any additional pressure or magma supply. This length is defined by the relationship (Dahm, 2000;Menand & Tait, 2002), The maximum thickness of this dike, H ef , is derived by a balance of P e and P f , such that, (Weertman, 1980). For G on the order of 30 GPa, K c of 0.3 MPa·m 1/2 and B of 1 km, this estimates a thickness of about 0.1 mm, much thinner than typically observed, on the order of 0.1-1 m for basaltic dikes (e.g., Geishi et al., 2012;Grasso & Bachèlery, 1995;Lénat et al., 2012).
The Weertman model traditionally assumes all the liquid remains in a discrete structure at the top of the fracture (the dike head) and, as the dike propagates, liquid leaves the lower, tail region. As Dahm (2000) points out, this is not physically sensible and would require an unreasonably large pressure gradient in the tail. This is resolved if some magma remains in the tail, causing the dike to gradually lose volume as it propagates. Volume loss can reduce the length of the dike head to less than the buoyancy length, so that propagation can only resume when magma migrates up into the head from the tail region.  found that under such conditions, the dike length grows proportionally to time, t, such that L ∼[ΔρgV 2 t/(4μB 2 )] 1/3 , assuming the volume, V, is V ∼LBH (their Equation 8). For constant magma rheology, volume and breadth, the propagation rate, ∂L/∂t, is also time dependent, so that ∂L/∂t ∼ΔρgV 2 /(4μB 2 )] 1/3 /(3t 2/3 ), assuming ∂B/∂t = 0. In length-dependent form, we can write this constant-volume (CV) velocity scale, ∂L/∂t CV , as, which is equivalent to buoyancy driven Poiseuille flow. The velocity therefore drops as the dike grows and thins.
Another model for a viscosity dominated dike is described by lubrication theory, in which a constant influx from the magma's source complements   (Rivalta et al., 2014;Roper & Lister, 2007). The propagation velocity of a constant-flux (CF) dike, ∂L/∂t CF , depends on the thickness of the dike tail, H bv , which in turn depends on the balance of buoyancy and influx to viscous pressure drop (H when P b ∼P v ; Roper and Lister (2007), their Equation 2.4), Note that Roper and Lister (2007) defined these per using the 2D flux (m 2 /s), which we have adapted assuming a negligible flux gradient along the dike Q 2D ≈ Q/B. The dike propagates at a constant velocity as influx at the base balances with growth at the leading tip. Volatile exsolution and fragmentation can be then incorporated to modify the dike's velocity and shape, in that exsolution and expansion cause the dike to accelerate and thin, whereas fragmentation causes the dike to expand and therefore decelerate .

Experimental Methods
We performed five scaled analog experiments (conditions shown in Table 2) using gelatin as the crustal analog and either vegetable oil or water as the magma analog. Gelatin is commonly used in such experiments, as it is brittle and elastic on the timescale of a propagating dike (Di Giuseppe et al., 2009;Kavanagh et al., 2013) and it is transparent, so we can observe the growth in real time.

Gelatin Preparation
We prepared the gelatins by dissolving 3-3.5 wt% of granular gelatin in a pot of deionized water heated to 60°C. We used Xiamen Yasin brand gelatin (250 bloom, 8-15 mesh), which does not guarantee a specific animal source, but is either pig-or bovine-skin. However, for our purposes, the source is inconsequential as Pang et al. (2014) measured the strength of bovine-skin gelatin and found it has similar storage and loss moduli to the pig-skin gelatin measured by Di Giuseppe et al. (2009), for similar concentrations and temperatures. We poured this gelatin solution into a 50 cm × 50 cm × 50 cm clear, acrylic tank (Figure 1), topped it up with more deionized water to a total volume ≥100 L (see Table 2 for volumes) and mixed well to ensure homogeneity. We then covered it with a layer of oil to prevent evaporation (and formation of a tough skin on top) and placed the tank into a cold room, to cure at 15°C for 3 days. This procedure results in a gelatin whose rheology scales to nature, with a dominantly brittle, elastic behavior (Kavanagh et al., 2013;Menand & Tait, 2002;Pansino & Taisne, 2020).
The gelatin's strength strongly controls dike propagation and defines how an experiment scales to nature. The resulting strength depends on a variety of factors, including gelatin concentration, curing time, salt content, pH, etc. (Brizzi et al., 2016;Kavanagh et al., 2013). For example, the values listed in Table 2 vary due to decreasing water quality from our water purifier, which was due for a filter change. We therefore estimated the Young's modulus of each gelatin by measuring the shear wave velocity, which is visible in gelatin using polarizing filters . Shear wave velocity, v s , is related to the shear modulus, G = ρv s 2 (Lee et al., 2017), where ρ is the gelatin density. We do not have density measurements of the gelatins, but we estimate the density of a similar sample of 3.3 wt% gelatin, by measuring the mass (to the nearest 0.1 g) and volume in a beaker. We used a photo to precisely measure volume, assuming the location of the solid surface meniscus is correctly identified within an error of 20 pixels, or approximately 10 mL of a 400 mL container. The corresponding density was 1,014 ± 5 kg/m 3 , but we assume all gelatins to have a density of 1,010 kg/m 3 for  simplicity. The Young's modulus and fracture toughness can then be estimated using Equations 6 and 7, assuming a ν of 0.5 and γ of 1 J/m 2 (Kavanagh et al., 2013;Pansino & Taisne, 2020).
The gelatin density also plays a role in the liquid-solid density contrast, which affects dike buoyancy. We approximate the vegetable oil and water to have respective densities of 910 and 1,000 kg/m 3 , so the corresponding density contrasts are 100 and 10 kg/m 3 . As noted, the actual solid density may be marginally higher, so the true density contrast could be similarly a bit higher.

Experimental Procedures
We performed two types of experiments, one driven at a constant flux by a peristaltic pump ( Figure 1a) and one driven by gravity ( Figure 1b). In either scenario, the supply of analog liquid was connected to the base of the tank using flexible silicone tubing. For pump-driven experiments, we drastically varied the flow rate between experiments to achieve a wide range of fluxes. For gravity driven experiments, we placed the liquid supply reservoir near the ceiling of the cold room, achieving a pressure head of 1.25 m, with a pressure decrease of less than 1% over the course of an experiment (<1 cm drop of the liquid level). The reservoir was placed on scale, which recorded mass against time, allowing us to determine the flux at any given time. We note that for experiment 4, with vegetable oil, the supply tubing did not have a sufficiently large diameter to mitigate viscous pressure drop within the tube. Gravitational forces balanced with the pressure drop, resulting in a constant influx condition in the experiment. In experiment 5, with water, we observed that the fluid could freely enter the dike with a highly variable rate, which we interpret is due to the evolving magnitudes of P f , P e , and P v (Equations 2-4). We recognize that this causes experiment 5 to different to the others in a number of ways, however as we will see, the results are consistent with how it should theoretically behave.
To initiate an experiment, we first used a retractable blade to make an initial cut in the gelatin to control the orientation of the ascending dike, so that it was vertical and perpendicular to our camera. We then injected the Re 10 −4 -10 −2 10 −3 -10 −1 10 −5 -10 −2 10 −4 -10 −2 10 −1 -10 1 10 −10 -10 −2 Note. Scaling parameters are given in italics. a Taisne and Jaupart (2009). b Wada (1994). c Meredith and Atkinson (1985). d Traversa et al. (2010). e Assuming a range of sizes. Lower estimates in experiments and nature given to estimate dimensionless numbers when the dike is small.  (Figure 1c). Note the retractable knife was contained in a small basal reservoir attached to the bottom of the experimental tank, visible in Figure 1. It is designed to be heated for temperature-dependent experiments, not presented in this study (see . The experiments were recorded using a video camera positioned at the side of the tank, recording at 10 Hz frequency. This setup allowed precise recording of the dikes' vertical and horizontal dimensions with time, which were picked automatically using MATLAB. We estimated the average dike thickness via the volume (as we know the influx and time of each video frame) and other dimensions, H ≈ 6V/(πLB) (from V∼ (4π/3) (L/2) (B/2) (H/2) = πLBH/6, for an ellipsoidal crack). Note this volume approximation follows the assumption of an elliptical cross section, which occurs during uniform loading (Pollard, 1976). It offers an order-of-magnitude relationship between volume and thickness.

Experimental Results
Experimental results are summarized in Figure 2. Measurements in the top row indicate the total length, maximum breadth, and calculated average thickness, such that H ≈ 6V/(πLB). The derivatives in the bottom row are estimated over a moving window of 10 data points. These are shown against time, normalized to the final measurement time, t/t f , which for experiments 1-4 is picked to be a few moments before eruption, such that the dike is a few centimeters below the surface and not undergoing an acceleration due to the free surface (as described by Rivalta and Dahm (2006)). For experiment 5, the influx was unintentionally restricted midway through the experiment, so the final time is the moment before this occurred. We observed that the growth rate in both the L and B dimensions was strongly dependent on the influx. Each dike initially grew radially, but dikes with vegetable oil transitioned to dominantly vertical growth as the dike approached the buoyancy length. Such transition has previously been observed and attributed to the liquid buoyancy (Menand & Tait, 2002). We note that the growth never became purely vertical, as each dike was subject to a continuous influx of liquid that could not be completely accommodated by vertical growth. As experiment 5, made with water, never approached the buoyancy length, it maintained radial growth throughout the propagation. We used a scale to record a time series of mass of the liquid in the supply beaker and estimate the volumetric flux. (c) We made an initial precut into the gelatin to control the orientation of the dike. The dike's growth, which we tracked with a video camera, is generally radial at first and then becomes more vertical when buoyant. (d) An image from experiment 2, which has a small en-echelon dike to the right. Image processing subtracts the background to make the dike stand out. Note particles embedded in the liquid are used for velocimetry, not discussed in this study. The edges of the dike are apparent and can be automatically tracked to record dimensions. (e) Definitions of dike dimensions and growth rates. The ∂B/∂t totals the growth on both sides of the dike.

Scaling
A dike's relative size can be characterized by the buoyancy length (Equation 5), which indicates if it is large enough to propagate without any additional source pressure. At this point, the propagation becomes dominantly vertical and the dike assumes a tabular shape (Menand & Tait, 2002). The dimensionless length, L bf *, can be defined as the length relative to the buoyancy length, * = . (10) When the buoyancy pressure is less than the fracture pressure (L bf * <1), the dike needs an additional driving pressure from the source to propagate. When this happens, the propagation is radial, so that the dike maintains a penny shape (Menand & Tait, 2002;. In our experiments, L bf * tends to have a value less than 2, whereas for nature it can be one to two orders of magnitude longer. In this sense, our experiments capture the initial stages of propagation for dikes with Δρ ∼100 kg/m 3 or the full propagation of dikes with Δρ ∼1 kg/m 3 (see Table 2 for approximate scaling parameters).
The influx dynamic can be quantified by the Reynolds number, Re, which is the ratio of inertial pressure, here defined by ρv 2 /L, to viscous pressure drop, defined by μv/H 2 (Schlichting & Gersten, 2017), Re = (ρv 2 /L)/(μv/H 2 ). For an influx Q ≈ vHB, this leads to In our experiments, and likely nature as well, propagating dikes are dominated by laminar flow due to their long, thin nature, so that Re << 1. In our experiments, the Reynolds number is relatively high at the start (10 −2 < Re < 10), as the viscous pressure is initially small in magnitude, and diminishes as it grows (10 −5 < Re < 10 −1 ). We would speculate that basaltic dikes behave similarly, assuming that they are not extremely narrow when they first initiate.
We will finally introduce a length scale to characterize dikes that are limited by viscous pressure drop versus those that are limited by fracture resistance, L fv . Menand and Tait (2002) show that dikes in gelatin tend to have low viscous pressure drop, which can increase as the dike grows, and much higher fracture pressure, which In the middle row, the growth rates are shown. Note that Experiment 2 has a small jump the ∂B/∂t measurements, as it had a secondary, en-echelon fracture due to the high flux and viscosity, which was eventually overcome by the main dike. In the bottom row are growth rates, relative to the current size (e.g., ∂L* = (∂L/∂t)/L).
decreases with dike growth. Since dikes are expected to be fracture dominated when small and viscously dominated when larger, the length scale comes from when these pressures are in balance with the elastic pressure, so that P f ∼ P e ∼ P v (respectively Equations 2-4). The H can be isolated by the balance of P e and P v , denoted H ev , Note that for a basaltic dike with parameters shown in The dimensionless length, relative to this scale, L fv *, is, * = .
For dikes in gelatin experiments (parameter values in Table 2), L fv is on the order of 30-3,000 m, depending on the viscosity and influx, much longer than our experiment size. This indicates they are effectively controlled by fracture pressure, in agreement with arguments made by Menand and Tait (2002).
Further evidence for these controlling pressures lies in the aspect ratios. For example, for the range of parameters given in Table 2, the H ef scale (Equation 6), predicts experimental dike aspect ratios H ef /B of 0.06-0.2, in strong agreement with our lab measurements, whereas H ev /L predicts 0.002 to 0.06. Note the reference dimension, L or B changes depending on the thickness scale being used. In contrast, for magmatic dikes, L fv <<1 m, indicating they immediately become viscously dominated. The H ev scale predicts an aspect ratio H ev /L of 7 × 10 −4 to 4 × 10 −6 , agreeing with typical field measurements of 10 −3 to 10 −4 (Kavanagh et al., 2013), whereas H ef /B predicts 2 × 10 −7 to 7 × 10 −7 .

Lubrication Theory Velocity Scale
We compare ∂L/∂t measurements against the propagation velocity scale ∂L/∂t CF defined by Roper and Lister (2007;Equation 9), for a constant influx dike, for benchmarking purposes. The other previously introduced velocity scale (Equation 7), presented by , was previously benchmarked and validated in their paper. The ∂L/∂t CF is derived for steady 1D growth and therefore overestimates ∂L/∂t measurements (Figure 3), which  (2007), ∂L/∂t CF , compared with measurements, ∂L/∂t. Symbols are as in Figure 2.
As dikes (a) become influenced by viscous pressure drop (P v /P f is large) and (b) grow primarily in one direction (Π V is small), the velocity scale converges with measurements. (c) The 3D view is oriented to show that these parameters are fairly confined to a log-planar surface (illustrated from side), and thus inter-related.
share influx with horizontal growth. The ∂L/∂t CF becomes more appropriate to describe experimental measurements as P v /P f approaches a value of 1 (i.e., as P f becomes negligible; Figure 3a) and as Π V becomes small (i.e., horizontal growth is negligible; Figure 3b). In such cases, dikes assume nearly vertical propagation and experimental conditions resemble the assumptions made to derive the velocity scale.

Experiment Boundary Conditions
We note here that the boundary conditions in our experimental setup are rigid tank walls and a free surface on top. Gelatin is an incompressible solid, meaning that the intrusion of a dike generates a small uplift of the surface, as well as a compressional stress field in the horizontal directions. It has been shown that for surface loading, boundary conditions have a non-negligible effect when the load diameter is greater than 5% of the container size (Kavanagh et al., 2013;Smittarello et al., 2021). However, dikes can generate their own stress field and ignore external stresses that are relatively small (Watanabe et al., 2002). We therefore quantify the far field horizontal compressional stresses, σ h , via the compressional strain, s, and the Young's modulus, so that σ h = Es. The strain is approximated via dike volume and the gelatin volume (100 L). For the vertical component, σ v , we estimate the stress via the uplift, h, and gelatin density, ρ, so that σ v = ρgh. The average uplift is approximately the dike volume divided by gelatin surface area (0.25 m 2 ). By quantifying the internal pressure as ΔP + P b , and comparing it to σ h and σ v , we can see that the external stresses are never significant ( Figure 4).  Models of ∂L/∂t, normalized by measurements, against dimensionless length, L bf *. Note colors and symbols are different from (a) as different models are being compared, using all experimental data for each model. All ∂L/∂t models perform best when the dike is small and there is significant horizontal growth. As the dike becomes long, the model underestimates ∂L/∂t. (c) Similar plots for ∂B/∂t.

Form of Growth and Numerical Model Description
We will define a numerical model based in part on observations from the analog experiments, in which the form of growth depends on buoyancy. We first introduce the ratio of horizontal to vertical growth, Π V , which approaches a value of 1 for radial growth and a value of 0 for vertical growth. For the buoyancy, we introduce the ratio of buoyancy to source pressure, Π P , We quantify the effective source pressure following Menand and Tait (2002, their equation 14), so that ΔP ∼[G(Δρg) 2 Q(∂L/∂t) −1 (1 − ν) −1 ] 1/3 . These parameters correlate to each other so that Π V approaches either 1 or 0, respectively for Π P << 1 and Π P >> 1 (Figure 5a). We describe the relationship using an error function of the form and solve for the constants C 1 and C 2 by maximizing the R 2 fit (values shown in Figure 5a). Note that the function has a maximum value of 1 to indicate purely radial growth, though data from experiment 5 somewhat exceeds this value.
We then define a few relationships to describe dike growth under different assumptions. Dikes are usually simplified to have either a constant-pressure or constant-influx source condition. However, in either case, it may be more helpful to think of the dikes as pressure-or influx-controlled, rather than with a constant source condition.
In the pressure-controlled scenario, one might imagine that the source reservoir is highly pressurized and loses negligible volume (and therefore pressure) as the dike propagates, and that the dike is the limiting structure for magma transport. The thickness and flow rate strongly depend on the elastic resistance of the country rock and both the thickness and velocity are expected to increase as the dike propagates.
In the influx-controlled scenario, one might expect that the shallow source reservoir is small in size and is continuously fed by more-deeply-sourced magma, ascending along a deep pathway (i.e., a dike or conduit). In this case, the driving pressure balances with the confining characteristics of the deeper path (presumably it would have to be narrower) and the shallow dike is relatively easy to inflate and propagate. The thickness and the velocity can therefore be defined via the balance of other pressures (i.e., P e vs. P f and P e vs. P v ), and the pressure in the dike can change as it grows. Similar to the pressure-controlled condition, the thickness increases as the dike grows, however, the constant influx restricts the flow velocity, causing it to decelerate.
To define these relationships, we first relate the dike dimensions to the volume, assuming an ellipsoidal shape V = πLBH/6. The influx is thus, The change in dike volume therefore equals the influx, assuming negligible magma compressibility. For convenience, we will write relative growth rates using a ∂* notation (e.g., ∂L* = (∂L/∂t)/L, ∂V* = Q/V), so that, * = * + * + * .
For an influx-controlled dike, influenced primarily by the balance of elastic and fracture pressures, fracture pressure limits dike growth and the thickness adjusts accordingly. Viscous pressure drop is negligible, allowing magma to flow relatively easily and the effective driving pressure in the dike decreases as the dike grows. The dominant forces remain in balance so that P e ∼ P f and ∂P e /∂t ∼ ∂P f /∂t. Taking the derivatives of Equations 2 and 3, we have, By substituting Equations 2 and 3 into Equation 20, and since P e ∼ P f , * ∼ 1 2 * .
which implies that we cannot yet fully resolve the dike growth in ∂H* and ∂B* dimensions without first knowing ∂L*. However as shown in Figure 5a, this can be resolved via the ratio of buoyancy force and driving pressure. We will label the empirical formula in Equation 17 as ∂B/∂L = F(Π P ) and convert it into relative growth rates, * = (Π ) * , and so, * ∼ * and * ∼ * We construct a similar set of equations for an influx-controlled dike, primarily under the influence of elastic pressure and viscous pressure drop. The pressures remain in balance so that P e ∼ P v and ∂P e /∂t ∼ ∂P v /∂t, Setting P e ∼ P v , this can be simplified to, * ∼ 1 4 ( * + * ) .
For a constant influx dike, ∂Q* = 0, however we will leave it in this general form. By substituting Equations 24 and 29 into Equation 19, the dike growth is, * ∼ Finally, for a dike controlled by driving pressure, where P e ∼ ΔP + P b and ∂P e /∂t ∼ ∂ΔP/∂t + ∂P b /∂t, Substituting Equations 1 and 3 into Equation 33, * ∼ Δ ∕ + * + * , where ∂ΔP/∂t = 0 for a constant pressure dike, though we will let this term remain. Substituting Equation 34 into Equation 19 gives * ∼ 1 2 By setting Equation 24 equal to Equation 35, * ∼ and * ∼ We compare estimates of ∂L/∂t and ∂B/∂t from Equations 25, 26, 30, 31, 36 and 37 with lab measurements (Figures 5b and 5c) and find that all the models produce similar results, since they contain similar terms of ∂V* and F(Π P ), although Equations 36 and 37 perform more poorly for dikes longer than the buoyancy length. Note the magnitude of ΔP (estimated following Menand and Tait (2002)) evolves as the dikes grows and so ∂ΔP/∂t ≠ 0.
Examples using each of these sets of equations are given for illustrative purposes (Figure 6). To do so, we begin with a dike that has initial dimensions L = B = 20 m. The magnitude of H is estimated using either P e ∼ P f , P e ∼ P v or P e ∼ ΔP + P b , depending on the model (we show all three), and then the volume estimated via V ∼ πLBH/6. For the influx-controlled models, we use a constant influx of 10 m 3 /s, whereas for the pressure-controlled models, influx is calculated assuming Poiseuille flow, Q ∼ (ΔP + P b )BH 3 /(12μL). The growth rates over a time step are calculated using Equations 25-27, 30-32, and 36-38, depending on the model. The model then repeats for the next time step, until a final time of 10 4 s. The values of other parameters used to calculate pressures (e.g., Δρ, μ, G, etc.) are given in the caption for Figure 6; these parameters are chosen for illustrative purposes, so that all curves can be displayed on a single plot and their shapes can be compared. For the P f dominated case, we estimate the source pressure via ΔP + P b = P f , whereas for the P v dominated case, it is estimated via ΔP + P b = P v , with H defined by Equation 12, so that ΔP + P b = B −1 [12μQLG 3 (1 − ν) −3 ] 1/4 . The models show that constant influx dikes begin with high growth rates, which then decelerate with time and vice versa for constant pressure dikes. As such, constant pressure dikes continuously accelerate as the dike grow in breadth and thickness and influx increases exponentially toward infinity. This obviously would not occur in nature as the magma supply depletes.

Case Study of Piton de La Fournaise (2000-2003)
We test the framework presented in Section 5 for a case study from Piton de la Fournaise, following a similar methodology to the examples given in Figure 6. The model is implemented for constant flux dike propagation, which can exhibit a constant seismicity rate (Traversa & Grasso, 2009;Traversa et al., 2010). It assumes each dike propagates through a homogenous medium and is controlled by fracture pressure or viscous pressure drop, depending on its size. By beginning with a dike of some initial geometry, it is possible to estimate how it will evolve in length and breadth, since at any moment, the driving forces determine the proportion of vertical growth to horizontal. We test our model using data presented by Peltier et al. (2005) for nine dikes at Piton de la Fournaise, which erupted from February 2000 to September 2003. The dikes had two phases of propagation, a relatively fast, deeper phase of vertical propagation, from a reservoir located at sea level, followed by a slower phase of shallow, horizontal propagation, leading to flank eruption. The vertical phase qualitatively well-matches our experiments, whereas the horizontal phase does not (we did not observe dominantly horizontal growth), which likely occurred because the magma reached a discontinuity within the edifice, such as the level of neutral buoyancy (Fukushima et al., 2010) or a rheological change (Geishi et al., 2012).
We forward model each dike incrementally using an initial length (L 0 ) and breadth (B 0 ) of 20 m and, together with magma buoyancy, viscosity, influx and shear modulus, which we vary between models (parameters shown in Table 3). We allow a wide range of magnitudes for each parameter to see which produces the best model. As a dike grows, it can be controlled by the balance of P f and P e , in which case we estimate the thickness via Equation 6 (H ∼K c (1 − ν) (2B/pi) 1/2 G −1 ), or the balance of P v and P e , in which case we estimate it via Equation 12 (H ∼ (12μQL(1 − ν)/G) 1/4 ). The transition between force balances occurs when L ≥ L fv (Equation 13), which in practice occurs when the dike is very small, usually less than the initial length. The magnitude of K c is estimated similarly to experiments, K c = (2Eγ) 1/2 , however the surface energy varies more widely in nature than experiments, from ∼1 to 1,000 J/m 2 (Atkinson, 1984). Mafic rocks commonly are commonly measured to have 1 < K c < 3 MPa·m 1/2 (Balme et al., 2004;Meredith & Atkinson, 1983, 1985, whereas studies of dikes in nature suggest they behave as though they have a higher effective value, K c ∼ 100 MPa·m 1/2 (Rivalta et al., 2015). This is not yet resolved, so we elect to use laboratory measurements and, assuming these had a shear modulus of 30 GPa, then the associated surface energy is 6 < γ < 60 J/m 2 . We assume a constant value of γ = 10 J/m 2 for simplicity. The ΔP is estimated following Menand and Tait (2002), when L < L fv and Δ = 1 ( ) when L ≥ L fv . Equation 39 is obtained by solving ΔP + P b ∼ P v , with H defined by Equation 12, and in cases where ΔP < 0 (when buoyancy is large), we set ΔP = 0. We then solve for ∂L/∂t using either Equation 25 or Equation 30 (depending on whether L > L fv ) and ∂B/∂t using Equation 24 (∂B/∂t = F(Π P )∂L/∂t). The length and breadth during the next time step are estimated using the growth rates over a time step of Δt, which depends on Figure 6. Examples of the dike growth numerical models. Solid lines represent each model, with Δρ = 0 kg/m 3 , and dashed lines represent the same models, with a buoyancy specified in the legend. The parameters for influx-controlled, P f dominated models are G = 3 GPa and K c = 3MPa·m 1/2 ; parameters for influx-controlled, P v dominated models are G = 30 GPa and μ = 30 Pa·s; parameters for pressure-controlled models are ΔP = 0.3 MPa, G = 0.1 GPa, and μ = 100 Pa·s. These values are chosen so that each model has a similar value of H and that ∂B/∂L decreases around 500 s, allowing the velocity curves to be easily compared.
the reported duration of the propagation (Δt = t t /10,000) and is generally about 1 s. The model runs until the time reaches the duration reported by Peltier et al. (2005;t t in Table 3) To assess the accuracy of each run, we fit the final dimensions in Table 3. We use length measurements, L t , derived from Peltier et al. (2005), for a dike length starting at 200 m above sea level and ending at the minimum elevation of the eruptive fissure (erupting on the edifice flank). They estimate the dike length associated with the vertical phase of propagation as such, and the shallower segment of the dike propagated laterally (which we do not model). We assume the target breadth, B t , is 1 km for all dikes, based on seismic measurements for an analog dike at Piton de la Fournaise, presented by Battaglia et al. (2005). There is uncertainty on the dike thickness at depth, so we run multiple models for thicknesses of 0.3-3 m, typical values for basaltic dikes and corresponding to the range reported by Grasso and Bachèlery (1995) for Piton de la Fournaise. The base error, ε 0 , is then calculated by adding the errors of each dimension (ε L = L(t = t t )/L t , ε B = B(t = t t )/B t , etc.), via, which we minimize. Note we give more weight to the ε L as it is based on individual measurements for each dike, rather than an assumed value. The resulting best fitting models match the dimensions and velocities with an average error of ε L , ε B , and ε H of 0.99 ± 0.07, where a perfect match has a value of 1 (Figure 7).
At this point, we know there is some combination of the model inputs Δρ, μ, Q, and G, which produce a dike that matches observations, but we do not know if the magnitude of those inputs are reasonable. We apply an adjustment to the error matrix based on a priori knowledge, so that errors associated with extreme magnitudes are penalized by multiplying them by a factor greater than 1 (as the minimum error indicates the best-fitting value). The form of this penalty is described after Equation 41. For Δρ, Gailler et al. (2009) found that the relevant section of the edifice has a density of 2,500-2,600 kg/m 3 (∼500 m to at least 1,500 m depth, relative to the summit; Gailler et al. (2009)), whereas the overlying rock has a lower density of 1,600-1,800 kg/m 3 . These dikes exhibit a phase of vertical growth in the lower, denser section (which we model), followed by lateral growth at the upper section, indicating that the magma likely has a bulk density between the two (i.e., it is buoyant in one, but not the other; Fukushima et al. (2010)). We expect there is at least a small density difference between the magma and host rock and set a constraint of 10 < Δρ < 1,000 kg/m 3 .
The other parameters can be similarly constrained. Typical viscosity measurements are on the order of 30 < μ < 300 Pa·s at eruption temperatures of ∼1400 K (Edwards et al., 2020;Villeneuve et al., 2008), which we set as a constraint. The time-averaged flux for each eruption is reported to be 0.5 < Q < 24.5 m 3 /s (Roult et al., 2012), however short-term variations are common, typically with a high magnitude at the onset (10-100 m 3 /s), that tapers off with time (Hibert et al., 2015). We cannot be sure that the influx during propagation is higher than the time-average flux, but it is unlikely to be lower, so we set a constraint of Q > 1 m 3 /s.

Table 3
Forward Model Parameters expected to be ∼20 GPa for intact rock (Got et al., 2013), although the effective magnitude may be lower if the edifice was substantially damaged during previous events (Carrier et al., 2014). We therefore set a constraint of 2 < G < 20 GPa.
The total error, ε, is: where ε p is the a priori penalty factor (tailored to error distributions, but similar in methodology for probability density functions, as described in Roy and Romanowicz (2017)). For Δρ, μ, and G, the ε p has the form of an inverse, log-normal distribution (ε p = 1/exp{−[ln(x) − b] 2 /a 2 }, where a and b are constants), with a minimum value of 1 and distribution that penalizes the ranges specified above by less than 5% (Figure 8). The penalty curve for Q has an exponential distribution (ε p = 1 + exp[−ln(x)]/20), achieving the same penalty that favors magnitudes over the specified value. The penalties are small in comparison with the base error distribution, and in practice only takes effect when there is more than one minimum error, thereby choosing the most reasonable model. For more information on these curves, see Supporting Information S1.
We finally attain a set of 4D error distributions, for each target thickness and each eruption. It can be a bit challenging to interpret multidimensional data, especially when there exist trade-offs between the different parameters. For example, in Equation 12 there is a tradeoff between μ, Q, and G. We find it helpful to view these grids in two complementary ways: (a) extract cross sections of the grid, located at the global minimum that show the best combination of parameters, which we will refer to as the conditional error (examples shown in Figure 8); and (b) for each dimension, find the mean values along the other dimensions (i.e., quasi-integrate the grid from 4D to 3D, then 2D, then 1D) to determine which parameters generally produce good models, which we will refer to as the marginal error. For more information on how these errors are quantified, see Supporting Information S1.
For a target final thickness of 0.3 m, the marginal errors shows that good models can be generated over fairly wide range of parameters (Figure 9 blue curves; 10 < Δρ < 300 kg/m 3 , 5 < μ < 3,000 Pa·s, 20 < Q < 2,300 m 3 /s, 10 9 <G < 10 11 Pa), which sometimes overlap with the best models in the conditional errors (Figure 9, red curves; 50 < Δρ < 300 kg/m 3 , 20 < μ < 200 Pa·s, 75 < Q < 750 m 3 /s, 5 < G < 18 GPa). A strong case can be made for parameter magnitudes in these overlapping regions. This tends to be the case for the viscosity and influx parameters, whereas the buoyancy and shear modulus show less agreement. For these mismatching parameters, the conditional probabilities better match the penalty curves shown in Figure 8 and we therefore consider them more accurate.
Note that for larger target thicknesses, the results look similar to those shown in Figure 9, but shifted to favor a higher viscosity and flux, and lower shear modulus (for H = 3 m, the conditional errors indicate 500 < μ < 5,000 Pa·s, 500 < Q < 5,000 m 3 /s, 0.3 < G < 1 GPa). They generally have a poorer match with a priori information than the models with thinner dikes. We will therefore focus our discussion on the results for H = 0.3 m.
The results from our model for 50 < Δρ < 300 kg/m 3 and G ∼ 10 10 Pa are reasonable, given the a priori constraints. For viscosity, the marginal and conditional errors overlap in the range of 30 < μ < 300 Pa·s, a good match with the a priori range (Edwards et al., 2020;Villeneuve et al., 2008). For the flux, the marginal and conditional errors  Table 3. Each parameter matches observations when equal to 1 and shown against different target thicknesses, H t .
overlap in the range of 50 < Q < 750 m 3 /s, which is generally higher than observed eruption fluxes, even at the onset, which can approach Q ∼ 100 m 3 /s (Hibert et al., 2015). It is possible however, that the flux is even higher during propagation. From the perspective of mass conservation, a dike with a given cross sectional area and propagation velocity would be associated with a certain flux. A simple, first-order approximation, Q ≈ HB(∂L/∂t), indicates that a 0.3 thick dike, with a 1 km breadth, propagating at 1 m/s, would be associated with a flux of 300 m 3 /s. Thinner or less-broad dikes could maintain the same velocity with a proportionally smaller influx.
The 1D conditional errors show relatively consistent estimates for Δ , μ, and G between eruptions, whereas the influx Q can be seen to vary more widely (Figure 9). Although Villemant et al. (2009) found that the magma composition became slightly more primitive with each event, we do not note any similar monotonic variations in the modeled magma rheology. Variations in propagation velocity observed by Peltier et al. (2005) appear to be primarily due to influx, with no significant variations due to the magma and host rock rheology.

Comments on the Case Study Results
The model we present generates good matches for the observed dike dimensions presented in Table 3. As we previously discussed, we have made assumptions of the horizontal breadth and the thickness, based on evidence from previous studies. The dike thickness, in particular, is known to be dependent on several parameters, including magma viscosity, driving pressure and host rock resistance (Geishi et al., 2012;Gudmundsson, 2000;Wada, 1994). For Piton de la Fournaise, typical thickness measurements range from 0.5 to 2.0 m, depending on the host rock, in which poorly consolidated layers tend to host thicker dikes than cohesive layers (Grasso & Bachèlery, 1995;Geishi et al., 2012). In the Dolomieu crater, dikes are generally less than 1 m thick (Lénat et al., 2012).
Our model generates thicknesses that are dependent on the viscosity, shear modulus and influx. Thin dikes (<1 m) are associated to viscosities in the range of previous estimates (30 < μ < 300 Pa·s; Villeneuve et al., 2008;Edwards et al., 2020) and thick dikes are above this range. Although there is some variation in composition for these eruptions, which became slightly more primitive with each event (Villemant et al., 2009), we assume they all had viscosities within this range. Our best-fitting models do not have a similar monotonic trend in the viscosity estimates.
There is uncertainty on the dike thickness at depth, as deep dikes do not necessarily have a clear surface deformation signal associated with a large thickness (Fukushima et al., 2010). If we estimate the thickness using Equation 12 (assuming 30 < μ < 300 Pa·s, Q ∼ 100 m 3 /s, G ∼ 10 GPa, and 870 < L < 2,050 m), we find that 0.29 < H < 0.64 m, which is on the low end of typical observations at the surface (0.5 < H < 2 m; Grasso and Bachèlery (1995)). This discrepancy is small, but could be explained by external stresses, in which dikes at depth are under greater compression and therefore thinner than at the surface (Fittipaldi et al., 2019;Geishi et al., 2010).
For the breadth dimension, the maximum size depends on the buoyancy and source pressure and can be approximated via Equation 16 (Π P ∼ 1, L ∼ B ∼ ΔP/(Δρg)). Prior to this size, L and B remain similar in magnitude, whereas afterward dikes grow predominately vertically. If we assume the density difference between magma and host rock to be Δρ ∼ 300 kg/m 3 and the source pressure to be ΔP ∼ 10 MPa, then B t ∼3 km. In our case study, we assumed B t ∼1 km, based on seismicity presented by Battaglia et al. (2005), which if correct indicates that the buoyancy is higher, the source pressure is lower or that the pressure diminished during the injection (i.e., the available magma volume was small). We estimate ΔP via Equation 39, which initially needs to be high (∼100 MPa) to supply the influx and propagate the fracture, but then drops as the dike grows (∼1 MPa).

Comments on the Model
The model we proposed in Section 5 has the potential to constrain the pressures acting on the dike through geophysical data. The breadth relates to the ratio of buoyancy to source pressure. For influx-controlled dikes, the thickness relates to either the balance between the elastic and fracture pressure or between elastic and viscous Figure 9. The error associated with each input variable, for each event. The marginal error for each parameter is obtained by averaging along the other dimension, whereas the conditional error is obtained by extracting a cross section of the error grid, at the global minimum. The 1D errors are shown along the main diagonal and 2D errors in the opposing corners. The curves in the 2D plots show the contour at the minimum error plus the median absolute deviation, times 0.05 (i.e., min(ε) + 0.05mad(ε)). Note that since the marginal error for each parameter takes the average over the other dimensions, the bounds may not perfectly coincide between plots (e.g., Δρ in panel (b) through panel (d)). For each 2D marginal error plot, the corresponding conditional errors are shown in gray for comparison, and vice versa for the conditional error plots.
pressure drop. For a pressure-controlled dike, it depends on the elastic and driving pressures. As shown in Section 6.1, the dike dimensions can ultimately be used to constrain valuable information, including the magma rheology and influx.
We note that our model has been constructed for a dike with some inflow at the source. Those that have a negligible influx act as discrete, CV structures which elongate and on average thin with time . Equation 19 (∂V* ∼ ∂L* + ∂B* + ∂H*) is conceptually the same, so that when ∂V* = ∂B* = 0, then ∂L* ∼ −∂H*. However, in their current form given in Section 5, propagation is modeled to stop when influx stops.
In reality,  show that CV, liquid-filled fractures, larger than the critical length L bf (Equation 5) theoretically grow proportionally to t 1/3 , decelerating with time. For magmatic dikes, this presumably continues as long as the magma temperature is above the solidus. Such cooling has a characteristic time scale dependent on the dike thickness and magma thermal diffusivity . Alternatively, it could be due to an effective fracture toughness due to the dike shape, which has a curved upper tip. Taisne and Tait (2009) argue that, since the tail of the dike never fully closes, the overall dike thickness decreases as it grows and the elastic pressure diminishes. Horizontal propagation eventually ceases along the upper edge of the dike and therefore vertical propagation ceases as well.
We additionally note that, while the model is designed for radial to vertical dike propagation in a homogeneous medium, it can potentially be implemented for lateral dike propagation due to a stress gradient (e.g., edifice loading). In this case, the horizontal stress gradient takes the place of buoyancy as a directional pressure. When this is relatively small in comparison to ΔP, we expect radial growth to dominate, and when it is relatively large, we expect horizontal growth to dominate. This is fundamentally the same argument made in Equation 17, where ∂B/∂L ∝ erf(ΔP/P b ), and we expect some form of the relationship could be applied for lateral dike propagation (perhaps with different empirical constants), although experiments similar to those by Urbani et al. (2018) are required to confirm this. In contrast to buoyancy, which grows proportionally with length, the average lateral stress gradient depends on the confining pressure at both ends of the dike, as well as the horizontal extent, so the evolution of this term depends on the topography. We expect that once a dike propagates significantly far from the edifice, its average lateral driving pressure drops, favoring a component of vertical growth.
The simplicity of our model implies it is not appropriate as is, in the context of layering or stress heterogeneities, as we assume a simple ellipsoidal shape in a homogeneous medium. Rigidity and/or density contrasts between layers are known to cause dikes to respond by assuming varying thickness. Furthermore, propagation may transition from dominantly vertical to dominantly horizontal when a dike encounters a layer that is either too rigid to fracture or not dense enough to ascend through. Discretizing the dike into smaller components is one way around this.

Discretization
Our model assumes that the entire dike can be approximated to an ellipsoidal geometry, albeit with one very thin dimension, which allows for efficient processing. However, there is the potential to adapt it in a way to allow for discretization along its length axis, similar to the model presented by Pinel et al. (2017). We will briefly describe here one method to do so.
For simplicity, we will assume a constant influx source condition, controlled by viscous pressure drop. The local thickness at distance from the base, z, depends on the driving pressures and the local rigidity of the rock, H(z,t) . For a constant influx dike, the effective ΔP drops with time and depends on the viscous pressure drop in the overall structure. We will characterize it in a way similar to Equation 39, but for a variable geometry for each discrete element. For an element with length, Δz, the viscous force, F v (z), on an element is the viscous pressure drop multiplied by the element surface area, F v (z,t) = [12μQΔzB(z,t) −1 H(z,t) −3 ] [ΔzB(z,t)] = 12μQΔz 2 H(z,t) −3 . The total viscous force equals the sum of each element, divided by the total surface area: Next, we will work out the amount of the influx that is used to propagate the dike vertically. We assume the vertical propagation rate depends in part on the influx and dike geometry, and in part dependent on the driving pressures, so that [∂L/∂t]/L ∼ [Q/V]/[1 + exp(−Π P (t))]. The term QL/V represents a quantification of velocity, where V/L ∼ HB, but allows for variable magnitudes of H(z,t) and B(z,t). The magnitude of QL/V is scaled by the term [1 + exp(−Π P (t))], which depends on buoyancy and source pressure. When buoyancy exceeds source pressure, it approaches a value of 1 (all the liquid is used to propagate vertically) and when source pressure is much greater, it approaches 0.5 (the liquid is shared evenly between vertical and horizontal propagation). For our constant-influx setup, the QL/V term stipulates that the velocity due to the influx decreases as the dike grows, as a function of time, so that Q/V ∼ t −1 . Combining the terms described above, we get, Note that Equation 43 is tested against experimental data and shows order of magnitude agreement (Figure 10a). The influx rate entering the tip is quantified by [V(t)/L(t)][∂L/∂t(t)], so change in volume of the rest of the dike is The inflation and horizontal propagation both depend on ∂V/∂t(t) in the rest of the dike. Since this is calculated on an element-by-element basis, we assume it is divided evenly by the number of elements. The relative growth rate, ∂V*(z,t) depends on the element volume, for number of elements, n(t). For an element of constant length, its volume adjusts by horizontal propagation and inflation, ∂V*(z,t) ∼ ∂B*(z,t) + ∂H*(z,t) (following Equation 19, where ∂L*(z,t) = 0 s −1 ). Recall our ∂* notation is the parameter growth rate over its magnitude (e.g., ∂B* = (∂B/∂t)/B, ∂V* = Q/V). The element volume also adjusts to maintain a balance between the local driving pressure and elastic pressure, so that ΔP(t) + P b (z) ∼ P e (z,t) and ∂ΔP/∂t(t) + P b (z)∂L*(z,t) ∼ P e (z,t)[∂H*(z,t) − ∂B*(z,t)]. Since the element length does not change, ∂L*(z,t) = 0 and ∂H*(z,t) ∼ ∂B*(z,t) + (∂ΔP/∂t(t)/P e (z,t) ∼ ∂V*(z,t) − ∂B*(z,t)). We can finally solve for horizontal propagation, * ( ) ∼ 1 2 Since a fracture cannot shrink, if ∂B* < 0, we set ∂B* = 0. Next, we find the inflation, * ( ) ∼ * ( ) − * ( ).
We show a simple test run of this discretized model to show how the dimensions evolve with time ( Figure 10b through Figure 10g). The model runs similar to the example in Figure 6, starting with an initial B and L of 20 m, finding the ∂L* and ∂B* as described, and then updating the dimensions over a time step of 5 s, until a final time of 5,000 s. We run the model for two layers, 1,000 m thick each, with the lower layer having G = 10 GPa and the upper layer having G = 1 GPa (Figures 10b and 10c), and then rerun the model reversing the layers (Figures 10e  and 10f). We choose the other parameters so that μ = 100 Pa·s, Q = 1 m 3 /s, and Δρ = 100 kg/m 3 . In both cases, there is a discontinuity in the geometry at the boundary between layers, with the dike in the harder layer being thinner and broader. In the case with the overlying hard layer, the dike velocity drops and causes the magma to accumulate in the underlying region. In the case with the overlying soft layer, the vertical velocity undergoes periods of deceleration associated with radial growth and periods of constant velocity associated with primarily vertical growth. We do not observe an acceleration when the dike enters the softer layer, which has been observed in analog experiments (e.g., Rivalta et al., 2005;Urbani et al., 2018). This is because the velocity is defined by the influx and dike geometry and the elastic pressure is uniform along the dike. In reality, we would expect an elastic pressure difference between the two layers to cause a temporary acceleration, as the harder layer squeezes magma into the softer one.

Comments on the Experiments
Scaling analysis in Section 4.1 suggests that dikes in our experiments are dominated by the fracture pressure (L fv * << 1), whereas dikes in nature are dominated by viscous pressure drop (L fv * >> 1), in agreement with arguments made by Menand and Tait (2002). For better similarity between experiments and nature, future work should aim to achieve higher magnitudes of viscous pressure drop, so that the fracture pressure is comparatively small. To the best of our knowledge, this has not yet been done in gelatin experiments.
We suggest the use of more-viscous liquids, together with powerful pumps capable of applying a sustained driving pressure. We can use Equation 13 to estimate the liquid viscosity and flux needed to observe this transition on the length scale of an experiment. For a typical tank size (a few tens of cm large) and a transition at a length scale L fv ∼ 15 cm, with other experimental parameters Δρ ∼ 10 kg/m 3 , G ∼ 10 4 Pa, and K c ∝ G 1/2 , one would need a combined viscosity and flux of μQ ≥ 6 × 10 −4 Pa·m 3 , which could be achieved with silicone oils with μ ≥ 100 Pa·s and pump flow rates Q ≥ 6 mL/s. Increasing the gelatin's density, by mixing in additives like sugar or salt, can reduce necessary viscosity. However, such additives can also affect the other rheological properties, by decreasing G and increasing K c , effectively the opposite of what is wanted (Brizzi et al., 2016;Choi et al., 2004;Smittarello et al., 2021). Due to the way Equation 13 is structured, as well as the relationship between K c and G, higher values of μQ are needed in softer gelatins. Performing experiments in a larger tank with a larger L fv reduces these values, so that μQ ∝ L fv −1 (i.e., for L fv ∼ 1.5 m, Q ∼ 6 mL/s, then μ ∼ 10 Pa·s).
Alternatively, a similar effect may be achieved by decreasing the fracture toughness. This could be accomplished by pre-cutting the entire gelatin block, or by thermally weakening it with prolonged exposure to a sheet laser.

Conclusion
We found through analog experiments that dikes transitioned from radial growth to vertical growth once the buoyancy pressure exceeds the source pressure. Based on these findings, we propose a numerical model of dike propagation, which simulates growth in both vertical length and horizontal breadth, so that the change in volume depends on the influx. The proportion of growth in the two dimensions is dependent on the ratio of source pressure to buoyancy, as seen in the experiments. The model is defined for various source conditions: influx-controlled dikes dominated by fracture pressure; those dominated by viscous pressure drop; and pressure-controlled dikes.
We then applied the model to a case study at Piton de la Fournaise, modeling the final geometry of a set of dikes. There is some uncertainty on the dike geometry at depth and, when we model thin dikes (∼0.3 m), the best-fitting magma and host rock properties (buoyancy, viscosity, and shear modulus) are plausible for Piton de la Fournaise. The best-fitting influxes for all models are higher than typically observed, as the influx during propagation may be higher than during eruption.