Numerical Investigation of Fracture Compressibility and Uncertainty on Water-Loss and Production Performance in Tight Oil Reservoirs

Multi-stage hydraulic fracturing along with horizontal wells are widely used to create complex fracture networks in tight oil reservoirs. Analysis of field flowback data shows that most of the fracturing fluids are contained in a complex fracture network, and fracture-closure is the main driving mechanism during early clean up. At present, the related fracture parameters cannot be accurately obtained, so it is necessary to study the impacts of fracture compressibility and uncertainty on water-loss and the subsequent production performance. A series of mechanistic models are established by considering stress-dependent porosity and permeability. The impacts of fracture uncertainties, such as natural fracture density, proppant distribution, and natural fracture heterogeneity on flowback and productivity are quantitatively assessed. Results indicate that considering fracture closure during flowback can promote water imbibition into the matrix and delay the oil breakthrough time compared with ignoring fracture closure. With the increase of natural fracture density, oil breakthrough time is advanced, and more water is retained underground. When natural fractures connected with hydraulic fractures are propped, well productivity will be enhanced, but proppant embedment can cause a loss of oil production. Additionally, the fracture network with more heterogeneity will lead to the lower flowback rate, which presents an insight in the role of fractures in water-loss.


Introduction
The development of conventional petroleum reservoirs has been unable to keep pace with the increasing global energy demand, which has shifted industrial attention to low permeability or tight reservoirs [1].Rapid development of horizontal wells and multi-stage hydraulic fracturing technology enables unconventional tight reservoirs to be exploited [2,3].In the process of stimulation, hydraulic fractures may connect and reactivate the already existing natural fractures or generate new induced fractures near the wellbore.Well productivity could be substantially increased if complex fracture networks are created by hydraulic fracturing [4,5].
Different from conventional fracturing, slickwater is the most commonly used fracturing fluid in unconventional reservoirs.Slickwater has a relatively low cost and low viscosity, which may help to create complex fracture networks [5][6][7].Horizontal wells that undergo multi-stage fracturing operations often require a large amount of water injected into the formation to create a large stimulated reservoir volume.After hydraulic fracturing is completed, the flowback treatment of fracturing fluid should be carried out, followed by long-term production [8][9][10].However, field data show that the Energies 2019, 12, 1189 2 of 20 recovery efficiency of tight oil or gas wells is generally low, and a remarkable fraction of fracturing fluid remains in the reservoir even after long-term production [11][12][13][14].Scholars have carried out some work to clarify the problems arising after hydraulic fracturing and to help improve the clean-up operation in low permeability formations [15,16].
If hydraulic fracturing fluid is trapped in the formation, then it must be in the rock matrix or fracture systems [17].An analytical model was presented for spontaneous imbibition of capillary force, which explained the fracturing fluid was imbibed into the shale matrix by strong capillary forces [18].However, Dutta et al. [19] showed that water-loss rate of low-permeability sands was far lower than that of high-permeability sands.Although the tight sands had strong capillary force, the low permeability of the formation also limited water imbibition.Similarly, experimental data showed that the water imbibition capacity of tight shale was very low, even for hydrophilic rock samples [20].Furthermore, it was difficult to explain the reason for low water recovery efficiency when the spontaneous capillary forces were absent in oil-wet tight rocks [21,22] or the wells were immediately cleaned up without shut-in.By analyzing the pressure transient data, natural fractures that re-opened during fracturing were apparently closed in subsequent production processes [23].In the process of fracture reopening and closing, water entering natural fractures continues to be imbibed into matrix, and only part of the water can flow back to the wellbore in the end.McClure [24] focused on the impact of fracture network complexity on water recovery efficiency.The results of his simulations indicated that the closure of unpropped fractures within fracture networks resulted in low water recovery.Fu et al. [25] collated and processed the early flowback data of seven multi-fractured horizontal wells in tight reservoirs.Data analysis suggested that most of the facture systems were un-propped after stimulation and held most of the fracturing water.
Whether the fracturing fluid remains in the matrix or fractures, it will affect the productivity of tight oil wells.Spontaneous imbibition of fracturing water into the matrix was considered as a possible mechanism to improve oil production [26,27].Extensive experimental and mathematical studies were devoted to explaining the oil displacement efficiency by water imbibition with various rock wettabilities, physical parameters, pore size distribution, and fracture characteristics [28][29][30].The results showed that water imbibition could drive oil out of the matrix pore, and the displacement effect would be better when the tight rocks contained fractures [29].However, Wang and Leung [31,32] established a series of numerical models to investigate the water-loss mechanisms during flowback operations.They reported that although prolonging shut-in time could enhance water imbibition into the matrix and increase initial oil rate, there was no benefit for long-term production.Low recovery efficiency means that a large amount of water remains in secondary fractures, which reduced oil mobility and had a negative impact on long-term production.Similar negative impacts were also found in studies of tight shale [33][34][35], where water retention might cause water blocking or damage to hydrocarbon phase relative permeability.Moreover, for tight formations, the invasion of fracturing fluid into the matrix did not require much depth to cause enough damage [15,16].
Although the mechanism of fracturing fluid retention underground and its impact on well productivity have been widely studied, there are still many problems to be further explored.Firstly, existing models failed to accurately characterize the distribution of slickwater in matrix and fracture systems during fracturing.These works either directly assign high water saturation value to the fractures or simulate the injection process of fracturing fluid through stress-dependent permeability models [36,37].Although the latter is slightly more reasonable than the former, its mechanism is unrealistic and overestimates the water absorption capacity of the matrix near fractures.The decrease of matrix permeability might be due to compaction, but it was hardly to reasonably increase matrix permeability [23].Secondly, after stimulation, secondary fractures and re-opened natural fractures held most of the fracturing water [38], and fracture depletion was observed during early flowback [39].The driving mechanism of fracture-closure is rarely considered in simulation and its impact on water-loss is unclear.Thirdly, the parameters of complex fracture networks are the key factors affecting the distribution of fracturing fluids in formations and subsequent productivity.Numerical Energies 2019, 12, 1189 3 of 20 studies mainly focus on the density and distribution of natural fractures [37,40], while the proppant distribution in complex fracture systems is seldom considered.The stress-dependent porosity and permeability correlations vary because of proppant concentrations, and the effect of proppant embedment on long-term production cannot be neglected [41,42].Finally, the impact of natural fracture heterogeneity on water retention has not been reported.
This paper reports on simulations to better understand water loss and production from complex fracture networks.A triple-porosity model is established by using a numerical reservoir simulator, where matrix and fracture systems are explicitly discretized, so that we can examine imbibition hysteresis and stress-dependent porosity and permeability.The adaptability of a mechanistic model is validated, and the influence of fracture-closure on fracturing water retention is investigated.Next, we quantitatively analyze the effect of uncertain fracture parameters, such as natural fracture density, proppant distribution, and natural fracture heterogeneity.The results of this work can provide a better understanding of the impact of fracture compressibility and uncertainty on water-loss and production performance in tight oil reservoirs.

Methodology
A series of numerical models are constructed by fully exploiting the functions of commercial reservoir simulation software (CMG IMEX).The whole reservoir consists of triple media: matrix (M), hydraulic fracture (HF), and natural fracture (NF).The model is essentially a single-porosity medium, while a fracture network is set up by adopting logarithmic local grid refinement.This approach solves the problem of huge differences between fracture width and matrix grid size, and ensures the stability of numerical calculation [43].Hydraulic fracture stages along horizontal wellbores are assumed to be uniformly spaced and symmetrical.Based on this symmetrical structure, only one hydraulic fracture segment is selected for simulation in this work, and its results can be mapped back to the whole horizontal well scale [44].The grid structure of a segment is shown in Figure 1, where the horizontal wellbore is along the x-direction and the hydraulic fracture is along the y-direction.proppant distribution in complex fracture systems is seldom considered.The stress-dependent porosity and permeability correlations vary because of proppant concentrations, and the effect of proppant embedment on long-term production cannot be neglected [41,42].Finally, the impact of natural fracture heterogeneity on water retention has not been reported.This paper reports on simulations to better understand water loss and production from complex fracture networks.A triple-porosity model is established by using a numerical reservoir simulator, where matrix and fracture systems are explicitly discretized, so that we can examine imbibition hysteresis and stress-dependent porosity and permeability.The adaptability of a mechanistic model is validated, and the influence of fracture-closure on fracturing water retention is investigated.Next, we quantitatively analyze the effect of uncertain fracture parameters, such as natural fracture density, proppant distribution, and natural fracture heterogeneity.The results of this work can provide a better understanding of the impact of fracture compressibility and uncertainty on water-loss and production performance in tight oil reservoirs.

Methodology
A series of numerical models are constructed by fully exploiting the functions of commercial reservoir simulation software (CMG IMEX).The whole reservoir consists of triple media: matrix (M), hydraulic fracture (HF), and natural fracture (NF).The model is essentially a single-porosity medium, while a fracture network is set up by adopting logarithmic local grid refinement.This approach solves the problem of huge differences between fracture width and matrix grid size, and ensures the stability of numerical calculation [43].Hydraulic fracture stages along horizontal wellbores are assumed to be uniformly spaced and symmetrical.Based on this symmetrical structure, only one hydraulic fracture segment is selected for simulation in this work, and its results can be mapped back to the whole horizontal well scale [44].The grid structure of a segment is shown in Figure 1, where the horizontal wellbore is along the x-direction and the hydraulic fracture is along the y-direction.

Model Description
For base case, a model with a size of 70 m × 300 m is established as shown in Figure 1, in which 70 m is the hydraulic fracture spacing and 300 m is the horizontal well spacing.Relevant known model parameters of Lucaogou Formation in Junggar Basin, obtained from field and literature reports [45,46], are summarized in Table 1.The parameters of complex fracture network formed after multi-stage fracturing along horizontal wells are unknown, so the typical values of other tight reservoirs [31,40] are assigned here as shown in Table 2. Hydraulic fractures are considered to be propped, while reactivated natural fractures with orthogonal distribution are considered to be un-propped [21].Since the wellbore flow pressure is always higher than the bubble-point pressure during production, the oil-water two-phase flow can be considered in the simulation process.The

Model Description
For base case, a model with a size of 70 m × 300 m is established as shown in Figure 1, in which 70 m is the hydraulic fracture spacing and 300 m is the horizontal well spacing.Relevant known model parameters of Lucaogou Formation in Junggar Basin, obtained from field and literature reports [45,46], are summarized in Table 1.The parameters of complex fracture network formed after multi-stage fracturing along horizontal wells are unknown, so the typical values of other tight reservoirs [31,40] are assigned here as shown in Table 2. Hydraulic fractures are considered to be propped, while reactivated natural fractures with orthogonal distribution are considered to be un-propped [21].Since the wellbore flow pressure is always higher than the bubble-point pressure during production, the oil-water two-phase flow can be considered in the simulation process.The two-phase relative permeability curves of matrix, or hydraulic fracture and natural fracture are assigned as presented in Figure 2a,b, respectively [37].Additionally, the clay content in the Lucaogou Formation is low [47], therefore, the clay swelling has not been considered in this work.two-phase relative permeability curves of matrix, or hydraulic fracture and natural fracture are assigned as presented in Figure 2a,b, respectively [37].Additionally, the clay content in the Lucaogou Formation is low [47], therefore, the clay swelling has not been considered in this work.

Capillary Pressure Curve
The impact of imbibition hysteresis is considered in this study, so the capillary pressure curves are defined as functions of normalized water saturation (Swn) according to Equations The normalized water saturation is defined as:

Capillary Pressure Curve
The impact of imbibition hysteresis is considered in this study, so the capillary pressure curves are defined as functions of normalized water saturation (S wn ) according to Equations (1)-(3) [36]: (1) The normalized water saturation is defined as: where P cD is drainage capillary pressure; MPa; P cI is imbibition capillary pressure; MPa; S w is water saturation; S wc is connate water saturation; and a 1 , a 2 , and a 3 are assumed to be equal to b 1 , b 2 , and b 3 , respectively.In shale gas reservoirs, the rocks are usually more strongly water-wet.Therefore, only the positive part of the capillary force curve was considered in previous studies, by assuming both b 4 and b 5 as equal to zero.However, in tight oil reservoirs, the rocks are not all hydrophilic, but also show as intermediately wet or oil-wet.Thus, in this work, the negative capillary pressures are not neglected and hysteresis effect in the matrix is considered.The capillary pressure curve of the matrix is shown in Figure 3a, in which b 1 , b 2 , b 3 , b 4 , and b 5 are assigned to 0, 4.883, 20, 86, and 5, respectively [37].
Energies 2019, 12, x FOR PEER REVIEW 5 of 22 where PcD is drainage capillary pressure; MPa; PcI is imbibition capillary pressure; MPa; Sw is water saturation; Swc is connate water saturation; and a1, a2, and a3 are assumed to be equal to b1, b2, and b3, respectively.In shale gas reservoirs, the rocks are usually more strongly water-wet.Therefore, only the positive part of the capillary force curve was considered in previous studies, by assuming both b4 and b5 as equal to zero.However, in tight oil reservoirs, the rocks are not all hydrophilic, but also show as intermediately wet or oil-wet.Thus, in this work, the negative capillary pressures are not neglected and hysteresis effect in the matrix is considered.The capillary pressure curve of the matrix is shown in Figure 3a, in which b1, b2, b3, b4, and b5 are assigned to 0, 4.883, 20, 86, and 5, respectively [37].For natural fractures, capillary pressure should not be zero, depending on wall roughness, fracture width and pore structure [48].Here the J-function correlation in Equation ( 4) [49] is used to calculate capillary pressure for natural fractures, as shown in Figure 3b.Hydraulic fractures assume that capillary pressure is neglected because of high conductivity: In Equation ( 4), Pc is capillary pressure in NF, KPa; σ is interfacial tension and is equal to be 72 dynes/cm; φ is porosity; k is absolutely permeability, md; and a1 ' , a2 ' , and a3 ' are equal to 1.86, 6.42, and 0.5, respectively.

Stress-Dependent Porosity and Permeability
During a hydraulic fracturing operation, a large amount of slickwater is pumped into a formation in a short time under high wellhead pressure, thereafter hydraulic fractures are created.Due to the limited fluid absorption capacity of the tight matrix, most of the high-pressure fluid leaks into natural fractures.With the increase of fluid pressure in natural fractures, the net closure stress within fractures decreases and the natural fracture width gradually increases until they are fully opened [50].During the well shut-in, clean-up, and production periods, with the fracturing fluid imbibing into the matrix or being recovered to the surface, the net closure stress within fracture networks increases and the fracture gradually closes.Additionally, the reduction of hydraulic fracture width is effectively restrained because of proppant support, while un-propped natural fractures will be closed or even completely locked [23].Previously, the matrix stress-dependent permeability was usually considered to enable rapid injection of so much fracturing fluid into the formation.However, it might not be realistic that the matrix permeability increases substantially For natural fractures, capillary pressure should not be zero, depending on wall roughness, fracture width and pore structure [48].Here the J-function correlation in Equation ( 4) [49] is used to calculate capillary pressure for natural fractures, as shown in Figure 3b.Hydraulic fractures assume that capillary pressure is neglected because of high conductivity: In Equation ( 4), P c is capillary pressure in NF, KPa; σ is interfacial tension and is equal to be 72 dynes/cm; ϕ is porosity; k is absolutely permeability, md; and a 1 ' , a 2 ' , and a 3 ' are equal to 1.86, 6.42, and 0.5, respectively.

Stress-Dependent Porosity and Permeability
During a hydraulic fracturing operation, a large amount of slickwater is pumped into a formation in a short time under high wellhead pressure, thereafter hydraulic fractures are created.Due to the limited fluid absorption capacity of the tight matrix, most of the high-pressure fluid leaks into natural fractures.With the increase of fluid pressure in natural fractures, the net closure stress within fractures decreases and the natural fracture width gradually increases until they are fully opened [50].During the well shut-in, clean-up, and production periods, with the fracturing fluid imbibing into the matrix or being recovered to the surface, the net closure stress within fracture networks increases and the fracture gradually closes.Additionally, the reduction of hydraulic fracture width is effectively restrained because of proppant support, while un-propped natural fractures will be closed or even completely locked [23].Previously, the matrix stress-dependent permeability was usually considered to enable rapid injection of so much fracturing fluid into the formation.However, it might not be realistic that the matrix permeability increases substantially with the increase of net pressure [23].Therefore, in this work, stress dependence of porosity and permeability in fractures is considered by Energies 2019, 12, 1189 6 of 20 fully excavating the software [51,52], as shown in Figure 4. Four paths in the figure correspond to two processes: (1) Elastic path, and (2) plastic/dilation path is used to simulate the increase of porosity and permeability in fracture system during the injection process, (3) unloading path, and (4) recompaction path are used to simulate the decrease of porosity and permeability in fracture system during the flowback or production process.with the increase of net pressure [23].Therefore, in this work, stress dependence of porosity and permeability in fractures is considered by fully excavating the software [51,52], as shown in Figure 4. Four paths in the figure correspond to two processes: 1) Elastic path, and 2) plastic/dilation path is used to simulate the increase of porosity and permeability in fracture system during the injection process, 3) unloading path, and 4) recompaction path are used to simulate the decrease of porosity and permeability in fracture system during the flowback or production process.• During the injection process The following empirical exponential formulas Equations ( 5) and ( 6) are used to describe the variation of porosity and permeability with the pressure in fracture system to realize the actual pumping rate and time [36,53]: In Equations ( 5) and ( 6), φf is current fracture porosity; φf0 is original fracture porosity; Cfi is the fracture compressibility during injection, KPa −1 ; Pnet is net pressure within fractures (difference between original pressure and current pressure), KPa; kf is current fracture permeability, md; kf0 is original fracture permeability, md; mi is the permeability changing factor, KPa −1 (exponent empirically determined).Tiab et al. [54] reported that the fracture compressibility was as much as two orders of magnitude higher than the matrix compressibility.Thus, Cfi is assigned a value of 5.97 × 10 −5 , which is within a reasonable range.The values of mi in HF and NF systems are assigned 6.67 × 10 −5 and 1.00 × 10 −4 [36,37].The calculated results of stress-dependence in fractures are presented in Figure 5.

• During the injection process
The following empirical exponential formulas Equations ( 5) and ( 6) are used to describe the variation of porosity and permeability with the pressure in fracture system to realize the actual pumping rate and time [36,53]: In Equations ( 5) and ( 6), ϕ f is current fracture porosity; ϕ f0 is original fracture porosity; C fi is the fracture compressibility during injection, KPa −1 ; P net is net pressure within fractures (difference between original pressure and current pressure), KPa; k f is current fracture permeability, md; k f0 is original fracture permeability, md; m i is the permeability changing factor, KPa −1 (exponent empirically determined).Tiab et al. [54] reported that the fracture compressibility was as much as two orders of magnitude higher than the matrix compressibility.Thus, C fi is assigned a value of 5.97 × 10 −5 , which is within a reasonable range.The values of m i in HF and NF systems are assigned 6.67 × 10 −5 and 1.00 × 10 −4 [36,37].The calculated results of stress-dependence in fractures are presented in Figure 5.  • During the production process In the process from well shut-in to production, fractures would gradually close, which is different from fracture dilation in the previous section.Many studies have been carried out to clarify the fracture-closure mechanism, among which the main difficulty lies in the determination of  • During the production process In the process from well shut-in to production, fractures would gradually close, which is different from fracture dilation in the previous section.Many studies have been carried out to clarify the fracture-closure mechanism, among which the main difficulty lies in the determination of fracture compressibility [55,56].For un-propped natural fractures, Equation (7) proposed by Jones [57] can be used: C fp = −1 P net ln(P net /P h ) In Equation ( 7), C fp is the fracture compressibility during production, KPa −1 ; P h is an apparent healing pressure and is equal to be 1.38 × 10 5 KPa.However, for propped fractures, the value of C fp can refer to the chart of fracture compressibility inferred by Aguilera [41] as shown in Figure 6.Aguilera [41] emphasizes that the graph is only an approximation, which is not to replace good laboratory work.The correlation is useful only if laboratory data are not available.Figure 6 shows the results of filling mineral ratio from 0% to 50%, which represents the percentage of minerals in fractures.Fu [25] and Williams-Kovacs [56] had proposed that proppants in fractures could be assumed as minerals in Aguilera model.Therefore, the mineral ratio in this work can be the percentage of proppants in fractures, so that the compressibility of propped fractures can be obtained.When the fracture compressibility is determined, the stress-dependent porosity and permeability in the production process can be calculated by Equations ( 8) and (9) [55,58].The calculation results of stress-dependence in fractures are presented in Figure 7a,b

Model Validation
Based on the established basic model, the processes of fracturing fluid injection and flowback are sequentially simulated, in which the imbibition hysteresis and stress-dependent porosity and permeability are considered.The model is verified by checking the bottom-hole pressure curve at the injection stage and the rate transient analysis at the flowback stage.

Bottom-Hole Pressure Curve
The injection of slickwater during hydraulic fracturing is simulated through three different mechanisms.The cumulative injection volume is 720 m 3 and the pumping rate is 6 m 3 /min, which are consistent with the field data reported by Chen et al. [59].The simulation results under different mechanisms are shown in Figure 8. Figure 8a presents the variation of bottom-hole pressure with pumping time.When stress-dependent porosity is considered, the injection can smoothly proceed at an actual field scale.Pressure rises sharply to the threshold at the early stage, then drops to the vicinity of fracture-closure pressure and slowly increases with the extension of pumping time.Although proppant is absent from the simulation, the pressure variation characteristics conform to the fracturing curve in the field [59], considering stress-dependent porosity.By contrast, the other two mechanisms reflect their weaknesses.When stress-dependent permeability and porosity are neglected, it is very difficult to inject water into the formation, for example, it takes nearly a week to inject 720 m 3 water, which may unrealistically overestimate the sweep range of pressure in the reservoir.When stress-dependent permeability is the sole mechanism, then pumping rate can barely reach 1 m 3 /min even though injection capacity is improved.
When stress-dependence is neglected, as shown in Figure 8b, more than half of the injected water leaks into the matrix.However, when stress-dependent porosity and permeability is considered, nearly 80% of water remains in the fracture systems.Fu et al. [24] analyzed flowback data from several wells in the field and predicted that an average of 77.7% of the injected water was in the fracture system.The distribution of fracturing water in our base case is very close to Fu's report, so the water-loss in the matrix should be overestimated when stress-dependent porosity is absent.

Rate Transient Analysis
Next, the flow-regime identification of fracturing fluid during early flowback is carried out for the base case.The water recovery is seldom at a constant rate, nor is it at a constant flowing pressure.The rate-normalized pressure (RNP) and its derivatives (RNP') are used to represent the flowback behavior that would be observed if the well were produced at a constant reference rate.The formulas are as follows [60]:

Rate Transient Analysis
Next, the flow-regime identification of fracturing fluid during early flowback is carried out for the base case.The water recovery is seldom at a constant rate, nor is it at a constant flowing pressure.The rate-normalized pressure (RNP) and its derivatives (RNP') are used to represent the flowback behavior that would be observed if the well were produced at a constant reference rate.The formulas are as follows [60]: In Equations ( 10)-( 12), P i is initial pressure in fractures, KPa; P wf is bottom-hole pressure, KPa; q w is water rate, m 3 /day; t MB represents the equivalent time, day; and Q w is the cumulative water production volume, m 3 .Figure 9 uses a log-log plot of RNP and RNP' vs. tMB during early flowback to show the sequence of flow regimes observed with: (1) Transient linear flow (half slope) within fractures occurs when the flow regime presents water flow from the hydraulic fractures to the perforations; (2) fracture depletion (unit slope): occurs when the pressure response reaches the boundary of the fractures; (3) transient linear flow (half slope) in the matrix occurs when the flow regime begins with oil or gas breakthrough from matrix to fractures.This transition might be very short depending on the fracture conductivity and fracture complexity; and, finally, (4) matrix boundary dominated flow (unit slope) occurs when the density of natural fractures is large and the spacing between hydraulic fractures is small, the pressure response quickly propagates to the matrix boundary.The flow regimes observed in Figure 9 are consistent with field observations and analytical models for flowback analysis [61][62][63].The first two flow regimes are well captured and last for more than two days.

Results and Discussion
Despite the good consistency of the base case results, an ideal historical match with actual flowback data is still difficult to achieve.There are two main reasons: (1) most of the flowback data obtained from the field are daily output, occasionally hourly data, which can hardly reflect the characteristics of early flowback [31]; and (2) properties of fracture systems are highly uncertain including fracture compressibility, density, and conductivity [37].These uncertainties not only make the quantitative validation of simulation results extremely challenging, but also make it difficult to understand the retention mechanism of water and its impact on production.Therefore, a series of models were established to quantify the impact of various factors on water-loss and production performance.The flow regimes observed in Figure 9 are consistent with field observations and analytical models for flowback analysis [61][62][63].The first two flow regimes are well captured and last for more than two days.

Results and Discussion
Despite the good consistency of the base case results, an ideal historical match with actual flowback data is still difficult to achieve.There are two main reasons: (1) most of the flowback data obtained from the field are daily output, occasionally hourly data, which can hardly reflect the characteristics of early flowback [31]; and (2) properties of fracture systems are highly uncertain including fracture compressibility, density, and conductivity [37].These uncertainties not only make the quantitative validation of simulation results extremely challenging, but also make it difficult to understand the retention mechanism of water and its impact on production.Therefore, a series of models were established to quantify the impact of various factors on water-loss and production performance.

Impact of Fracture Closure
Extensive analytical models have been used to reveal the fracture depletion process during early flowback.However, these models have some limitations, such as sequential flow and single-phase flow that make it difficult to explain the effect of fracture-closure on two-phase flow in the triple media after oil breakthrough.The flow region of fracture depletion was validated in the previous section, and the impact of fracture-closure on fracturing water retention is modeled here.Three models are built considering different mechanisms: (1) with stress-dependent porosity and permeability (base case), (2) ignoring stress-dependent porosity, and (3) ignoring stress dependence.The simulation result at the end of injection process in basic model is taken as the initial condition to better compare the impacts of three mechanisms on water-loss during flowback.The simulation process is as follows: shut-in for 22 h, then flowback for 10 days, followed by production for half a year.
Figure 10 shows the variation of bottom-hole pressure with time under different mechanisms during flowback.It can be seen from the figure that when considering stress-dependent porosity, the bottom-hole pressure is almost maintained near the fracture-closure pressure for the first half-day, then the pressure drop is gradually accelerated, even lower than the case without considering stress-dependent during the late flowback stage (from the insert figure in Figure 10).The drive mechanism of fracture-closure can be captured, and this mechanism dominates for nearly five days.The simulation results of production are shown in Table 3.The oil breakthrough time of the base case is the latest, and the single-phase flow duration is one day longer than that of the other two models, which further reflects the fractured storage effect.Moreover, the water recovery efficiency of the basic case is 55.3%, while ignoring stress-dependent porosity overestimates by 20.8% and ignoring stress dependence overestimates by 26.6%.
Energies 2019, 12, x FOR PEER REVIEW 11 of 22 case is the latest, and the single-phase flow duration is one day longer than that of the other two models, which further reflects the fractured storage effect.Moreover, the water recovery efficiency of the basic case is 55.3%, while ignoring stress-dependent porosity overestimates by 20.8% and ignoring stress dependence overestimates by 26.6%.The reasons for these differences are illustrated in Figure 11: (1) the Sw near HF (0.04 m to the fracture surface) in base case remains above 50% during early flowback, and Sw in the deeper matrix also substantially increased.While in the other two models, the Sw near HF has been reduced to 43% and almost no water leaks far away; and (2) the equilibrium point Sw = 0.4 shown in Figure 3a.The Sw near HF is much higher than 0.4 at the end of fracturing, so the hysteresis results in a capillary force that is resistant to water imbibition.With the decrease of water in HF during flowback, if the fracture-closure is ignored, the water near a fracture will flow back into the fracture, which leads to an increase in water recovery efficiency.On the contrary, when fracture compressibility is considered, the reduced water volume in fractures will be equal to the volume of fracture-closure, which provides an opportunity for water near fractures to infiltrate into the deeper matrix.The reasons for these differences are illustrated in Figure 11: (1) the S w near HF (0.04 m to the fracture surface) in base case remains above 50% during early flowback, and Sw in the deeper matrix also substantially increased.While in the other two models, the S w near HF has been reduced to 43% and almost no water leaks far away; and (2) the equilibrium point S w = 0.4 shown in Figure 3a.The S w near HF is much higher than 0.4 at the end of fracturing, so the hysteresis results in a capillary force that is resistant to water imbibition.With the decrease of water in HF during flowback, if the fracture-closure is ignored, the water near a fracture will flow back into the fracture, which leads to an increase in water recovery efficiency.On the contrary, when fracture compressibility is considered, the reduced water volume in fractures will be equal to the volume of fracture-closure, which provides an opportunity for water near fractures to infiltrate into the deeper matrix.From Figure 11, it can be found that the oil phase appears in NF after one day of flowback, when stress-dependent porosity is absent.This is also a good explanation for why the oil breakthrough time in the base case is delayed.By contrast, the cumulative oil production is the lowest for the basic case.Relatively high water saturation reduces the oil permeability because more water leaks into the matrix.The Sw within NF is relatively high (17%) for the base case, which limits the flow of oil in NF.However, the oil production gap between three models is not large because tight reservoirs do not require very high fracture conductivity, and proppant that is evenly distributed within a fracture should be more important for well productivity [64].From Figure 11, it can be found that the oil phase appears in NF after one day of flowback, when stress-dependent porosity is absent.This is also a good explanation for why the oil breakthrough time in the base case is delayed.By contrast, the cumulative oil production is the lowest for the basic case.Relatively high water saturation reduces the oil permeability because more water leaks into the matrix.The S w within NF is relatively high (17%) for the base case, which limits the flow of oil in NF.However, the oil production gap between three models is not large because tight reservoirs do not require very high fracture conductivity, and proppant that is evenly distributed within a fracture should be more important for well productivity [64].

Impact of Natural Fracture Density
Extensive studies of natural fracture density have been carried out [21,24,37], but the relationship between fracture complexity and fracture width has not been considered.Zou et al. [65] reported that the average width of the fracture network decreased as the complexity of a fracture system increased.If this relationship is neglected, the total volume of a complex fracture network should be overestimated with the increase of natural fracture density, which will inevitably interfere with the simulation of water injection and the flowback process.Here, a fracture complexity index (FCI) is introduced to quantitatively evaluate natural fracture density, as shown in Equation ( 13) [66]: where FCI is fracture complexity index; V nf is volume of natural fractures; and V hf is volume of hydraulic fractures.The value of FCI increases with the augment of fracture complexity.
A series of models are set up as follows: (1) considering that the half-length of HF is constant for all cases, increasing the number of natural fractures increases the value of FCI; and (2) by adjusting the width of HF and NF, the initial volume of fracture system in each case is equal, as shown in Table 4, and the width of un-propped NF is assigned less than 2.5 mm [17].The simulation process is as follows: (1) water injection for two hours and well shut-in for 22 h [32], then (2) clean-up at a constant rate for 10 days (in the field the choke size is gradually increasing reported by Clarkson [62], which is simplified here for simulation and analysis), followed by (3) production at constant pressure for one year.As can be seen from Table 4, with the increase of natural fracture density, the contact area between fracture and reservoir increases, which reduces the fluid efficiency during injection and more fracturing fluid leaks into the matrix.Figure 12a shows that oil breakthrough occurs in all cases during flowback, and the oil rate increases with the increase of NF density.On the first day of converting to constant pressure production, the oil rate of case 5 reached 47.5 m 3 /d, which was 4.5 times higher than that of case 1.However, Figure 12a illustrates that oil production rate decreases with the increase of NF density after half a year of production.The reasons for this phenomenon may be: (1) as the fractures become more complex, more water is imbibed into the deeper matrix, which reduces the oil permeability and interferes the long-term production; (2) due to stress-sensitivity in fractures, the longer production time is related to greater failure of un-propped NF conductivity; and (3) the simulated NFs here are completely connected, so the pressure response quickly reaches the matrix boundary, which limits the oil drainage area.As can be seen from Table 5, the cumulative oil of case 5 is 64% higher than that of case 1 in half a year of production, and the gap narrows to 40% after one year.
contact area, ×10    Figure 12b shows that with the increase of NF density, the oil breakthrough time is gradually shortened, and the water recovery efficiency decreases accordingly.For case 1, a single-phase flow lasts for five days under the drive mechanism of fracture closure, and more than 90% of water is recovered after long-term production.This is because of limited contact area between HF and the matrix, the effect of imbibition hysteresis, and high HF conductivity.With the increase of FCI, the drive of HF closure is limited, and transient linear flow in the matrix appears in advance, so the oil breakthrough time is advanced.In addition, even in the production process, spontaneous imbibition  Figure 12b shows that with the increase of NF density, the oil breakthrough time is gradually shortened, and the water recovery efficiency decreases accordingly.For case 1, a single-phase flow lasts for five days under the drive mechanism of fracture closure, and more than 90% of water is recovered after long-term production.This is because of limited contact area between HF and the matrix, the effect of imbibition hysteresis, and high HF conductivity.With the increase of FCI, the drive of HF closure is limited, and transient linear flow in the matrix appears in advance, so the oil breakthrough time is advanced.In addition, even in the production process, spontaneous imbibition is continuing, especially for un-propped NFs with low fracture conductivity.For Case 5, although there is no long-term shut-in, the flowback rate is still less than 50%.

Impact of Proppant Distribution
Although high fracture conductivity is not required for tight reservoirs, open fractures are better than collapsed ones.Proppants with various grain diameters (ranging from 0.104 mm to 0.838 mm) are commonly used to pack fracture systems during fracturing [67].As shown in Figure 6, the fracture-closure trend is related to proppant concentration, so it is necessary to assess the effect of proppant distribution on water and oil production.Moreover, proppant embedment depending on formation mechanical properties has been extensively studied [42,[68][69][70].The fracture width loss during production is associated with embedment, especially in soft or weakly consolidated formations, which may lead to a reduction of 60% in propped fracture width [68].Therefore, the effect of compaction and embedment on fracture conductivity are considered comprehensively in this work.
Due to the lithology of Lucaogou Formation in this paper is similar to the core samples used in Chen's model [42], we introduce his formulas, Equations ( 14)-( 16), for calculating proppant embedment that allows for modifying the stress-dependent porosity and permeability [42] with the results shown in Figure 13.
In Equations ( 14)-( 16), w f0 is fracture width before embedment, m; h is the proppant embedment, m; η = 2.1 × 10 −5 , λ = 2.8, when the proppant is 20/40-mesh; and η = 1.6 × 10 −5 , λ = 3.1, when the proppant is 30/60-mesh.A series of models are established and various proppant distributions are shown in Table 6. Figure 14a shows that the flowback rate is substantially increased when HNF is propped, which is 24.8% and 32.8% higher than un-propped.When NNF is propped, there is a slight increase of flowback rate because of the restraint of low conductivity and stress-dependence in NF.Water within HF and HNF is preferred to recover during early flowback, and a considerable portion of water in NNF is imbibed into the matrix.In addition, it can also be seen from the figure that the variation trends of oil production are similar to the results of water recovery, but there is an obvious difference at various proppant concentrations in HF.The oil production decreases by nearly 25% when proppant concentration in HF decreases, while water production decreases by less than 8% because water production is mainly at the flowback stage, which is mainly driven by fracture closure at the beginning.While oil production is mainly at the production stage, and stress dependence plays a significant role with the extension of time.
Figure 14b shows that when the proppant embedding is considered, the impact on flowback is negligible because water recovery mainly occurs in the early stage of production (flowback) when the main mechanism of fracture closure is compaction.Cumulative oil production is lower when the embedment occurs, and smaller proppant size results in greater loss of production.Furthermore, compared with the embedment in HF, there is greater impact of embedment in NF on oil production.Therefore, it is beneficial to use small size proppant to promote the propping of micro-fractures, but it is also necessary to examine its embedment in micro-fractures.
Figure 14b shows that when the proppant embedding is considered, the impact on flowback is negligible because water recovery mainly occurs in the early stage of production (flowback) when the main mechanism of fracture closure is compaction.Cumulative oil production is lower when the embedment occurs, and smaller proppant size results in greater loss of production.Furthermore, compared with the embedment in HF, there is greater impact of embedment in NF on oil production.Therefore, it is beneficial to use small size proppant to promote the propping of micro-fractures, but it is also necessary to examine its embedment in micro-fractures.

Impact of Natural Fracture Heterogeneity
Natural fractures opened during fracturing cannot be uniform, and in the subsequent flowback or production process, excessive closure of some positions may be a key mechanism of water retention.As shown in Figure 15, NF width is generally non-uniform [7,50].

Impact of Natural Fracture Heterogeneity
Natural fractures opened during fracturing cannot be uniform, and in the subsequent flowback or production process, excessive closure of some positions may be a key mechanism of water retention.As shown in Figure 15, NF width is generally non-uniform [7,50].Figure 14b shows that when the proppant embedding is considered, the impact on flowback is negligible because water recovery mainly occurs in the early stage of production (flowback) when the main mechanism of fracture closure is compaction.Cumulative oil production is lower when the embedment occurs, and smaller proppant size results in greater loss of production.Furthermore, compared with the embedment in HF, there is greater impact of embedment in NF on oil production.Therefore, it is beneficial to use small size proppant to promote the propping of micro-fractures, but it is also necessary to examine its embedment in micro-fractures.

Impact of Natural Fracture Heterogeneity
Natural fractures opened during fracturing cannot be uniform, and in the subsequent flowback or production process, excessive closure of some positions may be a key mechanism of water retention.As shown in Figure 15, NF width is generally non-uniform [7,50].this work, the heterogeneity of NF conductivity is studied using the Yu et al. [71] approach to describe heterogeneous distribution of matrix permeability.Three variation coefficients of permeability are used to represent weak heterogeneity, medium homogeneity, and strong heterogeneity, as shown in Figure 16.Other model parameters and simulation processes are the same as the basic case.Figure 17 illustrates that, after one year of production, the amount of fracturing fluid remaining in reservoirs increases with the increase of heterogeneity of NF permeability.Figure 18a shows that the flowback rate decreases with the increase of heterogeneity, in which the rate of strong heterogeneity is 11.2% lower than that of homogeneity.For weak homogeneity, the water recovery efficiency is close to that of homogeneous natural fracture.Figure 18b shows that the closed natural fracture area still maintains high water saturation during production when there is strong heterogeneity.In addition, the heterogeneity of NF conductivity is 0.1%, 1.5%, and 4.1% less effective for oil production, compared with homogeneity.The reason is that high conductivity is not required for natural fractures in tight reservoirs, and apart from over-closed areas, connected fracture networks are still contributing to productivity.
permeability are used to represent weak heterogeneity, medium homogeneity, and strong heterogeneity, as shown in Figure 16.Other model parameters and simulation processes are the same as the basic case.Figure 17 illustrates that, after one year of production, the amount of fracturing fluid remaining in reservoirs increases with the increase of heterogeneity of NF permeability.Figure 18a shows that the flowback rate decreases with the increase of heterogeneity, in which the rate of strong heterogeneity is 11.2% lower than that of homogeneity.For weak homogeneity, the water recovery efficiency is close to that of homogeneous natural fracture.Figure 18b shows that the closed natural fracture area still maintains high water saturation during production when there is strong heterogeneity.In addition, the heterogeneity of NF conductivity is 0.1%, 1.5%, and 4.1% less effective for oil production, compared with homogeneity.The reason is that high conductivity is not required for natural fractures in tight reservoirs, and apart from over-closed areas, connected fracture networks are still contributing to productivity.heterogeneity, as shown in Figure 16.Other model parameters and simulation processes are the same as the basic case.Figure 17 illustrates that, after one year of production, the amount of fracturing fluid remaining in reservoirs increases with the increase of heterogeneity of NF permeability.Figure 18a shows that the flowback rate decreases with the increase of heterogeneity, in which the rate of strong heterogeneity is 11.2% lower than that of homogeneity.For weak homogeneity, the water recovery efficiency is close to that of homogeneous natural fracture.Figure 18b shows that the closed natural fracture area still maintains high water saturation during production when there is strong heterogeneity.In addition, the heterogeneity of NF conductivity is 0.1%, 1.5%, and 4.1% less effective for oil production, compared with homogeneity.The reason is that high conductivity is not required for natural fractures in tight reservoirs, and apart from over-closed areas, connected fracture networks are still contributing to productivity.same as the basic case.Figure 17 illustrates that, after one year of production, the amount of fracturing fluid remaining in reservoirs increases with the increase of heterogeneity of NF permeability.Figure 18a shows that the flowback rate decreases with the increase of heterogeneity, in which the rate of strong heterogeneity is 11.2% lower than that of homogeneity.For weak homogeneity, the water recovery efficiency is close to that of homogeneous natural fracture.Figure 18b shows that the closed natural fracture area still maintains high water saturation during production when there is strong heterogeneity.In addition, the heterogeneity of NF conductivity is 0.1%, 1.5%, and 4.1% less effective for oil production, compared with homogeneity.The reason is that high conductivity is not required for natural fractures in tight reservoirs, and apart from over-closed areas, connected fracture networks are still contributing to productivity.

Conclusions
A stress-dependent porosity and permeability model is used to numerically simulate the opening and closing of fracture systems during fracturing and flowback operations.The simulation results indicate that fracture-closure is the main driving mechanism of fracturing water recovery in the early stage of clean-up, which promotes water imbibition into deeper matrix, resulting in increased water loss.
Uncertain parameters for the fracture network are assessed with a series of mechanistic models.With the increase of fracture complexity, more fracturing water remains underground, and the oil breakthrough time advances.Although tight reservoirs do not require high fracture conductivity, proppant packing in secondary fractures (opened natural fractures) can improve well productivity.The occurrence of proppant embedment will accelerate the decline of long-term production.In addition, the heterogeneity of natural fractures may be an important factor that results in low water recovery efficiency.Excessive closure of part positions results in water locking, while its impact on production is much less compared with other factors.This work provides a new insight into the mechanism of water retention in fracture systems.

Figure 1 .
Figure 1.Schematic diagram of grid structure of numerical model (zoom scale: X = 2 times).

Figure 1 .
Figure 1.Schematic diagram of grid structure of numerical model (zoom scale: X = 2 times).

Figure 2 .
Figure 2. (a) Relative permeability for matrix; and (b) relative permeability for HF and NF.

Figure 2 .
Figure 2. (a) Relative permeability for matrix; and (b) relative permeability for HF and NF.

Figure 3 .
Figure 3. (a) Capillary pressure for matrix; and (b) capillary pressure for NF.

Figure 3 .
Figure 3. (a) Capillary pressure for matrix; and (b) capillary pressure for NF.

Figure 4 .
Figure 4. Diagram of model for HF and NF.

Figure 4 .
Figure 4. Diagram of model for HF and NF.

RatioFigure 8 .
Figure 8.(a) Bottom-hole pressure curves during injection; and (b) the distribution ratio of injected water at the end of injection.

Figure 8 .
Figure 8.(a) Bottom-hole pressure curves during injection; and (b) the distribution ratio of injected water at the end of injection.

Figure 9 .
Figure 9. Water RNP and RNP' plot during early flowback for the base case.

Figure 9 .
Figure 9. Water RNP and RNP' plot during early flowback for the base case.

Figure 10 .
Figure 10.Bottom-hole pressure profiles with various mechanisms during flowback.Figure 10.Bottom-hole pressure profiles with various mechanisms during flowback.

Figure 10 .
Figure 10.Bottom-hole pressure profiles with various mechanisms during flowback.Figure 10.Bottom-hole pressure profiles with various mechanisms during flowback.

Figure 11 .
Figure 11.Water saturation distribution from HF or NF surface to deeper matrix under flowback and production periods.From top row to bottom row are: (a) basic case, (b) ignoring stress-dependent porosity, and (c) ignoring stress dependence.

Figure 11 .
Figure 11.Water saturation distribution from HF or NF surface to deeper matrix under flowback and production periods.From top row to bottom row are: (a) basic case, (b) ignoring stress-dependent porosity, and (c) ignoring stress dependence.
breakthrough time, days Recovery efficiency, % Energies 2019, 12, x FOR PEER REVIEW 14 of 22 breakthrough time, days Recovery efficiency, % Energies 2019, 12, x FOR PEER REVIEW 14 of 22 breakthrough time, days Recovery efficiency, % Energies 2019, 12, x FOR PEER REVIEW 14 of 22 breakthrough time, days Recovery efficiency, % Energies 2019, 12, x FOR PEER REVIEW 14 of 22 breakthrough time, days Recovery efficiency, %

Figure 12 .
Figure 12.(a) Oil production profiles at various NF densities; and (b) oil breakthrough time and flowback rate profiles at various NF densities.

Figure 12 .
Figure 12.(a) Oil production profiles at various NF densities; and (b) oil breakthrough time and flowback rate profiles at various NF densities.

Figure 13 .
Figure 13.(a) Stress-dependent porosity considering both compaction and embedment; and (b) stress-dependent permeability considering both compaction and embedment.

Figure 14 .
Figure 14.(a) Cumulative oil and flowback rate profiles at various propping patterns without embedment; and (b) cumulative oil and flowback rate profiles considering embedment.

Figure 14 .
Figure 14.(a) Cumulative oil and flowback rate profiles at various propping patterns without embedment; and (b) cumulative oil and flowback rate profiles considering embedment.

Figure 14 .
Figure 14.(a) Cumulative oil and flowback rate profiles at various propping patterns without embedment; and (b) cumulative oil and flowback rate profiles considering embedment.

Table 1 .
Known field data for the base case.

Table 2 .
Assumed fracture network parameters for the base case.

Table 1 .
Known field data for the base case.

Table 2 .
Assumed fracture network parameters for the base case. .

Table 3 .
Impact of various mechanisms on the half-year production process.

Table 3 .
Impact of various mechanisms on the half-year production process.

Table 4 .
The setting of fracture complexity for cases.

Table 4 .
The setting of fracture complexity for cases.

Table 4 .
The setting of fracture complexity for cases.

Table 4 .
The setting of fracture complexity for cases.

Table 4 .
The setting of fracture complexity for cases.

Table 4 .
The setting of fracture complexity for cases.

Table 5 .
Impact of various NF densities on Flowback and production performances.

Table 5 .
Impact of various NF densities on Flowback and production performances.