Finite element simulations of hydraulic fracturing: A comparison of algorithms for extracting the propagation velocity of the fracture

The finite element method is a powerful and general numerical method used to simulate subsurface processes. In this paper, we take recent hydraulic fracturing propagation algorithms and assess their performance when used within the finite element framework. In particular, we evaluate aperture and energy-based methodologies that are capable of extracting the propagation velocity of a hydraulic fracture propagating throughout the toughness and viscous regime. Such algorithms have the benefit of a quicker convergence on the fracture front. The aperture-based methodology consists of the multi-scale aperture asymptote that is yet to be applied with finite elements. On the other-hand, the energy-based methodology consists of a recently developed procedure for predicting the propagation velocity from the energy release rate, which is calculated using a 𝐽 -integral devised for hydraulic fracturing. A comparison of the accuracy and the number of iterations required to converge on the fracture length is undertaken, and found to produce similar results for both methods. Consequently, we conclude that the higher accuracy of energy-based methods in extracting stress intensity factors does not immediately translate to a higher accuracy in extracting propagation velocities, most notably in the toughness-dominated propagation regime. Given the similar performance of the methods, and the simplicity of the aperture-based approach, we then extend the evaluation of the multiscale aperture asymptote to the case of buoyancy-driven propagation. As a result, the aperture asymptote is shown to be a simple and efficient method for the simulation of subsurface processes using a finite element framework. for the fluid properties at 400 and 800 degrees, respectively. Hence, the ability of the HFM asymptote to simulate propagation in the mixed propagation regime is evaluated, while also incorporating the effects of the fluids buoyancy on propagation. HFM asymptote can be readily extended. Using linear shape functions, the HFM asymptote was shown to accurately reproduce the length and propagation velocity of a buoyancy-driven hydraulic fracture propagating within a mixed toughness/viscous propagation regime. The simplicity of the HFM asymptote, and its ease-of-use for buoyancy related propagation demonstrates the method to be highly suitable within finite element simulations.


Introduction
Deepening our understanding of subsurface processes is fundamental to securing a supply of sustainable geothermal energy [1][2][3], mineral resources [4,5] and safe carbon storage strategies [6,7]. In such cases, the hydraulic fracturing of rock via a pressurised fluid is a highly relevant mechanical process occurring across multiple scales and constraining the transportation of fluids, heat, and/or dissolved minerals [4,[8][9][10]. Since observations are typically scarce and costly, numerical simulations become an essential tool utilised to diminish the uncertainty around the possible fluid pathways, and to derive a fundamental understanding related to the exploitation of the resource.
The numerical simulation of hydraulic fractures often makes use of analytical solutions, termed asymptotes, to infer the extent of propagation of a fracture. Early asymptotes [11][12][13] were restricted to a particular regime of propagation; the toughness regime, Table 1 Methodological studies on hydraulic fracturing using FEM. The various propagation methodologies employed are the displacement correlation method using the toughness asymptote, domain/contour integral methods such as the J-integral or Interaction integral, the cohesive zone method, or the virtual crack closure technique and modifications thereof.
The accuracy and computational efficiency of the HFM asymptote in simulating hydraulic fracturing has not been evaluated within the finite element method (FEM). Indeed, the use of the HFM asymptote has remained confined to the boundary element method (BEM) or Displacement Discontinuity Method (DDM) [29,[31][32][33][34][35][36]. FEM has long been established as a powerful framework with which to conduct thermo-hydro-mechanical simulations of reservoir and crustal scale processes [25,[37][38][39][40][41][42][43][44][45][46]. Hence, using the HFM asymptote to predict the propagation velocity and track the fracture front would be a welcome development that augments the existing capabilities of finite element codes. An extensive, but not exhaustive, literature review of FEM studies with a focus on the verification of the propagation methodology are shown in Table 1. Fundamentally, the evaluation of the HFM asymptote within a finite element framework is warranted.
Recent advances in hydraulic fracturing have enabled energy-based methods to also extract the propagation velocity of the fracture [72]. Such an approach relies on the HFM asymptote to construct a semi-analytical evaluation of the -integral, which can be used to obtain the location of the fracture front using a numerically calculated energy release rate. The strength of such an approach has not been assessed relative to using the HFM asymptote directly to advance the fracture front. In particular, it remains unclear whether the superior accuracy of energy-based approaches in estimating the stress intensity factor [73] also confers benefits when predicting the location of a hydraulic fracture. Furthermore, since both algorithms require iteration to converge on to the location of the fracture front, an assessment of the relative speed of convergence between the two methods is warranted. Consequently, the aims of this manuscript are two-fold: (i) To evaluate a finite element implementation using the HFM asymptote to advance the fracture front. (ii) To compare the accuracy and speed of the HFM asymptote against energy-based approaches to locate the fracture front.
This manuscript is outlined as follows. In Section 2, the governing equations and finite element discretisations that describe the hydraulic fracturing problem are detailed. Then Section 3 introduces the asymptotics of a hydraulic fracture, and details how the stress intensity factor and propagation velocity can be extracted using the HFM asymptote. An overview of the energy-based approach of Pezzulli et al. [72] to estimate the propagation velocity is also given. In Section 4, the details of a propagation algorithm which avoids re-meshing by modifying the time step is discussed for the Khristianovich-Geertsma-de Klerk (KGD) problem. The results of the comparative analysis with the energy-based methodology ( −algorithm) developed by Pezzulli et al. [72] is then shown in Section 5. Finally, a further evaluation of the HFM asymptote to simulate buoyancy-driven propagation is provided in Section 6 to illustrate the capacity of the method for geomechanical finite element simulations.

Mathematical problem
In this work, the rock is assumed to undergo linear elastic deformation under plane-strain conditions at infinitesimal strains. We neglect the inertia in the rock and the fluid, and consider the rock as impermeable with isotropic and homogeneous mechanical properties. The fluid is incompressible with a constant Newtonian viscosity. We ignore the effects of leak-off, and do not simulate a fluid lag, whose presence is insignificant at sufficient depth [74,75].
The hydraulic fracturing problem amounts to finding the fluid pressure in the fracture , the vector displacements in the rock , and the geometry of the fracture . A geometrical conceptualisation of the problem is depicted in Fig. 1. A rock is bounded by a surface , subject to traction , or displacement , boundary conditions. The geometry of the fracture is characterised by the one-dimensional curve , that sits on the fractures mid-line. It has unit normal that points by convention from − to + , which are the location of the fracture surfaces and are delineated by the displacements + , − , respectively. Consequently, the discontinuity in the displacement field represents the fractures aperture = ( + − − ) ⋅ . The fluid is conceptualised to lie on , while exerting its pressure on − and + .

Governing equations
In view of the assumptions discussed in the preceding section, the displacements in the rock are constrained by the static force equilibrium and linear elasticity where = is the body force vector, which is related to the density of the rock , and the acceleration due to gravity written in vector form = ( , ) . The stress in the rock relates to displacements via the strains and the isotropic elasticity matrix where is the Young's modulus and denotes the Poisson's ratio. The Reynolds-lubrication equation describes the evolution of fluid pressure inside the fracture where the fluid flux in the fracture is given by the solution to the momentum equations for planar Poiseuille flow where is a coordinate in the direction aligned along the fracture , whose tangent intersects the horizontal at an angle . The fluid density is , while the lumped parameter ′ = 12 includes the Newtonian dynamic viscosity . The Dirac delta function is used to denote the location of the constant point volume source 0 entering the fracture. In the case of a vanishing fluid lag at the fracture tip, the fluid flux must be zero at the fracture tip to ensure finite fluid velocities [76]. In the vicinity of the fracture tip, the Reynolds lubrication can be shown to simplify to ′ = 2̂( 2) wherêis the inward distance from the moving fracture tip̂= − . Such an equation describes how the propagation velocity of the fracture is constrained in the absence of a fluid lag, and in the absence of leak-off [67,68,[76][77][78][79].

Fracture-tip asymptotics
Hydraulic fractures can be classified to propagate according to an array of propagation regimes [80]. The most fundamental propagation regimes are the toughness and viscous propagation regimes, where the dimensionless toughness  quantifies in which regime a plane-strain fracture is propagating. The lumped parameters are ′ = √ 32∕ and ′ = ∕(1 − 2 ), while is the fluid flux entering the fracture. Values of  ≲ 0.41 imply viscosity-dominated propagation, while  ≳ 4.1 implies toughness-dominated propagation. Based on the propagation regime, the aperture at the fracture tip displays a particular profile, herein termed the toughness asymptote and viscous asymptote [11,13,74,78,[80][81][82] wherêis again the inward distance from the moving fracture tip. The quantity is the propagation velocity and constant = 3.14735. Which asymptote describes the aperture profile depends on the distance from the fracture tip compared to the transition lengthscale where the toughness asymptote is valid at distanceŝ≪ , while the viscous asymptote describes the aperture for̂≫ , as shown in Fig. 2.
Which asymptote is required to advance the fracture front in numerical simulations depends on the resolution at the fracture tip ℎ in comparison to the transition lengthscale . Consequently, when the transition length-scale is very small (in the viscousdominated regime), apertures very close to the fracture tip (̂∼ 10 −3 ) would have to be numerically resolved if the toughness asymptote is to be applied to extract stress intensity factors [79]. On the other hand, a large transition length-scale (toughness-dominated regime) would render the viscous asymptote applicable at impractical distances from the fracture tip, where possibly asymptotic considerations are no longer valid. Such reasoning motivated the discovery of the multi-scale aperture asymptote (herein termed HFM asymptote) by Dontsov and Peirce [19,20] which approximates the aperture in the vicinity of the tip for all values of̂∕ , as shown by the black line in Fig. 2. Such a HFM asymptote initially follows the toughness asymptote for small values of̂∕ , before transitioning to the viscous asymptote at large values of̂∕ . As a result, the HFM asymptote provides a flexible and simple solution to extract the stress intensity factor and the location of the fracture front throughout the toughness-viscous propagation regimes. A finite element implementation using the HFM asymptote to advance the fracture front is elaborated in Section 3.  from the finite element solution to predict the distance to the fracture front * . The hydro-elastic equations are solved for a trial fracture front located at , for iteration . Using the HFM asymptote, an estimate for the location of the true fracture front , +1 can be obtained.

Finite element discretisation
Employing Galerkin finite elements for the spatial discretisation, while implicit Backward-Euler finite differences for the temporal derivatives leads to the following discretisation of the governing Eqs. (1) where hat superscripts denote nodal values which arise from the finite element approximation =̂, =̂. The shape functions and are 2-dimensional and 1-dimensional, respectively. The finite element matrices are and is the applied load vector, while is the 2-dimensional strain-displacement matrix [41] containing the derivatives of the shape functions .
The coupled and non-linear system of equations are solved for the fluid pressure and rock displacements using a Picard iterative procedure, whereby values for the aperture at previous iterations are used within the conductance matrix H and buoyancy vector G. The numerical model was developed within the Complex Systems Modelling Platform (CSMP++) [83] where the functionality of representing fractures as explicit discontinuities within the finite element mesh was developed. The linearised system is solved using the Algebraic Multigrid Methods for Systems (SAMG) Solver [84]. A minimum aperture = 1.5×10 −4 m is imposed in order to guarantee a minimum hydraulic conductivity of the fracture. Such an aperture was chosen since it delivers acceptable fluid pressures at the fracture tip (see fig 5a in Pezzulli et al. [72]), while still converging quickly with SAMG. Since the fracture is continuously open due to fluid injection, the minimum hydraulic conductivity only applies for the tip node where the aperture is zero. Smaller minimum apertures should be used if less-conductive 'closed' fractures are desired.

Algorithms for extracting the propagation velocity of a hydraulic fracture
This section is concerned with the prediction of the location of the front of a hydraulic fracture , +1 , and by extension its propagation velocity , +1 , within a finite element implementation. The iterative procedure requires the solution for the displacement , and fluid pressure , at a trial fracture front , for time step and trial iteration . An illustration of the fracture front geometry is given in Fig. 3, while a description of how the predicted fracture length , +1 is used to advance the fracture front is detailed in Section 4. Note, the existence of a fluid with non-zero viscosity theoretically constrains the velocity of propagation of a hydraulic fracture. This constraint enables the inference of a fracture front , +1 and propagation velocity , +1 which satisfies = . Such a feature can be understood by noting that the tip asymptotics depend on the propagation velocity as soon as viscosity is considered [78]. The dependence of the tip asymptotics on the propagation velocity allows one to estimate the location of the fracture front from numerical solutions for the aperture or energy release rate.

Finite element implementation of the HFM asymptote
The location of the fracture front is predicted with the HFM asymptote by finding the root of the following equation: which is derived by substituting the transition lengths-scale Eq. (5) into the HFM asymptote in Eq. (6), and rearranging to one side. The numerical aperture = (̂) is sampled at a distancêbehind the current trial fracture length , , while the unknown distance * describes the location to the predicted fracture front , +1 = , −̂+ * (see Fig. 3). The unknown propagation velocity , +1 = ( , +1 − −1 )∕ , depends on the location of the fracture front at the previous time step −1 , and can be expressed in terms of * , giving where , is the time step at the current iteration for time step . Note that ′ is known, as it depends only on the critical stress intensity factor . Hence, the distance to the fracture tip * can be found by finding the root of Eq. (8), giving also the propagation velocity , +1 via Eq. (9). Newtons method is used to converge on the root of equation Eq. (8), which typically converges in approximately 3 iterations.
The explicit representation of fractures with the finite element method enables a straightforward extraction of the numerical aperture . In this study, the aperture is sampled at the corner nodes of the first element along the fracture, regardless whether the mid-side nodes from quadratic interpolations exist (as shown in Fig. 3). Note that if the fracture is prescribed to propagate by one element's length each time step, then such a sampling procedure gives , +1 = * ∕ , . The aperture is calculated by taking the dot product of the unit normal of the fracture with the difference between the discontinuous displacement fields + , − defined on the top and bottom node, respectively. With the numerical aperture , the distance to the fracture tip * is predicted by inverting Eq. (8).

Finite element implementation of the energy-based -algorithm
The -integral is a powerful technique used in classical fracture mechanics to extract the energy-release rate of a fracture [85][86][87][88]. Recently, an adaptation of the -integral [72] has provided an energy-based methodology to calculate the location of the fracture front of a hydraulic fracture. The procedure is similar to the aperture-based approach above, but involves calculating the energy release rate for the current fracture front , instead of the aperture . Subsequently, the new location of the fracture front , +1 is predicted by finding the root of a semi-analytical -integral with the calculated energy release rate . The steps involved can be summarised as: 1 Calculate the energy release rate using the -integral (Eq. (19) in [72]).
2 Use a semi-analytical expression for the -integral (Eq. (27) in [72]) to solve for the location of the fracture front , +1 and the propagation velocity , +1 , given the calculated energy release rate in step 1.
For the details, the reader is referred to the −algorithm in Pezzulli et al. [72]. Two finite element rings surrounding the fracture tip are used herein for the integration involved in the -integral. The plateau function is also used, since it is a requirement of the method, and is constructed using bi-linear shape functions as in Pezzulli et al. [72].

Extracting the stress intensity factor
Given a velocity of propagation, the stress intensity factor can be calculated with the aperture or energy-based methodologies described above. For the HFM asymptote, rearranging Eq. (6) gives where the numerical aperture is at a known distancêfrom the tip, and , is the propagation velocity implied by the current fracture geometry , . When the fracture has not yet started propagating, can be estimated using Eq. (10) by setting , = 0, which is an identical procedure to the displacement correlation method in LEFM.
Similarly for the energy-based approach, the -integral, and -integral of Pezzulli et al. [72] can be used to find the energy release rate of the fracture given the current propagation velocity , . Setting , = 0 in their formulation returns the classical -integral for estimating the energy release rate appropriate for 'dry' fractures.

Hydraulic fracturing algorithm
The propagation algorithm for the energy and aperture based approaches follow the same structure. Both approaches require the computation of the stress intensity factor for an idle fracture, and a prediction and convergence on the propagation velocity  fig. 6 in [72]) is utilised, which details the use of the -integral, and the -algorithm to obtain the stress intensity factor and the propagation velocity, respectively. For the case of the aperture-based approach, an analogous propagation algorithm is used that makes use of the HFM asymptote instead of the various -integrals mentioned above. An overview of the propagation algorithm using the HFM asymptote is given herein, and shown in Fig. 4 for a time step .
An inner loop solves the governing hydro-elastic equations for the displacement of the solid , and the fluid pressure , , given a fracture length , and time-step , that are found by iteration over an outer loop. Note, that the modification of the global time step is not a procedure that can be generalised to handle multiple fracture tips, or 3-dimensions. In general settings, the underlying mesh would have to be modified, or the aperture at the tip would have to be prescribed according to the tip asymptotics [14,20,31].
Iteration on the correct time step , (and therefore fracture length , ) only occurs if the stress intensity factor estimated before fracturing ( = 0 and , = 0) is greater than the fracture toughness > . Such a stress intensity factor is calculated using Eq. (10), which amounts to using the displacement correlation method in LEFM. If the fracture is deemed propagating, its length is increased by one element ,1 = ,0 + ℎ of length ℎ, and iteration on the time-step , begins. For the propagating fracture, once the aperture is calculated, the location for the new position of the fracture-front , +1 is estimated by solving for * within Eq. (8). Consequently, the propagation velocity , +1 is obtained using Eq. (9). Given , +1 , the time step , +1 = ( , − −1 )∕ , +1 which would propagate the fracture to its current location , can be found. For both aperture and energy-based methods, convergence of the outer loop is achieved when whereby is a prescribed tolerance set to = 0.005. Re-meshing of the fracture tip is avoided by iterating on the time-step.

Khristianovich-Geertsma-de Klerk evaluation
In this section, we evaluate the accuracy and convergence speed of the HFM asymptote in predicting the location of the fracture front within a finite element framework. The results are compared to the energy-based -algorithm developed in by Pezzulli et al. [72]. The Khristianovich-Geertsma-de Klerk (KGD) problem [89,90] of a hydraulic fracture propagating in plane strain fed E. Pezzulli et al. The percent error in the fracture length = 100 ( − )∕ is shown in Fig. 6a and b for the toughness and viscous-dominated regimes, respectively. The analytical solution for the fracture length is taken from Garagash et al. [91] (toughness-dominated) and Garagash et al. [92] (viscous-dominated). For each time step, the fracture advances by at most one element of length ℎ, and the percent error is plotted over the relative discretisation of the fracture ℎ∕ . As the simulation progresses, some errors in the fracture length accumulate from the temporal discretisation, leading to results at small ℎ∕ (late times) having marginally higher accumulated errors from the temporal discretisation. Nonetheless, such an effect is small since the results in Fig. 6 remain unchanged for simulations with different mesh resolutions ℎ (and consequently different time steps).
In the toughness-dominated regime, the HFM asymptote is marginally less accurate than the -algorithm for most discretisations of the fracture. This can be seen in Fig. 6a for both linear and quadratic elements. More specifically, the percentage error for the -algorithm stays around ∼ −7% for linear elements, while converging to ≲ 2% when more than 20 quadratic elements (ℎ∕ < 0.05) discretise the fracture. The HFM asymptote achieves a marginally less accurate result with 20 quadratic elements, where the error ∼ −3% is one percentage point higher than the -algorithm at the same resolution. On the other hand, the HFM asymptote admits approximately an additional ∼ 2% error compared to the -algorithm when using linear elements in the toughness regime. The modest difference in accuracy between the two methods can be explained by the reliance of the HFM asymptote on values for the aperture at the fracture tip that exhibits a strong non-linear profile in the toughness-dominated regime which is difficult to resolve accurately, particularly with linear elements. On the other hand, the -algorithm relies on the displacements and fluid pressures among a ring of elements that do not encompass the fracture tip, such values are likely to be more accurate, hence the marginally higher accuracy of the -algorithm observed in Fig. 6a. In the viscous-dominated regime, both approaches achieve accurate results with less than 4% error in the fracture length reported for various fracture discretisations. In particular, both the HFM asymptote and the -algorithm have errors | | < 1% at discretisations ℎ∕ ∼ 0.03 and ℎ∕ ∼ 0.07 for linear and quadratic elements, respectively. On the other-hand, at courser resolutions, the HFM asymptote displays a superior accuracy relative to the -algorithm. Namely, the HFM asymptote has errors | | ≲ 2% for both linear and quadratic elements when less than 10 elements discretise the fracture (ℎ∕ ≳ 0.1), while the -algorithm admits errors | | ≲ 4% for both element types at the same resolution. The higher accuracy of the HFM asymptote relative to the -algorithm may be explained by the -algorithm's use of asymptotic relations outside of their domain of validity. Since the -algorithm makes use of several asymptotic relations in order to extract the propagation velocity, the method risks incurring errors earlier than the HFM asymptote since it applies the asymptotes at greater distances from the crack tip. This is the case if 2 or more element rings are used by the -algorithm, and may be the reason behind the increase in error relative to the HFM E. Pezzulli et al. asymptote at coarser discretisations. Overall, both methods achieve relatively low errors in the viscous-dominated regimes, with marginally better results for the HFM asymptote at coarser discretisations of the fracture.
Regarding the convergence speed of the methods, both the HFM asymptote and the -algorithm converge typically in around 3−5 iterations, with the HFM asymptote converging faster by one iteration on average. Such a theme can be observed for both linear and quadratic elements, and in both the toughness and viscous dominated regimes in Figs. 6c and 6d. For both methods, the number of iterations in the viscous-dominated regime show a consistently decreasing trend from 5 iterations at early times, to 2 iterations at late times. Furthermore, the difference between linear and quadratic elements is insignificant in the viscous regime, as the results for different element types essentially overlap each other. However, the HFM asymptote achieves its quickest convergence rate before the -algorithm, as it converges in 2 iterations by the time 30 time steps are surpassed, while the -algorithm begins to converge in 2 iterations at roughly 80 time steps. In the toughness-dominated regime, while the general trend is similar, the number of iterations display a higher variance. This is particularly present for the -algorithm, which continues requiring 4 iterations with quadratic elements till late times, while linear elements require between 2 − 15 iterations throughout the simulation (with also two time steps reaching ∼ 30 iterations, which have not been plotted to aid the visualisation of the remaining results). On the other hand, the HFM asymptote maintains a relatively small variance in the number of iterations, with few iterations observed above 4 after 50 time steps. Overall, while both methods converge on the fracture front with relatively few iterations, the HFM asymptote achieves more consistent results compared to the -algorithm, for both linear and quadratic elements.

Aperture and pressure snapshots
A comparison of the aperture and net pressure ( − 0 ) profiles predicted by the -algorithm and the HFM asymptote are shown in Fig. 7 for the toughness and viscosity-dominated propagation regimes. Quadratic finite elements have been used, and are compared to the semi-analytical solutions of Garagash [91] in the toughness-dominated regime, and Garagash & Detournay [92] and Perkowska et al. [93] in the viscous dominated regime. The solution for Perkowska is reproduced in Appendix. Both numerical models can be seen to give very low errors in the viscosity-dominated regime, while admitting some errors in the toughnessdominated regime. This is consistent with the errors reported in the fracture length in Fig. 6a,b. The largest errors are observed for the net pressure in the toughness-dominated regime, which can reach up to 10%. Such a feature results from the underestimation of the fracture length, which leads to elevated pressures in the fracture by conservation of mass. By the same mechanism, the use of linear elements gives larger errors in the net pressure (∼ 40%), although such results are not shown herein. One should note that the percent errors in the fluid pressure are significantly smaller.

Aperture vs energy-based methods
In the preceding section, only slight differences were observed in the accuracy of the HFM asymptote compared to the algorithm. However, energy-based methods like the -integral in LEFM are known to be significantly more accurate in predicting the stress intensity factor than aperture-based methods like the displacement correlation method [73]. Hence, the results in the preceding section beg the question as to why the HFM asymptote and the -algorithm share essentially similar levels of accuracy when estimating the propagation velocity to propagate the fracture front. The question is even more prominent within the toughness-dominated regime, where numerical estimates for the energy-release rate are generally, and significantly, more accurate than numerical solutions for the aperture. Yet, both methods display a ∼ −4% error in the fracture length even after 10 quadratic elements discretise the fracture (see Fig. 6a). A rigorous solution to such a question remains elusive, however some considerations are offered.
A concept worth noting is that small deviations in the estimated stress intensity factor lead to large changes in the estimated propagation velocity in the toughness-dominated regime. The reverse is true in the viscous-dominated regime. Such a functional relationship is illustrated in the schematic in Fig. 8, and a similar trend is observed in Fig. 4 of Shauer and Duarte [28] (plotted with respect to time instead of propagation velocity). The schematic in Fig. 8 is analogous to considering the change in stress intensity for a change in crack increment , since only linear operations are required to obtain the crack increment from the propagation velocity. The dominant feature to note regards the difference in gradient | ∕ | which is large in the viscous regime, and small in the toughness regime. Such a feature implies that small errors made in estimating the stress intensity factor may translate to larger errors for the propagation velocity in the toughness regime compared to the viscous regime ( > ). Either (i) errors are introduced via the extraction of the stress intensity factor from incorrect numerical values (which is likely to be the case for the HFM asymptote since it relies on numerical apertures, but is known to not be the case for the -integral in the toughnessdominated regime), or (ii) the conversion to velocities = −1 ( ) effectively introduces errors that get significantly magnified in the toughness regime (a feature which may explain the errors in the -algorithm). Overall, such a magnification or errors suggest it may be impractical to extract the velocity of propagation in the toughness-dominated regime if high levels of accuracy are warranted. Consequently, a promising strategy may involve using the −integral to converge on the critical stress intensity factor, without extracting propagation velocities using the −algorithm, since the −integral gives accurate estimations of the stress intensity factor [72].

Buoyancy-driven propagation
In the preceding section, the HFM asymptote was shown to provide an efficient methodology to simulate hydraulic fractures with the finite element method. Furthermore, its simplicity compared to the −algorithm motivates a further evaluation of the HFM asymptote method to more general hydraulic fracturing scenarios. Hence, in this section, we are interested to assess the performance of the HFM asymptote when applying the hydraulic fracturing algorithm to simulate buoyancy-driven propagation; a geomechanical process of significant relevance to the transportation of fluids, heat and dissolved minerals in many subsurface  processes. Consequently, we evaluate the HFM asymptote by simulating the ascent of a hot supercritical fluid at depth. While the -algorithm requires extension to include body forces, this is not the case for the HFM asymptote, pointing to an advantage of the method.
The vertical boundaries are fixed in the direction in order to simulate a uni-axial strain reference state of stress. Consequently, the horizontal stress increases with depth due to gravity according to ℎ = ∕(1 − ) [94], which depends on the vertical stress and the Poisson's ratio . The vertical stress is determined by the height of the body, and the overburden 0 = 54 MPa is applied at the top boundary (see Fig. 9 for schematic), while the bottom boundary is fixed in the direction. Gravity acts in the negative direction with = 0 and = −10 m∕s 2 and Linear elements are used with a discretisation ℎ = 1 m along the length of the fracture, that propagates up to a length of 300 m. The fracture toughness is set to   At late times, a hydraulic fracture containing a buoyant fluid develops a particular structure; a head region develops close to the tip of the fracture, which is responsible for generating the elastic pressures necessary to crack the rock [97,98]; while a tail region develops at the other end of the propagating tip, whose aperture determines the propagation velocity of the fracture [98,99]. For a fracture fed by a constant flux, the aperture of the tail is where the fluids buoyancy = − ℎ ∕ − is determined by the difference between the horizontal and hydrostatic stress gradients. Enacting volume balance calculations, along with Eq. (12), the propagation velocity of the fracture can be calculated to be where integrating Eq. (13) gives the analytical solution for the fracture length. Snapshots of the horizontal stress profile are plotted in Fig. 10 for times = 46 , 129 s and 201 s, where displacements in the horizontal dimension are magnified by a factor of 1200. The stress visualised at the fracture tip qualitatively reproduces the characteristic infinity shape, however, true tensile stresses are not observed at the fracture tip due to the interpolation of stress to the nodes, along with the use of linear elements. Resolving the stress singularity at the fracture tip is notoriously difficult, hence why aperture-based approaches to extract the stress intensity factor are more common compared to approaches which use numerical values of stress [100].
The HFM asymptote accurately captures the length and propagation velocity of the ascending buoyant fluid. This can be seen in Fig. 11 by the great match between the analytical and numerical solutions for the fracture length at late times during the simulation. At early times, a mismatch occurs since the head and tail of the fracture are not fully developed, therefore the constant-velocity analytical solution Eq. (13) is not valid at such a moment in time. The accuracy of the method is best assessed by looking at Fig. 11b. In particular, once the propagation velocity approaches that of the analytical solution at around = 200 s, it oscillates around the correct solution by less than approximately 5%. Such a feature is small enough to not significantly impact the solution for the fracture length.

Conclusion
The work herein has showcased the use of the HFM asymptote within the finite element method in order to predict the location of the front of a hydraulic fracture throughout the viscous-toughness regime. A finite element implementation of a front-tracking algorithm that uses the HFM asymptote is constructed in a similar manner common among implementations using the Boundary element method, or the Displacement Discontinuity Method [17,20,31]; the numerical solution for the aperture of the hydraulic fracture is inverted to predict the location of the fracture front, and therefore its propagation velocity. To evaluate the performance of the method without relying on costly re-meshing strategies, the time step is modified such that the fracture's length conforms with the underlying mesh.
Such a finite element implementation of the HFM asymptote is then evaluated against a recently-developed, energy-based algorithm used to predict the location of the fracture front [72]. The HFM asymptote was found to achieve accurate results in the viscous regime for both linear and quadratic elements even with coarse discretisations of the fracture, where the method marginally outperformed the -algorithm. On the other-hand, in the toughness-dominated regime quadratic elements were found to offer significantly higher levels of accuracy compared to linear elements, with marginally better results obtained by the -algorithm for most discretisations of the fracture. Both aperture and energy-based methods admit modest errors using linear elements in the toughness-dominated regime. An assessment on the sources of error suggest the use of the −integral without extracting propagation velocities holds promise for eliminating such errors. Overall, simulating hydraulic fracture propagation can be done at any scale as long as 10 − 30 elements discretise the hydraulic fracture.
Finally, the simplicity and generality of the HFM asymptote is taken advantage of to simulate buoyancy-driven propagation. While the -algorithm would require modification from its current form to include body-forces, the use of the HFM asymptote can be readily extended. Using linear shape functions, the HFM asymptote was shown to accurately reproduce the length and propagation velocity of a buoyancy-driven hydraulic fracture propagating within a mixed toughness/viscous propagation regime. The simplicity of the HFM asymptote, and its ease-of-use for buoyancy related propagation demonstrates the method to be highly suitable within finite element simulations.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.

Acknowledgements
This research received financial support from SNF grant number 200020_204499, 'Key problems of heat and mass transfer in the earth's crust'. The authors would also like to thank Brice Lecampion, Leonid Germanovic and Andreas Möri for insightful discussions related to hydraulic and gravity-driven fracturing.

Appendix. Semi-analytical solutions
We reproduced here the semi-analytical solutions derived by Wrobel & Mishuris [69] and later generalised by Perkowska [93] for the propagation of a hydraulic fracture with non-Newtonian fluid rheology. The solution is constructed for the viscosity-dominated regime, and is termed 'improved approximation' by the authors. We report the case of a Newtonian fluid herein, and the fracture is assumed to have zero initial length ( = 0 in the authors' formulation). The fracture half-length where represents the beta function and 2 1 denotes the Gauss hypergeometric function. The remaining constants arê0 = √ 3,