Critical Fluid Injection Volumes for Uncontrolled Fracture Ascent

Hydrofracturing is a routine industrial technique whose safety depends on fractures remaining confined within the target rock volume. Both observations and theoretical models show that, if the fluid volume is larger than a critical value, pockets of fluid can propagate large distances in the Earth's crust in a self‐sustained, uncontrolled manner. Existing models for such critical volumes are unsatisfactory; most are two‐dimensional and depend on poorly constrained parameters (typically the fracture length). Here we derive both analytically and numerically in three‐dimensional scale‐independent critical volumes as a function of only rock and fluid properties. We apply our model to gas, water, and magma injections in laboratory, industrial, and natural settings, showing that our critical volumes are consistent with observations and can be used as conservative estimates. We discuss competing mechanisms promoting fracture arrest, whose quantitative study could help to assess more comprehensively the safety of hydrofracturing operations.


Introduction
Official guidelines for hydraulic fracturing (e.g., EPA, 2016;Mair et al., 2012) outline safe operational practices for regulators. Such reports often state that during routine operations, fractures are unlikely to grow out of the target rock formation, as typical injection pressures are too low for this to occur. These claims are substantiated with empirical observations from closed access microseismic data of scarce vertical fracture growth following injection (Fisher & Warpinski, 2012). Evidence for unsafe vertical migration of such fluids remains ambiguous (Vidic et al., 2013).
Natural analogs of fluid migration by hydrofracturing include drainage crevasses in melting glaciers and magma transport by diking. Field and experimental observations provide some indication of typical rates of fracture ascent, in the order of mm/s to around half a m/s (Das et al., 2008;Tolstoy et al., 2006). For water-filled fractures in rock, this has not been observed; estimates from geochemical analysis supply similar rates of ∼0.01-0.1 m/s (1 km/day) (Okamoto & Tsuchiya, 2009). Theoretical arguments suggest that the migration velocity should have a dependency on volume (Dahm, 2000;Heimpel & Olson, 1994).
According to theory, tip-propagation occurs when a critical amount of fluid has accumulated, inducing enough stress to overcome the medium's fracture toughness, K c (Secor & Pollard, 1975). So far, critical "volumes" are given in terms of the fracture length, which is not directly observable and difficult to estimate from observations (Dahm, 2000;Secor & Pollard, 1975;Taisne et al., 2011); moreover, such analyses have been carried out in 2-D only, not capturing the fracture's 3-D shape and scaling of volume versus length.  Bell et al. (1990). Crack shown in red with length 2c. (b) Stress boundary conditions and 3-D crack wall displacement. (c) Cross sections of crack wall displacement, itp = interpenetration. (d) Shape of an ascending air-filled crack in gelatine from Rivalta and Dahm (2006). (e) Air-filled cracks tipline versus a penny-shaped tipline.
Here, after deriving a theoretical model and validating it with numerical simulations, we apply this to cracks filled with air, water, oil, and magma in solids of varying stiffness and toughness, across a wide range of length scales.

Hydrofracturing and Stress Gradients
We consider a pressurized penny-shaped crack of radius c and volume V in an elastic medium. The crack can only grow when the stress intensity K I at its tipline exceeds K c . The elastic parameters of the medium (shear modulus, , and Poisson's ratio, ) control the fracture's aperture. The internal pressure p 0 must overcome the stress normal to the crack walls (generally the minimum compressive stress, min ) by an amount accommodating the volume V against the elastic forces (Figures 1a and 1b).
When the crack is vertical, the gradient in the normal stress acting to close the crack and the gradient in the load due to the overlying fluid acting to open the crack, that is, r g and f g in Figure 1a, where r and f are the densities of the host rock and fluid, respectively, result in a net stress gradient Δ acting to push open the crack walls in an inverse "teardrop" shape (Figures 1b and 1c). When the crack is inclined, Δ needs to be adjusted by cos( ), where is the cracks' angle away from vertical. Quantitative formulations used to assess industrial fracture heights neglect stress gradients (e.g., Xu et al., 2019;Yue et al., 2019). This contrasts with routine observations of stress gradients from industry data ( Figure 1a) and the fact that these gradients are considered in the well design of industrial operations (Lecampion et al., 2013;Mair et al., 2012). When this gradient is included in formulations, stress intensity varies around the fracture's tipline ( Figure 2). Where K c is exceeded, the upper tipline advances. The contained fluid flows into this newly created fracture surface while the bottom edge of the fracture is pinched shut as the internal pressure drops. With a great enough volume, this fluid movement maintains a critically stressed upper tipline, and the fracture reaches a state of "self-sustaining propagation." Fluid viscosity will cause some fluid to stay trapped in the tail trailing behind the fracture; if fluid viscosity is low enough, the contained fluid is virtually all transported. Provided the fracture's shape and volume are maintained, no additional forces, such as pressure from injection, are required to aid this state of propagation.  Secor and Pollard (1975) define in 2-D the size and pressure inside a vertical fracture subject to Δ = ( r − f )g such that at the upper tip, K + I = K c , and at the lower tip, K − I = 0. We derive an analytical expression for the fluid volume needed for a 3-D crack to propagate in a self-sustained manner.

Analytical Formulation
K I for a Mode I penny-shaped fracture of radius c subject to a generic linear stress gradient can be expressed as the superposition of K I for a penny-shaped fracture subject to a uniform pressure p 0 : and that for a penny-shaped fracture subject to a linear pressure gradient Δ where pressure is equal to 0 at the fracture's midpoint (Tada et al., 2000, p. 355): where + refers to the propagating tip and -to the basal tip. Requiring K − I = 0 results in p 0 = 2Δ c∕3 and thus Requiring K + I = K c and rearranging for c yields We note that the 2-D plane strain critical length is ≈ 0.9c. The volume of the crack can be calculated based on the equation for a crack pressurized by uniform pressure p 0 , as the antisymmetric pressure contribution integrates to zero. Thus using (Tada et al., 2000), This equation requires validation in order to evaluate the bias due to approximating the shape of the propagating crack as circular (Figures 1d and 1e).

Numerical Model
To simulate propagation, we use a 3-D Boundary Element program where each element is a triangular dislocation with constant displacement (Figure 2). The program computes fracture opening and stress intensities, 10.1029/2020GL087774 based on fracture shape, rock and fluid parameters, and external stresses Davis et al. (2019). Our workflow during each iteration is as follows: 1. We invert for the uniform internal fluid pressure, P 0 , necessary to open the crack to match the required volume against all external and internal tractions. Nonlinear complementarity conditions are imposed such that the crack's faces cannot interpenetrate (Davis et al., 2019). 2. We calculate the crack opening and the stress intensity at the tipline using the method of Davis et al. (2019). In order to reduce artifacts, we smooth the stress intensity along the tipline by averaging the local K I with its two neighboring edges. 3. We calculate the advance or retreat of the tipline. At elements where K I exceeds K c , the tipline will advance proportionally to K I ∕K c . This approximation is akin to the "Paris fatigue law" (Lecampion et al., 2018). The maximum crack advance will occur at the triangle where K I is maximum; this advance is set equal to the mesh's average triangle size. The triangular elements that close are removed. The simulation assumes the fluid is inviscid, and, as such, we cannot retrieve time-dependent propagation rates. 4. Once the fracture's edge has been updated, it is remeshed and cleaned such that the triangles on the fractures tipline are approximately equal size and isosceles (Da & Cohen-Steiner, 2019).
For a description of the numerical methods accuracy, see Appendix A1. We start the simulation with a vertical penny-shaped crack. We fix the number of elements, K c , Δ , , , and the volume of fluid, V. We set the initial radius to 0.4c (Equation 4). In our 350+ simulations, we use variables spanning several orders of magnitude: G = 190-5·10 10 Pa, = 0.25-0.49, Δ = 7.8 · 10 2 − 2.2 · 10 4 Pa·m −1 , and K c = 1 − 1 · 10 8 Pa·m 0.5 . We state the fracture has reached self-sustaining ascent when its upper tip has traveled 4c upwards.
For all simulations, independent of mesh sampling, we find that if V = 0.7V an c , the numerical code returns a trapped fracture, and if V = 0.8V an c , the fracture always reaches self-sustaining propagation. Therefore, scaling Equation 6 by 0.75 supplies the numerical estimate of V c , independent of the scale we use: For all cracks that reached self-sustaining propagation, the horizontal and vertical lengths were greater than ∼0.6c and ∼1.14c, respectively.

Analog Gelatine Experiments
The analog study of Heimpel and Olson (1994) inspects critical volumes of fluids ascending in gelatine blocks of different stiffness and fracture toughness (Figure 3). The graph of volume versus speed from their experimental results shows a rapid increase in speed past a certain volume. The authors interpret that at velocities past ∼0.7 cm/s (crosses in Figure 3), the ascent transitions from a subcritical propagation regime (K I < K c ), where the fracture growth speed at the tip limits the velocity (Atkinson & Meredith, 1987) to a dynamic propagation regime. We test if our equation can predict volume of fluid that causes this transition in ascent speed. As Heimpel and Olson (1994) estimate K c directly from this change in velocity, to verify that we can use this estimation, we calculate K c differently, directly from the measured value of G. Strain energy release  increases with greater stiffness in gelatin solids: [N/m] ≈ 6.66·10 −4 G [Pa] (Czerner et al., 2016). This second independent estimate of K c lies within 5 Pa·m 1/2 of the original estimate. Using r = 1,000 kg·m −3 , = 0.5 and setting , f , and K c to match the experiments of Heimpel and Olson (1994), we find that our value of V num c independently captures the point at which the velocity transitions, described above, supporting the previous interpretation that this describes the transition to dynamic fracture propagation.

Magmatic Dikes
We consider magma propagation volumes at Piton de la Fournaise, La Réunion, to see how our equation matches observed dike volumes. The volumes of the dike intrusions observed between 1998 and 2016 range from 0.05-3.2·10 6 ·m 3 (Froger et al., 2004;Fukushima et al., 2005Fukushima et al., , 2010Smittarello et al., 2019). Using r − f = 100 kg·m −3 , = 5 GPa, = 0.25, and K c ranging from 29 to 112 MPa·m 1/2 (Delaney & Pollard, 1981;Fukushima et al., 2010), we retrieve V num c = 0.05 · 10 6 and 2·10 6 ·m 3 , respectively. The critical volumes we estimate are consistent with the observed dike sizes. As such, our approximation predicts the correct scale in natural settings, provided K c values estimated from field data are used, noting such field values appear to  Heimpel and Olson (1994). Equation 7 predictions shown as black lines. The thickness of the gray filled patches represents the velocity of the crack as the volume increases, normalized by maximum observed velocity. correct for a number of additional processes that we have disregarded, instead of being representative of the rock strength at the scale of a laboratory sample.

Water Injection Into Stiff Rock
The UK government defines hydraulic fracturing as operations that use over 1,000 m 3 of fluid per frack stage. During a hydrofracturing procedure, proppant is injected in the final phase to maintain an open fracture (e.g., spherical quartz grains). After the operation, not all injected fluid is recovered when the wellhead valve is opened: Vidic et al. (2013) report an average of only 10% fluid recovery in flowback waters, noting that this recovery volume decreases when shut-in times are longer. Using r = 2,700 kg·m −3 , f = 1,000 kg·m −3 , = 8.9 GPa, = 0.25, and K c in the range 0.36-4.05 to 7-25 MPa·m 1/2 , we obtain V num c = 6 · 10 −2 and 500 m 3 , respectively. These K c values are for laboratory-sized shale samples from 100 to 1.000 m confining pressure and effective K c values estimated for veins in the field, respectively (Gehne et al., 2020;Olson, 2003). Current operations use volumes around double our highest predicted limit. Few observations attest to the fact that industrial operations can cause ascent of fluids in fractures. One such example is the spectacular surface fissures created due to steam injection documented in Schultz (2016); additional examples can be found in Schultz et al. (2016). Geochemical data from aquifers above fracking operation sites have shown some evidence of the contamination of overlying units, which is attributed to poor well casing design, rather than fracture ascent (Vidic et al., 2013). Usually, microseismic monitoring of actual fracking operations show limited vertical extents of the fractures; however, these data are proprietary, and methodological descriptions are scarce (Fisher & Warpinski, 2012). Experimental fracturing data are of little help as volumes injected are typically below or close to our volumetric limit, with injected volumes of 2 to 20 m 3 (Pandurangan et al., 2016;Warpinski et al., 1982).
Natural degassing, such as CO 2 in the Cheb Basin, Czech Republic, has chemical signatures of fluids that have ascended over 20 km through the crust (Weinlich, 2014). Fracture-driven ascent can explain this phenomena without the requirement of permanent highly conductive fluid pathways at great depths. Supercritical CO 2 at depth has a similar density to water and as such may be a good natural analog for water-filled fracture ascent. We saw that in analog and magmatic examples, Equation 7 predicts the correct order of magnitude of critical volumes; at the same time, it appears that this equation is conservative for high volume water injection as fracture ascent in these settings has rarely been observed.

Discussion and Conclusions
In summary, Equation 7 provides an estimate of the minimum fluid volume for self-sustained propagation of fluid-filled fractures, ranging from cm to tens of km. V c is dependent on K 8∕3 c ; since K c is often poorly constrained, V c suffers from large uncertainties. Values of K c obtained in laboratory experiments show a strong dependency on pressure and temperature. Field estimates of effective K c from trapped fractures can be orders of magnitude larger. An effective way to estimate K c in Equation 7 incorporating all processes affecting the energy needed to extend the fracture at different scales would clearly be beneficial for any fracture-mechanics-based analysis of rock masses and the resultant interpretations.
In our derivation, we have neglected the effects of viscosity. Whether these effects will dominate over toughness in determining fracture growth can be assessed by evaluating the time scale needed for the fluid pressure to equilibrate within the crack, as this will mean that viscous dissipation is low and crack growth will be toughness dominated (Bunger & Detournay, 2007). The model of Bunger and Detournay (2007) assumes a constant injection rate with no stress gradients; we assume this still provides a rough estimate of the time scale until this transition. Typical industrial operations use fluid viscosity of 0.001-0.01 Pa·s, injection rates between 0.5 and 10 m 3 /min, and stiffnesses of 10-40 GPa. Using low values of K c from laboratory experiments in shale, 0.36 MPa·m 1/2 , this transition time ranges between 1 min and times exceeding the end of injection. On the other hand, setting K c higher from values for shale at depth, for example, 4 MPa·m 1/2 , significantly reduces this range from milliseconds to a maximum of 5 hr. This suggests that, depending on K c , Equation 7 can be a relevant estimate of V num c , independent of viscous forces.
While theory and experiments support Equation 7, this appears to be overly conservative in practice, as injections of quantities of fluid exceeding this do not result in significant ascent in most cases. In part, this discrepancy results from our simplification of the process as mass-conserving propagation in a homogeneous linear-elastic medium. Figure 4 shows a schematic of processes not quantified in relation to critical fluid volumes, which we review in detail below.
1. A series of mechanisms can reduce V during propagation and thus promote crack arrest. These include leak-off from the fractures faces, the fracture becoming multistranded/compartmentalized, fluid recovery (extraction), or fluid remaining in the tail of the fracture due to added proppant or viscous forces (Taisne & Tait, 2009). 2. Mechanisms that can lead to an effective increase of K c , and thus also promote crack arrest, include plastic tip processes, the fracture entering in a zone of damage of the host rock (Kaya & Erdogan, 1980;Sih et al., 1965), or seismicity surrounding the fracture, causing reduction in the system's energy/blunting the fracture's tip (Rivalta et al., 2015). Figure A1. Stress intensity factor approximation using the 3-D displacement discontinuity method. (a, c) Elliptical crack meshed with 650 and 1,500 triangles, respectively, (b, d) Numerical (dots) and analytical (solid line) results.
Results are normalized relative to the maximum analytical value of K I .

Geophysical Research Letters
10.1029/2020GL087774 1) (Atroshchenko et al., 2009). We note that under a stress gradient, K I for vertically aligned elliptical cracks is not maximal at its upper tip, due to the reduction in crack surface area proximal to this edge.
For a mesh with 650 triangles (Figures A1a and A1b), the greatest vertical separation between the numerical points and the analytical line is 0.09. For this test, we required the edge points of the mesh's triangles, not the midpoints of the triangles edge where K I is calculated, to lie on the tipline defined by the analytical solution. For a mesh with 1,500 triangles ( Figures A1c and A1d), the maximum vertical distance from the analytical solution is 0.06, noting that greater sampling does not necessarily converge to an improved accuracy; see appendix of Davis et al. (2019).
As a further test of numerical accuracy, we compare how well the numerical method approximates the opening volume of a penny-shaped crack subject to tension (Tada et al., 2000). We find that a sampling of 650 triangles overestimates the volume by 5.2%; by increasing the triangle count to 1,500, this drops to 3.5%.

Data Availability Statement
The code used for the numerical analysis of this study was the open-source code https://doi.org/10.5281/ zenodo.3694163 with an interface with the Computational Geometry Algorithms Library software (C++) for meshing.