The Impact of the Near-Tip Logic on the Accuracy and Convergence Rate of Hydraulic Fracture Simulators Compared to Reference Solutions

We benchmark a series of simulators against available reference solutions for propagating plane-strain and radial hydraulic fractures. In particular, we focus on the accuracy and convergence of the numerical solutions in the important practical case of viscosity dominated propagation. The simulators are based on different propagation criteria: linear elastic fracture mechanics (LEFM), cohesive zone models/tensile strength criteria, and algorithms accounting for the multi-scale nature of hydraulic fracture propagation in the near-tip region. All the simulators tested here are able to capture the analytical solutions of the different configurations tested, but at vastly different computational costs. Algorithms based on the classical LEFM propagation condition require a fine mesh in order to capture viscosity dominated hydraulic fracture evolution. Cohesive zone models, which model the fracture process zone, require even finer meshes to obtain the same accuracy. By contrast, when the algorithms use the appropriate multi-scale hydraulic fracture asymptote in the near-tip region, the exact solution can be matched accurately with a very coarse mesh. The different analytical reference solutions used in this paper provide a crucial series of benchmark tests that any successful hydraulic fracturing simulator should pass.


Introduction
The propagation of hydraulic fractures is a highly non-linear fluid-solid interaction problem involving a moving boundary (i.e., the propagating fracture front in the neighboroughood of which the governing equations degenerate). Simulating this class of problem numerically is challenging, especially properly tracking the evolving fracture front.

Propagation regimes
Plane-strain Radial Viscosity M (K = 0) [3] (with correction for small toughness) [6] Toughness K (K → ∞) [4] (with correction for small viscosity) [6,11]  In geoscience applications, hydraulic fractures propagate in a complex, often poorly characterized medium. Nevertheless, the description of the medium must be simplified in order to apply theoretical models. It is thus crucial that numerical implementations of such models for fracture growth be accurate such that differences from field observations can be attributed to model assumptions rather than poor numerical solution.
In the last ten years, a number of reference solutions (analytical and semi-analytical) have been obtained for propagating plane-strain [1][2][3][4][5] and radial hydraulic fractures [6,7] (see Table 1). These solutions provide invaluable benchmarks for numerical simulators. We compare a number of simulators (2D and 3D) that use different propagation algorithms against these reference solutions for hydraulic fractures driven by a Newtonian fluid under a constant injection rate. For the sake of clarity, we do not address fluid leak-off in our discussion. Of particular interest is the accuracy of the different simulators in tracking the moving fracture front, particularly in the so-called viscosity-dominated regime of propagation.
An outstanding question relates to the convergence and robustness of numerical simulators with respect to the multiscale near-tip behavior of hydraulic fractures. The coupled lubrication (fluid flow) and elasticity equations are known to degenerate near the fracture tip, such that the solution of a semi-infinite fracture propagating at a constant velocity is characterized by a multiscale singular behavior near the tip [8][9][10]. The nature of the dominant singularity depends on the relative importance of two dissipative processes (viscous forces and fracture energy), as well as the reference length scale. Such a multiscale behavior near the fracture tip in turn governs the evolution of the velocity of a finite hydraulic fracture during injection. We discuss the degree to which a numerical simulator needs to include and resolve the near-tip behavior in order to accurately match reference solutions in the light of different benchmarks.
where Q o is the volumetric injection rate per unit length in the out-of plane direction, E ′ denotes the plane-strain elastic modulus, µ ′ = 12µ f is an equivalent viscosity (with µ f the fluid viscosity), and K ′ = √ 32/π K Ic with K Ic the fracture toughness (see [14] for more details). Equivalently a dimensionless viscosity M = K −4 can be used. The complete solutions for the fracture length evolution, fracture width and net pressure have been obtained for the limiting cases of zero dimensionless toughness (equivalently infinite M) and zero dimensionless viscosity (infinite K) First order solutions for either small toughness or viscosity are also available [3,4]. Semi-analytical solutions for any finite values of dimensionless toughness or viscosity are also available [5].
We will restrict our comparisons to the case of relatively small dimensionless toughness (e.g. K < 1), which is known to be the more difficult condition to reproduce numerically. Therefore we express the solution in a so-called viscosity scaling. Because the solution is self-similar the time dependence can be obtained using dimensional analysis. We aim to compare the solutions provided by different numerical codes, which are typically developed in space-time. We thus introduce a dimensionless time τ = t/t c and a scaled coordinate ξ = x/ℓ, where t c is a characteristic time scale and ℓ is the fracture length. The fracture length, opening, and net pressure can be written as follows: where we have highlighted the correspondence between results obtained using a time-based algorithm (say γ, Ω, Π) to the self-similar solution dimensionless solution The dimensionless solution F m (K) = {γ m , Ω m , Π m } for small toughness developed in [3] will be compared with the numerical solutions from different simulators. More precisely, for three small values of dimensionless toughness K = 0.01, 0.1, 0.5, we will focus on the comparisons of the dimensionless fracture length γ m , opening profile Ω m (ξ) close to the fracture tip and error in the fracture volume etc. We are especially interested in the evolution of the error of the numerical solutions with respect to mesh sizes in the near-tip region of the fracture, where gradients are the largest. The solution, for small toughness (K < 1), is in fact governed by the hydraulic fracture viscosity tip asymptote: the opening behaves as w ∼ (ℓ − x) 2/3 and the net pressure as p ∼ (ℓ − x) −1/3 close to the fracture tip (see [10] for details). The tip region affected by the asymptote actually extends to about 10 to 20 percent of the plane-strain fracture length for dimensionless toughnesses below 0.5.

Radial hydraulic fracture
The growth of a radial hydraulic fracture spans both the viscosity and toughness regimes of propagation [6,11]. At early times, the perimeter and the opening of the fracture are small and most of the energy is spent in viscous flow, whereas at a later times, the fracture perimeter and opening are larger and the fracture energy required to extend the fracture dominates the energy required to drive the viscous fluid through the fracture. The radial solution is also dependent only on a dimensionless toughness which in this case is a function of time. Introducing the characteristic time , and the dimensionless time τ = t/t mk we have (see [6,14] for more details): Solutions for the case of zero and infinite dimensionless toughness (i.e. small and large dimensionless time) have been obtained semi-analytically [6]. The complete transient solution can be obtained only numerically. A reference algorithm [7] based on an explicit moving mesh algorithm with proper matching of the multiscale HF tip asymptotics [10] will provide the baseline for the comparisons for intermediate times. The fracture radius R, width w and net pressure p can be written as: where the dimensionless solution F = {γ, Ω, Π} depends only on dimensionless time τ = t/t mk and scaled position ρ = r/R along the fracture.
As before, we will pay particular attention to the case of small dimensionless toughness, i.e. early-time, which is the most challenging numerically. In the limit of zero-toughness/early-time (i.e. viscosity dominated propagation, here refereed as the M-vertex), the solution is self-similar and can be conveniently written in the following viscosity scaling: where again, we have highlighted the correspondence between the zero-toughness/ M-vertex self-similar solution F m0 = {γ m0 , Ω m0 , Π m0 }, which is independent of time and the dimensionless solution expressed as a function of the previously defined dimensionless time.
The zero-toughness/M-vertex solution F m0 has actually been found to correctly capture the propagation of hydraulic fractures up to a dimensionless toughness of K = 1, i.e. for dimensionless time τ ≤ 1 (t ≤ t mk ) [6]. We will thus investigate the convergence of different simulators to this zero-toughness/small time solution.
For dimensionless time above unity (τ > 1 (t > t mk )), the solution transitions from the viscosity dominated (early-time) to the toughness (large-time) dominated regime. The toughness dominated regime is reached for K ≈ 3.5 (τ ≈ 70000). Note that, for infinitely large dimensionless toughness (i.e. zero viscosity), the solution is also self-similar [6] and is also denoted as the K-vertex solution. We will also briefly investigate, for a subset of the simulators considered, the transition between these two regimes of propagation (M to K), focusing mostly on fracture length versus time.

Some practical numbers
Although rock properties and stimulation practices vary, it is interesting to compute the scales and dimensionless numbers previously introduced. Let's assume a "tight" rock with the following realistic properties: a plane-strain Young's modulus of 40 GPa and a fracture toughness of 1.5 MPa.m 1/2 . First, for a hydraulic fracturing treatment using a highly viscous fluid (e.g. gel-like with µ f = 100 cPoise) at a practical rate of 10 Barrels per minute, we obtain a transition time scale t mk for a radial fracture of 4.2 10 6 seconds! The propagation of such a hydraulic fracture for a realistic injection duration (i.e. less than two hours) will always be in the viscosity dominated regime of propagation. Remember that the dimensionless toughness evolves as (t/t mk ) 1/9 for radial fractures. For the plane-strain geometry (where the injection rate is per meter in the out-of-plane direction), using the same parameters we obtain a dimensionless toughness of 0.26 clearly indicating a viscosity dominated propagation.
For a slick-water treatment, popular in shale-gas reservoirs, the injection rates are usually much higher in order to compensate for the low viscosity of water and to obtain a sufficiently wide fracture to accomodate proppant (see the scales in front of the opening w in Equations (1) and (3)). Using a value of 20 Barrels per minute (a realistic value for a single perforation cluster) and the viscosity of water, the radial transition time scale t mk now reduces to two minutes. For radial fractures, the toughness dominated regime of propagation is obtained for dimensionless toughness above 3.5 (see [6]), which corresponds to t 75000 t mk , which translate for this particulate case to t 250 hours. This indicates that most of the duration of the treatment will take place in the transition from the viscosity to the toughness regimes of propagation for a radial fracture (assuming that no stress barriers affect the fracture geometry). For the plane-strain fracture geometry, we obtain a dimensionless toughness of 0.3 (assuming the same injection rate per meter in the out-of-plane direction), for which the propagation is still dominated by viscosity.
These examples show the importance of the viscosity dominated regime of propagation for oil and gas hydraulic fracturing applications. Numerical simulators therefore need to be able to capture this regime of propagation, which is a difficult task especially if the algorithm relies solely on the linear elastic fracture mechanics propagation condition that manifests itself at a length scale near the fracture tip that is much smaller than the modeling lengthscale in the viscosity dominated case [10,15].

Simulators tested
Two classes of simulators have been tested: codes simulating a two-dimensional configuration (plane-strain or/and axisymmetry) where the fracture is a one dimensional geometrical object, and three dimensional codes simulating planar fractures (which are two-dimensional objects in three dimensions). We now briefly describe the algorithms used by these different simulators.

MineHF2D
This simulator (see [16] for more details) handles the propagation of both straight and curved hydraulic fractures in plane-strain. The algorithm is based on a fixed grid. It uses the displacement discontinuity method to solve the elastic equations coupled to a finite difference scheme for the fluid flow within the fracture. This algorithm includes the presence of a fluid lag at the fracture tip; for the simulated case reported here, a large confining stress value was used to minimize the fluid lag (see [17,18] for more discussion on the effect of fluid lag). Explicit time-stepping is used and a volume of fluid method locates the fluid front. The fracture propagation criterion is based on the linear elastic fracture mechanics asymptote. The stress intensity factors are obtained using the displacement method with an adjusting factor of 0.88. A mesh with variable element sizes was used with refinement toward the fracture tip.

FEM_Cohesive
This code is based on a finite element model and the pore pressure cohesive element implemented in Abaqus (Abaqus 6-10.2, 2010). It can handle both plane-strain and axi-symmetric configurations (e.g. plane-strain and radial hydraulic fractures). In this model, a pre-defined surface made up of elements that support the cohesive zone traction-separation calculation is embedded in the rock and the hydraulic fracture grows along this pre-defined surface. The fracture process zone (unbroken cohesive zone) is defined within the separating surfaces where the surface tractions are nonzero. The fracture is fully filled with fluid in the fully damaged cohesive zone (where the cohesive traction is zero) and hence there will be no cohesive traction contribution, but fluid pressure is acting on the open fracture surfaces. So a coupled fluid pressure-traction-separation relationship exists between the cohesive zone defined by the traction-separation law and the pressurised fracture as found from solving the lubrication equation with the constraint that all tractions acting on the entire fracture and the cohesive zone must be in equilibrium. In this cohesive finite element model [19], the irreversible bilinear traction-separation cohesive law is adopted. An incompressible Newtonian fluid is injected at the centre of the fracture at constant injection rate. There is no fluid leak-off through the impermeable surfaces of the fracture, so only flow in the fracture radius direction is modelled. The cohesive elements at the injection point are defined as initially open to allow entry of fluid, and so that the initial flow and fracture growth is possible. Infinite elements surrounding the finite domain, which contains a hydraulic fracture, have been used to model the far-field boundary. Further details of the finite element model can be found in [19].

1DPlanarHF
This code, also based on a fixed mesh, simulates straight hydraulic fractures in two-dimensions (plane-strain and axisymmetric fractures) using a fully implicit scheme to solve for the coupling between the elasticity equation (discretized using the displacement discontinuity method), the fluid conservation (discretized using a finite volume scheme) and to locate the fracture front. An increment of fracture length is given and the corresponding time-step (to reach the new fracture length) is solved by satisfying the fracture propagation condition in the tip element in a weak form: i.e. the volume of the tip element is enforced to be equal to the LEFM square-root asympote (The algorithm is similar to the one described in [20], see also [21]). The HF viscosity tip asymptote [8] can also be used for the case of low fracture toughness, its performance will be compared to the LEFM asymptote. Results obtained using the LEFM and the viscosity asymptotes will be denoted as 1DPlanarHF_lefm, and 1DPlanarHF_m respectively. All the results presented here use a grid with a constant element size (i.e. without any refinements), a re-coarsening of the mesh during the simulation is possible.

EMMA
EMMA is an Explicit Moving Mesh Algorithm for radial geometry, which embeds the proper multiscale tip asymptotes of the hydraulic fracture depending on its velocity (see [7] for more details). It is extremely accurate and the moving mesh nature of the algorithm allows it to span more than ten orders of magnitude in dimensionless time. It notably provides a good solution for the transition between the viscosity and toughness dominated regime of propagation for a radial fracture.

Three dimensional codes
1. The Implicit Level Set Algorithm (ILSA) This algorithm [22] models the evolution a hydraulic fracture with an arbirarily shaped boundary that is assumed to propagate in a plane (a planar fracture in a 3D elastic medium), which is typically perpendicular to the minimum principal stress direction. The three dimensional elastic equilibrium equations are discretized using the dispacement discontinuity boundary integral method in which the fracture within the plane is represented by constant width rectangular elements that are collocated at element centres.
The Reynolds lubrication equation, expressing the conservations of mass of the viscous fluid contained within the crack surfaces, is discretized using a finite volume method also defined with respect to quantities sampled at the centres of the rectangular elements. At the periphery of the fracture, which may not conform to the structured rectangular mesh, the boundary is represented using a concept of partially filled tip elements that are used to define average fracture widths, which are also sampled at element centres. The distinguishing feature of this algorithm is its ability to locate the fracture free boundary using the asymptotic behavior of the hydraulic fracture width that is applicable at a particular point on the fracture perimeter. The free boundary is located by the following iterative process: given an initial guess for the fracture boundary ∂S, determine the corresponding trial fracture width w and fluid pressure field p f ; in the ribbon of elements that are completely filled with fluid and, which share at least one side with a partially filled tip element, use the trial width values to estimate the distance to the free boundary by inverting the applicable tip asymptotic behavior [10]; use these estimates of the distance to the free boundary as initial conditions for the eikonal equation |∇T(x, y)| = 1, whose level set curve T(x, y) = 0 is the free boundary. The fracture boundary is then moved to the curve T(x, y) = 0 and the iterative process is repeated until convergence is achieved. The algorithm uses the multi-scale hydraulic fracture tip asymptotics solution [10] and thus automatically captures the different type of propagation regimes with relatively coarse mesh. For this paper, a simplified version of the algorithm was also designed to only use the LEFM asymptote (hereafter denoted as ILSA_lefm) for comparisons with other algorithms (MineHF2D, 1DPlanarHF etc.). In this version, we adapted the ILSA code to damp the front advance by rescaling the level set function T(x, y) so that the the maximum distance between any point in the ribbon elements and the damped free boundary is no more than three element lengths. This sequence of damped front positions enables the trial widths to be relaxed until fracture width profile presents a close approximation to the viscosity dominated solution, in spite of the fact that the tip elements are, by the nature of the ILSA_lefm algorithm, locked into the LEFM asymptote.

HFLattice:
This code [23] simulates fracture propagation without limitation of shape, direction or number of fractures, as well as slip and opening along pre-existing joints. A 3D lattice formulation is used for simulation of deformation and fracturing. The lattice is a quasi-random assembly of nodes connected by non-linear shear and normal springs. The lattice resolution is given by the average node spacing. Newton's law of motion (for translation and rotation) is solved at the nodes using an explicit central difference scheme. The normal force in the spring is tested and a micro-crack is formed when breakage is detected (spring strength is adjusted to give the correct rock strength). A macro-fracture that develops in intact rock is thus characterized by an assembly of micro-cracks. Fluid flow and storage are based on a network of fluid nodes, located at broken springs or springs intersected by pre-existing joints, connected by pipes. The fluid network is updated continuously as new fracturing occurs. An explicit fluid pressure scheme is used to solve for fracture and matrix flow. The mechanical and flow models are fully coupled: fracture permeability depends on aperture (i.e. deformation of the mechanical components), fluid pressure affects deformation and strength of the solid model, and deformation of the solid model affects fluid pressures. In the algorithm, the lattice springs carry total forces, which affects force balance and motion. Also, effective stress is considered for joint slip or opening.

The plane strain small toughness benchmark
The solution for a plane-strain hydraulic fracture driven by the injection of a Newtonian fluid at a constant rate is self-similar (i.e. evolves as a power-law of time). Our comparisons here focus on the case of viscosity dominated fractures (K < 0.5, spefically K = 0.1, 0.001) which are the most difficult to simulate numerically. The simulators tested for that configuration are (Table 2): MineHF2D, FEM_Cohesive, 1DPlanarHF_lefm and 1DPlanarHF_m. Figure 1 displays, for the case K = 0.01 the dimensionless fracture opening profiles from the tip of the fracture obtained with the different simulators (at the last step of their simulations) as well as both the analytical solution and the viscosity HF tip asymptote [10]. Figure 2 is similar but for the case K = 0.1 . We first observe that the viscosity tip asymptote covers a region about 10 to 20 percent of the fracture from its tip. The different simulators provide width estimates that all correctly fall on the analytical solution "away" from the fracture tip. The distances from the tip at which the simulators recover the analytical solution appear to depend on both the mesh-size and the type of propagation condition used. For algorithms using the linear elastic fracture mechanics (LEFM) propagation condition (opening as a square root of the distance from the tip), this recovery distance from the tip is larger for coarser mesh sizes. The algorithm using a cohesive zone model appears to need significantly more refinement. In contrast, the algorithm using the viscosity HF tip asymptote (i.e. 1DPlanarHF_m) is able to capture the fracture opening exactly all the way to the tip and with a much coarser mesh than was used for the other computations.
This dependence of the convergence toward the exact solution on the mesh size and propagation condition can be further observed in Figures 3-5, which display the rate of convergence for the fracture length and fracture volume as a function of the ratio of the mesh size over fracture length (i.e. the inverse of the number of elements to discretize the frature for uniform mesh). All simulators converge correctly toward the analytical solution but at very different computational costs. We can see that for algorithms using the LEFM condition or a cohesive zone model, the mesh size required to reach the same level of accuracy is about 20 times smaller than for the algorithm that uses the correct HF viscosity tip asymptote.
In the case K = 0.01, 1DPlanarHF_lefm needs about 400 elements (h/ ∼ .0025) to obtain a relative error of about 1 and 3 percent in the fracture length and fracture volume respectively, while smaller relative errors are already obtained when using 20 elements (h/ ∼ .05) for 1DPlanarHF_m. The cost is even greater when a cohesive zone model is used: about 3000 elements (h/ ∼ 3 × 10 −4 ) are needed to reach a relative error of 1.5 and three percent in the fracture length and fracture volume respectively.

Viscosity dominated regime
Owing to its previously argued importance in practice, we focus on early-time, where the relevant analytical solution corresponds to the case of zero-toughness. The simulators tested for this geometry and regime of propagation are the two-dimensional codes under axi-symmetry; FEM_Cohesive and 1DPlanarHF_lefm, and the 3D codes ILSA, the simplified version ILSA_lefm and the HFLattice model ( Table 2).
The different simulators have been run for different ranges of dimensionless time all within the viscosity dominated regime (τ < 1), and are described in Table 3. The HFLattice Table 3. Range of dimensionless time of the simulations for the radial benchmark simulator has been run for the specific case of zero toughness (in fact zero tensile strength as per the formulation).
The convergence toward the zero-toughness solution as a function of the mesh size can be observed in Figure 7 (convergence of the fracture radius). Similar trends to the plane-strain case can be observed. The convergence requires a much finer mesh for the simulators using a cohesive zone model (axi-symmetric FEM_Cohesive) and the LEFM propagation condition (axisymmetric 1DPlanarHF and ILSA_lefm). ILSA, the only simulator using the appropriate HF tip asymptote, achieves the same accuracy with a much coarser mesh compared to all the other simulators. In particular, the version ILSA_lefm, which uses the LEFM asymptote, needs about an order of magnitude finer mesh compared to ILSA for the same relative error. The FEM_Cohesive algorithm requires a ratio h/R about two to three orders of magnitude smaller than ILSA for the same relative error. The HFLattice model, although less accurate, also exhibits convergence as h/R decreases.
The openings at the last time step of the simulation (refer to Table 3 for the corresponding dimensionless time) in the tip coordinate system are compared to the zero-toughness analytical solution in Figure 8. The FEM_Cohesive algorithm captures the solution away from the fracture tip well, i.e at distance larger than 3-5 % of the fracture radius from the tip. Closer to the tip, the opening from the cohesive zone algorithm appears to slightly overestimate the fracture opening. The results of the algorithms using the LEFM propagation condition (1DPlanarHF_lefm, ILSA_lefm) converge toward the analytical opening as their mesh gets finer. On the other hand, ILSA exactly matches the analytical solution for the opening with a relatively coarse mesh (the last element of ILSA corresponds to a partially fractured element for which the fracture width is also function of the fracture front location within the element).
A comparison of the net pressure profile obtained by the different simulators and the analytical solution is displayed on Figure 9. Similar trends to the fracture opening can be observed. One has to note that the HFLattice simulator approximates the point-source by a finite volume, hence the discrepancy on the net pressure with respect to the point-source solution close to the fracture inlet.  It is interesting to investigate more closely how an algorithm using the LEFM asymptote (e.g. ILSA_lefm) is able to converge to the zero-toughness analytical solution. Consider the case of modeling a radial hydraulic fracture starting at a very small initial time (τ = 10 −18 ), which corresponds to a dimensionless toughness K = 0.01, and is therefore very close to the M-Vertex (zero-toughness solution). If the LEFM asymptote is used in ILSA at this time, the asymptote would dictate that the fracture front needs to be advanced by roughly 10 4 element lengths! As already mentioned, to circumvent this problem, the ILSA_lefm code damps the front advance by rescaling the level set function T(x, y) so that the the maximum distance between any point in the ribbon elements and the damped free boundary is no more than three element lengths. This sequence of damped front positions enables the trial widths to be relaxed until the fracture width profile presents a close approximation to the viscosity dominated solution, in spite of the fact that the tip elements are locked into the LEFM. Indeed, the algorithm settles on a solution in which the LEFM tip widths are very small and contribute very little to the net volume of the fracture, while over the remainder of the fracture away from the leading edge, the widths are locked into the viscous asymptote as dictated by the conservation of volume. Thus the damped front advances continue until the conservation of fluid volume dictates that it should stop at which time it approximates the viscous asymptote. Results from the standard ILSA code with the correct viscous asymptote and damped ILSA_lefm code with the LEFM asymptote indicates that in order to obtain similar relative errors compared to the analytical solution, a mesh size that is an order of magnitude smaller is needed for ILSA_lefm compared to the standard ILSA code (see Figure 7).

Transition toward the toughness dominated regime
Finally, we investigate the performance of a subset of the algorithms on the transition of the solution toward the toughness dominated solution. Such a transition typically happens between τ = 1 (K = 1) and τ ∼ 70000 (K = 3.4). The Explicit Moving Mesh Algorithm (EMMA), the 1DPlanarHF_lefm and the ILSA codes are compared, focusing on the evolution of the dimensionless fracture radius with time. Figure 10 display the results. The fracture radius have been averaged over the fracture footprint for the results of ILSA (3D code), while the other codes are axi-symmetric by assumption.
We clearly see that the different algorithms capture the transition between the viscosity and the toughness propagation regimes extremely well. They are virtually indistinguishable.

Conclusions
In this paper, we have investigated the performance of a number of hydraulic fracture propagation algorithms based on different propagation conditions: LEFM, cohesive zone model/tensile strength and algorithms accounting for the multi-scale nature of hydraulic fracture propagation in the near-tip region. This exercise was made possible thanks to the existence of analytical solutions for both geometries of hydraulic fractures. All the simulators investigated here are able to capture the analytical solutions of the different configurations tested, but at widely different computational costs.