Ascent rates of 3-D fractures driven by a finite batch of buoyant fluid

Abstract

Propagation of fluid-filled fractures by fluid buoyancy is important in a variety of settings, from magmatic dykes and veins to water-filled crevasses in glaciers. Industrial hydro-fracturing utilises fluid-driven fractures to increase the permeability of rock formations, but few studies have quantified the effect of buoyancy on fracture pathways in this context. Analytical approximations for the buoyant ascent rate facilitate observation-based inference of buoyant effects in natural and engineered systems. Such analysis exists for two-dimensional fractures, but real fractures are three-dimensional (3-D). Here we present novel analysis to predict the buoyant ascent speed of 3-D fractures containing a fixed-volume batch of fluid. We provide two estimates of the ascent rate: an upper limit applicable at early time, and an asymptotic estimate (proportional to $t^{-2/3}$) describing how the speed decays at late time. We infer and verify these predictions by comparison with numerical experiments across a range of scales and analogue experiments on liquid oil in solid gelatine. We find the ascent speed is a function of the fluid volume, density, viscosity and the elastic parameters of the host medium. Our approximate solutions predict the ascent rate of fluid-driven fractures across a broad parameter space, including cases of water injection in shale and magmatic dykes. Our results demonstrate that in the absence of barriers or fluid loss, both dykes and industrial hydro-fractures can ascend by buoyancy over a kilometre within a day. We infer that barriers and fluid loss must cause the arrest of ascending fractures in industrial settings.

1. Introduction

1.1. Motivation

Fluid-filled fractures that propagate due to the buoyancy of the fluid are found in many geological settings. These fractures are discussed in the literature under various names including gravity-driven hydraulic fractures (e.g. Salimzadeh, Zimmerman & Khalili 2020), buoyant, liquid-filled cracks (e.g. Dahm 2000; Taisne & Tait 2009) and fluid-filled fractures under the influence of weight contrasts (e.g. Pollard & Townsend 2018). Examples of this phenomenon include dykes and veins in the crust, as well as melt-water-filled crevasses in glaciers. Understanding the rate and extent of propagation of such fractures is critical for making useful, quantitative predictions. For example, the speed with which a batch of magma rises through the crust determines the duration of the warning period between geophysical detection of a dyking event and the volcanic eruption. Fracture ascent speed and distance also have implications for industrial operations. This process will affect the size and shape of the final fracture footprint (Peirce 2016; Salimzadeh et al. 2020).

Field and laboratory observations of fluid-driven fracture provide some indication of typical rates of ascent; for example, in basaltic systems this is of the order of millimetres per second to ${\sim }0.5\ {\rm m}\ {\rm s}^{-1}$ (Tolstoy et al. 2006; Mutch et al. 2019). Observations of supraglacial lakes suddenly draining through growing crevasses indicate similar propagation rates (Das et al. 2008; Stevens et al. 2015). Although there is evidence that water- and gas-filled fractures can ascend through the crust (Schultz 2016; Cartwright et al. 2021), this process has not been observed directly. Estimates based on geochemical analysis give ascent rates of ${\sim }0.01\unicode{x2013}0.1\ {\rm m}\ {\rm s}^{-1}$ (one kilometre per day) (Okamoto & Tsuchiya 2009). Despite these estimates and their relevance for industrial operations, many published analytical and numerical models of industrial hydraulic fracturing do not account for the effect of buoyancy. Hence the influence of buoyancy on the ascent of fluid-filled fractures in this context may appear negligible or non-existent (e.g. Perkins & Kern 1961; Nordgren 1972; Detournay 2016). Here we show that buoyancy is an important control and should not be neglected. In industrial settings, ascent of hydro-fractures can be terminated at interfaces such as stiffness and stress contrasts introduced by layering (Bunger & Lecampion 2017). Furthermore, fluid can leak off into the surrounding porous rock formation, leading to closure of the fracture (Detournay 2016). We discuss the existing literature on buoyancy-driven hydro-fractures in the next section, highlighting the conditions under which these ascend, and the speed of propagation.

1.2. Previous work and current aims

Analytical predictions of the speed and height of a fracture have previously been obtained from two-dimensional (2-D) models. The 2-D models assume plane-strain elasticity, where there is zero strain out of the plane of interest ($x$–$z$) such that the fracture has an infinite extent with invariance along its third dimension ($y$). The 2-D models cannot capture effects of the finite lateral extent in $y$ of all real fractures. Nonetheless they provide useful physical insight and can be accurate along a central axis if the fracture is sufficiently broad in the $y$ direction.

Two classes of 2-D buoyant fractures have been characterised analytically. These are categorised according to the fluid source, either a constant area-injection rate or a fixed finite area (here, area refers to the volume per unit length in $y$). Early workers obtained solutions to the problem of a constant injection rate (e.g. Lister & Kerr 1991; Roper & Lister 2007; Rivalta et al. 2015). In these studies, injection is along a line-source of infinite length, at a rate that has units of volume per second per distance along the source (${\rm m}^2\ {\rm s}^{-1}$). At later times in these models, the fracture head ascends buoyantly above a constant-aperture tail, which transmits a steady supply of fluid upwards to the head.

In the other fluid-source class, Spence & Turcotte (1990) estimated the time-dependent ascent speed of a fixed, finite-area batch of fluid. Unlike the case of constant injection rate, this model predicts a diminishing aperture of the fracture tail, slowing fluid supply to the head, and hence causing the rate of propagation to diminish with time and ascent distance.

Analytical approximations for both of these classes have been shown to be consistent with 2-D numerical solutions (Lister 1990b; Roper & Lister 2007; Dontsov & Peirce 2015). The numerics confirm that the elastic constants and material toughness influence the shape of the propagating head but not the ascent speed. The latter is determined by the width of the tail.

Three-dimensional (3-D) fractures, with finite extent in the $y$ direction, pose additional challenges to analysis. Workers have begun to investigate the case of constant rate of volume injection in 3-D (Möri & Lecampion 2022), proving approximate analytical solutions for the height as a function of time, the width and the aperture of the fracture. In contrast, there has been no analysis of the ascent speed of a 3-D fracture containing a fixed volume of fluid. And yet such analysis would be valuable in various contexts. In volcanology, where the volume of intrusions is often well-constrained by geodetic data, a 3-D solution would help to estimate parameters that drive dyke ascent through the crust. In industrial operations, both the source rock properties and volumes of fluid injected are well constrained, but currently there is no simple way of predicting how fast an injected batch of fluid will ascend by buoyancy.

Here we develop an analytical solution that predicts the size and ascent speed of a 3-D fracture driven by a finite volume of buoyant fluid. In developing this solution, our first aim is to provide an upper bound on the ascent speed of fractures that are propagating due to the buoyancy of a finite batch of fluid. Our second aim is to understand how this ascent speed decays with time, as in the 2-D solution of Spence & Turcotte (1990).

Our strategy is to use a state-of-the-art numerical simulator (Zia & Lecampion 2020) to produce 3-D solutions to the full, nonlinear equations. We treat these as a benchmark for our novel analytical results, to show their applicability. The manuscript is organised as follows. In the next section we review existing analytical approximations and introduce the simulator of Zia & Lecampion (2020). We use this to simulate the ascent of a finite-volume fluid batch, reviewing this simulation in detail to describe the problem. In § 3, we use these insights and proceed to derive equations describing the ascent speed of a 3-D, finite batch of fluid. In § 4, we use 3-D simulations to test the validity of these results across a range of scales. We then test the ability of our approximations to predict the ascent speed of oil-filled cracks in gelatine solids. Lastly, in § 5, we discuss the implications of our results in relation to selected case studies, detailing limitations and potential avenues of future research.

2. Background

In this section we first review approximate analytical solutions for 2-D and 3-D fractures, then consider numerical simulators that solve the full, nonlinear problem.

2.1. Analytical solutions

2.1.1. Critical lengths and volumes

A fluid-filled fracture will either ascend or descend, depending on the density of the fluid relative to that of the surrounding solid medium. This has been well documented in analogue experiments by comparing, for example, the injection of air and mercury into gelatine solids as described in Heimpel & Olson (1994). For simplicity here, we consider only positively buoyant (less dense) fluids, and hence we model only fracture ascent, noting that there is no loss of generality in the solutions derived.

Fluid-filled fractures are driven to ascend by a weight contrast $\Delta \gamma$ between the fluid and surrounding rock. More specifically, the buoyancy force is defined as the difference between the vertical gradient of the horizontal stress in the rock column and the volume-specific weight of the magma contained within a vertical fracture. Assuming the horizontal and vertical stresses are equal, i.e. a lithostatic state of stress, then $\Delta \gamma =(\rho _{r}-\rho _{f})g$ (Secor & Pollard 1975), where $\rho _r$ and $\rho _f$ are rock and fluid density, respectively, and $g$ is the acceleration due to gravity.

Weertman (1971) and Secor & Pollard (1975) showed that fluid-filled fractures can ascend through their solid host by hydraulic fracturing, provided a critical areal extent $A_c$ is exceeded,

(2.1)\begin{equation} A_c=\frac{(1-\nu)}{2 \mu}\frac{K_c^2 }{ \Delta\gamma}. \end{equation}

Here $\nu$ and $\mu$ are Poisson's ratio and the shear modulus of the host rock, respectively. The stress intensity factor $K_I$ appears as $K_c$ to represent the fracture toughness of the host rock for mode-I fracture (Tada, Paris & Irwin 2000). This expression for $A_c$ is derived by requiring that the fracture walls close shut at the lower tip ($K_I(z=-c)=0$) and that the upper tip is critically stressed ($K_I(z=+c)=K_c$), where $c$ is the fracture's half-length and $z$ the vertical distance from its centre (Pollard & Townsend 2018). Working independently, Dahm (2000), Salimzadeh et al. (2020), Davis, Rivalta & Dahm (2020), Smittarello et al. (2021) and Möri & Lecampion (2022) each extended the 2-D analytical model of Secor & Pollard (1975) to quantify the critical volume of fluid for ascent of a 3-D fracture. The solutions derived in these studies have the same scaling with parameters but differ by around 10 % because of different numerical constants. The solution for critical volume by Davis et al. (2020) is

(2.2)\begin{equation} V_{c}= \frac{(1-\nu)}{16\mu}\left(\frac{9{\rm \pi}^4 K_{c}^8}{\Delta\gamma^5} \right)^{1/3}. \end{equation}

When the volume of fluid in a fracture exceeds $V_c$, the ascent of the fracture and contained fluid is self-sustaining. By this we mean that the ascent is due to buoyancy alone, and requires no additional forces such as a driving pressure (e.g. from a pressurised magma chamber or well bore). As the fracture ascends, the tail of the crack lengthens and thins, but the fracture head retains a volume sufficient to critically stress the medium ahead of it. Figure 1 shows numerical solutions, discussed below, which illustrate this.

Figure 1. Numerical simulation of an ascending, finite batch of fluid. In this simulation, $1.95\ {\rm m}^3$ of fluid is injected over 130 s. Two views of the fracture are shown at each time. On the left-hand side is a view looking onto the fracture's face, shaded by aperture ($w$) with the tip-line in black. On the right-hand side is a grey cross-section showing the profile of the aperture along the centre of the crack. The latter has a horizontal exaggeration of $2\times 10^4$. (a) Analytical solution that defines the radius $a$. The plotted solution contains fluid volume $V$ and has a radius such that $K_I(z=-a)=0$. (b) Numerical simulation of the fracture ascent using PyFrac (table 1). Simulation times are shown below the respective fractures. The time-dependant analytical approximation of the fracture's cross-section is shown as dashed lines (Roper & Lister 2007). The tail height $h$ from (3.10) is marked with a dot.

2.1.2. Ascent rate

Results cited above quantify the critical conditions for buoyant ascent, but not the speed of that ascent. The latter requires a consideration of the fluid flow within the fracture. In the tail of the fracture this can be usefully approximated as Poiseuille flow between parallel plates. At small Reynolds number, lubrication theory gives the mean flow speed $v$ and areal flux $Q_A$ as

(2.3a,b)\begin{equation} v=\frac{w^2}{12 \eta}\Delta\gamma,\quad Q_A=\frac{w^3}{12 \eta}\Delta\gamma, \end{equation}

where $w$ is the separation of the plates, $\eta$ is the dynamic viscosity and we have assumed that the gravity vector is parallel to the walls.

The ascent rate of a 2-D fracture with constant areal flux $Q_A$ is given by eliminating $w$ from (2.3a,b) to give (Lister & Kerr 1991)

(2.4)\begin{equation} v_Q=\left(\frac{Q_A^2\Delta\gamma}{12 \eta}\right)^{1/3}. \end{equation}

This result means that for $Q_A$ constant, the speed of the upper tip is controlled by viscous flow through a tail of constant width, which transports buoyant fluid upwards to the propagating head.

The ascent rate of a 2-D fracture with constant area $A$ is given by combining (2.3a,b) with a conservation of mass equation that relates the rate of change of aperture to the divergence of the vertical flux (Spence & Turcotte 1990). On this basis, Spence & Turcotte (1990) obtained a solution in terms of the velocity of flow through a half-ellipse with a fixed area. As this lengthens vertically its aperture decreases, which hinders fluid flow and slows the ascent of the fracture tip. The rate at which the half-ellipse lengthens therefore progressively decreases with time. Differentiating with respect to time the equation (22) of Spence & Turcotte (1990), which defines the $z$ distance between the upper tip and the initial centre point of the fracture (the injection point), the fracture's tip velocity is

(2.5)\begin{equation} v_{A}(t)=\left(\frac{A^2 \Delta\gamma}{48 \eta}\right)^{1/3} t^{{-}2/3}. \end{equation}

The ascent rate of a 3-D fracture with a constant rate of volume injection was considered by Lister (1990b). His equation (2.14) provides an approximate 3-D scaling for fracture height, breadth and aperture at a given injection rate (see also Germanovich et al. 2014). These scalings have been confirmed by Möri & Lecampion (2022) using 3-D numerical simulations computed with open-source codes (Salimzadeh et al. 2020; Zia & Lecampion 2020). Möri & Lecampion (2022) show that for a 3-D, constant injection-rate, ascending fracture, the shape of the tip-line depends on a dimensionless ratio that includes viscosity and fracture toughness. At early times, close to the injection point, viscous resistance limits tip-line propagation, resulting in a V-shaped tip line. Whereas later, sufficiently far from the injection point, fracture toughness limits tip-line propagation and the tip-line at either side of the fracture is vertical. In this study we are interested in how the ascent speed decays with time. In the case of buoyant ascent fed by a constant flux, the ascent speed in the viscosity-dominated regime is $v\propto t^{-1/6}$, whereas later, after transitioning to the toughness regime, the velocity is constant (Lister 1990b; Germanovich et al. 2014; Möri & Lecampion 2022). These limiting regimes correspond with the asymptotic solutions for an aperture near the tip-line (Detournay 2004).

The ascent rate of a 3-D fracture containing a constant volume of fluid is not known.

2.2. Numerical methods

We run 3-D simulations of an ascending fluid batch with a constant volume using a hydro-fracture simulator. This 3-D simulator generates numerical solutions to the nonlinear governing equations. These solutions guide the development of our approximate analysis for constant-volume fractures. Comparison of our analytical and numerical solutions then enables us to determine the validity of simplifying assumptions and assess the accuracy of our analytical predictions.

We simulate the propagation of buoyancy-driven fractures using the open source code PyFrac (Version 1.0, https://pyfrac.epfl.ch). The methods and implementation of PyFrac are extensively documented in Zia & Lecampion (2020). The calculations we present use the code as described in that original work, and hence we only summarise the algorithm here to give a sense of its aptitude for solving hydro-fracture problems.

The hydro-fracture problem comprises three coupled physical aspects: fluid flow and pressure inside the fracture; deformation of the surrounding medium due to the fracture; and propagation of the fracture's tip-line into the unfractured medium. These aspects are independently modelled as low-Reynolds-number fluid flow, quasistatic linear elasticity and linear-elastic fracture mechanics. Although each is linear independently, their coupling results in a nonlinear equation system. The time-dependent solution to this system is the location of the tip-line, the fracture aperture, and the fluid pressure and velocity. The PyFrac solver limits the complexity of this problem by restricting solutions to be planar fractures. This simplification enables precomputation of the Green's functions required for the displacement discontinuity method using a fixed-grid mesh (Peirce 2016). The displacement discontinuity method links fluid pressure inside the fracture to the opening of its walls and to deformation of the surrounding medium. On the same mesh, the lubrication equation is discretised using the finite volume method at element centres. The gradients in fluid pressure determine the rate of flow between cells. Crucially, to avoid the requirement of ultrahigh mesh resolution near the fracture tip, PyFrac uses near-tip asymptotic solutions for propagating, fluid-filled fractures. These asymptotic solutions resolve the subgrid-scale dynamics, and hence are used to assess the propagation criterion and rate of tip-line motion. By comparison of numerical solutions to existing analytical solutions (Peirce 2015, 2016; Zia, Lecampion & Zhang 2018; Moukhtari, Lecampion & Zia 2020; Zia & Lecampion 2020), it has been shown that PyFrac can accurately solve the coupled hydro-fracture growth problem.

For simulations in the present manuscript, the extent of the meshed domain in the horizontal direction is $2a$, where $a$ is the characteristic crack radius we define below. The vertical extent of the domain is sufficient to capture propagation of the fracture in the postinjection, buoyancy-driven regime without remeshing (a distance of $14a$). Along the horizontal direction we use a minimum grid size of 50 elements; the vertical grid is created such that the elements are square. We note that the meshed domain is only used to discretise the fracture itself. Displacements and stresses due to the fracture are described by Green's functions for an infinite space such that these hold anywhere in the body. The simulation presented here is initialised with a circular fracture of radius of $a/6$. The fracture aperture is a result of a uniform internal pressure such the critical stress intensity is attained around the tip-line ($K=K_c$). For this condition, an elapsed time since the start of injection is determined using an analytical solution, assuming the injection rate was constant (Detournay 2004). The remaining fluid volume is injected as part of the simulation, at a constant rate from a point source.

At each time step in the simulation, PyFrac returns the computed front velocity. This velocity fluctuates with a period of the order of the time step. The fluctuations, however, are small in amplitude and do not influence the general trend of the velocity. To remove them, we apply a moving average with a sampling width of 40.

The fluctuations arise from the numerical time-stepping scheme used to solve for the position of the moving fracture front, a scheme that is first-order accurate in time (Peirce & Detournay 2008). To compute the advance of the front, the user of PyFrac can choose between three time-stepping schemes: implicit; explicit; and predictor corrector (Zia & Lecampion 2019). We show a plot comparing filtered results for each of these schemes in Appendix B. Neither the scheme chosen nor the mesh size change the resulting trend of the front velocity. Following the example of Zia & Lecampion (2019), we have chosen the predictor–corrector front-advancing algorithm in this study.

2.3. Qualitative insights from numerical simulations of fracture ascent

Using the numerical method described above, we simulate the injection of a finite batch of buoyant fluid, closely following a case that was analysed by Salimzadeh et al. (2020). The parameter values for this simulation are provided in table 1.

Table 1. Parameter values used for the simulation shown in figures 1 and 2. These correspond to simulation-case II from the main text of Salimzadeh et al. (2020).

In figure 1 we plot snapshots of the fracture's shape as it evolves through time. Time increases from left to right in the plot. Looking onto the fracture's face, approximately eight minutes after the injection has ceased (10 minutes), the fracture's tip-line is circular. Note the vertically asymmetric aperture of the cross-section at this time, which indicates that the fluid is starting to flow towards the upper tip of the fracture. After six hours the top of the fracture has grown upwards whilst the tip-line in the lower portion has not changed shape; when looking at the face at this stage, the tip-line on either side of the fracture tapers progressively towards the top. The form of the upper part of the tip-line is a semicircle and remains so throughout the rest of the simulation. As time increases the aperture in the lower parts of the fracture progressively thins. By six hours, the fracture has reached a characteristic shape; in cross-section, the lower half of the fracture (tail) is thin and V-shaped while the upper half (head) has formed the teardrop profile of a Weertman crack (Roper & Lister 2007). Looking onto the face at 15 hours, the fracture's upper tip has continued to rise, resulting in a tip-line on either side of the fracture that is approximately vertical, where the tapering upwards is less evident. The cross-sectional aperture retains the typical structure predicted by 2-D analytical solutions, with a clearly defined teardrop-shaped head and tapered tail (Roper & Lister 2007). This structure dominates at late times ($t>15$ h), noting the thinning of the tail as the length increases, whilst the head's shape remains roughly constant.

In figure 2(a) we plot the speed of the upper-most point on the fracture's tip-line as a function of time (log scale). We refer to this speed as the ascent rate because it is identical to the speed at which the head of the fracture rises through the medium. A dashed line shows $v \propto t^{-2/3}$; this scaling appears in the 2-D solution (2.5). At late times greater than six hours, the speed from the 3-D numerical solution approaches this $t^{-2/3}$ asymptote.

Figure 2. Upper-tip ascent speed from the PyFrac simulation (solid line) and the predicted asymptotes and analytical approximations (broken lines). Times corresponding to steps in figure 1 are shown as asterisks. Parameter values are given in table 1. (a) Ascent rate versus time on a log–log plot. (b) Ascent rate versus height. Our prediction for the maximum ascent speed from (3.7) is shown as a horizontal dashed line. Also shown is the analytical approximation of the front speed produced using (3.10) and (3.11) (dash–dot line).

In figure 2(b) we plot the ascent rate as a function of the height of the fracture. Note the $y$-axis has a logarithmic scale and the times shown in figure 1 are plotted as asterisks. The plot shows two clear changes in front deceleration. The first occurs once the injection stops and the second occurs when the upper tip has propagated further than $2a$, a characteristic radius we examine later (Möri & Lecampion 2021). We attribute the second transition in the front speed to a shift in the process driving tip-line growth: from injection-driven to buoyancy-driven fracturing. This transition occurs once the upper tip of the fracture has propagated a given distance. If this interpretation is correct then at later times the decay in ascent speed should be related to the fluid flowing upwards with an increasingly long and thin tail (2.5) (Spence & Turcotte 1990). Such a thinning of the tail aperture can be seen in the snapshots of the ascending fracture in figure 1. In the following section we provide analytical expressions approximating: (i) the fracture's maximum lateral extent $a$ (figure 1); (ii) the maximal upwards propagation speed once injection has terminated (dashed line in figure 2b); and (iii) the deceleration of the ascending fracture with time (figure 2a).

3. Analysis of 3-D fracture ascent rates

In this section we first show that injection and buoyant ascent can occur sequentially, with minimal overlap in time. We then derive an approximation of the tip speed at the onset of buoyant rise. Next we derive a theory for deceleration of the tip over finite time. Finally, we compare our approximate theory with the simulation presented above.

3.1. Injection-driven front speed for constant injection rate

Previous studies of 2-D ascending finite batches of fluid placed the batch in the crust instantaneously and investigated how it subsequently evolves (Spence & Turcotte 1990; Roper & Lister 2007). Our 3-D approximate solution is informed by these 2-D solutions and also assumes that the finite fluid batch is fully emplaced before buoyant ascent occurs. In our 3-D numerical experiments, however, we have chosen to simulate injection of the fluid into the medium. We must therefore carefully choose a rate and duration of injection such that the fluid batch is emplaced before significant buoyant ascent has occurred, such that the approximation still holds. At early times, injection dominates the speed at which the fracture front advances; during this phase the fracture grows radially with a circular tip-line. In this section we consider how long injection of a given volume of fluid can last without overlapping the buoyant ascent phase. We therefore compare the front speed due to injection-driven growth and the front speed due to buoyant ascent.

The circular advance of the fracture tip driven by a constant injection rate in the absence of buoyancy is well understood (Detournay 2016). The growth of an injection-driven fracture can be separated into two regimes with distinct tip speeds: early times, when viscosity dominates the tip asymptote; and later times, where toughness dominates. Therefore, to compute the injection-driven front speed at a given time, we must determine which regime is dominant at that time. The transition time between viscous and toughness regimes occurs at

(3.1)\begin{equation} t_{mk}=\left(\frac{E'^{13} \eta'^5 Q_V^3}{K'^{18}}\right)^{1/2}, \end{equation}

where $K'=4(2/{\rm \pi} )^{1/2}K_c$, $E'=E/(1-\nu ^2)$ and $\eta '=12\eta$ and $Q_V$ is a constant volumetric fluid flux into the fracture (Detournay 2004).

We determine the injection-driven growth regime at the end of injection $t_I$ by comparing $t_I$ with $t_{mk}$. In the case that viscous forces dominate ($t< t_{mk}$), the speed of radial growth (${\rm d} r/{\rm d} t$) is

(3.2)\begin{equation} v_{rm}=b\frac{4}{9} \left(\frac{E' Q_V^{3}}{\eta' t^{5}}\right)^{1/9}, \end{equation}

where $b=0.6955$ (Detournay 2016). In the case that toughness controls growth ($t>t_{mk}$), Detournay (2016) describe its radial speed as

(3.3)\begin{equation} v_{rk}= \frac{2}{5}\left(\frac{9}{2{\rm \pi}^2 } \frac{E'^2 Q_V^2 }{K'^2 t^3}\right)^{1/5}. \end{equation}

Note, the volume injected at time $t$ is $Q_V t$. We use the relevant estimate of front speed at $t=t_I$ to make a comparison with the buoyant ascent speed. To this end, we next derive an approximate upper bound on the rate of buoyant ascent if the fracture is circular and filled with volume $V$.

3.2. Early-time ascent rate

Here we aim to approximate the ascent rate of a buoyancy-driven fracture after injection ceases, when the front speed becomes dominantly driven by buoyancy. This moment does not occur precisely at the end of injection; some growth beyond that time is driven by the release of elastic energy stored during the injection phase, which drives a radial motion of the tip line (Möri & Lecampion 2021). This energy is expended rapidly, after which buoyancy dominates the front speed. Our first calculation is an estimate of the speed of ascent at this point. We assume that a volume $V$ of fluid has been injected and resulted in a penny-shaped crack. The crack has extended radially to a size such that the deepest segment of the tip line, the bottom tip, has ceased its radial advance. The walls of the crack are subject to a downwards stress gradient of magnitude $\Delta \gamma$ that drives fluid upward. As this gradient becomes dominant post injection, the fracture above the bottom tip begins to drain fluid and pinch closed. Injection of a volume $V$ leads to a penny-shaped crack with mean internal pressure given by

(3.4)\begin{equation} p_0=\frac{3\mu}{8(1-\nu)}\frac{V}{a^3}, \end{equation}

where $a$ is the radius of the tip-line around the injection point (Tada et al. 2000). The mode-I stress intensity at the circular tip line arises from a combination of $p_0$ and the linear stress gradient (Tada et al. 2000),

(3.5)\begin{equation} K_I(\theta) = \frac{2}{\rm \pi}\sqrt{{\rm \pi} a}\left(p_0 + \frac{2}{3}\Delta\gamma a \cos(\theta)\right), \end{equation}

where $\theta$ is the angle in the $y$–$z$ plane away from vertical. If $z$ is the upwards distance from the injection point, eliminating the mean pressure from (3.4) and (3.5) and requiring that at the basal tip $K_I(z=-a)=0$ gives the radius of the crack as

(3.6)\begin{equation} a=\left(\frac{9 \mu }{16 (1-\nu) }\frac{V}{\Delta\gamma}\right)^{1/4}. \end{equation}

The mean aperture of this crack is simply given by $w = V/{\rm \pi} a^2$. Using this in (2.3a,b) to compute the characteristic fluid ascent speed $v$, then combining the result with (3.6) to eliminate $a$ we obtain

(3.7)\begin{equation} \tilde{v}_{V}=\frac{4(1-\nu)}{27{\rm \pi}^2\mu} \frac{\Delta\gamma^2 V}{\eta}\quad \text{for }V>V_c. \end{equation}

Here we have assumed that the fluid ascent speed $v$ is equal to the speed of the upper tip line $\tilde {v}_V$. The tilde indicates that this is an early-time solution and the subscript $V$ indicates that this is a fixed-volume, 3-D fracture. This speed is plotted as a horizontal dashed line in figure 2(b).

3.3. Radial growth versus buoyant ascent

Equation (3.7) is the speed of the upper tip while the crack is still penny-shaped, just after the phase of radial growth driven by the injection. As the upper tip propagates, driven by buoyancy of the fluid, the crack elongates vertically. The stress gradient drives fluid within the crack upwards such that part of the total fluid volume $V$ resides in the head of the crack, while the remainder forms a slowly thinning layer in the long tail (figure 1). Thus, (3.7), which describes the instant when the entire fluid volume is in the head and the lower tip reaches its maximum depth, gives an approximate upper bound on the buoyancy-driven ascent speed. Before proceeding to derive a late-time solution for the ascent speed, we use the result above to compute a maximum duration of injection $t_V$. If injection continues far beyond this, the assumption of rapid fluid emplacement in the 3-D analytical solution we introduce below will no longer hold.

The time at which buoyant ascent becomes dominant in determining the front speed can be approximated as the moment when the injection-driven front speed is equal to the buoyant ascent speed. We have approximated the latter with (3.7). We denote this crossover time $t_V$ and obtain values for each of the two injection regimes. For a viscosity-dominated circular crack, the speeds predicted by (3.7) and (3.2) are equal at

(3.8)\begin{equation} t_{Vm}=\left[ \left(\frac{b {\rm \pi}^2}{8 \Delta\gamma^2}\right)^9 \frac{E'^{10} \eta'^8}{ Q_V^6 } \right]^{1/14}. \end{equation}

For a toughness-dominated circular crack the front speeds described by (3.7) and (3.3) are equivalent at

(3.9)\begin{equation} t_{Vk}=\frac{9}{4}\left(\frac{3 {\rm \pi}^9}{5^5}\frac{ \eta^5 E'^7 }{ \Delta\gamma^{10} K_c^2 Q_V^3 }\right)^{1/8}. \end{equation}

We now show how to use these equations on an example case, the simulation in figure 1 (table 1). We first determine which regime the front was in at the end of injection (3.1). Because the duration of injection $t_I$ was 130 s and $t_{mk}$ is 400 s, we conclude that this was in a viscosity-dominated regime when injection ceased.

We then determine the time at which the injection- and buoyancy-driven front speeds would be equal. Equation (3.8) provides the time when (3.7) and (3.2) are equal, at $t_V = 720$ s. The injection time $t_I$ is much shorter than this, and hence the fracture is approximately circular when the injection stops. Figure 1 confirms this; the fracture is still circular around eight minutes after the injection stopped at $t_I$. If $t_I$ is below or of the same order of magnitude as $t_{V}$, the fracture will remain close to circular at the end of injection.

Having clarified the requirements to emplace the batch of fluid such that our assumptions hold, we next focus on the late time solution, long after injection has ended.

3.4. Finite propagation with deceleration

We now consider the vertical advance of the fracture over finite time. We define $h$ as the vertical distance between the fracture's upper tip and the initial centre of the crack (the injection point). Recalling that figure 2(a) shows the ascent rate of the 3-D simulation asymptotically scales with time as predicted by (2.5), this motivates us to approximate the 3-D ascent speed using this 2-D solution. To retrieve $h(t)$ and $v(t)$ for our 3-D fracture, we modify the similarity solution obtained by Spence & Turcotte (1990, (22)) for a 2-D fracture with constant area. Replacing that area $A$ with the initial cross-sectional area of the 3-D crack $2aw$ gives

(3.10)\begin{equation} h=\left(\frac{9 \left(2 a w\right)^2 \Delta\gamma t}{16 \eta}\right)^{1/3}. \end{equation}

Substituting (3.6) and $w = V/{\rm \pi} a^2$ into (3.10) and defining the speed $v_V(t) \equiv \textrm {d} h/\textrm {d} t$, we obtain

(3.11)\begin{equation} v_V(t)={\left(\frac{\left(1-\nu\right)}{9^2 {\rm \pi}^4 \mu } \frac{\Delta\gamma^3 V^3}{\eta^2 }\right)}^{1/6} t^{{-}2/3}. \end{equation}

This approximate solution might be expected to hold at later times, when buoyancy has driven the fracture to propagate away from the injection point, establishing distinct tail and head regions of the fracture.

Figure 2(b) shows the predicted initial and late-time speeds from (3.7) and (3.11), respectively. The figure also shows a numerical solution of a 3-D crack for the same parameter values. The ascent rate predicted by (3.11) at $h=2a$ is greater than the numerical result at the same position; at this point, the early-time approximation (3.7) fits better. At $h/2a = 3/2$, the two approximations are equal $\tilde {v}_V=v_V(t)$. As the fracture ascends beyond this height, the full numerical solution and the late-time solution converge. This convergence was also obtained in 2-D numerical solutions by Roper & Lister (2007, figure 9).

3.5. Comparison between the numerical results and analytical predictions

Using the initial cross-sectional area $2aw$ noted above, we plot as dashed lines in figure 1 the cross-section predicted by the 2-D similarity solutions of Roper & Lister (2007, § 6). See our Appendix A for the equations used in this calculation. The height of the tail in these cross-sections is given by (3.10). Focusing attention on the solutions at an elapsed time of 30 hours, the analytical solution provides a good fit to the width of the tail from the numerical solution. Furthermore, the simulated fracture's head length and shape are approximately that of the Weertman solution. The notable discrepancy is that the analytical solution predicts a greater propagation distance than produced by the simulation. This discrepancy is also evident in the comparison of fracture ascent speeds of figure 2(b) where, initially, the analytical prediction is for faster propagation, but as the simulation progresses and the fracture lengthens, the numerical and analytical speeds converge. Next we discuss the assumptions made in the derivation of (3.11) to explain why the numerical and analytical ascent rates and heights differ.

We begin with a reminder of three insights from 2-D results by Roper & Lister (2007) that also apply here. Firstly, as long as the cross-sectional area of the initial 2-D crack is above the critical head area given by (2.1), the fracture is predicted to ascend indefinitely, with a monotonically decreasing propagation rate. Secondly, in the 2-D, constant-area theory, elastic forces dominate in the head region, resulting in a head shape that remains constant over time and is described by the static Weertman solution (Weertman 1971; Rubin 1995). Thirdly, Roper & Lister (2007) show that the area of the head should be removed from the total area when computing the dynamics of the tail (e.g. (3.10)) and, in particular, in estimating the velocity and height of the fracture.

Keeping these insights in mind, we next evaluate assumptions of the 2-D, constant-area approximate solution in comparison with observations from our 3-D simulation from figure 1. In the 2-D solution, the areal extent of the head is constant during ascent, whereas in the 3-D simulation, the effective head volume decreases due to a tapering head breadth. Defining the head volume as open areas of the fracture from the upper tip down to a distance of $2c_h$ below this, we find that in simulations, the 3-D head volume decreases until it stabilises at a value around that of $V_c$. This reduction of the head volume is qualitatively shown by the lateral tip-lines in figure 1 that slightly converge upwards over time. Because of this feature, we have neglected to approximate the head volume, placing the entire cross-sectional area of the initial penny-shaped crack into the tail solution. Our analytical result is therefore based on an overestimate of the tail volume, and hence it should predict a faster ascent rate than the numerical simulation. This explains why initially, the predicted velocities are faster. Moreover, in the 2-D solution (3.10), the cross-sectional area of the fluid (tail plus head) is constant. This contrasts with observations in the simulation shown in figure 1 where along the centreline, the fracture's cross-sectional area increases by 36 % from the earliest (10 minutes) to the latest (60 hours) snapshots. This out-of-plane fluid flow is, by definition, not captured by the 2-D solution. Thickening of the cross-section in the 3-D simulation suggests that the ascent rate for a 3-D fracture will be faster than its 2-D analogue at later times, such that our solution will underestimate the ascent rate.

4. Validation of analytical approximations

4.1. Testing across length scales by comparison to numerical solutions

We now test our analytical approximations against numerical solutions across a wider range of physical parameters. Here, each parameter set represents a different physical context: analogue experiments with oil-filled cracks in gelatine; an industrial setting where water is injected into a shale sequence; and the ascent of a basaltic dyke through the crust (table 2). We compare our analytical prediction from (3.11) against rates obtained from PyFrac calculations using the three parametric combinations listed in table 2, each of which is simulated for three different injection volumes to give nine total models. Results are shown in figure 3, in comparison with analytically predicted rates. The comparison shows that $v_V(t)$ from (3.11) provides a reasonable match to the numerical results across all parameter sets considered. The analytical results predict the PyFrac velocities within a factor of two. The ratio of numerical to analytical ascent speed is approximately constant for a given simulation. In figure 4, we compare our analysis (3.10) with numerical results in terms of the time-dependent height of the upper tip above the injection point. As in the previous figures, we find the fracture must ascend by a length of more than $2a$ (box shown is $3a$) before buoyancy becomes the dominant driver of propagation. After that time, the results show the predicted and simulated heights are comparable. These comparisons show that our analytical approximations capture the leading-order fracture ascent speed across a broad parameter space. We further note that our analytical results for the extent and ascent speed provide a reliable means to forecast the required PyFrac domain size and simulation duration for a given set of parameters, such that the numerical model converges and completes the intended simulation.

Figure 3. Numerical versus analytical (3.11) speed estimates at times $t$ since injection. Plotted are nine simulations performed using the code of Zia & Lecampion (2020) and summarised in table 2. These comprise three parametric cases, each with three different injection volumes and rates. The plots stop at the point when the simulated fracture reached the edge of the meshed domain or Pyfrac terminated the simulation.

Figure 4. Numerical versus analytical (3.10) height estimates at times $t$ since injection. Height here is defined as the distance from the injection point to the upper tip. This plot uses the same numerical solutions as in figure 3.

Table 2. Properties used to numerically simulate fracture ascent for different physical processes. The injection rate for such cracks is set to $t_I=(V \, \tilde {v}_V)/a$, such that the injection is complete before the fracture ascends past $h=2a$. Note this results in $\tilde {v}_{V}/v_r$ to $\sim$2 for all results. Here $t_I/t_{mk}$ for $V_I/V_c=2$, 10 and 100 are $1.7\times 10^{4}$, 45 and $2\times 10^{-2}$, respectively; as such these span both viscous and toughness dominated radial growth regimes at the point at which the injection stops.

4.2. Comparison to analogue oil injection gelatine experiments

The gelatine experiments of Heimpel & Olson (1994) and Smittarello (2019) show that crack ascent rates are related to the injected volume. Our analysis can be used to understand what determines the speed of fractures in such experiments. Taisne & Tait (2009) noted that many fluids typically injected, such as air and some oils, are non-wetting when in contact with gelatine solids. When these flow through fractures in gelatine, the lubrication approximation breaks down; all fluid can be expelled from between the walls. Both effective closure of the fracture's tail and near-constant ascent rates are indications that non-wetting processes are active. As such, Taisne & Tait (2009) suggest that experiments that use non-wetting fluids may not be a good analogue for natural settings.

Silicon oils are wetting with respect to the gelatine. In figure 5 we plot the speed of silicon oils hydro-fracturing upwards through gelatine solids from Smittarello (2019). The material properties are given in table 3 and the experimental set-up is described in Smittarello et al. (2021). These experiments have well-constrained injection volumes and elastic parameters, and therefore provide a suitable basis for comparison with our equations. Note, the rate of injection was not constant in these experiments.

Figure 5. Upper tip height from analogue experiments where silicon oil is injected into gelatine solids. Parameter values are given in table 3. (a) Height versus time, showing that at late times, curves approach the $t^{1/3}$ asymptote (dash–dot). The triangles denote the end of injection ($t_I$). (b) The height from the experiments is plotted against the predicted speed using (3.11), where the time is that elapsed since the start of the injection. The predicted heights have been translated such that they are equal to the observed height at $t=2 t_I$ (squares).

Table 3. Properties of the analogue gelatine experiments of Smittarello (2019). Experiment reference numbers are shown in the first row. Note that the injection rates were not recorded. Note here the reported values of $K_c$ are estimated using an empirical formula based on the shear modulus, reported in § 3.1 of Davis et al. (2020).

In these experiments, a tank is filled with liquid gelatine that is allowed to solidify, then a buoyant fluid is injected at the base of the tank to form a fluid-filled fracture. This fracture ascends through the gelatine and its position is measured at regular time intervals. Figure 5(a) shows that after the postinjection transient ends, the length of the experimental fractures broadly follow a $t^{1/3}$ trend, as predicted by (3.10). Figure 5(b) shows the predicted late-time fracture lengths are around double that observed. Note that the predicted heights have been translated such that the fracture lengths are equal to the observed lengths following injection (at $t=2 t_I$). The results follow the $t^{1/3}$ trend at late times, but the magnitude of the prediction is too high. The critical volume $V_c$ required for ascent in table 3 is close to the injected fluid volume $V_I$, and our prediction assumes the entire fluid volume is in the fracture tail, as discussed earlier. When $V_I\approx V_c$, the majority of fluid is actually contained in the fracture head during ascent. Less fluid in the tail means slower draining out of the tail; this may explain why we predict ascent rates that are faster than observed. Moreover, model predictions neglect the presence of tank walls and free surface (Smittarello 2019). Experiment 1967 shows a deceleration as it approaches the surface where a load is placed on the gelatine to deflect the crack from vertical.

5. Discussion

5.1. Some insights into ascent speed

In our numerical simulations, we ensure that the injection rates are sufficiently high that the specified fluid volume has been injected into the fracture before significant ascent occurs. This means that at the end of the injection, the crack tip-line is approximately circular and $\tilde {v}_V$ should approximate the ascent speed at this time (§ 3.3). We expect that if the time scale of injection is much lower than the time scale of propagation, then once the crack has exceeded the critical radius of Davis et al. (2020), it will begin to elongate in the direction of the stress gradient. In these cases, the upwards ascent speed should lie somewhere between the two limiting regimes: that predicted by (2.5) for a batch of finite volume and the ascent rate determined for a constant rate of volume addition, (2.4) (Möri & Lecampion 2022). Even when the injection rate is low, once the entire fluid batch has been injected into the fracture and as $t$ increases, the general trend of the deceleration rate should follow a $t^{-2/3}$ asymptote.

Mutch et al. (2019) use geochemical techniques to retrieve magma ascent rates of the Borgarhraun eruption, northern Iceland. These results suggest the magma's ascent through the crust was rapid, in the range of 0.02 to $0.1\ {\rm m}\ {\rm s}^{-1}$. We now test whether our equations correctly predict such ascent speeds. The erupted lava volume of Borgarhraun was reported between $0.014\unicode{x2013}0.14\ {\rm km}^{3}$ and the magma density was around $2700\ {\rm kg}\ {\rm m}^{-3}$ (Maclennan et al. 2003; Hartley & Maclennan 2018). Assuming parameters value $\rho _r=2750\ {\rm kg}\ {\rm m}^{-3}$, $E=10\unicode{x2013}40\ {\rm GPa}$, $\nu =0.25$ and $\eta = 10\unicode{x2013}30\ {\rm Pa}\ {\rm s}$, we find the maximum ascent rate from (3.7) between 0.08 and $9.4 {\rm m}\ {\rm s}^{-1}$. Calculating the average speed from (3.11) between $2a$ and the reported distance traversed by the batch of 24 km, we find that it is between 0.06 and $9.6\ {\rm m}\ {\rm s}^{-1}$. Hence, we observe that by using typical crustal parameters, the model returns propagation velocities of the order of those observed in nature. In particular, the model provides a physical explanation for how a relatively small batch of magma can traverse the crust within a week. We note the large uncertainties in the predicted rates of ascent that are a consequence of uncertainties in input parameters.

Our results are also important in the context of hydro-fracturing operations. They can be used to quantify the time it would take for a fluid to pass into overlying formations. Furthermore, they provide an estimate of the area of rock exposed at the crack surfaces. We aim here to give an indication of how this formula can be applied to industrial operations such as hydro-fracturing. We envision the case of injecting a fluid volume of $25\ {\rm m}^3$, where the fluid viscosity ranges between $10^{-3}\unicode{x2013}10^{-2}\ {\rm Pa}\ {\rm s}$, the rock stiffness is 10–40 GPa and assuming the rock and fluid density are $2700\ {\rm kg} {\rm m}^{-3}$ and $1000\ {\rm kg}\ {\rm m}^{-3}$, respectively. For this range of properties, using (3.10), it should take between 15 minutes and five hours for the fracture to propagate 1000 m vertically and 2–40 hours to propagate 2000 m. These ascent rates and distances suggest this is an efficient way to transport fluid through the crust, and they raise the question: Which processes act to slow or stop this ascent?

We recall that in layered formations, both stress and stiffness interfaces can trap fractures (Bunger & Lecampion 2017). Additionally, low-viscosity fluids in porous formations can leak into the surrounding rock. This process is known as fluid leak-off; it would reduce the effective volume driving the fracture upwards and could, in some cases, change the dynamics of flow-driven fracture growth within the tip (Dontsov & Peirce 2015). Note that in the context of a dyke, solidification of the magma along the dyke walls is mathematically equivalent to leak-off. This process is well understood for constant-inflow, 2-D buoyant cracks (Lister 1994b). It is of interest to quantify how the upwards propagation speed and trapping height for a finite volume of fluid is changed by fluid leak-off/solidification during ascent (Detournay 2004). We leave this extended analysis to future studies where more experimental data is available, but note that Lister (1994a, (5.1)) provides a 2-D approximation of the maximal distance a fixed, finite-area dyke can propagate when solidification is considered. Recall that in § 3.4 we have shown how to convert similar equations to 3-D fractures. Note that 2-D solutions including solidification for constant flux boundary conditions also exist (Lister 1990a; Dontsov 2016). The effect of layering and rock porosity are specific to a given site. To understand how far a given fracture can ascend before termination in a given setting, it is important to quantify these properties and model them accordingly.

6. Conclusions

We have provided an analytical approximation of the maximum ascent speed of a 3-D, buoyant fluid-driven fracture containing a finite fluid volume. We verified this by comparison with outputs from a hydro-fracture simulator. We showed that the ascent speed decays away from this maximum due to viscous drag in the growing tail region, at a rate that is asymptotic to $t^{-2/3}$. Our quantitative approximations help to explain why a dyke can traverse the crust in a time that is of order of days. Although the same is true for typical injection volumes in hydro-fracturing operations, other processes will determine when the fracture is trapped and its final extent.