Modeling Hydraulic Fracturing Using Natural Gas Foam as Fracturing Fluids

: Stranded gas emission from the ﬁeld production because of the limitations in the pipeline infrastructure has become one of the major contributors to the greenhouse effects. How to handle the stranded gas is a troublesome problem under the background of global “net-zero” emission efforts. On the other hand, the cost of water for hydraulic fracturing is high and water is not accessible in some areas. The idea of using stranded gas in replace of the water-based fracturing ﬂuid can reduce the gas emission and the cost. This paper presents some novel numerical studies on the feasibility of using stranded natural gas as fracturing ﬂuids. Differences in the fracture creating, proppant placement, and oil/gas/water ﬂowback are compared between natural gas fracturing ﬂuids and water-based fracturing ﬂuids. A fully integrated equation of state compositional hydraulic fracturing and reservoir simulator is used in this paper. Public datasets for the Permian Basin rock and ﬂuid properties and natural gas foam properties are collected to set up simulation cases. The reservoir hydrocarbon ﬂuid and natural gas fracturing ﬂuids phase behavior is modeled using the Peng-Robinson equation of state. The evolving of created fracture geometry, conductivity and ﬂowback performance during the lifecycle of the well (injection, shut-in, and production) are analyzed for the gas and water fracturing ﬂuids. Simulation results show that natural gas and foam fracturing ﬂuids are better than water-based fracturing ﬂuids in terms of lower breakdown pressure, lower water leakoff into the reservoir, and higher cluster efﬁciency. NG foams tend to create better propped fractures with shorter length and larger width, because of their high viscosity. NG foam is also found to create better stimulated rock volume (SRV) permeability, better fracturing ﬂuid ﬂowback with a large water usage reduction, and high natural gas consumption. The simulation results presented in this paper are helpful to the operators in reducing natural gas emission while reducing the cost of hydraulic fracturing operation.


Introduction
The booming of shale oil development has rapidly increased the amount of flared gas in the past decade in the United States ( Figure 1). Due to the lack of pipelines in the major shale oil production sites, typically operators have no way to utilize the produced natural gas and so it must be flared or reinjected. The CO 2 emissions and toxins produced by the flaring gas threaten human health and the environment. Solutions to the gas flaring Different hydraulic fracturing fluids have been developed for different purposes. Different fluids can be used in the periods of hydraulic fracturing such as fracture initiation, fracture propagation, proppant placement, and flushing. Slick water, linear gel, crosslinked gel, gas, and foam are the most common fracturing fluids and slick water is the most popular fracturing fluid at present. Different hydraulic fracturing fluids have very different properties (viscosity, density, leakoff, the ability to carry sand) depending on their base fluid and the additive properties. The viscosities of several fracturing fluids are listed in Figure 2 and Table 1. Slick water is a Newtonian fluid while the others are shearthinning power-law fluids. It is not uncommon to use foam as a fracturing fluid. The advantages of foam include (1) well production enhancement, (2) less water usage, and (3) better proppant transport [4]. Foam or gas (e.g., CO2, N2, NG) has the potential of being used as the fracturing fluids, which can improve oil recovery starting from the time of well completion while reducing the CO2 or flare emission. This is very important as a step towards carbon-neutral or net-zero emissions. Many researchers have performed experimental studies on foam rheology. [5] studied the viscosity of guar-based nitrogen foam at various foam qualities, temperatures, and shear rates and developed empirical correlations for foam viscosity. [6] measured foam leakoff coefficients and return permeability in a closed flow-loop connected to a dynamic filtration cell under simulated hydraulic fracturing conditions. [7,8] conducted laboratory-scale and field-scale experiments to measure NG foam fluid rheology. Table 1 list their published experimental data for 70% quality foam with different base fluids. The advantages of using NG instead of other gases include: 1. NG is often free. Operators just need to collect and compress the NG from the field while they need to buy other gases (such as CO2, N2) from other sites and transport them to the asset. 2. Stranded natural gas (NG) or flared gas is richer than the gas fed into the pipeline for sale. So, the use of NG or NG foam as a fracturing fluid has advantages over other gases in terms of the phase behavior for some reservoir fluid types. [9] systematically studied the effects of injection fluid composition on the hydrocarbon recovery in a Different hydraulic fracturing fluids have been developed for different purposes. Different fluids can be used in the periods of hydraulic fracturing such as fracture initiation, fracture propagation, proppant placement, and flushing. Slick water, linear gel, crosslinked gel, gas, and foam are the most common fracturing fluids and slick water is the most popular fracturing fluid at present. Different hydraulic fracturing fluids have very different properties (viscosity, density, leakoff, the ability to carry sand) depending on their base fluid and the additive properties. The viscosities of several fracturing fluids are listed in Figure 2 and Table 1. Slick water is a Newtonian fluid while the others are shear-thinning power-law fluids. It is not uncommon to use foam as a fracturing fluid. The advantages of foam include (1) well production enhancement, (2) less water usage, and (3) better proppant transport [4]. Foam or gas (e.g., CO 2 , N 2 , NG) has the potential of being used as the fracturing fluids, which can improve oil recovery starting from the time of well completion while reducing the CO 2 or flare emission. This is very important as a step towards carbon-neutral or net-zero emissions. Many researchers have performed experimental studies on foam rheology. [5] studied the viscosity of guar-based nitrogen foam at various foam qualities, temperatures, and shear rates and developed empirical correlations for foam viscosity. [6] measured foam leakoff coefficients and return permeability in a closed flow-loop connected to a dynamic filtration cell under simulated hydraulic fracturing conditions. [7,8] conducted laboratory-scale and field-scale experiments to measure NG foam fluid rheology. Table 1 list their published experimental data for 70% quality foam with different base fluids. The advantages of using NG instead of other gases include: 1.
NG is often free. Operators just need to collect and compress the NG from the field while they need to buy other gases (such as CO 2 , N 2 ) from other sites and transport them to the asset.

2.
Stranded natural gas (NG) or flared gas is richer than the gas fed into the pipeline for sale. So, the use of NG or NG foam as a fracturing fluid has advantages over other gases in terms of the phase behavior for some reservoir fluid types. [9] systematically studied the effects of injection fluid composition on the hydrocarbon recovery in a volatile oil reservoir and found that richer gas leads to a higher oil and gas recovery. This indicates that using NG from a gas flare as the fracturing fluid is extremely beneficial to the volatile oil reservoir. 3.
The use of NG as a fracturing fluid can also contribute to a reduction in carbon emissions and the carbon-neutral/net-zero goal. Several methods have been proposed to model hydraulic fracture propagation [10]. However, majorities of existing fracture propagation models assume single-phase (usually water phase) flow and can only model one fracturing fluid within one simulation [11][12][13][14][15][16][17][18][19]. Ref. [20] presented a hydraulic fracturing model which incorporated thermal and compositional effects in addition to fracture mechanics and proppant transport. [21] presented a two-phase fracturing model with one fluid phase and one solid phase. The main limitation of these models came from the simplified fracture mechanics. Their models assume the fracture geometry follows the Perkins-Kern-Nordgren analytical results and thus is not appropriate in heterogeneous reservoirs. One simulator, EFrac-3D, is capable of modeling fracture propagation using the finite element methods considering waterbased fracturing fluids (slick water, gel) and energized fracturing fluids (gas, foam) [4]. The compositional 3-D fracturing model is coupled with a simplified wellbore productivity model. However, this model also has some limitations. It only models a single halfwing fracture propagation using far-field stress (without computing the evolving stress field). A single leakoff coefficient value is used in the simulations as it is hard to estimate  Several methods have been proposed to model hydraulic fracture propagation [10]. However, majorities of existing fracture propagation models assume single-phase (usually water phase) flow and can only model one fracturing fluid within one simulation [11][12][13][14][15][16][17][18][19]. Ref. [20] presented a hydraulic fracturing model which incorporated thermal and compositional effects in addition to fracture mechanics and proppant transport. [21] presented a two-phase fracturing model with one fluid phase and one solid phase. The main limitation of these models came from the simplified fracture mechanics. Their models assume the fracture geometry follows the Perkins-Kern-Nordgren analytical results and thus is not appropriate in heterogeneous reservoirs. One simulator, EFrac-3D, is capable of modeling fracture propagation using the finite element methods considering water-based fracturing fluids (slick water, gel) and energized fracturing fluids (gas, foam) [4]. The compositional 3-D fracturing model is coupled with a simplified wellbore productivity model. However, this model also has some limitations. It only models a single half-wing fracture propagation using far-field stress (without computing the evolving stress field). A single leakoff coefficient value is used in the simulations as it is hard to estimate the leakoff coefficients when different phases co-exist. What is more, the equilibrium values method is used for hydrocarbon phase behavior modeling instead of using an equation of state. For the well productivity model, infinite reservoir size, ultra-low reservoir permeability, and infinite fracture conductivity are assumed.
An in-house simulator for integrated hydraulic fracturing and reservoir modeling has been developed with well lifecycle modeling capability [9,[22][23][24][25][26][27]. This simulator solves within three domains (reservoir, fracture, wellbore) with fully coupled physics including porous flow, solid mechanics, fracture mechanics, energy balance, fracture slurry flow, etc. (Figure 3). Options for single-phase, multi-phase black-oil, multi-phase equation-of- state compositional, or water-steam two-phase geothermal flow can be selected for flow modeling. Different solution algorithms are also built-in, and Figure 4 shows an implicit pressure, explicit concentration algorithm that will be used in the paper. ues method is used for hydrocarbon phase behavior modeling instead of using an equation of state. For the well productivity model, infinite reservoir size, ultra-low reservoir permeability, and infinite fracture conductivity are assumed.
An in-house simulator for integrated hydraulic fracturing and reservoir modeling has been developed with well lifecycle modeling capability [9,[22][23][24][25][26][27]. This simulator solves within three domains (reservoir, fracture, wellbore) with fully coupled physics including porous flow, solid mechanics, fracture mechanics, energy balance, fracture slurry flow, etc. (Figure 3). Options for single-phase, multi-phase black-oil, multi-phase equation-of-state compositional, or water-steam two-phase geothermal flow can be selected for flow modeling. Different solution algorithms are also built-in, and Figure 4 shows an implicit pressure, explicit concentration algorithm that will be used in the paper.     In this paper, we extended the capabilities of the model to include the flow of non-Newtonian fracturing fluids (foams and gels) for hydraulic fracturing. We then apply this new capability to investigate the use of natural gas-based foam fracturing fluids. A reservoir layer model and a fluid model were built based on public datasets in the Permian Basin and a well lifecycle analysis was performed. We simulated hydraulic fracture propagation using different fracturing fluids, shut-in and fracture closure, and long-term production. The rest of this paper is structured as follows: The technical details and new In this paper, we extended the capabilities of the model to include the flow of non-Newtonian fracturing fluids (foams and gels) for hydraulic fracturing. We then apply this new capability to investigate the use of natural gas-based foam fracturing fluids. A reservoir layer model and a fluid model were built based on public datasets in the Permian Basin and a well lifecycle analysis was performed. We simulated hydraulic fracture propagation using different fracturing fluids, shut-in and fracture closure, and long-term production. The rest of this paper is structured as follows: The technical details and new model developments in this paper are discussed in the Model Description section. Then the input parameters and settings for the illustrative case are shown in the Case Setup section. Simulation results using different fracturing fluids are presented and analyzed in the Results and Discussion section. Key findings and conclusions are presented in the Conclusion section.

Model Description
Main governing equations embedded in the integrated compositional fracturing and reservoir simulator are introduced first in Section 2.1. Then the details in modeling the foam phase behavior and modeling the stress dependent permeability during fracturing and production are presented in Sections 2.2 and 2.3, respectively.

Main Governing Equations
The simulator used in this paper integrates the reservoir, fracture, and wellbore domains. The main governing equations related to this paper are summarized as follows.
where λ is the Lame first parameter, µ is the Lame second parameter, n is the surface normal vector, u n is the normal displacement, u t is the tangential displacement, σ 0 is the in-situ stress, α is the Biot's coefficient, p is the pore pressure, Γ is the surface area, N i is the moles for component i, V b is the bulk volume, n p is the number of phases, x ij is the mole fraction of component i in phase j, ξ j . k rj , µ j , p crj , ρ j , S j , u pj are the molar density, relative permeability, viscosity, capillary pressure, mass density, saturation, and proppant velocity of phase j. k is the matrix permeability, D is the depth, φ is the porosity, K is the diffusion coefficient, q i stands for the source/sink term which may come from the leak-off from the fracture or a well model, V t is the total fluid volume, V 0 p is the original pore volume, c f is the formation compressibility, w is the fracture width, S is the fracture surface area, q Ini is the source term due to fluid injection/production from/into the wellbore domain, q Leaki is the source term due to fluid leak-off into/from the reservoir domain, V pi , q pi are the volume and volumetric injection rate of proppant i, γ ij is the phase partitioning factor of proppant i in phase j, ρ f is the fluid density, v f is the fluid velocity, c is the proppant concentration, and k f is the propped fracture permeability. Equation (1) (7) and (8) describe the slurry distribution among clusters. Equation (9) formulates the permeability of propped fractures. All the governing equations in our model are discretized in space, using the finite volume and finite area methods and in time using the backward Euler method. The coupled systems of equations are solved using the Newton-type formulations. Readers are encouraged to refer to the previous publications for more details in the governing equations, solution algorithms, model implementation, and model validation [9,[22][23][24][25][26][27].

Modeling of Foam Rheology
In the previous work, the water-based fracturing fluid is modeled as a Newtonian fluid with constant viscosity and the viscosity of gas fracturing fluid is modeled using the Lohrenz-Bray-Clark correlation [28]. However, many hydraulic fracturing fluids such as gel and foam are non-Newtonian fluids, and their viscosity is a strong function of shear rate. Typically, they can be modeled as a shear thinning power low fluid. In this work, we extend the model to include the power-law fluids. We have implemented a general formulation (10) to describe the fracturing fluid viscosity as a function of the pressure gradient and width inside the fracture as well as fluid rheology parameters by [29].
where µ f is the viscosity of fracturing fluid, ∇p is the pressure gradient within the fracture, k and n are the fluid consistency index and behavior index. Regarding the effect of proppant on the viscosity increase of the proppant laden, [29] adopt the method by [30] to describe k and n as functions of proppant concentration. However, the parameters in the functions are different for different fracturing fluids and need to be measured by experiments. In this paper, we keep our correlation (11) to be the same as [31].
where µ s is the viscosity of the slurry and c max is the packing concentration of proppants.
When the simulator is used to model foam, in the fracture domain, there are two phases: a water phase and a hydrocarbon phase, and the two phases are considered to flow as one pseudo-phase. The viscosity of the two phases is calculated using the above viscosity formulation using the input values for k and n. The density of the gas phase is calculated using the EOS. The slurry density is calculated using the volume average of water, hydrocarbon, and proppant densities. Since one pseudo-phase is considered, we assume the foam is stable and there is no velocity difference between the water and gas phase. Thus, the two phases are flowing at the same velocity and proppants are considered to flow in the pseudo-phase.

Stress Dependent Permeability
Changes in permeability and porosity because of the reservoir depletion or recharging are important in stress-sensitive reservoirs [32]. During hydraulic fracturing, in addition to Energies 2021, 14, 7645 7 of 28 the main hydraulic fracture growth, there are numerous micro-cracks created, forming an SRV region, and there is an enhanced permeability in the SRV region. During production, the hydraulic fractures and micro-cracks close, and the enhanced permeability decreases as the pressure decreases and the effective stress become more compressive. In the past, linear relationships, exponential relationships, and fractal models have been proposed to model the dynamically changing SRV permeability [33]. In the exponential relationship, the SRV permeability is usually described as an exponential function of pore pressure or effective stress. In this paper, we adopt the exponential relationship and model the SRV permeability as an exponential function of the effective stress as shown in Equation (12).
where k is the matrix permeability, k 0 is the initial matrix permeability, σ is the effective stress, σ 0 is the initial effective stress, and γ is the permeability-stress exponent that can be obtained from experiments or history matching the field data. An illustration of stress-dependent permeability changing with effective stress during hydraulic fracturing and production is shown in Figure 5. During hydraulic fracturing, the permeability increases because of the formation of micro-fractures around the main hydraulic fractures and a decrease in the effective stress (it becomes more tensile). A larger change in the effective stress indicates a larger stress shadow and more micro-cracks, leading to a larger SRV permeability enhancement. The stress and permeability at the beginning of the production simulation are saved and used as initial stress and permeability for the production simulation. During production, pore pressure decreases, causing the effective stress to become more compressive, and the SRV permeability gradually decreases, representing the closure of the micro-cracks. Different permeability-stress exponents may be used for the hydraulic fracturing and production periods. the SRV permeability as an exponential function of the effective stress as shown in Equation (12).
where is the matrix permeability, 0 is the initial matrix permeability, is the effective stress, 0 is the initial effective stress, and is the permeability-stress exponent that can be obtained from experiments or history matching the field data. An illustration of stress-dependent permeability changing with effective stress during hydraulic fracturing and production is shown in Figure 5. During hydraulic fracturing, the permeability increases because of the formation of micro-fractures around the main hydraulic fractures and a decrease in the effective stress (it becomes more tensile). A larger change in the effective stress indicates a larger stress shadow and more microcracks, leading to a larger SRV permeability enhancement. The stress and permeability at the beginning of the production simulation are saved and used as initial stress and permeability for the production simulation. During production, pore pressure decreases, causing the effective stress to become more compressive, and the SRV permeability gradually decreases, representing the closure of the micro-cracks. Different permeability-stress exponents may be used for the hydraulic fracturing and production periods.

Case Setup
Per estimates made by the U.S. Geology Survey (USGS), the Wolfcamp shale contains 20 billion barrels of oil and 16 trillion cubic feet of gas. The recoverable oil in the Wolfcamp shale is also estimated to be three times larger than the Bakken-Three Forks shale [34].
Given the huge amount of associated natural gas in the Wolfcamp formation in the Permian Basin and the fast growth of shale oil development and flared gas in the State of Texas (Figure 1), we conducted this simulation study using different kinds of fracturing fluids in the Permian Basin. Layer properties and fluid properties are the two key sets of inputs for our simulator. Because of data confidentiality, very limited field data in the Permian Basin is publicly available. Ref. [35] presented fracture geometry characterization in the Permian Basin using integrated microseismic, completion, and production data. Ref. [36] studied the impact of parent well depletion on stress changes and infill well completion. In this paper, we inherit the published layer data from Patterson and Sangnim-

Case Setup
Per estimates made by the U.S. Geology Survey (USGS), the Wolfcamp shale contains 20 billion barrels of oil and 16 trillion cubic feet of gas. The recoverable oil in the Wolfcamp shale is also estimated to be three times larger than the Bakken-Three Forks shale [34].
Given the huge amount of associated natural gas in the Wolfcamp formation in the Permian Basin and the fast growth of shale oil development and flared gas in the State of Texas (Figure 1), we conducted this simulation study using different kinds of fracturing fluids in the Permian Basin. Layer properties and fluid properties are the two key sets of inputs for our simulator. Because of data confidentiality, very limited field data in the Permian Basin is publicly available. Ref. [35] presented fracture geometry characterization in the Permian Basin using integrated microseismic, completion, and production data. Ref. [36] studied the impact of parent well depletion on stress changes and infill well completion. In this paper, we inherit the published layer data from Patterson and Sangnimnuan to build up our geological model. We focus on the Wolfcamp shale formation (A1, A2, A3, B1, B2, B3). Porous properties (porosity and permeability), mechanical properties (Young's modulus and Poisson's ratio), and initial conditions (water saturation and temperature) vary between different Wolfcamp shale layers as shown in Table 2. It is worth noting that in this paper, an isothermal condition is assumed which means the temperature in Table 2 will be used directly. This is because the total injection volume is small and thermo-elasticity is only important when injecting a large amount of cold fluid into the reservoir (e.g., during the water flooding). Our equation-of-state compositional simulation requires a PVT model for the reservoir fluids. [37] studied the effect of fluid phase behavior on well production, in which they reviewed several PVT reports in the Wolfcamp, but only black-oil properties were released in their paper. [38] presented a study on the effect of injection gas composition for huffn-puff improved oil recovery in the Permian reservoirs in which they presented a Peng-Robinson EOS model of a Wolfcamp live oil sample. Component mole fractions (z), critical pressure (P c ), critical temperature (T c ), molar weight (M w ), binary interaction parameters (BIP) of this oil sample are presented in their paper, but other data such as acentric factor (α), critical volume (V c ), and volume shift parameters were not included. In this paper, we adopt the fluid data presented by [38] and added reasonable estimates for acentric factor and critical volume for the Wolfcamp oil sample. The component mole fraction and critical properties are listed in Table 3. The binary interaction parameters between different hydrocarbon components are listed in Table 4.
Stone-II relative permeability model is used to model the relative permeability of different mobile phases and the water and oil relative permeability as a function of the water saturation is shown in Figure 6. Table 3. PVT critical properties of the reservoir hydrocarbon fluid (partially from [38]).   Stone-II relative permeability model is used to model the relative permeability of different mobile phases and the water and oil relative permeability as a function of the water saturation is shown in Figure 6. Other important input parameters are listed in Table 5. In the simulation, the well landing depth is 2577 m, which is in the Wolfcamp B2 formation. We use the same completion design for different hydraulic fracturing fluids, which is 10 perforation clusters per stage and 3 perforations per cluster. The cluster spacing is set to be 15 m. The initial perforation diameter is set to be 1. 27 and perforation erosion and perforation diameter increase are modeled using the methods proposed by [39,40].  Other important input parameters are listed in Table 5. In the simulation, the well landing depth is 2577 m, which is in the Wolfcamp B2 formation. We use the same completion design for different hydraulic fracturing fluids, which is 10 perforation clusters per stage and 3 perforations per cluster. The cluster spacing is set to be 15 m. The initial perforation diameter is set to be 1.27 cm and perforation erosion and perforation diameter increase are modeled using the methods proposed by [39,40]. A computation domain with a size of 960 × 960 × 360 m 3 is set up for fracture propagation and shut-in simulation as shown in Figure 7. The well is landed in the middle of the domain (true vertical depth = 2577 m). Layer properties summarized in Table 2 are assigned to each cell according to the depth. We simulate the fracture propagation and shut-in of one stage. In the production simulation, a reservoir simulation box with a size of 150 × 400 × 280 m 3 is selected within the computation domain and cells outside the box are set as null blocks, which is a representation of the control reservoir volume of this particular stage. Permeability-stress exponent for shut-in and production 5 × 10 −7 Pa −1 A computation domain with a size of 960 × 960 × 360 m 3 is set up for fracture propagation and shut-in simulation as shown in Figure 7. The well is landed in the middle of the domain (true vertical depth = 2577 m). Layer properties summarized in Table 2 are assigned to each cell according to the depth. We simulate the fracture propagation and shut-in of one stage. In the production simulation, a reservoir simulation box with a size of 150 × 400 × 280 m 3 is selected within the computation domain and cells outside the box are set as null blocks, which is a representation of the control reservoir volume of this particular stage. In the simulation, we first simulate fracture propagation for a one-hour slurry injection schedule. Three separate cases are run using slick water, 0.7 quality NG foam 30 ppt J580, and hybrid NG (first 10 min), and 0.7 quality NG foam 30 ppt J580 (last 50 min). The viscous properties of the three fracturing fluids are shown in Figure 2 and Table 1. The slurry pumping rate is set to 0.15 m 3 /s. Proppant laden fluid is injected starting from 10 min. A 100-mesh proppant is injected from 10 min to 30 min and a 30/50-mesh proppant is injected from 30 min to 60 min. The proppant volume fraction ramps up gradually every 10 min and reaches 20% in the last 10 min (Figure 8). The composition of natural gas used in the simulation is assumed to be 70% C1, 20% C2, 5% C3, and 5% C4. After one hour of In the simulation, we first simulate fracture propagation for a one-hour slurry injection schedule. Three separate cases are run using slick water, 0.7 quality NG foam 30 ppt J580, and hybrid NG (first 10 min), and 0.7 quality NG foam 30 ppt J580 (last 50 min). The viscous properties of the three fracturing fluids are shown in Figure 2 and Table 1. The slurry pumping rate is set to 0.15 m 3 /s. Proppant laden fluid is injected starting from 10 min. A 100-mesh proppant is injected from 10 min to 30 min and a 30/50-mesh proppant is injected from 30 min to 60 min. The proppant volume fraction ramps up gradually every 10 min and reaches 20% in the last 10 min (Figure 8). The composition of natural gas used in the simulation is assumed to be 70% C 1 , 20% C 2 , 5% C 3 , and 5% C 4 . After one hour of injection, the well is shut-in for 5 days. During the shut-in time, the fractures are allowed to grow further (tip extension) and close. Fracturing fluid (water and gas) continues to leak-off into the reservoir and proppant settles inside the fractures. After the shut-in period, the well is produced for 5 years with a decreasing downhole pressure from 32.8 MPa (initial pore pressure) to 10 MPa for the duration of the production time. Cumulative surface oil, gas, and water production rates are computed. injection, the well is shut-in for 5 days. During the shut-in time, the fractures are allowed to grow further (tip extension) and close. Fracturing fluid (water and gas) continues to leak-off into the reservoir and proppant settles inside the fractures. After the shut-in period, the well is produced for 5 years with a decreasing downhole pressure from 32.8 MPa (initial pore pressure) to 10 MPa for the duration of the production time. Cumulative surface oil, gas, and water production rates are computed.

Results and Discussion
Simulation results during the injection, shut-in, and production periods are compared and analyzed in this section to demonstrate the differences in the performance of the different fracturing fluids. Figure 9 shows the fracture width distribution at the end of the one-hour slurry injection period. One trend seen in all three sets of simulation is the non-uniform growth of the hydraulic fractures. Because of the stress shadow between all the growing fractures (which has been considered implicitly in our model by solving the solid deformation equation), adjacent fractures propagate in a way to achieve the least overlapping area. It is also observed that the fracture width varies more in the slick water case than in the two NG foam cases. By comparing the two NG foam cases, one can find that the case using hybrid NG and NG foam gives a smaller fracture surface area than the case using pure NG foam, as the injection of low-viscous pure NG causes larger leakoff and less fracture propagation. Since we do not know the real stress profile in the Wolfcamp formation, a single minimum horizontal stress value (41.34 MPa) is applied to all the layers. The fracture geometry may be different if a different stress profile is put into the simulator.

Results and Discussion
Simulation results during the injection, shut-in, and production periods are compared and analyzed in this section to demonstrate the differences in the performance of the different fracturing fluids. Figure 9 shows the fracture width distribution at the end of the one-hour slurry injection period. One trend seen in all three sets of simulation is the non-uniform growth of the hydraulic fractures. Because of the stress shadow between all the growing fractures (which has been considered implicitly in our model by solving the solid deformation equation), adjacent fractures propagate in a way to achieve the least overlapping area. It is also observed that the fracture width varies more in the slick water case than in the two NG foam cases. By comparing the two NG foam cases, one can find that the case using hybrid NG and NG foam gives a smaller fracture surface area than the case using pure NG foam, as the injection of low-viscous pure NG causes larger leakoff and less fracture propagation. Since we do not know the real stress profile in the Wolfcamp formation, a single minimum horizontal stress value (41.34 MPa) is applied to all the layers. The fracture geometry may be different if a different stress profile is put into the simulator. Figure 10 shows the proppant mass per unit area plots at the end of the one-hour treatment in the three cases. It is found that hydraulic fracturing using foam (Figure 10b,c) leads to better propped hydraulic fractures than slick water (Figure 10a). The created fracture surface and the propped fracture surface area are also calculated as shown in Table 6. We define the propped fracture to be the part of the fracture with an average of at least one layer of proppant. It can be clearly seen that slick water is able to create a larger fracture surface area than the foam cases but the propped surface area in the slick water case is smaller than the foam cases. The propped area in the foam cases is around 10% larger than the slick water case. Two reasons can explain this fact. The first reason is that foam typically has a higher viscosity than slick water, even under high temperature and high shear rate conditions inside the fractures. Thus, it has less leakoff and creates fractures with larger apertures, which makes it easier for proppant to flow and reduces the pinch point effects. The second reason is that foam with higher viscosity has better sand carrying capability than lower viscosity slick water, especially for complex height growth scenarios.   Figure 10 shows the proppant mass per unit area plots at the end of the one-hour treatment in the three cases. It is found that hydraulic fracturing using foam (Figure 10b,c) leads to better propped hydraulic fractures than slick water (Figure 10a). The created fracture surface and the propped fracture surface area are also calculated as shown in Table  6. We define the propped fracture to be the part of the fracture with an average of at least one layer of proppant. It can be clearly seen that slick water is able to create a larger fracture surface area than the foam cases but the propped surface area in the slick water case  The shear rate for the fracture flow at the end of the injection period is shown in Figure 11. We notice that the shear rate is within the experimental measurement range reported by [7]. The shear rate for the fracture flow at the end of the injection period is shown in Figure 11. We notice that the shear rate is within the experimental measurement range reported by [7]. Other properties for fractures from each cluster can also be visualized for detailed analysis for the three simulation cases. In Appendix A, we include the fluid efficiency, fracture half-length, fracture height, created fracture surface area, created fracture volume, cumulative proppant mass, and propped fracture surface area for each cluster during hydraulic fracturing. Figure 12 contains four sub-plots showing the individual properties changing with treating time for the three cases. Figure 12a shows the bottom hole pressure versus treating time. Slick water shows the highest breakdown pressure whereas cases using foam have a lower breakdown pressure. This is mainly because of the high compressibility of the NG foam (gas phase in foam has high compressibility). For the hybrid case, one can observe a big bump after 10 min of injection when NG foam is injected into the wellbore following the previous injected NG and the friction suddenly increases. Figure 12b shows the created fracture surface area versus time. As expected, the slick water case gives the largest created fracture surface area while the hybrid case gives the lowest fracture surface area, but the difference is not significant (around 6 percent). This is because slick water has a lower viscosity and lower compressibility and tends to create longer fractures. Foam has a higher viscosity and higher compressibility and tends to create shorter fractures. Natural gas has the lowest viscosity and highest compressibility and tends to create the shortest fractures. Figure 12c shows the average fracture width versus treating time. As expected, foam gives the largest average fracture width and gas gives the smallest fracture width. Figure 12d plots the calculated standard deviation of fracture surface area for each cluster. A larger standard deviation means more non-uniform fracture growth and a smaller standard deviation implies that the fractures from each cluster are growing more uniformly. It is found that foam cases yield more non-uniform fracture growth compared to the slick-water cases, which may cause potential parent-child well interaction problems in a newly developed field [36,41]. Other properties for fractures from each cluster can also be visualized for detailed analysis for the three simulation cases. In Appendix A, we include the fluid efficiency, fracture half-length, fracture height, created fracture surface area, created fracture volume, cumulative proppant mass, and propped fracture surface area for each cluster during hydraulic fracturing. Figure 12 contains four sub-plots showing the individual properties changing with treating time for the three cases. Figure 12a shows the bottom hole pressure versus treating time. Slick water shows the highest breakdown pressure whereas cases using foam have a lower breakdown pressure. This is mainly because of the high compressibility of the NG foam (gas phase in foam has high compressibility). For the hybrid case, one can observe a big bump after 10 min of injection when NG foam is injected into the wellbore following the previous injected NG and the friction suddenly increases. Figure 12b shows the created fracture surface area versus time. As expected, the slick water case gives the largest created fracture surface area while the hybrid case gives the lowest fracture surface area, but the difference is not significant (around 6 percent). This is because slick water has a lower viscosity and lower compressibility and tends to create longer fractures. Foam has a higher viscosity and higher compressibility and tends to create shorter fractures. Natural gas has the lowest viscosity and highest compressibility and tends to create the shortest fractures. Figure 12c shows the average fracture width versus treating time. As expected, foam gives the largest average fracture width and gas gives the smallest fracture width. Figure 12d plots the calculated standard deviation of fracture surface area for each cluster. A larger standard deviation means more non-uniform fracture growth and a smaller standard deviation implies that the fractures from each cluster are growing more uniformly. It is found that foam cases yield more non-uniform fracture growth compared to the slick-water cases, which may cause potential parent-child well interaction problems in a newly developed field [36,41].

Results for the Injection Period
From the simulation results, we can calculate the total amount of fracturing fluid that is used in the treatment. The total mass of water in this stage for the slick water case is 486 tons. The total mass of the natural gas and water in this stage for the NG foam case is 112.68 and 145 tons, respectively (using the pumping schedule shown in Figure 8). Supposing this well has 30 stages, then in the slick water case, 14,580 tons of water is used. In the NG foam case, the total mass of natural gas and water used in the fracturing job are 3380.4 and 4350 tons. By comparison, 10,230 tons of water can be saved using the 0.7 quality NG foam as the fracturing fluid instead of water. From the simulation results, we can calculate the total amount of fracturing fluid that is used in the treatment. The total mass of water in this stage for the slick water case is 486 tons. The total mass of the natural gas and water in this stage for the NG foam case is 112.68 and 145 tons, respectively (using the pumping schedule shown in Figure 8). Supposing this well has 30 stages, then in the slick water case, 14,580 tons of water is used. In the NG foam case, the total mass of natural gas and water used in the fracturing job are 3380.4 and 4350 tons. By comparison, 10,230 tons of water can be saved using the 0.7 quality NG foam as the fracturing fluid instead of water.

Results for the Shut-In Period
The well is then shut-in for 5 days after the one-hour slurry injection. At an early time, immediately after shut-in, water hammer effects are visible during which there is complex slurry redistribution among clusters and wellbore. As Figure 13 illustrates, during the shut-in time, some fractures obtain more fluids and fracture tip extension happens in some of the clusters and the fractures from these clusters grow further. Slurry comes out of the other fractures and flows back into the wellbore because of fracture closure caused by the complex stress shadow effects from the adjacent clusters. Figure 13 also indicates that the water hammer effect mainly happens at the first several hours of injection and gradually dies away. Moreover, Figure 13 demonstrates the strong capability of our simulator in lifecycle modeling.

Results for the Shut-In Period
The well is then shut-in for 5 days after the one-hour slurry injection. At an early time, immediately after shut-in, water hammer effects are visible during which there is complex slurry redistribution among clusters and wellbore. As Figure 13 illustrates, during the shut-in time, some fractures obtain more fluids and fracture tip extension happens in some of the clusters and the fractures from these clusters grow further. Slurry comes out of the other fractures and flows back into the wellbore because of fracture closure caused by the complex stress shadow effects from the adjacent clusters. Figure 13 also indicates that the water hammer effect mainly happens at the first several hours of injection and gradually dies away. Moreover, Figure 13 demonstrates the strong capability of our simulator in lifecycle modeling.
Fracture closure and fluid leakoff happen throughout the shut-in time. Figure 14 shows the pore pressure changes on the horizontal plane that crosses the wellbore after the shut-in time for the three cases. Comparing Figure 14a,b, one can observe that the slick water case has more fluid (water) leaking into the reservoir than fluid (NG and water) leaking into the reservoir in the NG foam case because of the high viscosity of foam. However, more water leak-off into the reservoir is potentially harmful to the reservoir and may cause severe damages including clay swelling, relative permeability reduction, and water blocking which can be extremely harmful to gas wells. The leaked water in the reservoir is also hard to flow back because of micro-crack closure and capillary pressure effects. On the contrary, the use of foam can minimize the water invasion into the reservoir and more natural gas, which is miscible with the reservoir hydrocarbon fluid can flow into the reservoir and mitigate the impact of water leak-off. A large pore pressure increase can be seen in the near-wellbore region in Figure 14c for the hybrid case, which is the result of pure NG injection. NG flow into the reservoir can potentially benefit oil production by lowering the hydrocarbon viscosity, increasing oil relative permeability, and building up pore pressure. Fracture closure and fluid leakoff happen throughout the shut-in time. Figure 14 shows the pore pressure changes on the horizontal plane that crosses the wellbore after the shut-in time for the three cases. Comparing Figures 14a,b, one can observe that the fects. On the contrary, the use of foam can minimize the water invasion into the reservoir and more natural gas, which is miscible with the reservoir hydrocarbon fluid can flow into the reservoir and mitigate the impact of water leak-off. A large pore pressure increase can be seen in the near-wellbore region in Figure 14c for the hybrid case, which is the result of pure NG injection. NG flow into the reservoir can potentially benefit oil production by lowering the hydrocarbon viscosity, increasing oil relative permeability, and building up pore pressure. We plot the cumulative fluid efficiency as a function of time during the hydraulic fracturing time period (first 1 h) and the shut-in period (from 1 h to 5 days) for the three fracturing fluid designs ( Figure 15). We notice that the fluid efficiency of the slick water and NG foam is very high during the hydraulic fracturing period, which indicates that the leakoff is very low for both cases. For the hybrid NG and NG foam case, during the NG injection (first 10 min), the fluid efficiency is very low, which is due to the fact that NG viscosity is very low, leading to larger leakoff. During the NG foam injection (last 50 min), the fluid efficiency increases due to the increasing fluid viscosity and decreasing leakoff by the NG foam. During the shut-in period, the fluid efficiency decreases in all three cases because no more slurry is injected into the wellbore and more fluid leaks off into the reservoir during the 5-day shut-in. It can also be observed that 90% of slick water leaks into the reservoir whereas only 50% of the NG foam leaks into the reservoir because of the difference in the fluid viscosity and reservoir permeability enhancement. This indicates the advantage of easier flowback for the NG foam fracturing fluids. We plot the cumulative fluid efficiency as a function of time during the hydraulic fracturing time period (first 1 h) and the shut-in period (from 1 h to 5 days) for the three fracturing fluid designs ( Figure 15). We notice that the fluid efficiency of the slick water and NG foam is very high during the hydraulic fracturing period, which indicates that the leakoff is very low for both cases. For the hybrid NG and NG foam case, during the NG injection (first 10 min), the fluid efficiency is very low, which is due to the fact that NG viscosity is very low, leading to larger leakoff. During the NG foam injection (last 50 min), the fluid efficiency increases due to the increasing fluid viscosity and decreasing leakoff by the NG foam. During the shut-in period, the fluid efficiency decreases in all three cases because no more slurry is injected into the wellbore and more fluid leaks off into the reservoir during the 5-day shut-in. It can also be observed that 90% of slick water leaks into the reservoir whereas only 50% of the NG foam leaks into the reservoir because of the difference in the fluid viscosity and reservoir permeability enhancement. This indicates the advantage of easier flowback for the NG foam fracturing fluids.

Results for the Production Period
After the shut-in period, the well is produced for 5 years. The downhole pressure is decreased from 32.8 MPa to 10 MPa in the first 100 days and is fixed at 10 MPa in the following production period. Figure 16 shows the pore pressure distribution within the horizontal plane across the wellbore at the end of the production period. Only the region within the stage of focus is being depleted because of the fact that the grid blocks outside the stage of focus are set to be null blocks. One can observe that the pore pressure in the matrix is larger in the slick water case than in the two foam cases. This is because of the

Results for the Production Period
After the shut-in period, the well is produced for 5 years. The downhole pressure is decreased from 32.8 MPa to 10 MPa in the first 100 days and is fixed at 10 MPa in the following production period. Figure 16 shows the pore pressure distribution within the horizontal plane across the wellbore at the end of the production period. Only the region within the stage of focus is being depleted because of the fact that the grid blocks outside the stage of focus are set to be null blocks. One can observe that the pore pressure in the matrix is larger in the slick water case than in the two foam cases. This is because of the low residual fracture width and poor sand placement near the wellbore regions after the shut-in as indicated by Figures 9 and 10. Poor proppant placement in some clusters near the wellbore leads to low propped fracture conductivity. The unpropped fracture part will close and form pinch points after the shut-in, which prevents the reservoir fluid from flowing through those fractures into the wellbore. The fractures from these clusters may contribute less to the production of the full well. This results in low cluster efficiency and lower well productivity. Thus, the fractured reservoir volumes with unpropped fractures will produce less and the pore pressure will remain high. The fractured reservoir volumes with well-propped fractures will produce more and the pore pressure will be low. Comparing the three depleted reservoir pressure profiles in Figure 16, we found that the reservoir is better depleted in the NG foam case and the hybrid NG and NG foam case. This is because of the better sand placement in the two foam cases as shown in Figure 10 and Table 6. Some final important caveats to the application of our results: First, since this is a simulation study and many assumptions have been made in the simulation, the field implementation of this strategy may yield less favorable results when using the NG foam. More precise data input, model calibration with existing field data, and field pilot tests are still essential to obtain more predictive results of NG foam fracturing. Second, safety concerns associated with the use of NG and NG foam may prevent the use of such fluids Figure 16. Pore pressure distribution on the horizontal plane through the wellbore at the end of the production period. (a) Case using slick water, (b) Case using NG foam, (c) Case using hybrid NG and NG foam.
The cumulative oil and gas production versus production time at surface conditions is listed in Table 5 and plotted in Figure 17. We observe that the cases using NG foam or hybrid NG and NG foam give much higher production than the case using slick water. An increase of 95% and 65% in oil production is observed from the simulation. A similar trend is also observed in cumulative gas production. As discussed in the previous parts of this paper, several reasons are contributing to the production of foam fractured wells surpassing the slick water fractured well. (1). More microcracks are created (reflected by a larger change in the effective stress) using NG foam and NG, and thus higher SRV permeability enhancement. (2). Created fractures are better propped using NG foam, and thus higher cluster efficiency. (3). NG enters the reservoir and reduces the amount of water entering the reservoir, benefiting the oil production by relative permeability improvement, viscosity reduction, and energizing the reservoir. Some final important caveats to the application of our results: First, since this is a simulation study and many assumptions have been made in the simulation, the field implementation of this strategy may yield less favorable results when using the NG foam. More precise data input, model calibration with existing field data, and field pilot tests are still essential to obtain more predictive results of NG foam fracturing. Second, safety concerns associated with the use of NG and NG foam may prevent the use of such fluids in many instances in the field. These concerns can also potentially increase the cost of pumping such fluids (even when the gas is available for free). These safety concerns are offset by the benefits of a substantial reduction in flaring of produced natural gas.

Conclusions
An equation of state compositional hydraulic fracturing and reservoir model has been extended to include the modeling capability of non-Newtonian fracturing fluid for hydraulic fracturing in this paper. We built a precise layer model and fluid model using published data from the Wolfcamp formations in the Permian Basin and performed lifecycle analysis to the reservoir. We simulated hydraulic fracture propagation using slick water, NG foam, and hybrid NG and NG foam, shut-in and fracture closure, and long-term Some final important caveats to the application of our results: First, since this is a simulation study and many assumptions have been made in the simulation, the field implementation of this strategy may yield less favorable results when using the NG foam. More precise data input, model calibration with existing field data, and field pilot tests are still essential to obtain more predictive results of NG foam fracturing. Second, safety concerns associated with the use of NG and NG foam may prevent the use of such fluids in many instances in the field. These concerns can also potentially increase the cost of pumping such fluids (even when the gas is available for free). These safety concerns are offset by the benefits of a substantial reduction in flaring of produced natural gas.

Conclusions
An equation of state compositional hydraulic fracturing and reservoir model has been extended to include the modeling capability of non-Newtonian fracturing fluid for hydraulic fracturing in this paper. We built a precise layer model and fluid model using published data from the Wolfcamp formations in the Permian Basin and performed lifecycle analysis to the reservoir. We simulated hydraulic fracture propagation using slick water, NG foam, and hybrid NG and NG foam, shut-in and fracture closure, and long-term production. Results from different periods of the well lifecycle are visualized and analyzed.
The key findings are summarized below: 1.
The use of NG foam reduces both gas flaring and water usage in well completions potentially saving both capital and reducing the environmental impact. The simulation results presented here can be used to conduct an economic analysis that includes the drilling cost, completion cost, additional cost for gas storage and foam generation, oil/gas/water daily production, oil price to evaluate the cash flow performance of different fracturing fluids and completion designs.

2.
The use of NG foam gives lower breakdown pressure, smaller created fracture surface area, and larger average fracture width during the hydraulic fracturing than slick water. 3.
NG foam causes higher effective stress (more compressive) and a larger stress shadow effect than slick water during hydraulic fracturing. More microcracks and higher SRV permeability are expected using NG foam. 4.
NG foam has a significantly better sand transport capacity and can thus achieve better proppant placement inside the fractures.

5.
Higher cluster efficiency and higher oil and gas productivity are observed using NG foam than slick water. 6.
NG foam leads to more non-uniform fracture growth which may leave unstimulated regions in a newly developed reservoir. 7.
Precise data input, model calibration, and field pilot tests are essential to obtain more accurately predictable results for NG foam fracturing.

Appendix A
In this Appendix A, the properties of each cluster during the fracturing simulation are plotted for the three simulation cases, respectively. Figure A1 shows the fluid efficiency of each cluster. Figure A2 shows the fracture half-length of each cluster. Figure A3 shows the fracture height of each cluster. Figure A4 shows the created fracture surface area of each cluster. Figure A5 shows the fracture volume of each cluster. Figure A6 shows the cumulative proppant mass into each cluster. Figure A7 shows the propped fracture surface area of each cluster.