Multiphase Multicomponent Numerical Modeling for Hydraulic Fracturing with N-Heptane for Efﬁcient Stimulation in a Tight Gas Reservoir of Germany

: Conventionally, high-pressure water-based ﬂuids have been injected for hydraulic stimulation of unconventional petroleum resources such as tight gas reservoirs. Apart from improving productivity, water-based frac-ﬂuids have caused environmental and technical issues. As a result, much of the interest has shifted towards alternative frac-ﬂuids. In this regard, n-heptane, as an alternative frac-ﬂuid, is proposed. It necessitates the development of a multi-phase and multicomponent (MM) numerical simulator for hydraulic fracturing. Therefore fracture, MM ﬂuid ﬂow, and proppant transport models are implemented in a thermo-hydro-mechanical (THM) coupled FLAC3D-TMVOCMP framework. After veriﬁcation, the model is applied to a real ﬁeld case study for optimization of wellbore x in a tight gas reservoir using n-heptane as the frac-ﬂuid. Sensitivity analysis is carried out to investigate the effect of important parameters, such as ﬂuid viscosity, injection rate, reservoir permeability etc., on fracture geometry with the proposed ﬂuid. The quicker fracture closure and ﬂowback of n-heptane compared to water-based ﬂuid is advantageous for better proppant placement, especially in the upper half of the fracture and the early start of natural gas production in tight reservoirs. Finally, fracture designs with a minimum dimensionless conductivity of 30 are proposed.


Introduction
Continued development of new technologies has resulted in better exploitation of unconventional reservoirs, leading to efficient and expeditious production of petroleum resources. The increased petroleum production is attributed to advanced technologies of horizontal well drilling and hydraulic fracturing in tight gas reservoirs. Hydraulic fracturing has accounted for producing 50% of the natural gas and 33% of petroleum production in the U.S. [1][2][3][4][5][6][7].
Since inception, hundreds of thousands of fracturing operations have been performed [2,3,[8][9][10][11]. Billions of liters of water have been utilized in these operations yearly In this regard, a frac-fluid is proposed which consists of n-alkanes: pentane to decane (C n H 2n+2 :C 5 -C 10 ) [26,27]. Through this scheme, the use of water, as well as the environmental issues associated with it, can be largely minimized. The proposed fluid can be handled as a liquid at surface conditions and is compatible with the formation fluid as it is a petroleum product and hydrocarbon in nature. The density (ranging between 626-730 kg/m 3 ) and viscosity (0.2-0.702 mPa·s) is also considerably less compared with water, therefore, the fracture closure and flowback can be significantly faster than water-based fluid. A single component, such as n-heptane, or a combination of different components according to different reservoir conditions can be used. As it is considerably less dense than water, the hybrid fluid concept, whereby initiating fracturing with 30-40 m 3 water is followed by light alkanes injection, can help reduce large surface injection pressure requirements. The phase behavior of the individual components of the proposed fluid are shown in Figure 2.
The objective of this research is two-fold. Firstly, the development of a full 3D thermohydro-mechanical (THM) coupled numerical model to completely simulate the hydraulic fracturing process for a multi-phase multi-component (MM) fluid system. Secondly, the application of the numerical model to simulate hydraulic fracturing for the optimization of well x in a tight gas reservoir in Germany with alternative frac-fluid (n-heptane) through sensitivity analysis of the most important parameters. However, it should be noted that hydraulic fracturing methods are not used industrially in Germany.

N-Heptane as Frac-Fluid
The tight gas reservoir considered for optimization in this study (discussed in Section 4) has a temperature range of more than 415 K, therefore n-heptane is selected as injection fluid as it has a boiling point of 371.6 K. This temperature difference between reservoir and injection fluid causes fluid viscosity and density reduction due to the rise in fluid temperature following fluid-rock heat exchange. Some properties of n-heptane are presented in Table 1 and its phase behavior is shown in Figure 3. The density of n-heptane is 689.5 kg/m 3 which is much lower than water. Therefore, to avoid excessive surface injection pressure requirements, a hybrid fluid concept will be utilized in this work, in which fracturing is initiated by injecting a small volume of water followed by an n-heptane injection. Table 1. Important properties of n-heptane [29].

Numerical Modeling
Numerical modelling is a powerful tool to deal with the complex phenomenon of hydraulic fracturing. FLAC3D (fast Lagrangian analysis of continua in 3D) is a numerical modeling software which utilizes an explicit finite volume method to solve a wide variety of geotechnical problems [30]. Since FLAC3D deals with continuum mechanics, a new numerical 3D-model built on tensile fracture (discontinuous media) criterion in the presence of 3D stress state and hydromechanical coupling between fracture and matrix was developed by Zhou and Hou [31]. The effects of stress redistribution after the tensile failure and fluid leakoff to the matrix were all numerically modelled. To simulate proppant transport with gelled fluid, the solid-liquid two phase flow in the fracture was integrated, including proppant concentration, shear rate, fluid viscosity, proppant, and fluid densities, etc. [32]. These models were integrated into FLAC3D. The effect of reservoir heterogeneity on fracture orientation was also investigated in tight gas reservoirs using XFEM and FVM approaches [33]. To numerically study the heat transport in the fracture and heat exchange between the fracture and formation, a new thermal module was added to FLAC3D [34]. To model the THM processes in an MM environment, two well established numerical simulators, FLAC3D and TOUGH2, were coupled to exchange data with non-linear coupling functions [35][36][37]. The in-house upgraded version of FLAC3D, called FLAC3D plus , was then coupled with TOUGH2MP to model hydraulic fracturing in different reservoirs such as oil, gas, and geothermal [37][38][39][40].
To simulate the hydraulic fracturing with n-heptane, it is important to model its behavior in the presence of formation fluids, i.e., water and gas. In addition, the determination of properties such as viscosity, density, and enthalpy, etc., of an alternative fluid at reservoir conditions is also imperative to properly understand its performance. Therefore, there is a need for a three-phase multi-component model which can simulate the flow of n-heptane as a separate phase in a non-isothermal environment. TMVOC (belonging to the family of TOUGH2) can model the MM flow of volatile organic chemicals in a non-isothermal manner [41]. Parameters such as saturation, relative permeability, viscosity, density, specific enthalpy, capillary pressure, and diffusion for every phase are considered and updated at each successful Newton-Raphson iteration. Hence, the appearance and disappearance of phases and components in different phases is modeled with reasonable accuracy.
Thus, numerical modeling of hydraulic fracturing with a variety of fluids (alternative and conventional) can be performed by coupling FLAC3D plus (full 3D rock mechanical simulator) and TMVOC (reservoir simulator). In this approach, a fracture model, MM fluid flow model, and proppant transport model are implemented in the THM coupled FLAC3D plus -TMVOC framework [38,39]. In addition, fracture elements residing in the host matrix elements in a pre-defined path perpendicular to least principal stress are considered.

Mechanical Deformation
The main formulation utilized for mechanical calculations is based on FLAC3D. However, as the fracture is a discontinuous media, the discontinuous displacement due to tensile failure needed to be modelled. The flow in the fracture can be considered between two parallel plates. Continued injection of fluid after fracture initiation leads to the leakoff of fluid into the formation matrix from the fracture. Consequently, the pressure in the fracture changes. A simple bi-wing fracture can be observed in Figure 4. In order to perform the mechanical calculations, FLAC3D plus formulation, which relies upon the elasto-plasticity theory, is used. In this regard, the displacement increment in a time interval is determined by the solution of the equation of motion (Equation (1)). The strain and stress increments in a time step can be determined using the continuum and constitutive equations (Equations (2) and (3)) [30].
The tensile failure of the rock occurs when the injected fluid pressure exceeds the combined effect of the smallest principal stress and tensile strength. The fracture propagates in the path of least resistance, which is perpendicular to the least principal stress. Considering compressive stress to be positive, the failure criteria can be described as [31,32,42]: σ min : minimum principal stress (Pa); σ tens : tensile strength of rock (Pa); and P f p : pressure in fracture (Pa). Due to the fluid injection, the pressure inside the fracture changes, effectuating the deformation of fracture elements. Hence, strain change perpendicular to fracture is induced. This strain increment is dependent upon the pressure inside the fracture working against the normal compressive stress, which is mathematically expressed by Equation (5): The fracture width increment due to strain change, considering the small width of the host element, is calculated as follows: where ∆ε f : strain increment (-); α 1 = K + 4G/3 (rock toughness); P f p (t + ∆t): pressure at new time step (Pa); σ n (t): normal stress at previous time step {Pa); K: bulk modulus of rock (Pa); G: shear modulus (Pa); ∆w: fracture width change (m); and l c : width of the zone perpendicular to fracture (m). Injection pressure that is higher than the closure stress gives rise to width augmentation. However, as the fluid injection is stopped, the pressure inside the fracture decreases and the fluid leaks off to the formation, resulting in width reduction. To avoid complete fracture closure, proppants are injected which keep the fracture open. This requires the inclusion of contact stress in the strain increment Equation (5), which is written as: where, σ ct : contact stress σ ct (t + 1) = 0 i f C ≤ 0.65 and w ≥ w resd σ ct (t + 1) = σ con (t) + α 1 .∆ε o i f C > 0.65 and w < w resd [MPa]; C: proppant concentration (m 3 /m 3 ); ∆ε o : over reduced strain (-); and w resd : residual width (m). The fracture propagation in this model utilizes the concept of subdividing the zones into fractured, partially fractured, and unfractured elements. The zone containing the fracture tip is crucial and is therefore further divided into sub-elements to enhance the precision of numerical modeling ( Figure 5) [43]. Once enough sub-elements are fractured in an element, its status is changed to fractured element and the next element becomes the tip or partially fractured element. The zones beyond the fracture tip are considered unfractured.

Multi-phase Multi-Component (MM) Fluid Flow and Proppant Transport
For the MM flow of fluid in the fracture and matrix, TMVOC formulation is utilized [41]. The mass conservation for different components and phases is given as: M κ : mass accumulation term for component k (mol/m 3 ); F β : mass flux (mol/m 2 /s); β: liquid/gas/NAPL (non-aqueous phase liquids) phase; x κ β : component k mole fraction in phase β (-); and q κ : source or sink (mol/m 3 /s).
Solving the MM flow problem requires space discretization. A model is divided into small discrete blocks and the integral finite difference method is used for averaging. A discrete volume element (Vn) in the system is shown in Figure 6, bounded by the closed surface Γn. The interface area between blocks is represented by A ( Figure 6). Then, surface integrals for mass or energy can be approximated by taking the integral of Equation (8).
The mass of a component k includes its share from all the phases in which it is present. Mass accumulation for the formation (matrix) for different phases in the reservoir is based upon the following relation: where φ: porosity (-); S β : saturation of phase β (-); and ρ β : density of phase β (mol/m 3 ).
The flow simulation can be characterized into three categories, which include flow in fracture, flow between fracture and matrix, and matrix flow. The discrete elements are subdivided into matrix or reservoir elements and potential fracture elements. Darcy law governs the fluid flow in these elements, which can be written for a fluid phase β as: u β : Darcy velocity (m/s); k: permeability (m 2 ); k rβ : relative permeability of phase β; µ β : viscosity of phase β (Pa·s); ∇P β : pressure gradient (Pa/m); and g: gravity (m/s 2 ). Due to hydraulic fracturing, a fracture with a certain width is present. Therefore, the mass conservation equation for a fracture involves the inclusion of fracture width to Equation (8), given by: where, F: fraction of fractured sub-element: no. o f f ractured sub elements total no. o f sub elements (-); and w: fracture width (m).
Fluid viscosity determination is of prime importance for determining the fluid flow behavior in porous media. During the hydraulic fracturing operation, the fluid viscosity is affected by gelling, proppant transport, shear and, especially, temperature variations. For mixture viscosity, the following relation is utilized in TMVOC: For proppant transport and avoiding unnecessary fluid loss, fluids are normally gelled using guar gum. During the fracturing and proppant transport process, the fluid viscosity varies, accordingly the term apparent fluid viscosity is used which depends upon the shear rate, proppant concentration, and temperature differentials, etc. The following relations modified from Torres et al. and Zhang et al. [44,45] are utilized for apparent viscosity to include the effect of shear rate, gelling agent concentration, and temperature variations: where, µ ap : apparent viscosity (Pa·s); µ 0 : zero shear viscosity (Pa·s); γ: apparent shear rate (1/s); n: flow index; a 1 , b 1 , k: constants (-); C g : guar concentration (g/l) E a : activation energy of viscous flow (J/mol); R: ideal gas constant (J/mol/K) T: temperature (K); and T 0 : reference temperature (K). Once the proppant injection starts, the viscosity of the proppant-carrying slurry changes. Barre and Conway [46] presented the following correlation to include the effect of proppant concentration on fluid viscosity: C p : proppant concentration in the slurry (-); C pmax : maximum proppant concentration (-); and a: correlation coefficient.
To model the fluid flow including the effect of shear, temperature variation, gelling agent concentration, and proppant concentration on viscosity, the above equations are combined to give the new model as given by Equation (17) [41,44,45]: where µ ap, µ The proppant-carrying slurry velocity and proppant settling velocity is found by the following set of equations [32]: where The injection of NAPL (n-heptane or other hydrocarbons) into a gas reservoir in the presence of connate water requires three-phase fluid modeling. To model the multi-phase flow, a modified version of Stone's three-phase relative permeability function (available in TMVOC) is applied as defined by the following equations [41,47]: The injected fluid saturation is found as S a = 1 − S g − S w and k rg : gas phase relative permeability (-); k rw : water phase relative permeability (-); k ra : NAPL (alkane: n-heptane) phase relative permeability (-); S g : gas phase saturation (-); S w : water phase saturation (-); S a : NAPL (n-heptane) saturation (-); S gr : irreducible gas saturation (-); S wr : irreducible water saturation (-); and S ar : irreducible NAPL (n-heptane) saturation (-).
Moreover, the three-phase capillary function from Parker et al. [48] is used in the TMVOC for measuring the capillary pressure given by the following relations [41]: where , m = 1 − 1 n , P cga : capillary pressure gas-alkane (n-heptane), and P cgw : capillary pressure gas-water. The fracturing, fluid flow and proppant transport models are implemented in the FLAC3D plus -TMVOC coupling. The data sharing between software and specific calculations performed at a time step can be observed in detail in Figure 7. The computing time of the simulation depends upon the size of 3D model, number of grid blocks, isothermal or non-isothermal process, number of components for injection and production, etc., and therefore may range from a few hours to a few days.

Verification of Developed Numerical Model
The verification is divided into the following three parts:

Fracturing Ability
Hydraulic fracturing was performed in the GeneSys (Generated Geothermal Energy Systems) project in which two wells were drilled. For the purpose of verification, one well, Gross Buchholz Gt1 drilled in 2009, was considered. Hydraulic fracturing was performed to generate high permeability flow paths for heat production. Massive hydraulic fracturing treatment, injecting 20,000 m 3 of fresh water with five injection-shutin cycles, was carried out to create a fracture area of more than 0.5 km 2 [49].
To validate the hydraulic fracturing ability of the numerical model, the first hour of injection was simulated and matched with published results, although 106 h of injection were carried out to complete the massive fracturing operation. The generated 3D quarter model in FLAC3D, as per the available geological data, can be observed in Figure 8a. The generated model lies between the depth of −3287 m to −3850 m and has dimensions of 1700 m × 350 m in the x-y plane, respectively. The reservoir pressure and stress state can be observed in Figure 8b [39,49]. The simulated bottom-hole fracture pressure was matched with published data [38]. The pressure rose initially for fracture initiation by overcoming the minimum horizontal stress and tensile strength of the rock and then normalized during the propagation stage. However, it remained slightly higher than 80 MPa during the simulation period of 1 h (Figure 9).  [39], sim: simulated, meas: measured, and inj_rate: injection rate.

Non-Isothermal Flow
To validate the model for its ability to simulate multi-phase multi-component nonisothermal flow of fluids in porous media, displacement of the Benzene-Toluene mixture in a column with steam was carried out [41]. It involved the injection of benzene and toluene in a column, followed by water flood and subsequent steam flood. The sand column properties can be found in Table 2.  Figure 10 presents the variations in saturations of gas, water, and NAPL phases (threephase). It can be observed that the initial high saturation of NAPL changed and shifted downward in the column due to water injection. In addition, steam injection vaporized the fluid. The NAPL bank was reduced from grid block 38 to 42. Additionally, due to high solubility, benzene was initially depleted. The effect on the temperature of the laboratory column after 5000 s of steam injection could also be observed. The temperature in grid block 36 reached 373.15 K, whereas the temperature in grid block 45 was still at 295.15 K, which was the initial temperature of the column.  [41,50]. Sg: gas saturation, Sw: water saturation, So: NAPL saturation, and Temp: temperature.
Isothermal flow was considered in the first and second stage of injection. The first stage involved the injection of NAPLs into the 15th grid block from the top in the laboratory column. In the second stage, water was injected from the top for displacement. As can be observed in Figure 11, the maximum saturation of NAPLs reached about 40% before the water-flood. Soon after the water-flood, the saturation of NAPLs continued decreasing until the end of the injection to below 5%.

Proppant Transport
The modelling of proppant transport and settling was verified by matching it with published data [32,46,51], through modelling the proppant transport and settling in a large slot. Table 3 summarizes the utilized properties of fluid, proppant, and injection rate, etc.   The profile of the simulated proppant transport and settling was recorded for different time points and the contours are presented in Figure 13. Analogous settling profiles were obtained while comparing them with numerical results [51]. The contours at 5, 10 and 150 s showed a very good match, while a slightly different proppant settling profile was observed at 100 s. In addition, a comparison was also made with the experimental results ( Figure 14). It can be observed that simulated proppant transport and settling was in reasonable agreement with the experimental and numerical published data [32,46,51].

Case Study
In this section, a case study for a wellbore x in a tight gas reservoir in Germany is presented. According to the P/Z analysis, the reservoir had initial reserves of 645 million m 3 . A bottom-hole temperature and pressure of 423.15 K and 67 MPa, respectively, were recorded. Hydraulic fracturing was performed in the year 2000 targeting the Wustrow, Dethlingen and Havel formations. However, there was no significant productivity increase, and the production continued only for a few months with intermittent shutin periods recovering only around 2 million m 3 gas. Considering homogenous conditions, a quarter (1/4) 3D geometric model of the reservoir was generated. Figure 15a,b present the 3D model and pressure and stress profile of the reservoir. It is a multilayer reservoir; the perforations are located at depths of 4671.4-4676.2 m. From breakdown and minifrac testing, the breakdown and frac gradient were determined to be 0.213 MPa/10 m and 0.175 MPa/10 m, respectively, and the closure stress was measured as 80.82 MPa. The formation properties, obtained from core analysis and well logging, are summarized in Table 4. The geometric model was first verified by history matching with previous frac job. According to the pumping schedule, a total of~424.592 m 3 fluid and 120-ton 3465 kg/m 3 density proppants were injected. The injection continued for a period of 150 min. The injection rate remained low until 63 min and then gradually increased to 4.89 m 3 /min. A proppant slug at 5 ppg (pounds/gallon) (=599 kg/m 3 ) was initially placed to remove the existing tortuosity. Afterwards, the proppant injection began at 96 min and was increased to a maximum proppant concentration of 7 ppg (=839 kg/m 3 ) at the end. As per the main frac data, the simulation of the fracking operation was performed with the developed model. Figure 16 presents the bottom-hole pressure (BHP) match of the simulated frac job with measured data. After pumping only 23.2 m 3 of injection fluid, the radiator shaft on the blender sheared off. Thus, the operation had to be stopped and commenced the following day. Due to the presence of still x-linked gel in the tubing, high tubing pressures were recorded and the first 2 m 3 of fluid could only be injected at 0.2 m 3 /min. The crosslinker was switched off and displacement was carried out with linear gel, however, the injection rates remained below 0.6 m 3 /min until the first 60 min. Consequently, a pressure match was not feasible due the abnormal pressure recorded, as depicted in section A of Figure 16. Afterwards, the pressure finally dropped, and the injection rate was increased. It can be observed that a reasonable pressure match was obtained during the main course of injection ( Figure 16, section B).
The simulation was continued after shutin until 10 h to observe complete fracture closure. Figure 17 shows that the maximum fracture half-width reduced from the initial value of 0.007 m to 0.0024 m at closure. The fracture aperture decreased until the fracture walls came in contact with the proppants. Due to large closure time, the proppants settled to the bottom of fracture, depriving the upper half of the fracture of proppants. Consequently, the fracture width was very low at the level of perforations and the propped fracture height was even below the injection level away from the wellbore. Therefore, the fracture could not contribute efficiently towards production and may be a major reason behind the lower production from wellbore x.
In addition, the viscosity contour of the injected fluid can be observed in Figure 17. The injected fluid viscosity reduced from 0.205 Pa·s to about 0.02 Pa·s due to the gel break. The temperature of the fluid in the fracture remained low, especially near the injection zone, as fluid was continuously injected at surface temperature conditions. Thus, the highest viscosity contour at shutin was in the region of the lowest temperature. However, due to the temperature difference between fluid and formation, heat exchange took place, which raised its temperature. The temperature profile showed that the temperature increased from about 40 • C (313K) to the reservoir's initial temperature range, which lead to gel beak. After shutin, the fracture closure, due to leakoff, reduced the fracture volume from an initial volume of 280 m 3 to 56 m 3 (Figure 18).

Sensitivity Analysis with N-Heptane as Frac-Fluid
To understand the effect of important parameters such as injection rate, fluid viscosity, and injection time on hydraulic fracturing with n-heptane, simulation tests were performed with the help of a developed numerical model. For the case of proposed fluid injection, a hybrid fluid concept was utilized and fracturing was initiated with a small volume of water.
The utilized parameters of relative permeability and capillary pressure function are presented in Table 5. Three different injection rates of 4, 6 and 8 m 3 /min were considered to analyze their effect on fracture propagation. The injection program is discussed in Table 6. The increase in injection rate increased the stimulated reservoir volume (Figure 19a,  b). By looking at the fracture width, volume, half-length and height, it becomes clear that doubling the injection rate from 4 to 8 m 3 /min almost tripled the fracture volume with more fracture half-length, height and width. As it was not an ultra-tight gas reservoir, the leakoff of injected fluid was higher. At lower injection rates, fluid had more time to percolate into the surrounding matrix, resulting in lower fracture volume. Therefore, increasing the injection rate increased the fracturing rate. In addition, injecting gelled fluid from the beginning will also reduce excessive fluid leakoff.

Fluid Viscosity
Fluid viscosity is a major parameter which affects the hydraulic fracturing job. To investigate this, fracturing was performed for 0.15 Pa·s, 0.25 Pa·s, 0.35 Pa·s and 0.45 Pa·s fluids at injection rates of 6 m 3 /min and 8 m 3 /min ( Table 7). As was expected, increasing viscosity increased the simulated reservoir volume, and fractures with maximum widths of 0.9-1.4 cm can be generated with higher injection rates and fluid viscosities ( Figure 20).

Reservoir Permeability
Reservoir permeability is an uncontrollable factor which significantly affects fracture propagation and geometry. In ultra-tight reservoirs, the leakoff of the fluid, especially those like the proposed fluid, will be very low due to formation tightness. Resultantly, larger stimulated reservoir volumes, compared to tight or less tight reservoirs, can be generated. Low viscosity fluid at lower injection rates can be utilized in such reservoirs. Three cases are presented in this regard (Table 8). Two cases, case 1 and case 2, bear the original reservoir properties with the only difference being in the fluid injection schedule. In case 2, all gelled fluid was injected, but in case 1 pure fluids and then gelled fluid were injected. In case 3, reservoir permeability was reduced by a factor of 0.01, while maintaining the injection program of case 1. Table 8. Scenarios for the effect of reservoir permeability on fracture growth.

Original-1
Water + n-heptane + gelled fluid n-heptane (0.15 Pa·s) Original-2 Gelled water + gelled n-heptane (0.15 Pa·s) Water + n-heptane + gelled n-heptane (0.15 Pa·s) A significant increase to a fracture half-length of 149 m in case 3, from 51 m and 55 m in case 1 and 2, respectively, was found. Similarly, an almost 500% increase in fracture volume for case 3 was observed (Figure 22a,b).  Figure 23 provides the trends of the fracture volume for the different cases. For case 1, when pure fluids were injected initially, very little fracture volume was created as most of the injected fluid was lost to the formation. In case 2, only high viscosity fluid was injected, therefore, a slightly higher fracture volume was generated. Above all, the largest fracture volume of 200 m 3 was generated with lowest leakoff in case 3 due to the lower leakoff. Therefore, pure fluid at lower injection rates can be utilized in such reservoirs for fracturing.

Fluid Flowback
To perform the fluid flowback analysis, proppants were injected so that the fracture did not completely close upon shutin. An injection rate of 4.89 m 3 /min was maintained for 77 min, during which 375 m 3 of injection fluid and 67 tons of proppant were injected. The details of fluid injection schedule are provided in Table 9. Table 9. Injection schedule for conventional water-based and proposed hybrid fluid fracking.   Figure 25 shows the fracture and injection volume of the fracking operation. The fracture closure for the water-based fracture took place at about 400 min, whereas the fracture closed at 150 min for n-heptane. Due to the difference in heat capacities and volatility, the hybrid fluid caused quicker fracture closure, which could lead to better proppant placement, especially in the upper half of the fracture.
After the fracture closure, the flowback was carried out for a period of 7 days. The reservoir gas and injected frac-fluid saturation on day one and seven days after the hydraulic fracturing operation can be observed in Figure 26. As can be seen from saturation contours, most of the injected conventional fluid (water-based) remained in the reservoir even after seven days of flow back. In contrast, most of the alternative fluid (n-heptane) had flowed back from the propped fracture zone when the reservoir gas had reached maximum saturation (red contour: 95%). Since n-heptane is hydrocarbon and a petroleum product, the phase trapping was minimized and, due to the lower viscosity and density, the flowback was quick, resulting in an efficient fracking job. The speedy flowback of the proposed fluid ensures early production of reservoir gas as compared to conventional fluids [22].

Design Proposals
Generally, stimulation treatments are designed according to fracture conductivity. The higher the fracture conductivity, the better the wellbore productivity due to the ease of fluid flow. The dimensionless fracture conductivity is the ratio of the ability of the fracture to transmit the fluid to the ability of reservoir formation to feed the fracture (Equation (26)) [53].
Dimensionless fracture conductivity has been utilized in the petroleum industry for designing hydraulic fracturing operations. However, they cannot identify the impact of hydraulic connection between propped fracture and perforations on conductivity. Therefore, a better technique would be to include the position and concentration of proppants with reference to injection level (perforations) for predicting the fracture conductivity. Recently, Hou et al. presented the weighted fracture conductivity (Equation (27)) which includes the proppant placement with reference to injection level [54]. However, it can be considered a conservative approach.
where, w : fracture width (m); k f : fracture permeability (m 2 ); d: distance between fracture element and perforation (m); A: fracture element area (m 2 ); c p : proppant concentration (-); c max : maximum proppant concentration; k: reservoir permeability (m 2 ); x f : fracture half-length (m) and n: total fracture elements (-). From the outcome of the sensitivity analysis, hydraulic fracturing with gelled nheptane (0.15 Pa·s) and injection rates of 6-8 m 3 /min for a duration of 105 min were proposed. Two fracture designs are presented in this paper, in which fractures with weighted dimensionless conductivity of 30 and 44 are created with n-heptane. Due to quick fracture closure and leakoff, better proppant placement and improved borehole-fracture hydraulic connection was established in comparison with previous frac job. The fracture geometry and design parameters are presented in Figure 27. At this point, it is important to address the economic considerations regarding the field application of n-heptane (or light alkanes) as a frac-fluid. Studies for the application of hydrocarbon-based fluids have demonstrated economically attractive prospects and shown that comparable or even lower costs will be incurred in comparison with waterbased fluids due to lower well cleanup requirements, better flowback ability, less additives, injection fluid flowback, reuse or sale, no special disposal requirements, and increased productivity [22,[55][56][57][58][59][60]. Moreover, unlike gas-based hydrocarbon fluids such as LPG, there is no need for compression for liquification for the fluid reuse [61]. The boiling point of n-heptane is 371.6 K, therefore, it can be separated from the gas production stream as a liquid at the surface.

Discussion
In this research, a new MM numerical model to perform hydraulic fracturing operation with n-heptane in a 3D stress state was implemented and verified. Three-phase relative permeability and capillary pressure functions were utilized. The model was then applied to a wellbore x in a tight gas field in Germany. After the verification of the model through pressure history match with a previous fracture job, it was found that inadequate boreholefracture hydraulic connection could be created, which can be a major reason behind less recovery. Sensitivity analysis was performed for hydraulic stimulation with n-heptane and the effect of different parameters, such as fluid viscosity, injection rate and reservoir permeability, were analyzed. To avoid excessive leakoff, gelled fluid injection from the beginning of the operation is beneficial. However, higher fluid leakoff is beneficial for quick fracture closure after shutin to hold support agents (proppants) in the upper half of fracture to maintain adequate borehole-fracture connection.
A higher viscosity fluid will generate more fracture width, which may end in less propped fracture heights as a larger width fracture will require more time to close, giving a longer time for proppants to settle at the bottom of fracture due to gravity. Therefore, a reasonable viscosity fluid (such as 0.15 Pa·s), to avoid excessive leakoff and create sufficient width for better proppant placement, is desirable.
Increasing the injection rate increases the fracturing rate compared to leakoff rate. Thus, fractures with half-lengths more than 100 m can be created with an 8 m 3 /min injection rate. The leakoff is low in ultra-tight reservoirs, so pure n-heptane can be utilized for fracture propagation and the gelled fluid may be required for the proppant transport period.
The comparison between flowback of the proposed and conventional water-based fluids has shown encouraging results, where it is possible to recover almost 100% of the alternative frac-fluid within a few days of flowback from the propped zones. This means that the reservoir gas production can start as early as the first day of flowback. Another reason for early reservoir gas production is because of the hydrocarbon nature of n-heptane which is an alkane, like natural gas. Therefore, phase trapping of injected fluid can be largely minimized due to its compatibility.
The design proposals put together from sensitivity analysis can help optimize wellbore x's performance. Proposal 1, with an injection rate of 8 m 3 /min, can generate a fracture with weighted dimensionless conductivity of 44 (Fcd: 180) and more than 100 m half-length fracture. While proposal 2, with a lower injection rate (i.e., 6 m 3 /min), can create weighted dimensionless fracture conductivity of 30 (Fcd: 123). In both designs, excellent hydraulic connection between borehole and fracture is predicted.

Conclusions
With attributes of quick fracture closure, better propped fracture height, minimizing phase trapping and quick flowback, the proposed fluid is recommended for hydraulic fracturing operations in tight and ultra-tight gas reservoirs. However, based upon sensitivity analysis, higher injection rates, such as 6-8 m 3 /min, and injecting gelled fluid will result in frac jobs with better fracture geometries in wellbore x. The alternative frac-fluid can be chosen according to different reservoir conditions and desired fracture geometry and conductivity. The developed numerical model with the capability to simulate hydraulic fracturing process in full 3D approach with MM fluid flow can be utilized to perform hydraulic stimulation with a variety of injection fluids. By promoting the use of alternative fluids, environmental hazards associated with water-based fluids can be minimized, in addition to efficient hydraulic fracturing with better fracture conductivities and borehole-fracture connection. Fluid recovery and reuse for subsequent fracturing operations can make its application economically attractive and even comparable with conventional water-based fluid fracking.
The use of rod-shaped proppants in comparison with spherical proppants will be considered for improved fracture conductivity in the future work.