Efficient Modeling of Hydraulic Fractures with Diamond Boundary in Shale Gas Reservoirs

The e ﬀ ects of hydraulic fractures with complex boundaries on the well performance of a shale gas well, considering a more realistic corner point geological model, have rarely been studied previously. In this study, the nonintrusive embedded discrete fracture model (EDFM) method was employed to investigate the e ﬀ ects of di ﬀ erent boundary shapes of hydraulic fracture, coupling a network of sophisticated natural fractures discrete fracture network (DFN) on well performance. First, by implementing the powerful EDFM technology, concepts of two categories (rectangle and diamond) of hydraulic fracture with di ﬀ erent boundaries were designed. Next, the geometric equations de ﬁ ning vertices of multiple rectangular- or diamond-shaped hydraulic fractures in arbitrary coordinate systems were derived. Subsequently, the horizontal well with multistaged hydraulic fractures and sophistically oriented 3D natural fractures was inputted into the reservoir model to perform history matching. After history matching, the results were further analyzed to compare the production forecast from the two categories. The results show that 20-year cumulative gas productions for rectangle- and diamond-shaped fractures are approximately 1 : 237 × 10 8 m 3 and 1 : 486 × 10 8 m 3 , respectively. In other words, the diamond category can produce 20.1% more gas than the rectangle category. For cumulative water production, the diamond category produces 3 : 8 × 10 4 m 3 , as against the 3 : 0 × 10 4 m 3 produced by the rectangle category (or 26.7% more). This implies that the diamond-shaped fractures have the potential to reach the far ﬁ eld region of the reservoir away from the wellbore. This means that more intersections with natural fractures DFN can be achieved, and more drainage area is unlocked. The visualization of pressure distributions and drainage volume was easily shown, and these results further con ﬁ rm that the extent of ﬂ uid drainage by the diamond fracture is larger compared to that by the rectangular fracture given the same total surface area.


Introduction
In the last half decade, shale gas is becoming progressively important for the production of natural gas resources. According to an earlier report [1], estimated shale gas reserves are 49.3, 40.6, 25.7, and 38.9 trillion cubic meters from major countries in North America, South America, Europe, and Asia. The tremendous amount of worldwide reserves and large-scale exploitation of shale gas is of profound impact to alleviate the problems of global resource shortage and environmental pollution. The recent technology of horizontal drilling and hydraulic fracturing unlocks the economic potentials of the shale gas reservoirs.
Through recent years, extensive studies have been conducted to research different aspects of shale reservoirs. From one perspective, extensive experimental and mathematical works have been dedicated to study a fine-detailed physical mechanism observed in the shale reservoir, such as the water flowback phenomenon [2,3], and molecular dynamics, such as adsorption and phase behavior [4][5][6][7][8]. From another perspective, great efforts have also been focused on the macroscale analysis of different modelling workflows and various data sources such as microseismic and geomechanics [9][10][11][12]. Besides these analyses, various aspects of optimization strategies have also been investigated, such as fracturing parameter optimization [13][14][15], well spacing optimization [16], and choke management [17]. Through the tremendous amount of qualitative and quantitative examinations of shale gas reservoirs, the role of shale gas development elevates into a strategical level for various global companies and multiple international governments.
Despite numerous general studies on shale gas reservoir, the role of fractures is cardinal in the shale gas reservoir development. Without accurate and efficient characterizations of fracture systems, it would be impossible to evaluate shale reservoirs holistically. Numerical fracture modeling in unconventional reservoir coupling with a natural fracture system is still a painstaking issue in the industry. Various methods have been proposed to solve this complicated problem. First, there is the method of dual porosity dual permeability model (DPDK). The mathematically difficult principles associated with this method makes it impractical for reservoir simulation. It is also exclusively limited to simulating simple fractures and easy reservoir models [18,19]. The second method is the unstructured grid. Although it is more accurate than DPDK and can handle a more complex reservoir model, the unacceptable computational costs [20] also render it as an impractical approach. The next method is called the local grid refinement (LGR), which is even more accurate than the unstructured grid method. Nonetheless, it is immensely tedious and time consuming to consider a more complex fracture and model [21,22]. As an alternative, the    [26]. 2 Lithosphere embedded discrete fracture model (EDFM) comes to address all the above issues [23,24]. Nevertheless, no work has investigated the EDFM's capability to model complex fracture boundary shapes under the more realistic corner grid model and match field production data. Furthermore, the comparison of history match results and long-term production forecast between a regular fracture with a rectangular boundary and an irregular fracture with a diamond boundary is unprecedented.
In this study, a two-phase (gas/water) flow numerical model is built for fractured horizontal wells of the shale gas reservoir that considers key physical mechanisms including adsorption/desorption and pressure-dependent matrix and fracture permeability. And then, complex hydraulic fracture geometries are created. Specifically, rectangle and diamond boundary-shaped hydraulic fractures are established for fracture simulation. On top of that, comparative results of those two categories of the hydraulic fracture are studied through history matching, production forecasting, and pressure visualizations using actual field data. Through this investigation, the unprecedented corner-point shale gas reservoir modelling through the EDFM method is proven to be highly dependable. More importantly, it is found out that the characterization of the lateral reach of hydraulic fractures and the formation of complex fracture networks with natural fractures are very important in evaluation of long-term shale gas well productivity.

Modeling Complex Hydraulic Fracture Boundary and Natural Fractures
In order to simulate gas transport mechanisms in shale formation, the coupling of seepage and desorption must be considered. The study of the entire flow process is based on the following assumptions [25]: (1) It is assumed that the shale reservoir is a primitive structure, and the shale reservoir is characterized by ultralow water saturation, that is, the original  3 Lithosphere water saturation of the matrix is far less than the irreducible water saturation (2) When the water saturation of the shale reservoir matrix reaches its maximum irreducible water saturation, the imbibition process wherein the fracturing fluid invades into the reservoir is over The flowback water saturation is the initial water saturation during the production process (4) Considering the nanopore structure of shale reservoir, it is assumed that water does not flow in the matrix and only exists in the single-phase flow (gas phase) (5) Adsorption-desorption follows the Langmuir isothermal adsorption model The following governing equations control the general simulation physics in the shale gas reservoir. Transport mechanism function for the matrix system: where γ is a conversion factor, K m is the apparent permeability, μ m is the gas viscosity in the matrix, B g is the volume factor of the gas phase, P m is the matrix pressure, t is the time, ϕ m is the matrix porosity, ρ gsc is gas density in the standard condition, V L is the Langmuir volume, and P L is the Langmuir pressure. Transport mechanism function for fracture system: Gas phase flowing equation: where γ is the conversion factor, K rg is the relative permeability of the gas phase in the fracture system, K f is the fracture permeability, μ f g is the gas viscosity in the fracture, P f g is the gas phase pressure in the fracture, ϕ f is the fracture porosity, S g is the gas saturation, ξ is the shape factor, and q g is the gas production. Water phase flowing equation: where K rw is the relative permeability of the water phase in the fracture system, μ w is the water viscosity in the fracture, B w is the water saturation, P w is the water phase pressure in the fracture, S w is the water saturation, and q w is the water production.
In the previous decade, complex fracture networks were modeled successfully by using EDFM which was proven to have superior efficiency and thus popularity. The EDFM method discretizes each fracture plane into a bunch of smaller portions and then embeds those fracture segments, which are regarded as the additional grids, into the matrix grids. The connection between the matrix mesh and the crack mesh is treated by introducing a nonneighboring connection (NNC). The accuracy of EDFM is better than that of the equivalent continuum model, and the complicated unstructured mesh division process of the discrete fracture network (DFN) is avoided, which greatly reduces the computational complexity.
The founding principle of the EDFM method will be introduced in the following part, which is illustrated graphically in Figure 1. The general concept of NNC transmissibility factor (T a ) is introduced to describe the first three kinds of connections [26]: where T a is an arbitrary NNC from the three kinds of connections specified in Xu et al. [26], k a is the permeability of the arbitrary medium permeability, A con refers to the common surface area between the NNC pairs, and the equivalent distance between the NNC pair is define by d e . All the types of nonneighboring connections are based on this general equation, but different forms of equation variation are maintained by different connection mechanisms.
In addition, the modified Peaceman equation is applied to model the well index for the connection between the wellbore and fractures [26], shown as r e = 0:14 where k f is the fracture permeability, w f is the fracture aperture, r e is the effective radius, r w is the wellbore radius, L is length of the fracture segment, and W corresponds to the height of the fracture segment.
With the powerful EDFM technology, two categories of the hydraulic fracture with different boundaries can be simulated with ease. The first category has a rectangular boundary, and its simplified, projected vertices can be defined as where x rect is the possible x coordinates of the four vertices of the rectangle, x perf is the projected x coordinate of the   Lithosphere perforation, and xf rect is the half-length of the rectangular fracture.
where z rect is the possible z coordinates of the four vertices of the rectangle, z perf is the z coordinate of the perforation, and hf rect is the height of the rectangular fracture. The illustrative diagram of the rectangular hydraulic fracture is provided in Figure 2.
The second category has a diamond boundary, and its simplified, projected vertices can be defined as where x diam is the projected x coordinates of the side vertices of the diamond, x°d iam is the projected x coordinates of the middle vertices of the diamond, x perf is the projected x coordinate of the perforation, and xf diam is the half-length of the diamond fracture.
where z diam is the projected z coordinates of the middle vertices of the diamond, z°d iam is the projected z coordinates of the side vertices of the diamond, z perf is the projected z coordinate of the perforation, and hf diam is the height of the diamond fracture. The illustrative diagram of the diamond hydraulic fracture is provided in Figure 3. After a comprehensive list of simplified rectangle/diamond boundary coordinates is defined, it is crucial to define the surface areas of these categories of fractures, so that equivalent surface area can be guaranteed for both rectangle fractures and diamond fractures. The area of the rectangle fracture is displayed as follows first: where A rect is the surface area of the rectangular fracture, xf rect is the half-length of the rectangular fracture, and hf rect is the height of the rectangular fracture. The following equation shows the area of the diamond fracture: where A diam is the surface area of diamond fracture, xf diam is the half-length of the diamond fracture, and hf diam is the height of the diamond fracture. It is easily deduced that if A rect and A diam are set equal to each other, the condition of xf diam = 2xf rect must be satisfied if the heights of both rectan-gular and diamond fractures are controlled by layer laminations, or in other words, hf rect = hf diam . After the simplified and projected formulations are derived, the derivation of vertex coordinates of multiple rectangular and diamond fractures with arbitrary spatial orientation and arbitrary 3D configurations is introduced. The formulation is based on the work of Liu [27] but extended in this work. The first step is to find the coordinates of two endpoints of the middle axis coordinates of the rectangular fracture, as shown in the following: where x mid , y mid , z mid are the x, y, and z coordinates of the middle axis endpoint coordinates of the fracture; x u , y u , z u are the x, y, and z coordinates of an adjacent well trajectory point upstream (having a lower measured depth than) the perforation of this fracture; x d , y d , z d are the x, y, and z coordinates of an adjacent well trajectory point downstream (having a higher measured depth than) the perforation of this fracture; MD f is the perforation of the fracture; MD u is the measured depth of the upstream well trajectory point; MD d is the measured depth of the downstream well trajectory point; xf is the half length of the rectangular fracture; γ is the fracture rotation angle; and α is the fracture strike angle. The concept of the rotation angle is illustrated in Figure 4(b). The purpose of rotation angle is to fit hydraulic fractures into the corner point model, since the corner point model's depth and thickness are nonuniform.
Only after endpoint coordinates of the fracture middle axis fare obtained, coordinates of all fracture vertices can then be calculated as follows: where x rect , y rect , z rect are the four x, y, and z coordinates of the rectangular fracture vertices; x mid , y mid , z mid are the obtained middle axis endpoint coordinates of the fracture; hf is the height of the rectangular fracture; β is the fracture dip angle; α is the fracture strike angle; and γ is the fracture rotation angle. It is important to note that for different fractures along the entire well, different values of xf and hf can be assigned to consider the engineering variation and underground heterogeneity. Figure 4 explains the concept of generating multiple rectangular fractures and the idea of fracture rotation.
As the vertices of the rectangular fracture are defined, the vertices of the diamond fracture can be easily derived. The coordinates of the diamond fracture vertices are given as where x diam , y diam , z diam are the coordinates of four vertices of the diamond fractures; x i rect , y i rect , z i rect are the coordinates of one arbitrary vertex (i) of the rectangular fractures; and x i+1 rect , y i+1 rect , z i+1 rect are the coordinates of the adjacent vertex (next to vertex i) of the rectangular fractures. Figure 5 explains the concept of generating multiple diamond fractures and the idea of fracture rotation. Since the EDFM supports the format of inputting multiple vertices of any complex fracture with a sophisticated boundary, results from our formulation can be easily modelled by EDFM.

Basic Shale Gas Reservoir Model
In this study, a multistage fractured horizontal well with real trajectory coordinates and two categories of fracture boundaries (rectangular and diamond) are modelled by the EDFM method to study their effects on well productivity. Based on the information from the geology team working on an actual field, Figures 6(a) Figure 7(a) as the rectangular fracture and in Figure 7(b) as the diamond fracture. As can be observed, the two outer fractures are designed to be longer than the inner fracture within each stage. The reason behind this design is to simulate the heel bias effect typically observed in the field. Furthermore, the effective flowing surface area is set equal to each other for  Figure 8 shows the comprehensive rendering of the reservoir model, hydraulic fracture model, and complex natural fracture DFN model in 2D and 3D perspectives. Through the 2D view, it is clearly observed that the diamond fractures can better reach the region further away from the wellbore. Through the 3D view, the variation of the diamond fractures' height is better visualized.   Table 1 lists the simulation-related properties that were utilized in this study. The simulated reservoir size is 1000 m × 2200 m × 76 m. The grid cell size for the model is 15 m × 15 m in the horizontal directions and 10 layers with different thickness. Thus, the total number of grid blocks is 92,600 (64 × 144 × 10). The initial reservoir pressure is 70 MPa, with a reservoir temperature of 150°C and a reservoir depth of around 3600 m. The matrix water saturation is assumed to be given and fixed at 0.3. The water compressibility is 6:18 × 10 −4 MPa −1 , with the reference pressure at 0.101 MPa. Gas adsorption is also modelled using a Langmuir adsorption volume of 0.002 m 3 /kg, Langmuir adsorption inverse pressure of 0.00029 kPa -1 , and a bulk shale density of 3560 kg/m 3 . The relative permeability curves used in this study are defined as where k rw is the water phase relative permeability, k rg is the gas phase relative permeability, k°r w is the water phase endpoint relative permeability at residual gas saturation, k°r g is the gas phase endpoint relative permeability at irreducible water saturation, S w is the arbitrary water saturation within the fluid-flowing saturation range, S wirr is the irreducible water saturation, and S gr is the residual gas saturation. These 10 Lithosphere relative permeability parameters are fixed in this study and not included in the history matching process. Their values are tabulated in Table 1. Next, a set of a compaction table is also defined for hydraulic fractures as follows: where ϕ P /k P are the reduced hydraulic fracture porosity/permeability at arbitrary pressure P below the initial reservoir pressure P init in units of MPa, ϕ P init /k P init is the original permeability at the initial reservoir pressure, and α/β are the compaction coefficients for porosity and permeability, respectively. ϕ P init /α used in this study are 1.0 and 0.0015, respectively. k P init /β used in this study are tabulated in Table 2.

History Matching.
Gas flow rate was used as the well constraint to match the bottomhole pressure (BHP) and water flow rate, considering a production history of 334 days. The adjusted parameters were matrix permeability, compaction coefficient, fracture water saturation, hydraulic fracture conductivity, natural fracture conductivity, inner fracture half-length, outer fracture half-length, and fracture height. The tuned values of these parameters for both scenarios (rectangle and diamond) are listed in Table 2. The added grid block numbers generated by EDFM to calculate nonneighboring connections are 18,720 for rectangular hydraulic fractures and 20,160 for diamond hydraulic fractures. The elapsed simulation time for a single instance of a history matching trial is around 5 minutes for the rectangular scenario and 3 minutes for the diamond scenario. Figure 9 shows the result of the history matching. As stated above, Figure 9(a) is the gas flow rate constraint for both the rectangular and diamond fracture scenarios. Figures 9(b) and 9(c) display the BHP and water flow rate matching results. Both scenarios produce an excellent match between simulation and the history data.

Production
Forecasting. Gas and water EUR forecasting are the most crucial variables for evaluating well productivity or degree of engineering success of any shale gas well.

Lithosphere
A 20-year production forecasting is studied for both rectangular and diamond fracture scenarios in this section. It is assumed that for both scenarios, the BHP after the history period is directly dropped to the minimum operating constraint BHP of 800 kPa, as shown in Figure 10(a). EUR prediction curves for the gas and water phases of two different scenarios are shown in Figures 10(b) and 10(c). The 20year cumulative gas production for the rectangular fracture is about 1:237 × 10 8 m 3 while the 20-year cumulative gas production for the diamond shape of the fracture is approximately 1:486 × 10 8 m 3 . Apparently, the cumulative gas production for the diamond shape is 20.1% more than that of the rectangle shape. The cumulative water production of rectangular fractures is about 3 × 10 4 m 3 , and that of the diamond fracture is roughly 3:8 × 10 4 m 3 , which is 26.7% more than that of the rectangle fracture shape. These results further substantiate our previous statements regarding the diamond fractures' superior ability to form more complex DFN with far-field natural fractures. In addition, the cumulative water production is slightly higher for diamond fractures than rectangular fractures. This is due to higher rates of water flowback as the well is shut-in for some period and then is put back online because of diamond fractures' stronger support from the far field of the reservoir. Based on these outcomes of the production forecast, it is paramount to characterize the hydraulic fracture's lateral reach and the number of intersections between hydraulic and natural fractures, thus the degree of general DFN complexity in shale gas reservoir. These variables are very sensitive to evaluate the long-term performances of shale gas well.

Visualization of Pressure Distributions and Drainage
Volume. In order to further understand the detailed differences of rectangular and diamond fractures on the production potential, the most straightforward method is to visualize the changes of the gas drainage volume and the fracture pressure drawdown profile. In Figures 11(a), 11(c), and 11(e), drainage volumes of rectangular fractures are visualized for three distinctive time periods at the end of the production history, after 10 years of production, and after 20 years of production, respectively. Moreover, the drainage volumes of diamond fractures shown in Figures 11(b), 11(d), and 11(f) are also simulated for the same time periods which are specified above. It is evident that the drainage volumes of diamond fractures are larger than those of the rectangular fractures at all time periods, especially in the lateral direction. This observation again supports the conclusion of the last section. Figure 12 shows the visualization of the fracture pressure drawdown profile. Consistently, the pressure drawdown is prevalent even in the far-field natural fractures for the diamond fracture scenario but is negligible for the rectangular fracture scenario. Once again, the power visualization tools further reinforce the previous judgements.

Conclusions
This robust flexibility and practicality of the nonintrusive EDFM method, along with an in-house reservoir simulator, is proven reliable through this work to properly simulate a complicated corner point geology model with porosity heterogeneity, complex boundary categories (rectangle/diamond) of the hydraulic fracture with arbitrary spatial configurations, and sophisticated 3D natural fracture DFN. In particular, each fracture can be assigned a rotation angle that can adjust flexibly in accordance with the geological structure of the ups and downs. Meanwhile, the complex geometries of hydraulic fractures can be modeled extremely efficiently, when compared to LGR, which requires a tremendous amount of work if the reservoir model is in the corner point gridding system. Moreover, simulations of two categories of hydraulic fractures are performed on the basis of the history-matching-calibrated model. Based on the obtained results, the following conclusions can be drawn from this study: (1) The investigation of long-term well performance through EDFM, coupling corner point porosity heterogeneity model, real shale gas field production data, rectangular and diamond hydraulic fractures with model-conforming rotation angles and varying half-lengths honoring the heel bias effect, and complex 3D-oriented natural fracture DFN are pioneered in this work (2) History matching and production forecasting of diamond fractures are analyzed, and the resulting comparisons with rectangular fractures are thoroughly studied in this work, which is a first-time research in the domain of shale gas reservoir simulation. A similar match can be achieved by only tuning hydraulic fracture conductivities for both scenarios with 90.82 md-m for rectangular fractures and 20 md-m for diamond fractures (78% lower), because of the enhanced intersections of diamond fractures with the natural fracture DFN (3) The diamond fracture is discovered to have better well productivity than the rectangular fractures based on a 20-year production forecasting, given that the surface areas are the same. Specifically, the cumulative gas productions of diamond fractures and rectangular fractures can reach 1:486 × 10 8 m 3 and 1:237 × 10 8 m 3 (20.1% higher), respectively. Cumulative water productions of diamond fractures and rectangular fractures are 3:8 × 10 4 m 3 and 3 × 10 4 m 3 (26.7% higher), respectively (4) This first-time explorative work reveals that the characterization of the fracture boundary in the real field has a profound influence on the degree of intersections with natural fractures, scope of gas drainage volume, and long-term EUR

Data Availability
All data used to support the findings of this study are available from the corresponding author upon request.