**3. Methodology**

#### *3.1. 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 *Vinj*/*SRV*, where *Vinj* 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.

#### *3.2. 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).

$$\eta(\varepsilon) = a \ast (\tanh((\frac{\varepsilon}{\varepsilon\_{thresh}})^a))^a \tag{1}$$

where *η* is permeability, *a* is a maximum permeability, *e* is the void ratio, *ethresh* 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),

$$\begin{aligned} X\_{\text{trans}} &= X \ast \frac{L\_{\text{scal}}}{L\_{\text{max}}}\\ Y\_{\text{trans}} &= Y \ast \frac{L\_{\text{scal}}}{L\_{\text{min}}}\\ Z\_{\text{trans}} &= Z \ast \frac{L\_{\text{scal}}}{L\_{\text{int}}} \end{aligned} \tag{2}$$

where *X*, *Y*, and *Z* are coordinates of detected events in the microseismic cloud along the maximum, minimum, and intermediate stress directions, respectively; *Lmax*, *Lmin*, and *Lint* specify the SRV length parameters along the maximum, minimum, and intermediate stress directions, respectively;

*Lscal* = √3 *LmaxLminLint*, and *Xtrans*; *Ytrans*, and *Ztrans* 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.

**Figure 1.** Applying the normalisation method by [44] on a notional stimulated rock volume (SRV) microseismic cloud. The notional microseismic cloud is assumed to be within anisotropic 2-dimensional ellipse. After normalisation, the microseismic cloud can be bounded with a circular shape reducing the analysis to a 1-dimensional scale.

The permeability function that is obtained using the resultant normalised data are denormalised using Equation (3),

$$\begin{aligned} k\_x &= (\frac{L\_{max}}{L\_{scal}})^2 \ast k\_0\\ k\_y &= (\frac{L\_{min}}{L\_{scal}})^2 \ast k\_0\\ k\_z &= (\frac{L\_{int}}{L\_{scal}})^2 \ast k\_0 \end{aligned} \tag{3}$$

where, *k*0 is the permeability estimated using one-dimensional analysis and *kx*, *ky*, and *kz* 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.

#### *3.3. 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.

**Figure 2.** Hoadley field well orientation (**a**) and simplified stratigraphy (**b**). Depth histogram of microseismic events. Gray bars indicate total number of microseismic events recorded at a given depth interval during the hydraulic fracturing [18].

**Figure 3.** Field recorded injection pressure and rate for stage 12 hydraulic fracturing.

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–48]. The trend shown in Figure 4 suggests unique volumetric strain when fractures open, which allows relating the fractured rock permeability to volumetric strain.

**Figure 4.** Volumetric strain comparison using the injected volume and minimum in-situ stress approaches. The plotted relation suggests that fracture propagates when strain due to injection equals the strain due to minimum in-situ stress.

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.


**Table 1.** The summary of Glauconite data.

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 *Lmax*, *Lmin*, and *Lint* are 150 m, 50 m, and 75 m, respectively. The normalised microseismic cloud for stage 12 is presented in Figure 5b.

**Figure 5.** The planar microseismic event coordinates for stage 12 rotated to align with maximum in-situ stress direction and translated to be relative to injection port coordinates (**a**). The normalised stage 12 microseismic event coordinates using *Lmax* = 150 m, *Lmin* = 50 m, and *Lint* = 75 m (**b**).

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.

**Figure 6.** The best matching void ratio dependent permeability relation for normalised and isotropic microseismic cloud (**a**). The denormalised void ratio dependent permeabilities along maximum, minimum, and intermediate in-situ stress directions for stage 12 hydraulic fracturing (**b**).

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 m3/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.

**Figure 7.** The boundary conditions and injection port location. The single headed arrows triangles indicate loading directions and symmetry surfaces, respectively. Red dot represents the location of injection port.

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.

**Figure 8.** The injection pressure comparison between Equivalent Continuum Methods (ECM), Perkins–Kern–Nordgren (PKN), and Khristianovic–Geertsma–de Klerk (KGD) models against the field observation. As it is observed, the analytical solutions slightly overestimate the injection pressure while ECM shows a good match.

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*) = √<sup>4</sup>*πDt* where *r*(*t*) is the fracture length as a function of time in *m*, *D* is diffusivity in m2/s and time in *s*. The estimated diffusivity using the microseismicity along maximum stress direction is 4 m2/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.

**Figure 9.** Comparison of fracture length increment along maximum in-situ stress direction using diffusivity, PKN, KGD, and ECM models. The orange dots indicate recorded microseismicity along the maximum in-situ stress direction.

**Figure 10.** Comparison of fracture length increment along minimum in-situ stress direction using diffusivity and ECM models. The green dots indicate recorded microseismicity along minimum in-situ stress direction.
