A semi‐analytical model for capturing dynamic behavior of hydraulic fractures during flowback period in tight oil reservoir

Hydraulic fracturing has been successfully employed for unconventional oil and gas recovery for decades. During flowback, the closure of the fracture may exhibit with the pressure drop of fracturing fluid dewatering. However, fracture closure always is ignored or treated as stress‐dependent fracture properties in previous flowback models. This paper presented a dynamic fracture model, which can comprehensively capture the dynamic behavior of hydraulic fractures during the flowback. A nonlinear relationship between fracture aperture and contact stress acting on the fracture surfaces is adopted to simulate fracture closure. The fracture aperture calculated by the displacement discontinuity method (DDM) is used to characterize the fracture pore volume and fracture conductivity, which will be dynamically updated in the flow model. Then, the pressure and saturation of each phase, along with the displacement on the fracture surface, are calculated by solving flow equations and geomechanics equations with iterative coupling approach. The new semi‐analytical model is validated by comparing it with a fully coupled stress‐porosity pressure numerical simulation model setup by ABAQUS® and CMG. Then, the dynamic behaviors of hydraulic fractures are investigated in detail by several cases. Results show that fracture closure is an important reason for the decline in production during the flowback and early production. And it is more important to enhance the properties of the stimulated reservoir volume (SRV) than to only create a fracture with high conductivity. Lastly, the key parameters (the fracture effective length and fracture conductivity under variable contact stress) can be interpreted by history‐matching the field flowback data.

els. This paper presented a dynamic fracture model, which can comprehensively capture the dynamic behavior of hydraulic fractures during the flowback. A nonlinear relationship between fracture aperture and contact stress acting on the fracture surfaces is adopted to simulate fracture closure. The fracture aperture calculated by the displacement discontinuity method (DDM) is used to characterize the fracture pore volume and fracture conductivity, which will be dynamically updated in the flow model. Then, the pressure and saturation of each phase, along with the displacement on the fracture surface, are calculated by solving flow equations and geomechanics equations with iterative coupling approach. The new semi-analytical model is validated by comparing it with a fully coupled stress-porosity pressure numerical simulation model setup by ABAQUS ® and CMG. Then, the dynamic behaviors of hydraulic fractures are investigated in detail by several cases. Results show that fracture closure is an important reason for the decline in production during the flowback and early production. And it is more important to enhance the properties of the stimulated reservoir volume (SRV) than to only create a fracture with high conductivity. Lastly, the key parameters (the fracture effective length and fracture conductivity under variable contact stress) can be interpreted by history-matching the field flowback data.

K E Y W O R D S
flowback, fracture closure, fracture dynamic behavior, semi-analytical model of the fracture network is of vital importance for production evaluation and prediction. Monitoring the hydrocarbon production data is an important method to determine the characteristic of the fracture networks. However, during the flowback of the fracturing fluid, only a small amount of hydrocarbon production data can be recorded, causing it difficult to extract the fracture parameters. In order to quickly evaluate the effectiveness of fracturing treatment, accurately predict horizontal well production, and make up suitable treatments particularly in the early time of well life, analyzing frequently monitored fracturing fluid flowback data to characterize fracture networks is all important.
There have been three methods to quantitative analyze the flowback data: (a) The analytical model 3 was used to analyze the flowback data in water-dominated production stage to extract the hydraulic fracture parameters according to the linear relationship between rate normalized pressure (RNP) and material balance time (MBT). (b) The most accurate method for quantitatively analyzing the fracturing fluid flowback data is the numerical simulation method, 4,40 which can rigorously account for multiphase flow, gravity segregation, and pressure-dependent permeability. Many numerical flowback models 5,6 have been established to comprehensively describe the flow regimes and rigorously characterize the reservoir flowback. (c) The semi-analytical model describing the multiphase flow during flowback has received more attention because of its convenience and efficiency. Many two-phase (oil + water) flow semi-analytical models [7][8][9] have been presented to capture the fracture key parameters (eg, effective fracture length and fracture permeability) by history-matching the flowback data. The initial condition of the flowback can be accurately considered by coupling hydraulic fracture model, 10 and the semi-analytical model also can be used in complex fracture network by the star-delta transformation used to solve the interplay flow of fracture intersection. 11 However, current flowback models always ignored the fracture closure or simply treated the process of fracture closing as stress-sensitive fracture permeability and porosity, which cannot accurately describe the fracture closure because of the inaccurate effective stress calculation.
Laboratory experiments indicated that the fracture will gradually close with pressure drop within the fracture 12 and the permeability of the fracture will also decay with effective stress increase during a reduction of fracture pore pressure caused by production. 13 During the flowback period, there are many researches 6,7 found that the parameters of the fracture calculated by two-phase flow data after hydrocarbon breakthrough are less than that calculated by single-phase flow data before hydrocarbon breakthrough. In addition, the permeability loss of fracture during the production was also observed by well-testing analysis 14 and production decline analysis 15 of field data. These evidences show that the closure of the fracture, causing the fracture properties (eg, fracture permeability) deteriorate significantly, occurs with the reduction of fluid pressure within the fracture throughout the flowback and production. The process of fracture closure, therefore, has an enormous impact on well production performance and should be considered during the flowback and production.
The impact of fracture closure on well production has been considered and investigated in many previous works. Moinfer et al 16 incorporated the dynamic behavior of fracture into an embedded discrete fracture model (EDFM) in tight oil reservoir. The empirical joint model 6 or the fracture conductivity table was used to represent the dynamic behavior of fracture, and the impact of fracture closure on production was discussed by their research. Similarly, the correlation between fracture permeability and effective stress proposed by Raghavan and Chin is used in Ayber's model to evaluate the long-term effect of natural fractures closure on gas production from unconventional reservoirs. 17 However, these models had not fully coupled geomechanics and flow, and the stress caused by fracture deformation cannot be calculated, causing the inaccurate calculation of effective stress. Seth et al 18 developed a model that couples the flow and width in the fracture. The stress variations introduced by fracture deformation during the fracture closure can be captured by modeling hydraulic fractures as explicit discontinuities. Though there are many methods to evaluate the impact on well production, where the fracture properties are sensitive to stress, the stress field calculation is not coupled into the two-phase flow model during flowback and early-time production. In addition, the inputs of the relationship between fracture parameters and contact stress are more dependent on laboratory data.
The objectives of the current work are to develop a model coupling two-phase flow and geomechanics to accurately describe the dynamic behavior of the hydraulic fracture during flowback and early production in tight oil reservoirs and capture the dynamic parameter of fractures by quantitative analysis of flowback data. Considering the stress caused by fracture deformation, the dynamic behavior of the fracture aperture is accurately calculated by combining the DDM with the relationship between fracture aperture and contact normal stress. The flow equation and the stress equation are coupled by continuously updating the pore volume and fracture conductivity, which can be represented by the fracture aperture in fracture flow equation. The fracture aperture associating with fracture pressure and fracture saturation can be solved by the coupled model. The advantage of this model is that the stress acting on the fracture surfaces associated with the dynamic behavior of the hydraulic fracture aperture can be accurately calculated. When the closure of the hydraulic fracture is taken into account, the dynamic parameters of fracture can be obtained by history-matching the frequently monitored flowback and early-production data.  Figure 1A shows a horizontal well treated by three-stage hydraulic fracturing. Due to the injection of a large amount of fracturing fluid, complex fracture systems are formed with hydraulic fractures connecting the natural fractures. In order to model this system, it is reasonable to assume that all the fracturing stages are located uniformly along the wellbore with the same properties because of the common field practice to design equally spaced hydraulic fracturing with the same properties. As Figure 1B illustrates, each fracturing stage is idealized as two regions: a primary hydraulic fracture (PHF) and an adjacent reservoir region with enhanced permeability due to nature fractures reactivation. 9,19 The PHF is assumed to be initially filled with treatment water and majority of proppants. Although treatment water will leak off to the formation during stimulation, it is difficult to quantify the distribution of water in the matrix. It is therefore assumed that the mobile water is only in the hydraulic fracture. This assumption has been adopted in several models to describe two-phase flowback. 7,11,20,21 The distribution of the water in the matrix will be investigated in our future work. Because of the symmetry of the system, a representative fracturing stage is considered, and thus, the system shown in Figure 2A can be used for modeling.
Due to ultra-law permeability of unconventional tight oil reservoir and short duration of flowback, the contribution of the reservoir beyond the stimulated volume is usually negligible. 22,23,25 Thus, as Figure 2A illustrates, only the stimulated reservoir volume is considered in this model. And the fracture is divided into several fracture segments, as Figure 2B illustrates. Considering the drainage volume is limited to a very small region surround the fractures, Jia et al 4 indicated that the matrix linear flow can be equal to the case where the matrix is divided into several isolated segments with no-flow boundary and the matrix flow in the isolated region of each segment is independently modeled by transient linear flow, as illustrated in Figure 2C. And the out boundary specified by the blue dashed lines in Figure 2C can be assumed as an infinite boundary because of short duration of flowback.

| Fracture stress model
The high-pressure fluid is injected into the formation to open the fractures during the hydraulic fracturing and then shut-in. The fluid pressure in the fracture gradually decreases with fracturing fluid leak-off during the shut-in period. However, due to the ultra-low permeability and the insufficient shutin time, the initial fluid pressure within the fracture may be higher than the minimum horizontal principal stress. The fracture, which called open fracture, is supported by the internal high-pressurized fluid and maintains its maximum aperture. At this time, the proppant is in suspension state and cannot form an effective supporting effect, as illustrated in Figure 3A. As a result, the total normal pressure acting on the fracture surfaces only is the internal fluid pressure, and the shear stress on the fracture surfaces is zero because the liquid cannot withstand shear stress, as Equations 1 and 2 show: where n is the normal stress acting on the surface of the fracture. p F is the fluid pressure in the fracture. is the shear stress.
The fracture aperture will decrease during the reduction of the fracture pressure caused by dewatering, and the fracture will gradually close on the rough surfaces or proppant. When the fracture closes on the rough surfaces or proppant, as illustrated in Figure 3B, the total normal pressure acting on the surfaces of fracture is fluid pressure plus a contact (1) n = p F (2) = 0 F I G U R E 1 A basic representation of a multistage fractured system. A, Schematic of a fractured horizontal well. B, Idealized multifractured model. Each fracturing stage is idealized as two regions: a primary hydraulic fracture (PHF) and an adjacent reservoir region with enhanced permeability due to nature fracture reactivation during the stimulating treatment pressure. Coulomb's law is suitable for the shear stress on the surfaces of the fracture, as Equations 3 and 4 illustrate: where c is the contact stress between the fracture surface and the rough fracture surface or proppant. is the fraction coefficient. S 0 is the cohesion.
A homogeneous, isotropic linearly elastic formation with preexist hydraulic fracture is assumed in this paper. Figure 3 shows the stress model used in our coupled model, the fracture is divided into several elements, which is the same as the flow model, and the fracture surfaces stress state depends on whether the fracture is closed or not.

| Characterization of fracture aperture and its relationship with fracture flow parameters
The fracture closure gradually exhibits as the fluid pressure drop during the flowback. In this process, the change of the fracture aperture controlled by the fluid pressure and the stress acting on the fracture surfaces has a significant influence on the fracture pore volume and the fracture conductivity. The fracture mechanical aperture refers to the distance between the two surfaces of the fracture, reflecting the geometrical information of the fracture under formation conditions. While the fracture hydraulic aperture reflects the fracture conductivity, get by fracture conductivity test. These two fracture apertures are the key parameters to capture the fracture dynamic behaviors and couple seepage field and stress field.  12 present that there is still non-negligible residual aperture (compared to the maximum aperture during fracture propagation) after fracture closure. Wang and Sharma analyzed field cases of diagnostic residual fracture tests and showed the existence of stressdependent residual fracture width after fracture closure. So the closure of fractures does not mean that fractures do not have conductivity. 25 Even though in the absence of proppant, the roughness of the fractures surfaces still leads to a certain conductivity after fracture closure.
When fracture closes on proppant pack, an empirical formula is introduced to relate contact stress to fracture mechanical aperture 18,25-28 : where c,prop,Wref is the contact reference stress for the proppant pack, which denotes the contact stress at which the fracture aperture is reduced by 90%, and W c,prop is the local propped fracture aperture when the contact stress is equal to 0.
However, when considering the roughness of the fracture surfaces and the proppant in the fracture, the fracture mechanical aperture cannot be directly used to indicate the conductivity of the fracture. Therefore, the fracture hydraulic aperture is induced to indicate fracture conductivity. 29 The closed fracture hydraulic aperture has an exponential relationship with the contact normal stress, and the shear displacement also has an influence on the fracture hydraulic aperture. 13 Equation 5 also can be used to describe the relationship between the fracture hydraulic aperture and contact normal stress 27 : where w c is the local propped fracture hydraulic aperture when the contact stress is equal to 0, c,Wref is the contact reference stress for the proppant pack, which denotes the contact stress at which the fracture hydraulic aperture is reduced by 90%, and wdil is aperture dilation angle. D w,eff is the effective shear displacement. Lee and Cho pointed out that when the shear displacement reaches a certain level, it does not cause an increase in fracture aperture. 13 Thus, D w,eff is defined as equal to D if D ≤ D w,eff,max , and equal to D w,eff,max otherwise.

| Method to distinguish open and closed fracture conditions
When the contact stress is equal to 0, the fracture mechanical aperture,W c , will be given as the critical fracture aperture to distinguish the closed and open fracture. After each iterative calculation, the calculated fracture aperture of each element will be compared with the critical fracture aperture to determine the closed state of the fracture elements.
The following equation is adapted to describe the opening fracture mechanical aperture in order to ensure the continuity of the fracture aperture during fracture closure 27 : where W and w are fracture mechanical aperture and fracture hydraulic aperture separately, and W open is the difference between the opening fracture aperture W and the critical fracture aperture W c .

| The quantitative relationship between fracture apertures and flow parameters
The fracture mechanical aperture can be used to describe the change of the fracture pore volume during the closure. It assumes that the proppant is incompressible; however, the fracture pore volume is compressible. The distribution of proppant in the fracture during the closure is shown in Figure 4. At the beginning of the flowback, the proppant in the fracture may be in a suspended state ( Figure 4A) when the fracture surfaces are not exposed to contact stress. With the fracture pressure declining, the fracture will gradually close on the proppant, and Figure 4B is a critical state when the contact stress acting on the fracture surface is exactly zero. Because of the shear stress acting on the proppant, the proppant particles will slip causing the fracture pore volume to be compressed. Figure 4C shows the state in which the pore volume of the fracture pore volume is compressed.
The critical proppant volume concentration, c c , is defined as: where V p is the proppant volume and V Fc is the critical fracture volume when the contact stress is equal to 0.
It is assumed that the proppant is incompressible; therefore, the fracture porosity can be defined as: where V F is the fracture volume, V por is the fracture pore volume, F is the porosity of the fracture.
Cubic law was found to hold whether the fracture is open or closed by Witherspoon et al. 30 The fracture hydraulic aperture is induced to indicate fracture conductivity. And the fracture conductivity is shown as: where w is the fracture hydraulic aperture. k F is the permeability of the fracture.
However, Equations 6 and 11 only have been valid and used for the unpropped fractures. In order to investigate whether Equations 6 and 11 can also characterize the conductivity of fracture closed on proppants during the flowback, we compare proppant pack conductivity calculated by our model and Barree's general model, 31 as illustrated in Figure 5. Barree's general model is a set of correlations that can predict the proppant pack conductivity based on a long history of proppant testing. The conductivity of the different particle diameter of brown sand at the uniform proppant concentration under different contact normal stresses can be calculated by the general model. The conductivity of the fracture can be calculated by the following equation 31 : k1, m1, s1, m2, w1, m3, m4, Sharp and Exp are regression coefficients, which is different for different proppant types. K 0 is the unstressed pack permeability. d m is the media diameter of proppant. C is proppant pack concentration. K s is pack permeability at a given stress.
By adjusting the corresponding parameters (w c , c,wref ), the conductivity of the fracture calculated by the Equations 6 and 10 can be matched with the calculation result of Barree's model with the contact normal stress <8000 psi. Although there will be a larger offset when the contact normal stress is larger, 0-8000 psi is completely suitable for the contact stress on the fracture surface of the general unconventional reservoir. The type of proppant, the size of proppant particles, the proppant concentration, and the formation mechanical properties all are upscaled into three parameters (w c , c,wref , and wdil ), so this model cannot describe the influence of the proppant properties in details. This influence will be considered in our next work, and the general model proposed by Barree will be used in our next model to accurately describe the relation between fracture aperture and contact normal stress with considering the different proppant properties.

| Mathematical model
This model includes two domains: matrix and hydraulic fractures, and solves two systems of equations: flow equation system and geomechanics equation system. The fracture is divided into N F segments, and 5N F unknowns (pressure, saturation, matrix flux, normal displacement, and shear displacement for each segment) are solved in this model. During the flowback, formation fluid flows from the matrix into the fracture and then along the fracture into the wellbore. The matrix single-phase (oil) flow is simulated by the analytical transient linear flow equation and dynamically coupled with the fully numerical fracture multiphase (oil + water) flow equation by imposing continuity of pressure and flux on the surfaces of the fractures. At the same time, a nonlinear formula is used to describe the relationship between fracture aperture and the stress acting on the surfaces of the fracture, and the fracture aperture calculated by the displacement discontinuity method (DDM) is used to characterize fracture pore volume and fracture conductivity, which is dynamically updated in flow equations, at every time step.

| Numerical model for two-phase (oil + water) fracture flow
The volume leaking into the surrounding matrix can be ignored due to the extremely low permeability of the matrix and the short duration of flowback period. Darcy flow is assumed. Therefore, under the constraints of S w = 1 − S o and p w = p o − P cow S o , the one-dimensional diffusion equation of the two-phase flow in the fracture can be expressed as follows 23 : for the oil phase, and for the water phase, where S Fo and S Fw are oil saturation and water saturation in the fracture, respectively. t is the time. q SCFo is a source term representing the mass flow of oil from the matrix to the fracture per surface area of the fracture. q SCWo and q SCWw represent the oil and water production rater per surface area of the fracture, respectively. B o and B w are the formation volume factor of oil and water, respectively. k F is the fracture permeability. k rl is the fracture relative permeability of the fracture fluid. l is the viscosity of the fluid. p Fl is the pressure of the fluid in the fracture. From the finite volume method, the diffusion equations (Equations 17 and 18) are discretized for each individual element. The two-phase flow equation for each fracture segment can be expressed using the following nonlinear matrix equation: where A and b are the coefficient matrix for each fracture segment, as defined in Appendix A. X is a vector composed of each fracture segment unknown, where X = X 1 ⋯X N F T and q is a vector consisting of the flow rate from the matrix to the fracture segment. q = q 1 ⋯q N F T and

| Linear flow model for matrix flow
As Figure 2C illustrates, the matrix is divided into isolated percolation regions that are consistent with the fracture segments, and the flow of the matrix to the fracture within each region can be considered as an independent linear flow. The flow in each matrix region can be expressed using the transient linear flow equation without considering the effects of other fracture segments. 4 Since the flowback process lasts for a short period of time and the drainage volume is confined around the fractures, the outer boundary condition can be considered as an infinite boundary for each matrix region. The flow rate of the matrix into the fracture changes with time, so the internal boundary condition of the matrix region is a variable-rate inner boundary. The planar linear flow equation proposed by Wattenbarger et al 32 for the infinite boundary with constant production can be used to solve the pressure of the fracture segment i as follows: where p I is the initial formation pressure. c2 and c2 are the unit conversion factor. h is the thickness of the formation. ΔL F i is the length of the fracture segment i. M is matrix porosity. c tM is total compressibility of matrix. k M is the matrix permeability. In addition, the fracture segment can be thought of as the source within the matrix unit, so the pressure of the fracture segment p Foi is equal to the pressure of the source terms P Mi of the matrix unit. The superposition principle is used to calculate the pressure of the fracture segment with the variable-rate inner boundary: With the introduction of H, which is defined in Appendix B, Equation 21 applied to each fracture segment can be simplified and rearranged as follows.
where H and h are the coefficient matrix for each fracture segment, as defined in Appendix B. q SCFo is a vector consisting of the flow rate from the matrix to the fracture segment.p Fo is a vector consisting of each fracture segment pressure. P I is a vector consisting of initial formation pressure.

| Geomechanics Model for fracture aperture
Displacement discontinuity method, 31 which is a special boundary element method, is used to calculate the fracture aperture. DDM has been widely used in many simulators to model fracture growth because of its efficiency of computing comparing to the finite element method. 26,33,34 Besides, it is also an outstanding method for matching the semi-analysis solution of the flowback, because both problems are solved only on the fracture domain and the interaction of the matrix domain can be simplified. When the fracture is deformed, the normal stress acting on the fracture surfaces also changes. The deformation of a fracture segment causes a change in the entire stress field, so the stress caused by the deformation of any fracture segment is equal to the superposition of the stress caused by the deformation of each fracture segment at its position. Therefore, the stress caused by deformation in fracture segment i at n + 1 time step can be calculated by the following formulation 31 : A planar fracture with a height equal to the thickness of the formation is assumed in this model. Therefore, a correction factor G adj,i,j derived by Olson is used to correct the DDM. This factor is also be adapted in many hydraulic models. [34][35][36] where d is the distance between the centers of the two segments i and j. h is the thickness of the formation. is an empirical coefficient equal to one, and is an empirical coefficient equal to 2.3. 27,34,36 The normal stress acting on the surfaces of the fracture segment i is equal to the component of the in situ stress in the direction perpendicular to the fracture surface plus the normal stress caused by the fracture deformation. And the shear stress is equal to the shear stress component of in situ stress on fracture surface plus the shear stress caused by the fracture deformation.
For the open fracture segments: For the closed fracture segments, as mentioned above, the aperture of the fracture has a relationship with the contact normal stress, and the corresponding solution equation can be obtained.
where n+1 c,i is the contact stress of the segment i and can be calculated by the following equation: The shear stress of the closed fracture segments: The slip state of the fracture needs to be judged at each iteration. If the left side of the equation is larger than the right side, we define that the fracture slips, and the less than or equal sign of the equation becomes an equal sign. On the contrary, we define that the fracture is locked, and the shear discontinuous displacement of the previous time step is assigned to be the same with the last time step.
Whether the fracture segment is closed or open, it will have a normal stress equation and a shear stress equation, respectively. Therefore, there are 2N F geomechanics equations for N F fracture segments.

| Solving the coupled system
In the coupled model, there are 5N F unknowns (pressure, saturation, matrix flux, normal displacement, and shear displacement for each fracture segment) need to be solved with 2N F equations (Equations 19, 22, and 28-32 for each fracture segment). For efficiency and convenience, the iterative solving method adopted by McClure is used in this model. 27,34 At n + 1th time step, (a) the shear stress equations (Equations 29 and 32) are solved with shear displacement D n+1 as unknowns holding pressure, saturation, and normal displacement constant. The slip state of the fracture segment will be checked after each subloop, and the corresponding shear stress equation will be used. and D′ n+1 , if the difference is less than the tolerance, the solution process converges; otherwise, let D′ n+1 equals D n+1 and repeats step (c) until convergence occurs. The flow chart in Figure 6 describes the structure of the program.

| Validation
We simulated two cases using the new semi-analytical model, and the simulation results compared to ABAQUS and CMG, respectively.   Figure 7) is established using ABAQUS. The global grid system is with the grid configuration of 64 × 250 × 1. The grid size in the y direction is 1 ft, and a transition mesh (grid size from 1 to 10 ft) is used in the x direction. Single-phase water is injected for 100 s to create a hydraulic fracture at a constant injection rate of 0.001 m 2 /s, followed by a constant rate of 0.0001 m 2 /s to product. The setting of the model input parameters refers to Chen et al 37 and Zanganeh and Clarkson 38 as shown in Table 1. The fracture length and internal fluid pressure at the end of the injecting obtained by the ABAQUS hydraulic fracturing model are input into our coupling model as initial conditions. The fracture is divided into 21 segments with 6 ft length, and the initial pressure of each segment is shown in Table 2. The water saturation in the fracture is set to 1, and the other parameters are consistent with the finite element model established by ABAQUS, as Table 1 shows. The flow parameters and mechanical parameters of fracture under constant dewatering rate calculated by our model and ABAQUS are compared. The curve of the bottom hole pressure, the fracture aperture and the fluid pressure internal of the fracture calculated by the two models at the initial time (the beginning of the production), 125, 305, 485, 665, and 864 seconds are plotted in Figures 8,  9 and 10, respectively. Figure 8 illustrates that without considering leak-off, the bottom hole flow pressure under constant dewatering rate is linear with the time until the bottom hole flow pressure is close to the minimum horizontal principal stress (4350 psi). When the fracture aperture is large, the fracture has a large conductivity and then the fluid pressure distribution is almost uniform, as shown in Figure 10 with a production time <864 seconds. However, when the fluid pressure in the fracture is close to the closed stress, the fracture can no longer be regarded as the infinite conductivity fracture because of its small aperture, and an increasing pressure drop funnel will be formed near the wellbore, as shown in Figure 9 with a production time equal to 864 seconds. At this time, the fracture shape will no longer be the arch with the widest width at the bottom of the well, as shown in Figure 9 with production time <864, but forming an "M" shape in the symmetric fracture, as shown in Figure 9 with production time equal to 864 seconds. When the crack is closed, this phenomenon is more pronounced because the conductivity of the crack is greatly reduced. Figure 9 shows that the fracture aperture gradually decreases with the pressure drop within the fracture under a uniform production rate. Usually, we think that the fracture does not completely close, but there is a residual width. Figure 10 shows the global grid system with the grid configuration of 101 × 95 × 1. The grid size for the reference solution by the CMG is 10 ft × 10 ft. In addition, the grids around the fractures are logarithmically refined (grid size from 0.01 to 10 ft) as shown in Figure 10. In total, 9595 grids are used in the CMG model. The inputs of reservoir and fracture properties are shown in Table 3. The relationship between contact stress and fracture permeability is calculated by Barton's formula and given by a lookup table, as shown in Table 4. Figure 11 shows the corresponding PVT properties for the oil. And the oil and water relative permeability is model with the following equations 4,8 : where n and m are both assumed to be 1.

| Case 2: Comparison with CMG model
The water saturation of fracture grids is 100% initially, while the matrix is assigned 100% oil saturation (single-phase flow). The well is produced at a constant bottom hole pressure (3500 psi). Figure 12 provides a comparison of the results of the fully numerical model and our new semi-analytical model. Excellent agreement between the oil and water production profiles from the commercial simulator and the new model suggests our new semi-analytical model is reasonable and accurate.

| Simulation case
In this section, three cases are simulated by this coupled model to describe the fracture dynamic behaviors during the flowback and to investigate the effect of initial fracture hydraulic aperture and contact reference stress of fracture closure on production performance, respectively.
Case 1 is designed to characterize the dynamic variation of fracture aperture and conductivity during the flowback and early production. The inputs of reservoir and fracture properties are shown in Table 5. The initial pressure within the fracture is above the minimum horizontal principal stress to ensure that the initial fracture is open. The initial oil saturation in the fracture is set to 0 to ensure single-phase flow until the hydrocarbon breakthrough. The well is produced at a constant rate of 500 bbl/d, and the minimal BHP is 3500 psi. The oil PVT properties and the oil and water relative permeability are the same with that used in Section 3.1.2.
Since the fracture pore volume does not change abruptly when the fracture is closed, the fracture hydraulic aperture characterizing the fracture conductivity is more important for the oil/water production profile and performance during flowback and early-time production than the fracture mechanical aperture characterizing the fracture pore volume. The relationship between the fracture hydraulic aperture and the contact normal stress during the flowback is expressed by Equation 6, where the effects of the hardness of the rock, the type of proppant, the particle size of the proppant and the sanding concentration are upscaled to three parameters (initial fracture hydraulic aperture, contact reference stress for fracture aperture, and friction coefficient). Because the symmetrical hydraulic fracture perpendicular to the direction of the minimum horizontal principal stress is assumed in the model, the shear displacement does not occur during the fracture closing process, and the influence of the friction coefficient on the closed fracture aperture is negligible. In order to characterize the production performance and determine the appropriate range of parameters (initial fracture hydraulic aperture and contract reference stress) values during the flowback within 100 days, case 2 and case 3 are designed to examine the effects of two parameters: (a) initial fracture hydraulic  Table 3 except contact reference stress and initial fracture hydraulic aperture. And the well is produced under a constant bottom hole pressure (3500 psi). The parameters investigated are listed in Table 6.

| Case 1: Dynamic behaviors of fracture during flowback and early production
In this section, the dynamic variation of fracture aperture and conductivity is characterized by the coupled model. Figure 13 illustrates the fracture aperture and fracture conductivity at 0, 0.1, 0.3, 0.8, 5, 10, and 100 days during the flowback. Analysis of Figure 13 shows that based on the stress state of the fracture, the fracture closure may be divided into three stages. (a) The fracture, only supported by high-pressure fluid within it, is in the open state with a large mechanical aperture and conductivity. The shape of the fracture is wedge-shaped as shown in Figure 13 at 0 days. The fracture aperture and the fluid pressure within the fracture decrease rapidly as a result of fluid withdrawal during flowback. In this stage, the portion of flowback caused by effective fracture pore volume closure is the main drive mechanisms. However, this state always occurs in the shut-in period because the fracturing fluid leak off into the formation causing the fracture pressure to be below the minimum horizontal principal stress. This state occurs in the flowback only when the shut-in time is insufficient or "forced closure" measures are taken. (b) The pressure gradient within the fracture is small because of the high conductivity; therefore, the fracture aperture decreases with the similar trend until the fracture closes on the proppant, as shown in Figure 11 from 0 to 0.8 days. The fracture gradually closes on the proppant along the fracture tip to the heel. In this process, the proppant gradually provides contact stress to support the fracture. (c) Then, the fracture is completely closed on the proppant, and the pressure inside the fracture still drops rapidly due to the absence of formation fluid supply resulting in a decrease in the fracture aperture and conductivity before hydrocarbon breakthrough. The pressure in the fracture will be lower than the breakthrough pressure during dewatering, and the formation fluid flowing into the fracture causes the pressure in the fracture to decline at a slower rate resulting in a slow decrease in the fracture aperture and conductivity as Figure 13 shows from 5 to 100 days. It is illustrated in Figure 13 that the fracture aperture associated with the fracture pore volume will decrease rapidly until the hydrocarbon breakthrough. At the initial conditions (the fracture pressure is equal to 5150 psi), the average aperture of the fracture is 0.038 ft. Without considering the leak-off, it can be calculated that the average of the fracture is 0.011 ft, at the beginning of the hydrocarbon breakthrough (the fracture pressure is approximately 4000 psi). Based on the material balance, 525 STB fracturing fluid will be produced because of the reduction of the fracture pore volume, while only 0.72 STB fracturing fluid will be produced because of the fracturing fluid volume expansion due to the fracture pressure drop. This means that the portion of flowback recovery caused by effective fracture pore volume closure as pressure falls is the main drive mechanisms before hydrocarbon breakthrough. After the hydrocarbon breakthrough, the formation fluid flux plays a more and more important role in drive mechanisms.

| Case 2: Effect of contact reference stress
The contact reference normal pressure is related to clay content, 24 and it can be set the same value for fracture hydraulic aperture and fracture mechanical aperture. When clay content is high, the corresponding contact reference stress is small, and the proppant particles are more easily embedded in the fracture surface causing a rapid decline in closed fracture conductivity and fracture pore volume. The contact reference stress varying from 5000 psi to 50 000 psi is considered.
Contact reference stress and initial fracture hydraulic fracture mainly affect the properties of fracture when they close on the proppant pack. As for the open fracture, this effect can be ignored due to the pretty large fracture aperture and fracture conductivity. The initial open fracture, therefore, is almost the same under different contact reference and different initial fracture hydraulic aperture, as shown in Figures 13  and 14. Figure 14 shows the dynamic variation of fracture aperture and conductivity under different contact reference stress. The progressive closure is not the same as case 1 when the well is produced at a constant rate. Due to the high flow rate at the bottom of the well, a non-negligible pressure gradient within the fracture causes the heel of fracture to close before the middle of the fracture. As demonstrated in Figure 14, the residual aperture and residual conductivity of the fracture at 100 days will increase with elevated contact reference stress. That is because the larger contact reference stress corresponds to a larger fracture aperture and conductivity with the same contact stress. After the breakthrough, because of formation fluid flux, the pressure in the fracture and the stress caused by deformation change gently. The difference in contact normal stress acting on the fracture surfaces under different contact reference stress is small. Hence, increasing contact reference stress, the residual fracture mechanical aperture and the conductivity loss can be reduced. Figure 15A illustrates that water flow rate declines faster for higher contact reference stress and it takes less time to dewater for higher contact reference stress. Besides, progressive closure caused a clear break in the water flow rate as Figure 15A shows. Figure 15B demonstrates that the oil production will reach a maximum value after hydrocarbon breakthrough, and then decline. The oil production decline rate with higher contact reference stress is significantly higher than the lower contact reference stress. However, the oil production at 100 days still increases significantly with elevated contact reference stress, indicating that the higher contact reference stress can cause a better production performance. This case tells us that for reservoirs with low clay content, the initial high-production wells have a better production performance, but may experience rapid production decline due to fracture closure caused by a dramatic decrease in pressure within the fracture.

| Case 3: Effect of initial hydraulic fracture aperture
Initial hydraulic fracture aperture is related to the type of proppant, the particle size of the proppant, and the sanding  Figure 16 indicates that with an increase in initial fracture hydraulic aperture, the fracture residual conductivity at 100 days is increasing. The primary reason is that with the similar contact stress acting on the fracture surface after hydrocarbon breakthrough, the fracture with larger initial hydraulic aperture has greater conductivity. Figure 17A illustrates that stimulated well with larger initial hydraulic aperture results in the faster clean-up process. And because of the sharp pressure drop within the fracture caused by the higher conductivity, the fracture with a larger initial hydraulic aperture closure is faster.
It can be seen from Figure 17B that when the initial fracture hydraulic aperture is 0.21 mm, the oil production will reach a maximum value of 28.45 STB/D and then decrease to 10.75 STB/D at 100 days, which is reduced by 62%, and when the initial fracture hydraulic aperture is 0.15 mm, the oil production increases to the peak value of 13.59 STB/D, and then declines to 8.56 STB/D, which reduced by 37%, while, when the initial fracture hydraulic aperture is 0.09 mm, the decrease after the maximum value of 3.95 STB/D is negligible. This is because the fluid pressure within the fracture with higher conductivity decreases rapidly, and the fracture conductivity also decreases rapidly leading to a dramatic decline in initial capacity until the fluid pressure drops slowly due to the balance of well production and matrix flux. This case tells us that the early productivity of the well is a comprehensive result of the dynamic closure of the fracture and the fluid supply of the matrix formation. When the fracture has a high conductivity due to effective support, improving the reservoir parameters in the SRV is beneficial to maintain the fracture conductivity along with high early productivity.

| Field case
Two multistage fracturing horizontal wells from Canada and China Changqing Oilfield were analyzed.

| Case 1
Flowback data obtained from a multifractured horizontal well completed in a tight oil reservoir are analyzed by the new coupled model. This field case was previously analyzed by Clarkson and Williams-Kovacs, Clarkson, Zhang and Clarkson and Jia et al. 4,[7][8][9]38 However, only static parameters such as fracture permeability were obtained by analyzing the flowback data. Therefore, the object to revisit this example is to capture reliable fracture dynamic parameters by historical matching the field data.
The horizontal well was fractured in 18 stages with one perforation cluster per stage, and a total of ~ 24 500 bbls of fracturing fluid were injected into the well. Microseismic monitoring data indicated that hydraulic fracturing results in mostly planar fractures with a consistent NE-SW orientation. 9 It is reasonable to assume that a multistage fracturing horizontal well can be represented as 18 rectangular stimulated zones (include hydraulic fracture and SRV) with the same properties due to the same fracturing measures, and one of them is analyzed. After hydraulic fracturing, the well was shut in 12 days, followed by flowback, and the flowback data were recorded every 15 minutes. Since the well shut-in time is long enough, it is reasonable to assume that the initial statement of fracture is closed.
The field case model is controlled by historical well bottom hole pressure, and we adjusted the fracture parameters (fracture half-length, initial fracture hydraulic aperture, and contact reference stress) to match the actual flow rate data. The input parameters are shown in Table 7. Noted that the subject well flow is always higher than the bubble point pressure, the flow in the fracture is considered to be a two-phase flow, and the gas production is equal to the oil production multiplied by the gas-oil ratio. The oil PVT property and relative permeability to oil and water are the same as the above case, and the effect of capillary force is assumed to be negligible.
The main adjustment parameters to match the flowback data are the fracture half-length, the contact reference stress, and the initial hydraulic fracture width. The historical match results of the new coupled model are shown in Figure 18. Good consistency of the three-phase production indicates that the historical match result is credible.
The derived contact reference stress is 80 000 psi, and the initial hydraulic fracture aperture is 0.0853 mm (fracture conductivity is ~167 mD ft). The relationship between fracture conductivity and contact stress can be calculated by Equations 6 and 8 with shear displacement equal to zero. Ghanizadeh et al 38 measured the unpropped/propped fracture conductivity on Montney core obtained from the same interval within the Upper Montney Formation with the subject multifractured horizontal well. The experimental data measured by Ghanizadeh are compared with the result derived by the new coupled model, as Figure 19 shows. The result indicates that when the contact stress is small, the difference between the data captured by history-matching and the experimental data is significant; however, with the increase in contact stress the difference is becoming smaller. This relationship between contact stress and conductivity, therefore, can be used to forecast well productivity.
The interpreted fracture half-length is 1010 ft, which is greater than the 450 ft derived by Clarkson and 725 ft derived by Jia et al, 4,8 but comparable with 1043 ft derived by Zhang. 10 The reason may be that the constant fracture conductivity without consideration of fracture closure is used in Clarkson and Jia's flowback model. The fracture conductivity derived by Clarkson is ~292 mD ft, and the fracture conductivity derived by Jia is ~161 mD ft. 4,8 Both are bigger than the fracture conductivity extracted by the new coupled model considering the fracture closure, as shown in Figure 19. Therefore, in order to match the water production, a larger fracture conductivity is used resulting in a shorter fracture half-length captured by the flowback model without considering fracture closure, while our result is similar to the fracture half-length derived by flowback model proposed by Zhang et al, 10 where laboratory-based petrophysical and geomechanical inputs were used to constrain the flowback history-matching process.

| Case 2
Flowback data obtained from a multifractured horizontal well completed in Changqing tight oil reservoir are also analyzed by the new coupled model. The horizontal well was fractured in 31 stages with one perforation cluster per stage. It is reasonable to assume that a multistage fracturing horizontal well can be represented as 31 rectangular stimulated zones (include hydraulic fracture and SRV) with the same properties due to the same fracturing measures, and one of them is analyzed. The field case model is controlled by historical well bottom hole pressure, and we adjusted the fracture parameters (fracture half-length, initial fracture hydraulic aperture, and contact reference stress) to match the actual flow rate data. The input parameters are shown in Table 8. The oil PVT property is shown in Figure 20. The relative permeability to oil and water is the same as the above case, and the effect of capillary force is assumed to be negligible.
The main adjustment parameters to match the flowback data are the fracture half-length, the contact reference stress, and the initial hydraulic fracture width. The historical match results of the new coupled model are shown in Figure 21. Good consistency of the three-phase production indicates that the historical match result is credible.
The interpreted fracture half-length is 420 ft. The derived contact reference stress is 50 000 psi, and the initial hydraulic fracture aperture is 0.06706 mm (fracture conductivity is ~82 mD ft). The relationship between fracture conductivity and contact stress can be calculated with shear displacement equal to zero. The relationship between contact stress and conductivity is shown in Figure 22.

| DISCUSSION
In this work, a hybrid fracture model coupling flow and geomechanics is developed for flowback and early-time  Oil viscosity Oil FVF production simulation. The dynamic behaviors of fracture during the flowback can be described in detail by this new coupled model. Through history-matching the field flowback data, it can extract the key parameters of fracture and the relationship between fracture conductivity and contact stress, which can be used to forecast the well productivity during the production period. However, there are several improvements that can be made, as discussed in this section.
First, due to the use of the nonlinear relationship between fracture aperture and contact stress, the type of proppant, the size of proppant particles, the proppant concentration, and the formation mechanical properties are upscaled into three parameters (w c c,Wref , and wdil ). Therefore, it is difficult to describe the effect of each proppant parameters on well early production. Different models to calculate the fracture aperture and fracture conductivity under various contact stress have been proposed. 31,39,40 Hence, a more accurate model to calculate the fracture aperture and fracture conductivity will be combined with our new coupled model to evaluate the effect of different proppant parameters on early production and to guide the fracturing treatment in our future work.
Second, it is demonstrated by Zhang that the initial inputs are generally unknown at the start of the flowback, creating significant uncertainty in the quantitative analysis of flowback data. Therefore, in order to recreate initial condition, a new coupled model considering the fracturing fluid leak-off for the shut-in and flowback period will be proposed in our future work. 10 The initial fracture half-length and the fluid  pressure within the fracture at the end of hydraulic fracturing will be set as the initial condition of the new coupled model.

| CONCLUSIONS
In this paper, a flow-geomechanics coupled model considering the fracture closure during flowback period is presented. From the result of this work, the following important conclusion can be obtained: 1. The two-phase flow model is fully coupled with the stress model by the relationships between fracture aperture and fracture flow parameters (fracture pore volume and fracture conductivity). The dynamic fracture behaviors during the flowback and early-time production, therefore, can be comprehensively described by the coupled model. 2. The fracture closure may exhibit three stages. (1) The fracture, only supported by high-pressure fluid within it, is in the open state with a large mechanical aperture and conductivity. (2) The fracture gradually closes on the proppant pack during fracture pressure reduction caused by dewatering. And the proppant gradually provides contact stress to support the fracture. (3) The fracture is completely closed on the proppant pack. The fracture aperture is a relationship of contact stress acting on the fracture surfaces. 3. Larger contact reference stress and initial hydraulic fracture aperture can lead to a better production performance during flowback and early production. A stimulated well with a large initial fracture aperture will have a better production performance at the beginning of the flowback; however, it will decline rapidly until the flux from the matrix can effectively supply the energy within the fracture. Therefore, it is more important to effectively increase the permeability of the SRV than to only create a fracture with high conductivity.
4. The fracture effective length and the relationship between fracture conductivity and contact stress can be captured to forecast production by history-matching the field flowback data.

Fracture f low model
We bring the second-order difference for Each spatial derivative term and the first-order backward difference for the time derivative term into Equations 12 and 13, to obtain the implicit finite difference form of the oil and water flow equation: In Equations A.1 and A.2, the terms (q n+1 SCWoi , q n+1 SCWwi P ′ cowi,j , and i ) and the coefficients (T and C of the pressure and saturation terms) are defined as follows: where T i,j is the conductivity between adjacent segments i and j. C is a coefficient used to indicate that the fracture volume changes with pressure and saturation. a i is the half-length of the fracture segment i. c1 is the unit conversion factor.
In Equations A.3 and A.4, G i,j is the geometric transmissibility between segments i and j, calculated using the harmonic average method.
where i is the transmissibility of fracture segment i: and, w i is the hydraulic aperture of fracture segment i. c1 is the unit conversion factor.
The horizontal well production equation is given by 17 : where W is a collection of all segments connected to the wellbore. T n+1 oi,W and T n+1 wi,W are the oil and water transmissibility between wellbore and fracture segment i at n + 1 time step, respectively. p W is wellbore pressure. T n+1 oi,W and T n+1 wi,W are the oil and water transmissibility between wellbore and fracture segment i at time level n + 1, respectively. We assume that the effective fracturing aperture of the wellbore is infinite 24 ; therefore, the conductivity between the wellbore and the segment i can be expressed as: where the submatrix A i,j is a 2 × 2 matrix, which can be expressed as: The term b i is given by:

Matrix f low model
This method has been used to describe the flow of fluid in the matrix during the flowback. 2 The planar linear flow equation for infinite boundary with constant production in each isolated region can be expressed as follows: The analytical solution proposed by Wattenbarger et al 32 and the superposition principle can be used to solve the above equations under a variable-rate inner boundary to obtain the Equation 19. 29 Let