**1. Introduction**

Hydraulic fracturing is a fluid injection process intended to create fractures in order to increase formation permeability [1]. This process is applicable in both vertical and horizontal wells, but is now performed mainly in extended horizontal wells [2] to allows multiple hydraulically fractured zones to be completed in a single well [3]. Each hydraulically fractured zone is then referred to as a stage, and a single extended horizontal well may have more than 20 stages. Geomechanical analysis is crucial in the pre-assessment of the geological formation to ensure the safety and effectiveness of the practice [4]. Moreover, considering the high cost of drilling and completion programs, modelling is a key assessment tool for hydraulic fracturing. There are a multitude of techniques and methods available for the modelling of hydraulic fracturing [5–11].

Analytical approaches are common for modelling hydraulic fractures. These methods are primarily based on the well established Perkins–Kern–Nordgren (PKN) and Khristianovic– Geertsma–de Klerk (KGD) models [12–17]. These models remain highly useful for a quick hydraulic fracturing assessment at a site. The PKN model is better suited for cases of fractures with length equal to or greater than twice the fracture height. These conditions are generally applied to hydraulic fracturing in shale reservoirs [18]. The KGD method tends to bea better representative of cases where the fracture length is equal to or less than its height. Such fracture conditions are mainly applicable to brownfields applications where fracturing is used to revive production [19,20].

Based on these analytical solutions, there are numerous numerical methods available for hydraulic fracture modelling. Among the most common numerical modelling techniques are the cohesive zone method (CZM), extended finite element method (XFEM), and discrete element method (DEM) [9,21–26].

Hybrid techniques such as the finite element method-discrete element method (FEM-DEM) have also started to emerge. Such techniques attempt to combine the advantages of two methods to improve their performance [27,28].

Another type of modelling approach that has received less attention is a class of methods called equivalent continuum techniques. These methods attempt to represent fractured media through a generalised approach that provides descriptive results for a sufficiently large scale. Previous applications of equivalent continuum approaches have tended to largely focus on relatively shallow depths of analysis [29–32].

#### **2. Literature Review on Hydraulic Fracturing Modelling Methods**

The practice of hydraulic fracturing has resulted in an extensive body of literature focused on modelling techniques. Some of these techniques are briefly described below.

## *2.1. Analytical Models*

Analytical solutions tend to be highly generalised and simplified through a set of assumptions. For example, assumptions that underlie the PKN and KGD models are: Homogeneous, isotropic, and linear elastic rock, laminar flow, no effects of additives such as proppant, and simple fracture geometry. The main distinction between PKN and KGD models is the assumed cross-sectional shape of the hydraulic fracture. In the PKN model, the hydraulic fracture is assumed to have an elliptical cross-section, while the KGD model assumes a rectangular shape. The KGD model can easily be converted to radial coordinates, which is useful to represent penny-shaped cracks. General solutions for pressure, width, and length for PKN, KGD, and penny-shaped hydraulic fracture models are widely known [17,33].

The use of analytical solutions can be extended to the formations with natural fractures using the leak-off coefficient. The introduction of leak-off coefficient to the analytical solutions improves the agreemen<sup>t</sup> of planar hydraulic fracture geometry with field data. However, analytical solutions are constrained by the simplified definition of approximate hydraulic fracture geometries. On the other hand, numerical models are capable of simulating complex hydraulic fracture geometries and their effects on host medium.

#### *2.2. Cohesive Zone Method (CZM)*

The cohesive zone method utilises cohesive elements within a continuum model to represent fractures in the rock. Cohesive elements use traction-separation relation for mechanical strength weakening [34]. In the cohesive zone, this relationship defines the energy release equation, wherein if a specific threshold strain/stress is reached, then the cohesive element is relaxed and allowed to deform more freely. In addition, the permeability relation for the cohesive zone uses a parallel plate theory, which is an effective representation for fracture permeability. As outlined below, a recent enhancement of the CZM is known as XFEM.

#### *2.3. Extended Finite Element Modelling (XFEM)*

Extended finite element modelling (XFEM) uses the same fundamental approach as the cohesive zone method, except that the fractured zone is defined differently. In particular, for the XFEM, heaviside, junction, and branch enrichment functions are introduced into the domain. These functions are used to calculate where the fracture tip will propagate and how the deformation of the fracture will be imposed onto surrounding rock [35].

Both CZM and XFEM attempt to model fractures individually, which complicates the computation process due to the increased number of equations in the model. Additionally, the large dimensional difference between the fractures and general modelling domain imposes challenges for model convergence. Another challenge for XFEM approaches is mesh dependency due to the use of enrichment functions, the dimensions of the element specifying the orientation and displacement of fractures become highly mesh dependent. Although the CZM approach does not suffer from the same degree of mesh dependency, it requires a pre-definition of the fracture geometry, which imposes limitations on how a fracture will propagate.

#### *2.4. Discrete Element Method (DEM)*

The discrete element method (DEM) is a modelling approach wherein each mesh element is considered as a separate continuum particle and the intervening space between is treated as fractures [36,37]. These fractures use the same type of physics as in the case of CZM, with respect to strength degradation and fracture initiation. A superiority of the DEM over the CZM is an ability for fractures to propagate in any direction in 3-dimensional space. The convergence rate is improved, compared to XFEM, due to a similar dimensionality of the fractures and continuum elements.

The improvements provided by the DEM approach are based on the treatment of fractures as contact points between the elements. However, this method also has high computational requirements, as it generally requires a large number of model elements in order for the model to be realistic and representative.
