A Two-Phase Flowback Model for Multiscale Diffusion and Flow in Fractured Shale Gas Reservoirs

A shale gas reservoir is usually hydraulically fractured to enhance its gas production. When the injection of water-based fracturing fluid is stopped, a two-phase flowback is observed at the wellbore of the shale gas reservoir. So far, how this water production affects the long-term gas recovery of this fractured shale gas reservoir has not been clear. In this paper, a two-phase flowback model is developed with multiscale diffusion mechanisms. First, a fractured gas reservoir is divided into three zones: naturally fractured zone or matrix (zone 1), stimulated reservoir volume (SRV) or fractured zone (zone 2), and hydraulic fractures (zone 3). Second, a dual-porosity model is applied to zones 1 and 2, and themacroscale two-phase flow flowback is formulated in the fracture network in zones 2 and 3. Third, the gas exchange between fractures (fracture network) and matrix in zones 1 and 2 is described by a diffusion process. The interactions between microscale gas diffusion in matrix and macroscale flow in fracture network are incorporated in zones 1 and 2. This model is validated by two sets of field data. Finally, parametric study is conducted to explore key parameters which affect the short-term and long-term gas productions. It is found that the two-phase flowback and the flow consistency betweenmatrix and fracture network have significant influences on cumulative gas production.Themultiscale diffusion mechanisms in different zones should be carefully considered in the flowback model.


Introduction
Shale gas reservoirs have gradually become a major source of hydrocarbon fuel in the global energy market [1][2][3].However, the extraction of shale gas confronts difficulty because shale gas reservoirs have extremely low permeability and high formation pressure [4,5].Horizontal drilling and multistaged hydraulic fracturing techniques have made shale gas production economically feasible [6,7].A horizontal well is hydraulically fractured through the injection of fracturing fluids.Generally, multistaged hydraulic fracturing techniques are used to fracture the shale gas reservoirs, forming a stimulated reservoir volume or fractured zone (SRV).The permeability of this zone is largely enhanced due to the hydraulic fractures and the microfractures in the SRV zone.The shale gas diffuses from matrix into the fractured zone and finally flows into the production well along hydraulic fractures.When the bottomhole pressure is drawn down or injection is shut in for gas production, the injected fracturing fluid (assumed water) in the hydraulic fractures and the fractured zone will flow back.This poststimulation flow period is called "flowback" [8].Usually, higher total injectedwater volume leads to larger effective fracture pore volume [9].In most cases, these flowback data are ignored if singlephase flow is used for the prediction of shale gas production [10][11][12].This ignorance of two-phase flowback makes the predicted gas production curve inconsistent with the actual field observation in the early period [13].In this paper, the two-phase flow in the fractured zone and hydraulic fractures is studied through our fully coupled model.The relationship between shale gas production rate and water flowback rate is explored.It is found that the two-phase flowback has significant impacts on the gas production rate, thus being not ignorable.
The fracture properties of fractured shale gas reservoirs have been characterized by many methods.For example, Crafton and Gunderson [14] used the high frequency singlephase flowback data to characterize fracture properties.Further, Williams-Kovacs et al. [15] developed a flowing material balance equation with single-phase flow.However, these single-phase flow models ignore the analysis for water flowback data.Their prediction accuracy thus was not well evaluated in the two-phase flow period and later stages.The immediate two-phase flow after the shut-in period was observed in the Barnett shale and the Marcellus shale [16].A transient flowing model considering the effect of twophase flow after the hydraulic fracturing was also developed for gas production [13].This model further analyzes the response of transient pressure.Yang et al. [17] developed a semianalytical model to simulate the two-phase flow during the flowback period with complex fracture networks.They found that increasing the fracture network complexity is favorable to gas production enhancement.Ezulike and Dehghanpour [6] applied a dynamic relative permeability on the two-phase flow in the hydraulic fractures.Their results showed that the relative permeability varies with reservoir parameters in the early flowback period.Xu et al. [18] analyzed the mechanisms for flowback behaviors and put forward a material balance approach to estimate the effective fracture volume in the Horn River Basin.Ezulike et al. [19] developed a two-phase flowback tank model for estimating fracture pore volume independent of fracture geometry.Their results indicated that the effective fracture pore volume is the most sensitive to fracture pore volume compressibility.Alkouh et al. [20] provided an effective method to estimate fracture volume through the water flowback and gas production data.Therefore, different methods have been used to characterize the fracture properties but the interaction between fracture flow and shale gas diffusion has not been considered in the two-phase flowback stage.
The above-mentioned two-phase flowback or singlephase gas models have an obvious limitation: the microdiffusion and macroflow mechanisms are not consistently described by a fully coupled model for the prediction of gas production rate.This inconsistency may affect the prediction of gas production rate in the two-phase flowback period and even the full gas production period.Therefore, it is necessary to analyze the two-phase flow of the fractured shale reservoir after hydrofracturing treatment.In this paper, the full interaction between the microscale gas diffusion in matrix and the macroscale two-phase flow in fractures is considered through our fully coupled model.

A Conceptual Model
In order to enhance the permeability of shale matrix, a big amount of water-based fracturing fluid is injected through the horizontal well, forming a fractured zone or SRV [21].The water is retained in the fractures and increases the water saturation.With gas production, 10-40% of the water can be recovered along the hydraulic fractures [22].In this period, the capillary pressure is directly linked with the process of two-phase flowback, which is necessary to be discussed.Hydraulic fractures are the main channels for the flow from the fractured zone to the well, whose properties are important to gas production.In the fractured zone, macroscale watergas two-phase flow in the fractures and microscale diffusion in the shale matrix occur simultaneously.With the process of flowback or continuous gas extraction, the gas pressure in the matrix began to decline and gas desorbs from the surface of matrix.This free gas diffuses from the matrix to the fractured zone and further flows into the wellbore along the hydraulic fractures.Therefore, the flow consistency in the matrix, fractured zone, and hydraulic fractures should be the focus.
Based on different flow regimes, a given computational domain is roughly divided into three zones as shown in Figure 1: shale matrix or naturally fractured zone (zone 1), (artificially and naturally) fractured zone (zone 2), and hydraulic fractures (zone 3).In zone 1, shale matrix is usually regarded as a single-porosity model in the traditional models.However, experimental measurements (such as SEM and CT) have observed a mass of micropores even microfractures in the shale [23].The gas diffusion and flow in the microfractures of shale matrix are considerably complicated [24,25].Thus, the traditional models neglecting the microdiffusion in the matrix block are not applicable.In this study, a dualporosity model is proposed to describe the gas diffusion and flow mechanisms in the zone 1.A new apparent permeability model for matrix is introduced into the dual-porosity model to consider three important physical processes: viscous flow, free molecular diffusion, and surface diffusion.In zone 2, the fractured zone is also defined as a dual-porosity model, but they have different properties.The two-phase flow is considered in the fracture and relative permeability varies with water saturation.The microfractures in this zone are different from those in zone 1, which are also induced by water fracturing.It is noted that the effective connectivity of different kinds of microfractures (in zones 1 and 2) is not described in this paper.In zone 3, the two-phase flow is still considered along the hydraulic fractures.The capillary pressure and relative permeability models are same with the fractured zone.
Overall, a series of interactions are included in the shale gas production within the reservoir.The three zones interact with each other and have a consistency impact of flow and pressure on short-term and long-term shale gas productions.In the following section, we will formulate this conceptual model by a set of governing equations in the three zones.

Governing Equations for Gas and Water Flows in Multiscale Porous Media
The three-zone model is shown in Figure 1.The flows in each zone have different diffusion and flow mechanisms.Figure 2 presents the details of the diffusion and flow in each zone.Their governing equations are stated below.

Gas exchange
Figure 2: Flow and diffusion processes in multiscale zones.

Shale Gas Flow in the Microfractures of Matrix (Zone 1).
The mass conservation law for the gas flow in the microfractures of matrix (zone 1) is governed by where  gf is the gas density in the microfractures and   is the porosity of microfractures.It is noted that only free gas is taken into consideration in microfractures of this matrix (zone 1). nw and   represent the gas viscosity and pressure, respectively. mapp is the apparent permeability for microfractures in the matrix (zone 1).This apparent permeability is directly expressed as [26] Here   is the aperture of microfractures in this zone and can be measured by many experiments such as SEM, AFM, and CT tests. ℎ reflects the degree of tortuosity in porous medium. is a Klinkenberg factor determined by the pore structure of the matrix and the reservoir temperature.
where  is the gas molecular weight. represents the universal gas constant. is the reservoir temperature.The righthand term of (1), the  mf , is the source of gas.It expresses the gas exchange between micropores and microfractures in the matrix.
where   is the gas pressure in the matrix.  is the apparent permeability of the matrix which will be described in the next section. gm is the gas density in the shale matrix. is a shape factor of matrix block.For isotropic matrix, the shape factor is [27] where   and   are the fracture spacing in the -direction and the -direction, respectively.

Knudsen Number.
In molecular dynamics of gas, the probability of gas collisions or diffusion is directly related to an important index, the molecular mean free path.This free path is the mean path for a molecular free motion in an interval of two consecutive collisions [28]: The gas molecular flow usually obeys different mechanisms, depending on the scale of flow channel.Knudsen number is defined as the ratio of the mean free path of gas flow to the characteristic size of gas flow channel:

Diffusion in Matrix.
The mass conservation law for the gas diffusion in matrix can be given as The gas mass (  ) in the matrix consists of both free phase and absorbed phase gases. gk is the gas density in the matrix.− mf means the gas exchange between matrix and microfractures of the matrix.  is the permeability of matrix which is obtained by where  is the effective diffusion coefficient of the matrix and has three components: The diffusion coefficients of   ,   ,  es are determined as follows.
For a given temperature, the gas mass in matrix transports mainly in three forms: (1) viscous flow; (2) molecular and Knudsen diffusion; (3) surface diffusion in the absorbed layer.
The viscous flow induced diffusion coefficient is [24] where   is the porosity of matrix.  is the diameter of pores in matrix.Further, the effective diffusion coefficient varies from the limit of molecular diffusion (  → 0) to fully developed Knudsen diffusion (  → ∞) [29].Therefore, the effective molecular-Knudsen diffusion coefficient is where  is the ratio of molecular size to average pore diameter.
means the fractal dimension of the pores surface and is assumed to be 2.5 in this study.  represents the Knudsen diffusion coefficient which is obtained as In addition to the gas desorption from the surface of matrix pores, gas can migrate within the absorbed layer [30,31], forming a surface diffusion.This diffusion is usually significant when the temperature is high relative to the normal boiling point of the adsorbed gas [32].Its effective surface diffusion coefficient is given as where   and  ga are the matrix density and the gas density under standard conditions, respectively.  is the Langmuir volume constant and   the Langmuir pressure.  is the surface diffusion coefficient which is expressed as [33]   = 8.29 × 10 −7 ×  0.5 exp (− 20  ) .
Combining ( 11), (12), and ( 14), the final effective diffusion coefficient for the whole matrix is given as 3.3.Two-Phase Flow in the Fractured Zone.The two-phase flow of gas and water occurs in the microfractures of the fractured zone, where the wetting phase is water and the nonwetting phase is gas.Their governing equations are obtained from the mass conservation laws of water and gas, respectively [34].
For the water phase: For the gas phase: where  is the absolute permeability of the fractured zone.
is the porosity of the fractured zone.  and  nw are the viscosities of water and gas under in situ conditions, respectively.  is the pore pressure of water in the fractures and  nw is the pore pressure of gas in the fractures.  denotes the saturation of water and  nw denotes the saturation of gas. rw and  rnw are the relative permeabilities of water and gas, respectively.Saturation, capillary pressure and relative permeability satisfy capillary pressure model and relative permeability model as given in the appendix.At the right term,   and  nw are the sources of water and gas, respectively. nw has two sources: one is the gas exchange between fractures and matrix in the fractured zone (zone 2); the other is the generated gas source at site.The gas from zone 1 is input into zone 2 as a boundary source.

Two-Phase Flow in Hydraulic
Fractures.In the twodimensional computational model of this study, the hydraulic fractures are expressed by one-dimensional lines.The twophase flow in hydraulic fractures still follows the Darcy law.The porosity of hydraulic fractures is usually enhanced or limited by the effects of pressure, swelling, and geochemical reaction [35].Xu et al. [36] indicated that 30% of the effective fracture volume is lost due to pressure depletion during early-time water flowback.In this study, these effects are not considered for hydraulic fractures and thus the porosity of each hydraulic fracture is constant.Therefore, the mass conservation equations along a hydraulic fracture are as follows.
For water For gas where  hf is the average width of a hydraulic fracture. hf and  hf are the porosity and permeability of the hydraulic fracture, respectively.  represents the source or sink of water and  nw represents the source or sink of gas.

Model Validations
The above fully coupled model describes the multiphysical processes and interactions in different scales, including twophase flow, capillary pressure, relative permeability, and multiscale diffusion.This formulates a fully coupled flowback model for the two-phase flow simulations of the shale gas production in a fractured shale reservoir.COMSOL Multiphysics, a commercial partial differential equation solver, is used for its implementation.This fully coupled model is validated through following two examples.

Example 1: Comparison of Single-Phase and Two-Phase
Flows.The gas production is predicted and compared with single-phase flow and two-phase flow models.Yu and Sepehrnoori [16] performed the history matching with two sets of field gas production data from the Barnett shale and the Marcellus shale.They analyzed the contribution of gas desorption and geomechanics on gas production through single-phase flow.The immediate two-phase flow was also observed in the Barnett shale [37].Therefore, our fully coupled model is applied to the Barnett shale.A typical computational model is set up in Figure 3.This model is 550 m long and 145 m high, and the initial reservoir pressure is 20.34 MPa.This horizontal well is located at the bottom boundary with a fixed pressure (3.45 MPa), and no-flow is assumed at its left and right boundaries.The simulation parameters and hydraulic fractures for the Barnett shale are listed in Table 1.The parameters used in this simulation are taken from the paper of Yu and Sepehrnoori [16].Most of the parameters are the same as those they used.This example has no water flowback data available, thus we take longer gas production time than two-phase flowback period for comparison.The numerical results are used to validate the prediction accuracy of long-term gas production with two-phase flowback model.The comparison of the simulation results and field data is presented in Figure 4.This figure observes that the two-phase flowback model better matches with actual gas production data than single gas flow model.After 200 days, the gas production rate declines by 80%, reaching the value of 5 × 10 4 m 3 /d.The flow rate of simulation is little lower than that from the field data during the later production.That is because the permeability in the actual condition would be higher than the uniform hydraulic fracture in the simulation model.The hydraulic fractures are usually multiscale and contact with natural fractures to form a complicated fracture network.Thus, uniform fractures are not sufficient to describe this problem.The effects of fracture properties should be further discussed in four aspects: fracture spacing, fracture width, fracture uniformity, and fracture geometry.This will be done in Section 5.2.

Example 2: Two-Phase Flowback Data of a China Shale
Gas Well.Further, both gas and water productions from an actual shale gas well [17] are used to verify this model.This computational model is 1380 m long and 250 m high and has 11 hydraulic fractures.The initial reservoir pressure is 27.4 MPa and the bottomhole pressure is 19.67 MPa.Its left and right boundaries are assumed to be no-flow.The parameters in simulations are taken from Yang et al. [17] and listed in Table 2.The gas production rate is presented in Figure 5(a) and the water flow rate is compared in Figure 5(b).It is noted that gas flow rate rapidly increases and reaches the maximum of 2 × 10 5 m 3 /d at the 10th day and then follows a slight decline.The water flowback rate has a rapid decrease from 22 m 3 /d to 2.8 m 3 /d after the first 50 days and then keeps stable during the flowback period up to 200 days.Our simulations are in a good match with both production rates of gas and water.
These two examples verify the applicability of this fully coupled model to different time periods.Example 1 focuses on the prediction accuracy of long-term gas production when the two-phase flowback is considered.Example 2 compares the production data of both water and gas in the early production period.These two examples suggest that this fully coupled model is feasible to evaluate the gas production in short-term and long-term period.

Impact of Two-Phase Flowback on Gas Production in the Early
Period.The impact of two-phase flowback on gas production in the early period is investigated here.The twophase flowback usually occurs in zones 2 and 3, because the water-based fracturing fluid can only penetrate into such a range [38].Figure 6 compares the cumulative gas productions with and without the consideration of two-phase flowback.
Generally, the cumulative gas production increases rapidly in the early 100 days.This increase then starts to slow down until 230 days.At this time, the cumulative gas production reaches 1.3 × 10 7 m 3 for the case of two-phase flowback and is 2.07 × 10 7 m 3 for the case of single-phase gas flow.Obviously, the cumulative gas production is much lower if the two-phase flowback is considered.On the other hand, the two-phase flowback would affect the gas production in the early period.This may be the reason why the previous models based on single-phase gas flow usually have a higher production rate than actual field data.

Impacts of Fracture Properties on Gas Production.
As a fine tune, the impacts of fracture properties on the longterm gas production are investigated.The following four parameters are chosen: the fracture spacing, the fracture width, the fracture uniformity, and the fracture geometry.They describe the fracture properties in the two-phase flow process.This parametric study is to explore the favorable fracture properties for the enhancement of long-term gas production.

Effect of Fracture Spacing.
The fracture spacing is an index to express the density of fractures generated in a fixed domain which has a definite volume.In this computation, each fracture has the same properties, but 14 hydraulic fractures are distributed in the hydraulic fractured zone with different spacings.Two fracture spacings of 36.2 m and 25 m are assumed.Figure 7 compares the gas production rate under these two spacings.In the same region, more sparse distribution of hydraulic fractures means larger stimulated reservoir volume (zone 2) and less matrix volume (zone 1).Thus, the pressure nearby hydraulic fracturing zone decreases more quickly and the gas in the fractured zone can flow into the well more easily.With the process of gas production, the difference of gas production rate under different fracture spacings gradually vanishes due to the completion of pressure.The gas production rate with fracture spacing of 25 m is higher than that of the original model at later production time because the gas in the zone 1 begins to flow into the well.

Effect of Fracture Width.
The fracture properties such as fracture conductivity heavily change with fracture width.In this subsection, the effect of fracture width on gas production is examined.The fracture width is taken as 0.03, 0.003, and 0.0003 m, respectively.As shown in Figure 8, the gas production rate declines rapidly.However, the history curves are almost identical regardless of fracture width.This implies that the gas production rate is not controlled by the conductivity of fracture.This is reasonable because fracture only provides a channel for gas flow.When its capacity is larger than the demand for gas flow, the fracture does not affect the gas flow any more.

Effect of Fracture Uneven
Length.After the treatment of hydraulic fracturing in shale gas reservoirs, a complex fracture network can be generated in the fractured zone.In the previous discussion, the hydraulic fractures are assumed to have equal length, called uniform fracture length.The fracture length may vary in the actual field.In order to evaluate the impact of fracture length on gas production rate, the total length of fractures is kept constant but each fracture has different length.This case is called uneven fracture length.The gas production rates of these two cases are compared in Figure 9.The gas production rate with the uniform fracture length decreases quickly and soon falls to a low point.After the production of 230 days, the gas production rate decreases to approximate 1.06 × 10 4 m 3 /d.The gas production rate with the uneven fracture length only declines to about 7 × 10 3 m 3 /d at 230 days.This may be a good explanation why the simulated gas production rate with uniform fracture length is lower than the actual filed data.

Effect of Fracture Geometry.
In the simulation model, the hydraulic fracture is usually assumed to be a straight line.This may deviate from the real case where the complicated fracture network is formed.The assumption of straight line may affect the distributions of gas and water pressure, further changing the gas production rate.In order to reveal the effect of fracture geometry on gas and water production, the gas drainage maps are presented in Figure 10.After the 36.5 days, the fractures in both geometry 1 and geometry 2 have their own drainage area.The hydraulic fractures would have an interference with neighboring fracture network to form an integrated drainage area.The drainage area varies with the geometric shape of fractures.If the height of fracture is the same, the length of fracture in geometry 2 is longer because of tortuosity.This means that the hydraulic fracture has a larger contact area with the surrounding fractured zone.Figure 11 compares the gas production rates in geometries 1 and 2. Because geometry 2 adds the contact area, its gas production rate is always higher.This effect is similar to that of fracture uniformity.
As a summary, both fracture uniformity and geometry shape add the complexity of fracture network.Their effects on gas production rate are similar.Denser distribution of hydraulic fractures means higher level of hydraulic fracturing and can enhance the gas production.However, the increase of hydraulic fracture width has less influence on gas production, if the hydraulic fracture width is larger than some value.

Evolution of Capillary
Pressure.Capillary pressure (  =  nw −   ) is the difference between gas pressure  nw (the nonwetting phase) and water pressure   (the wetting phase).
It is an important parameter in two-phase flow and is affected by many factors like wettability, contact angle, pore size distribution, and permeability.When the capillary pressure is higher, the gas production rate is greater at the beginning  of flowback but declines later [39].The water saturation is a function of capillary pressure as This equation shows that the capillary pressure is directly related to water saturation.The point A (36,20) near the hydraulic fracture is selected for observation.The evolution of water saturation with capillary pressure is plotted in Figure 12(a).The changes of water saturation and capillary pressure with time at this point A are plotted in Figure 12(b).During the gas production, the water originally in the fracture is extracted with the two-phase flowback.The capillary pressure changes from 2.4 MPa to 5 MPa.This reflects the process of gas-water displacement.As a sequence, the water saturation declines from 0.78 to 0.66 after almost 100 days.In computation, the irreducible water saturation is specified as 0.6.The free water in the fractured zone is too little to form water flow when the water saturation at point A decreases to 0.66.The remaining water adsorbs on the surface of matrix and may form a thin film of water [40].As the gas production continues, the water would be extracted with gas flow.

Impact of Entry Capillary
Pressure.Entry capillary pressure describes the occurrence of gas-water displacement and directly depends on the aperture of fracture [35].This entry capillary pressure varies with rock type and a typical range of 0.1-48.3MPa was observed [41].In this study, the effect of entry capillary pressure on gas production is investigated when entry capillary pressure is taken as 2, 4, and 6 MPa, respectively.The effect is presented in Figure 13(a) on gas production rate and in Figure 13(b) on water production rate.Both gas and water production rates decrease with the increase of entry capillary pressure.At production of 100 days, the water flowback rate declines to almost 98.8 m 3 /d for   = 2 MPa, which is higher 35.5% than that at   = 6 MPa.This difference of gas production rate due to entry capillary pressure is larger at 100 days.These results show that the increase of entry capillary pressure has much more increasing impacts on gas production rate than on water flowback rate.

Contribution of Multiscale Diffusion Mechanism to Gas
Production.The shale gas production rate has been investigated by single-scale models, in which the fractured zone has hydraulic fractures with high permeability and matrix with low permeability.Such models ignore the microfractures in matrix.This may lose the prediction accuracy of gas production rate.If the matrix is only treated as a porous medium with extremely low permeability, the microscale gas diffusion cannot be described.This produces lower gas production than the actual field data.This section will use our fully coupled model (called multiscale model) to discuss the evolution of permeability in the fractured zone and the apparent permeability in microfractures of matrix.The gas pressure in matrix and the gas exchange rate as the gas source between microfractures and matrix are also discussed in a long-term gas production.5.5.1.Fracture Permeability.When a large amount of waterbased fracturing fluid is extracted with two-phase flow, the gas pressure in the fractured zone is lower than its initial value.At this time, the gas absorbed on the surface of matrix begins to flow into the fractured zone through the microscale diffusion and flow mechanism.Figure 14(a) gives the comparison of gas production rates in single-scale and multiscale models.The gas production rate with multiscale diffusion is always higher than that with single-scale diffusion.Because of microscale diffusion, the gas in microfractures would offer a stable source to gas production.The cumulative gas production is shown in Figure 14(b).Before the production of 100 days, the cumulative gas productions in single-scale and multiscale models are almost the same.The gap between   single-scale model and multiscale model is 21.9% after 580 days.These results indicate that gas flows in zones 2 and 3 dominate the early gas production.During this period, the contribution of microscale gas diffusion in the matrix can be ignored.As the gas production continues, most of gas in the fractures is extracted and the gas pressure begins to decline.Due to the pressure difference, the gas in matrix flows into the fractured zone.This process becomes significant in gas production and the effect of microscale diffusion and flow on gas production is presented.Therefore, the microscale diffusion and flow mechanism in matrix should be carefully considered, especially for the gas production in the later period.
The permeabilities of microfractures in matrix and fractured zone are in different orders of magnitude.The evolution of permeabilities in matrix and fractured zone at point B (275, 80) are presented in Figure 15.Due to the decline of pore pressure induced by gas extraction, the pore volume in the fractured zone would decrease accordingly.Therefore, this permeability declines with production time.Such a decline is consistent with the curve of gas production rate.Finally, a 0.7% decrease of permeability in the fractured zone is observed.However, the apparent permeability of microfractures in matrix has the opposite variation.Because of the gas desorption, the matrix shrinks and the permeability of microfractures in matrix gradually increases, approaching to the permeability in the fractured zone.Though the permeability of microfractures in matrix is a few orders of magnitude lower than that in the fractured zone, it increases by 10.4% compared to its initial value.If the microfracture size is bigger, the contribution of microscale diffusion and flow to gas production is greater.

Apparent Permeability.
The above results indicate that the aperture of microfractures in matrix becomes a critical parameter to the gas exchange between fractured zone (zone 2) and matrix (zone 1).The aperture was observed to vary from 10 nm to 50 m [42].The value of aperture is used to divide the matrix system into microfractures and micropores in this paper.Here, the apertures of microfractures are assumed to be 100 nm and 900 nm, respectively.Figure 16 presents the variation of apparent permeability with aperture at point B. The black line represents the base case that was presented in Figure 15.The red dashed line represents the apparent permeability when the microfracture aperture is 900 nm.This figure shows that the aperture of microfracture has important impacts on apparent permeability.When the aperture is 900 nm, the apparent permeability is nearly two orders of magnitude higher than the base case.This value would get closer to the permeability of fractures (zone 2) as the aperture increases.Previous study shows that with the pressure decline in the fractured zone, the gas in matrix flows into the fractured zone and dominates the gas production.The aperture of microfractures has a salient influence on the apparent permeability.Therefore, the gas in matrix through the microscale diffusion and flow is the key contribution to gas production, especially in the later production period.

Gas Exchange between Microfractures and Matrix.
Pore size is an important parameter to gas production because it is the channel of gas diffusion from the matrix into the microfractures.When the diameters of pores in matrix are assumed to be 1, 3, and 9 nm, respectively, Figure 17 presents the variation of gas pressure with diameter.It is found that the pore size is greater and the pressure drops faster.More gas desorbs from the surface of matrix when the diameter of pore in matrix is 9 nm.The impact of pore diameter on the gas exchange rate between matrix and microfractures is presented.These variation curves all present the shapes of mountain in Figure 18.Due to faster gas pressure drop, the gas exchange rate firstly increases when the diameter of pore is 9 nm.With the decrease of pore diameter, this increase of gas exchange rate would delay correspondingly.Further, the diameter of pore is bigger, the diffusion coefficient is bigger, and the gas exchange rate reaches higher peak.That is because bigger diameter of pore is directly linked to higher matrix permeability.Finally, the gas exchange rate becomes zero.It means that the pressure in matrix and the pressure in microfractures reach a balance.More gas desorbs from the surface of matrix when the diameter of pore is bigger.This is the gas supply to gas production after the gas in the fractured zone is extracted.

Conclusions
This study numerically investigated the impact of water flowback in the water-based fracturing fluid on the shortterm and long-term shale gas production.A fully coupled multiscale two-phase flowback model was proposed for the three zones model and the multiscale diffusion mechanisms were incorporated.Both water and gas production rates in the early period and the gas production rate in the later production period were studied under different fracture properties like fracture spacing, fracture width, fracture uniformity, and fracture geometry.The evolution curves of capillary pressure and water saturation, the contribution of multiscale diffusion to gas production, and the gas exchange between microfractures and matrix in different zones were explored.From these studies, the following conclusions can be made.
First, this fully coupled multiscale two-phase flowback model can not only describe the process of two-phase flowback in the early production period, but also well describe the impacts of fracture properties and multiscale gas diffusion in the later production period.It thus provides a useful tool for the assessment of gas production in the waterbased fracturing process.The two-phase flow affects the gas production in the early and long-term periods.Our example shows that the cumulative gas production at 230 days has a 58.2% decline after considering the effect of twophase flowback.Thus, the two-phase flowback should be carefully included when the gas production in fractured shale reservoirs is evaluated.
Second, the hydraulic fracture properties of fracture spacing, fracture width, fracture uniformity, and fracture geometry have variable influences on gas production rate.The increase of fracture density has an obvious positive enhancement on gas production rate.The uniformity and geometry of fractures have similar effects on gas production rate, particularly at the early production period.The fracture width may have no effect when it is larger than some value.The hydraulic fracture network in the actual field is interconnected and not uniform.Their impacts on gas production rate should be further studied.
Finally, the gas production rate with multiscale twophase flow model is always higher than that with singlescale model.In microscale diffusion and flow mechanism, the gas in microfractures of matrix is a stable source to gas production.With gas extraction, the permeability of microfractures in matrix gradually increases and approaches to the permeability in the fractured zone.This phenomenon reflects the flow consistency in the macroscale flow and the microscale diffusion.The contribution from microfractures is

Figure 5 :
Figure 5: Comparison of numerical simulation with field production data at some China shale gas well.

3 )Figure 6 :
Figure 6: Comparisons of cumulative gas production for singlephase gas flow and two-phase flow.

Figure 9 :
Figure 9: Effect of fracture uniformity on gas production rate.

Figure 10 :Figure 11 :
Figure 10: Effect of fracture geometry on gas pressure at the early and late production periods.
Variation of water saturation and capillary with time

Figure 12 :
Figure 12: Change of water saturation and capillary with time at point A (36, 20) near the hydraulic fracture.

Figure 13 :
Figure 13: Impact of entry capillary pressure on gas and water flowback rates.

Figure 14 :
Figure 14: Comparison of gas production prediction by multiscale and single-scale models.

Figure 17 :Figure 18 :
Figure 17: Gas pressures in matrix at different pore sizes.

Table 1 :
Model parameters for Barnett shales.Figure 4: Comparison of numerical simulation and site production in Barnett shale reservoir.

Table 2 :
Simulation parameters for field data.
Figure 16: Evolution of apparent permeability with time under different aperture.