A Novel Equivalent Continuum Approach for Modelling Hydraulic Fractures

: Hydraulic fracturing has transitioned into widespread use over the last few decades. There are a variety of numerical methods available to simulate hydraulic fracturing. However, most current methods require a large number of input parameters, of which the values of some parameters are poorly constrained. This paper proposes a new method of modelling the hydraulically fractured region using void-ratio dependent relation to deﬁne the permeability of the fractured region. This approach is computationally efﬁcient and reduces the number of input parameters. By implementing this method with an equivalent continuum representation, uncertainties are reduced arising from heterogeneity and anisotropy of earth materials. The computational efﬁciency improves modelling performance in stress sensitive zones such as in the vicinity of the injection well or near faults.


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][6][7][8][9][10][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][13][14][15][16][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][22][23][24][25][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][30][31][32].

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.

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 agreement 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.

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.

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.

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.

FEM-DEM
FEM-DEM is a hybrid approach that combines the advantages of DEM with finite element modelling (FEM) to reduce the computational load. The DEM component is applied within model micro-scale sections, and the stresses resultant from the DEM are then imposed within the FEM, which models the macro-scale domain [38]. This method is most applicable if the domain can be discretised into high-and low-sensitivity subdomains. It should also be mentioned that FEM-DEM codes are available that combine these two methods for different purposes [27,28].

Equivalent Continuum Methods (ECM)
Equivalent continuum methods assume that within a representative element volume (REV) the medium can be treated as homogeneous. Instead of defining fractures individually and modelling the strength degradation, this means that the properties of the domain are generalised and modelled as a continuum. The region of interest can then be discretised based on the regions of high and low fracture density, allowing for the modelling of the fractured and virgin zones of the reservoir [18,39]. The use of an ECM representation makes it simpler to apply well established numerical approaches such as the finite element method (FEM) to solve governing equations that describe the physical processes and behaviour of the medium.
Most of previous applications of ECM are focused on soils and near-surface materials [29][30][31][32], rather than for modelling hydraulic fracturing in rocks at a greater depth. This may be due to the general availability of a wide range of empirical methods to define the properties of fractured materials at shallow depths [40,41].
A novel ECM method introduced in this paper attempts to address the necessity of optimised hydraulic fracturing modelling. The introduced method in this paper uses the dependency of the permeability on void ratio to define the hydraulic fracture propagation. This method allows for an easier convergence due to uniform mesh dimensions and the exploitation of an already present coupled fluid flow and continuum relations.

Average Volumetric Strain in Fractured Zone
During hydraulic fracturing, it is a common practice to record microseismicity, injection pressure, and rates. These data are then used to analyse the effectiveness of the hydraulic fracturing program. One type of analysis used in this paper is a comparison of estimates of volumetric strains based on microseismicity and injection pressure. In estimating the volumteric strain using the injection pressure, rock is assumed to be cohesion-less and fractures to open when change in injection pressure equals minimum in-situ stress. Thus, the volumetric strain using injection pressure is given by σ min /K, where σ min is the minimum in-situ effective stress and K denotes bulk modulus of the fractured domain. The volumetric strain using microseismicity is estimated using V inj /SRV, where V inj is the injected volume and SRV represents stimulated rock volume (SRV). SRV is defined as a three dimensional rock volume where new fractures and opened natural fractures during hydraulic fracturing are contained [42]. The uniqueness of the strain when fracturing occurs suggests that both of the volumteric strains should exhibit same solution. Thus, this uniqueness allows one to relate the hydraulic fracturing permeability to volumteric strain which then can be related to the void ratio.

Void Ratio Dependent Permeability
Permeability and porosity of a medium are often closely related. Relationships between these two quantities are widely used in reservoir engineering through semi-log porosity-permeability plots [43]. In our ECM method, a similar approach is used to estimate formation permeability. Unlike conventional porosity-permeability relations, our proposed relation uses a step-like function. This provides a good representation for hydraulic fracturing, due to steep spatial gradients of permeability in the vicinity of a hydraulic fracture [17]. Practicality and robustness of this approach arise from the availability of built-in relationships in ABAQUS FEM between permeability and void ratio [34]. This functionality is described by Equation (1).
where η is permeability, a is a maximum permeability, e is the void ratio, e thresh is the minimum void ratio after which fracture permeability starts to dominate, and α is a factor related to rock brittleness and fracture density. Void ratio is a scalar quantity while permeability is a tensor. Therefore, the anisotropy of permeability, i.e., in which it varies with direction, must be handled through transformation. In our implementation, permeability anisotropy is introduced using a normalisation and denormalisation method proposed by [44]. In this method, first the microseismic data are normalised using the dimensions of the SRV using Equation (2), where X, Y, and Z are coordinates of detected events in the microseismic cloud along the maximum, minimum, and intermediate stress directions, respectively; L max , L min , and L int specify the SRV length parameters along the maximum, minimum, and intermediate stress directions, respectively; L scal = 3 √ L max L min L int , and X trans ; Y trans , and Z trans are transformed coordinates of detected events in microsesimic cloud.
To illustrate this transformation, assume a notional SRV characterised by an elliptical microseismic cloud in a 2D plane. The elliptical shape implies that the resultant fracture propagation rate and formation permeability are both anisotropic. After normalising the microseismic data in this canonical example, the resultant microseismic cloud will be circumscribed by a circular boundary. These steps are graphically depicted in Figure 1. This approach allows for a simplification of the analysis into a more straightforward one-dimensional analysis. The permeability function that is obtained using the resultant normalised data are denormalised using Equation (3), where, k 0 is the permeability estimated using one-dimensional analysis and k x , k y , and k z are permeabilities in X, Y, and Z directions, respectively. To provide a case study of our ECM approach, the proposed relation was applied to the Hoadley field data to demonstrate the applicability of this method using real data.

Case Study-Hoadley Tight Sand Reservoir
The Hoadley gas field is located northwest of Red Deer and southwest of Edmonton, Alberta, Canada. The pay formation is the Glauconite tight sand member of the Cretaceous Mannville Group, which is overlain by the coal-bearing Medicine River Formation. Open-hole, multi-stage hydraulic fracturing was carried out in two horizontal injection wells and monitored using a 12-level geophone array from a nearby vertical well [45]. The configuration of the wells and directions of maximum and minimum in-situ stresses are shown in Figure 2a. In Figure 2a blue and red line and, the big blue dot indicate two injection and vertical observation wells, respectively. The scattered dots are the observed microseismicity during hydraulic fracturing. The different colour is assigned for microseismicity observed at each stage. Simplified stratigraphy of the Hoadley field, as well as a depth histogram showing the number of microseismic events at given depths, is provided in Figure 2b. This paper is focused on a single stage of the hydraulic fracturing program in well 1. The injection pressure and rate profiles for stage 12 hydraulic fracturing in well 1 are plotted in Figure 3.  Furthermore, the volumteric strain using the injected volume and injection pressure methods are plotted in Figure 4. The other field data are also included for comparison [45][46][47][48]. The trend shown in Figure 4 suggests unique volumetric strain when fractures open, which allows relating the fractured rock permeability to volumetric strain. Despite most of the microseismicity apparently concentrating above the target unit, this paper models only the Glauconite member. We focus on the in-zone microseismic response due to the abundance of geomechanical property data for Glauconite member, whereas other units are depicted through correlations only. The input parameters of the Glauconite member for ABAQUS FEM are provided in Table 1. Poisson's ratio 0.23 [18] The first step in ECM simulation of the hydraulic fracturing is to normalise the microseismic data to reduce the analysis to one-dimensional form. Prior to normalisation, the microseismic data was rotated to align with maximum horizontal stress. Then the coordinates were translated to be relative to the injection port of stage 12. The planar view of rotated and translated microseismic event coordinates for stage 12 are illustrated in Figure 5a. The values for L max , L min , and L int are 150 m, 50 m, and 75 m, respectively. The normalised microseismic cloud for stage 12 is presented in Figure 5b. The void ratio dependent permeability relation is constructed to match the normalised microseismic cloud. The resulting relation is shown in Figure 6a. As outlined above, this relationship is then denormalised to obtain anisotropic permeability. The resultant permeabilities along the maximum, minimum, and intermediate stress directions, respectively are plotted in Figure 6b. The anisotropic permeability relation is then used in ABAQUS to combine with the continuity equation, including the stress change effect. The constructed model uses C3D8P three-dimensional pore pressure solid element with an approximate mesh size of 20 m. The modelled region dimensions are 1000 m by 250 m by 225 m at the maximum, minimum and vertical stress directions, respectively. The modelled region dimensions are much larger than the SRV to avoid boundary effects. At the outer and top surfaces of the model, a constant pressure boundary is applied to ensure a constant far-field effective stress effect. This injection rate is only 12.5% of the average injection rate since the constructed model uses symmetry boundary along three surfaces. The inner and bottom surfaces use a symmetry boundary. The injection is defined using concentrated fluid flow load on a single node at a constant rate of 0.015 m 3 /s. In-situ stresses and pore pressure are introduced through pre-defined fields in ABAQUS. A constant value at 1900 m is applied for the in-situ stresses and pore pressure. The boundary conditions and injection port locations are indicated in Figure 7 using a simplified sketch. The results from the ECM are compared to field data, as well as the PKN and KGD models. The injection pressure profile obtained from three models and field observation are provided in Figure 8. The PKN and KGD models account for the intact rock cohesion of 11.55 MPa. Furthermore, to account for the natural fractures in the formation, the leak-off coefficient is introduced into analytical solutions. The leak-off coefficient is estimated based on fracture half-length of 150 m at the end of 38 min injection period. As indicated in Figure 8, the analytical models slightly overestimate injection pressure while ECM shows a good match with reported data. The analytical and ECM estimates for fracture length are compared against observed microseismicity. Fracture length along maximum stress comparison is provided in Figure 9. Besides the previously mentioned methods a diffusivity model estimate is included in Figure 9. The diffusivity model estimates the fracture length growth according to r(t) = √ 4πDt where r(t) is the fracture length as a function of time in m, D is diffusivity in m 2 /s and time in s. The estimated diffusivity using the microseismicity along maximum stress direction is 4 m 2 /s. This indicates the case for fracture length increment with assumed constant permeability. As observed from Figure 9, both analytical solutions show the same fracture length elongation profile as the estimation of fluid leak-off coefficient based on the assumed final length of the fracture. For the ECM it is assumed that microseismicity is observed at a given node when pore pressure reaches 12.12 MPa. This satisfies the Mohr Coulomb failure criterion for Glauconite formation [18]. The ECM estimate for fracture length is lower compared to other methods. This is possible due to some of the microseismicity being not due to pore pressure change but due to total stress change around the increased pore pressure region.
Furthermore, the fracture length estimates along minimum in-situ stress direction are compared in Figure 10. As the PKN and KGD models are one-dimensional solutions, they do not show any fracture extend along minimum in-situ stress direction. In addition, the diffusivity model assumes a radial flow. Therefore, it overestimates the fracture extend along the minimum in-situ stress. The ECM approach underestimates the fracture length for the minimum in-situ stress direction, too. This further confirms the possibility of some of the microseismicity due to total stress change effects. The model can be calibrated to match better the microseismic cloud however, the injection pressure profile would differ from the field observation. Considering the higher reliability of injection data compared to microseismicity, it is more preferential to match the injection data.
Another advantage of the proposed ECM has optimised modelling efficiency. The ECM and XFEM models were run using a personal computer with Core i7 2.4 GHz quad-core processor. The ECM required 8768.1 s of CPU time to model 3600 s (or 60 min) of injection period with no divergence problem. However, XFEM required about 28,800 s CPU time before running into convergence problem after 300 s (or 5 min) of the injection period. Therefore, for given generalised data set ECM approach poses as a more efficient modelling method. The XFEM might become a more preferred option if more detailed geomechanical data will allow for a better optimisation of this method.

Conclusions
Based on increasing interest in the use of optimised modelling techniques to simulate hydraulic fracturing, we developed a robust and computationally efficient equivalent continuum method. Our ECM method implemented using ABAQUS applies a step-like relationship between void ratio and permeability and incorporates microseismic data to define the stimulated volume of a reservoir. Anisotropy is accounted for by applying a transformation that enables the use of a standard one-dimensional modelling approach with a uniform mesh. This method is particularly useful for cases where there is insufficient data available to define parameters for other numerical approaches, including stress/strain sensitive zones in the target formation. In addition, the use of representative element volume allows the homogenisation of rock elastic and fracture properties for a modelled region. The proposed ECM is more representative compared to PKN, KGD, and diffusivity models and more efficient compared to the XFEM approach. A case study of hydraulic fracturing of tight sand in the Hoadley gas field, Alberta, illustrates the utility of our method for characterising hydraulic fracturing.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: