Numerical Study on the E ﬀ ects of Imbibition on Gas Production and Shut-In Time Optimization in Woodford Shale Formation

: In shale gas formations, imbibition is signiﬁcant since the tight pore structure causes a strong capillary suction pressure. After hydraulic fracturing, imbibition during the period of shut-in a ﬀ ects the water recovery of ﬂowback. Although there have been many studies investigating imbibition in shale formations, few papers have studied the relationship between gas production and shut-in time under the inﬂuence of imbibition. This paper developed a numerical model to investigate the e ﬀ ect of imbibition on gas production to optimize the shut-in time after hydraulic fracturing. This numerical model is a 2-D two-phase (gas and water) imbibition model for simulating an imbibed ﬂuid ﬂow and its e ﬀ ect on permeability, ﬂowback, and water recovery. The experimental and ﬁeld data from the Woodford shale formation were matched by the model to properly conﬁgure and calibrate the model parameters. The experimental data consisted of the relationship between the imbibed ﬂuid volume and permeability change, the relative permeability, and the capillary pressure for the Woodford shale samples. The Woodford ﬁeld data included the gas production and ﬂowback volume. The modeling results indicate that imbibition can be a beneﬁcial factor for gas production, since it can increase rock permeability. However, the gas production would be reduced when excessive ﬂuid is imbibed by the shale matrix. Therefore, the shut-in time after hydraulic fracturing, when the imbibition happens in shale, could be optimized to maximize the gas production. Author Contributions: Conceptualization, Z.Z.; Methodology, Z.Z. and S.W.; Formal analysis, S.W.; Data Curation, S.W.; Investigation, Z.Z. and S.W.; Resources, R.L.; Validation, R.L. and X.L.; Writing–Original Draft, Z.Z. and S.W. research Natural Science Foundation China (Grant


Introduction
Hydraulic fracturing is a method to improve gas production in shale formations. Well shut-in is one of the necessary steps during the hydraulic fracturing operation to calm down the formation and prepare the flow back. According to field data, the shut-in time is related to the production of shale gas formations. Almulhim et al. found that the gas production peak and cumulative gas production could be higher when the shut-in time was longer in some shale gas formations, which had strong water-wet characteristics [1]. Bertoncello et al. also indicated that up to two months of shut-in could bring a better performance in the cumulative production of the shale gas formation compared with the most current shut-in times, which were very short in the field [2]. The long shut-in time leads to a small water recovery. The common opinion is that as much water recovery as possible is important to prevent formation damage. However, based on the field reports in shale formations from Alkouh and Wattenbarger, less than 50% of the injection fluid was recovered when the shut-in took longer to achieve a better production [3].
The reason for the small water recovery in shale gas formations is that a longer shut-in time can result in fluid flow from fracture faces to deep inside the rock matrix [4]. This fluid flow is called imbibition. Imbibition is defined as a flow of immiscible fluids in a porous rock. It can cause a larger stimulation reservoir volume in shale formations with complex natural fractures so that the percentage of fracture fluid recovery is small and takes a much longer time to complete flow back [5]. During flowing, most of the remaining fluid is imbibed by the rock matrix. Hence, imbibition during the shut-in is considered to result in a small amount of water recovery [1,[6][7][8].
Roychaudhuri et al., Makhanov et al., and Zhou et al., have proved that a large volume of the fracturing fluid can be imbibed by the shale samples [9][10][11][12]. However, few studies have investigated the relationship between imbibition, shut-in time, water recovery, and gas production. Therefore, this paper developed a numerical model to simulate an imbibed fluid flow and its effect on permeability, production, and water recovery. Additionally, the model could be applied to optimize the shut-in time based on this relationship to achieve the best gas production performance. The highlight of this study is to quantitatively combine the imbibition with the gas-water flow model.
In this paper, the experimental and field data of the Woodford shale formation were considered by the model to study the effect of imbibition on gas production and shut-in time optimization. The data included the imbibition volume, rock permeability change, relative permeability, capillary pressure, daily flowback rate and volume, and gas production, which were matched by the model to achieve the parameters for simulating. The contribution of this model is to understand the relationship of the imbibition to the gas production in shale formations and close the gap between experimental research and field operations.

Formation Characterization
The field data used in this study are from the Cana Field of western Oklahoma. The field is located in the Anadarko Basin with the initial development program started in Canadian County, as shown in the map of Figure 1. Cana Field produces from the Upper Devonian Woodford Shale, and the production is mainly composed of gas, condensate, and natural gas liquids. The stratigraphic chart of the target formation is shown in Table 1.
XRD testing is summarized in Table 2 for the Woodford shale samples. It indicates that the major mineral is quartz and the second is illite. Hence, the formation is clay-rich.
Other geological and rock properties are shown in Table 3. Cana Field produces from the Upper Devonian Woodford Shale, and the production is mainly composed of gas, condensate, and natural gas liquids. The stratigraphic chart of the target formation is shown in Table 1. XRD testing is summarized in Table 2 for the Woodford shale samples. It indicates that the major mineral is quartz and the second is illite. Hence, the formation is clay-rich. Other geological and rock properties are shown in Table 3.

Numerical Model
In this paper, we introduce a numerical model which considers imbibition to simulate the fluid flow from porous media into hydraulic fractures. The experimental and field data from the Woodford Shale formation are collected [15][16][17]. The data are used to show the relationship between imbibition, permeability, capillary pressure, and complex transport mechanisms. There are two main assumptions for the model when studying the shut-in time. One is that the length of the hydraulic fracture is assumed constant after shut-in. Although some further fracture propagation can still take place, there is not a sufficient amount of proppants to support those lengths to make a significant contribution to the flow. Another assumption is that the fluid pressure in hydraulic fractures keeps constant. This assumption is because the permeability in shale formation is too small, so it is difficult to impact the pressure in fractures by formation. In addition, it neglects some effects from the proppant and gel filter cake on the fracture conductivity during the shut-in period.

Fundamental Equations
The mass conservation equations for gas and fracturing fluids in porous media are: where φ is the rock porosity. ρ w and ρ g are, respectively, the water density and gas density. s w and s g are, respectively, the water and gas saturation. K is the absolute permeability. k rw and k rg are, respectively, the phase permeability of water and gas. k is the apparent permeability of gas, which is the enhanced gas mobility considering different transport mechanisms [18]. The water phase is assumed to follow Darcy's law, thus its apparent permeability is 1. µ w and µ g are, respectively, the water and gas viscosity. p w and p g are, respectively, the water and gas pressure. p c is the capillary pressure, which is the function of water saturation, and the specific function of capillary pressure can be obtained through mercury injection experiments. The detailed solving procedure is included in the Appendix A.
The equation of state for the real gas can be expressed as: where R is the universal gas constant. T is the formation temperature. C is the gas concentration. Z is the ratio of gas volume in the real condition to the gas volume in the ideal condition, and it can be expressed as [19,20]: The reduced density can be defined as: where T cr and p cr are, respectively, the critical temperature and critical pressure, and T r = T/T cr stands for the reduced temperature. The eleven coefficients of Equation (6) are given below [20]:

Imbibition Effect during the Shut-In Period
During the shut-in period, the fracturing fluid has two aspects of influence on the shale matrix. Its first is to flow into the matrix to reduce the relative gas permeability. Second is that the imbibed fracturing fluid induces micro-fractures in the shale matrix to increase the absolute permeability of shale reservoirs. Those two influences were proved in the experiment [17]. The model in this paper can combine the influences to optimize the shut-in time in shale gas formations. The gas-water relative permeability and capillary curves are essential to the calculation of the gas-water flow in porous mediums. The experimental data of gas-water relative permeability and capillary pressure can usually be fitted with the power functions given below [21].
where s e is the effective water saturation, and s e = (s w − s wr )/(1 − s wr − s or ). λ is related to the pore size distribution [22]. The original data of the relative permeability curves and capillary curves are shown in Figures 2 and 3. They were originally from the experimental results from Charoenwongsa (2011Charoenwongsa ( , 2013 and Lu (2014) [15,16,23]. The shale cores in the imbibition experiment and relative permeability experiment were also from the Woodford Shale formation. Due to the fact that the Brooks-Corey relation originated from the experimental data of sand cores, these data from shale cores cannot be well fitted by Equations (9) and (10). The specific power functions used to fit these data are given in Equations (11)- (13), and the value of the determination coefficients of Equations (11)-(13) are all above 0.98, which shows the perfect agreement between these fitting equations and the experimental data [24].
where dry K is the permeability of the dry shale samples. K is the permeability of the wet shale samples. Because Equation (14) is obtained by fitting the experimental data and is optimized with minimum error, K is not equal to dry K when the water saturation is 0, which contradicts the physical meaning of the value. However, when the initial water saturation of the shale gas reservoir is 0.46, and the irreducible water saturation is 0.45, the unphysical value of / dry K K will not appear.
where dry K is the permeability of the dry shale samples. K is the permeability of the wet shale samples. Because Equation (14) is obtained by fitting the experimental data and is optimized with minimum error, K is not equal to dry K when the water saturation is 0, which contradicts the physical meaning of the value. However, when the initial water saturation of the shale gas reservoir is 0.46, and the irreducible water saturation is 0.45, the unphysical value of / dry K K will not appear.

Imbibition Effect
The traditional opinion is that the shut-in time during hydraulic fracturing should be as short as possible to prevent permeability decline due to water blockage. In shale formations, the imbibed fluids can induce microfractures, which improve the absolute matrix permeability. According to the experiment, the shale permeability can increase by 2 orders of magnitude in the imbibition experiments [17]. The overall permeability is proposed to describe the absolute permeability of the wet shale sample. The pressure build-up method was used for the wet shale sample permeability measurement. This is one of the most efficient methods to measure the permeability of tight rocks [12,17].
This paper considers the effect of imbibition on permeability based on Zhou's experiments (2016). The experimental data of the imbibed water saturation and overall permeability are given in Figure 4, and the data are from Zhou et al. [17]. Enhanced permeability is the permeability ratio of wet shale samples to dry shale samples. Figure 4 indicates that the overall permeability of the shale sample increases when the imbibition starts. When the imbibition saturation reaches 0.45, the enhanced permeability of the shale sample reaches up to 36.2. The imbibition experimental data can be well fitted with the exponential function: where K dry is the permeability of the dry shale samples. K is the permeability of the wet shale samples. Because Equation (14) is obtained by fitting the experimental data and is optimized with minimum error, K is not equal to K dry when the water saturation is 0, which contradicts the physical meaning of the value. However, when the initial water saturation of the shale gas reservoir is 0.46, and the irreducible water saturation is 0.45, the unphysical value of K/K dry will not appear.
physical meaning of the value. However, when the initial water saturation of the shale gas reservoir is 0.46, and the irreducible water saturation is 0.45, the unphysical value of / dry K K will not appear.     [17]. K dry is the permeability of the dry shale core.

Shut-In Time
Shut-in time influences the imbibition region, as shown in Figure 5. When the shut-in time is short, the imbibition region is small, so that in matrix gas relative permeability is not affected by water invasion. At the same time, based on the experiment, the imbibition would induce microfractures to increase the permeability of the formation (see Figure 5a). When the shut-in time is longer, the imbibition region is bigger in the formation to reduce the gas relative permeability in the rock matrix (see Figure 5b). Meanwhile, the induced microfractures would not change too much. Hence, that is the reason that it is necessary to optimize the shut-in time. Figure 6 shows the relationship between the enhanced gas permeability (i.e., k·k rg /k 0 ) and water saturation in the shale matrix, and the enhanced gas permeability is 1 in the initial water saturation. It can be seen from Figure 6 that when we consider the influence Energies 2020, 13, 3222 7 of 17 of imbibition, the higher the water saturation, the larger the gas permeability. The enhanced gas permeability goes up sharply when the water saturation is over 0.8. gas permeability is 1 in the initial water saturation. It can be seen from Figure 6 that when we consider the influence of imbibition, the higher the water saturation, the larger the gas permeability. The enhanced gas permeability goes up sharply when the water saturation is over 0.8.
During the shut-in time, micro fractures do not occur until the shale matrix has been in contact with the fluids for some time, over 1 day in Zhou's experiments (2016a). Therefore, the enhanced gas permeability is given at the end of the shut-in time. Once the microfractures appear, the matrix permeability will not decrease with the decline of water saturation. Thus, the process of permeability enhancement is not reversible in the production process. As shown in Figures 5 and 6, coupling the imbibition effect with the gas-water flow model optimizes the shut-in time in shale gas reservoirs.

Complex Transport Mechanisms
Due to the co-existence of nano-pores and micro-fractures in shale, the shale gas flow is controlled by different transport mechanisms. Based on the ratio of the molecular-free path to the pore characteristic length, the Knudsen number is defined to determine the flow mechanism in pores of different scales [25]. The range of Knudsen numbers and the corresponding flow mechanisms are shown in Figure 7. the influence of imbibition, the higher the water saturation, the larger the gas permeability. The enhanced gas permeability goes up sharply when the water saturation is over 0.8.
During the shut-in time, micro fractures do not occur until the shale matrix has been in contact with the fluids for some time, over 1 day in Zhou's experiments (2016a). Therefore, the enhanced gas permeability is given at the end of the shut-in time. Once the microfractures appear, the matrix permeability will not decrease with the decline of water saturation. Thus, the process of permeability enhancement is not reversible in the production process. As shown in Figures 5 and 6, coupling the imbibition effect with the gas-water flow model optimizes the shut-in time in shale gas reservoirs.

Complex Transport Mechanisms
Due to the co-existence of nano-pores and micro-fractures in shale, the shale gas flow is controlled by different transport mechanisms. Based on the ratio of the molecular-free path to the pore characteristic length, the Knudsen number is defined to determine the flow mechanism in pores of different scales [25]. The range of Knudsen numbers and the corresponding flow mechanisms are shown in Figure 7. During the shut-in time, micro fractures do not occur until the shale matrix has been in contact with the fluids for some time, over 1 day in Zhou's experiments (2016a). Therefore, the enhanced gas permeability is given at the end of the shut-in time. Once the microfractures appear, the matrix permeability will not decrease with the decline of water saturation. Thus, the process of permeability enhancement is not reversible in the production process. As shown in Figures 5 and 6, coupling the imbibition effect with the gas-water flow model optimizes the shut-in time in shale gas reservoirs.

Complex Transport Mechanisms
Due to the co-existence of nano-pores and micro-fractures in shale, the shale gas flow is controlled by different transport mechanisms. Based on the ratio of the molecular-free path to the pore characteristic length, the Knudsen number is defined to determine the flow mechanism in pores of different scales [25]. The range of Knudsen numbers and the corresponding flow mechanisms are shown in Figure 7.
where λ and r are respectively the mean-free path of molecules and the pore characteristic length. Figure 7 gives the range of Knudsen numbers and the corresponding flow mechanisms.
Energies 2020, 13, 3222 8 of 17 where λ and r are respectively the mean-free path of molecules and the pore characteristic length.
where B K is the Boltzmann constant. c λ is the collision diameter of gas molecules. ζ is the dimensionless rarefaction coefficient of gas, which is also the function of the Knudsen number [26]: where A, B, and 0 ζ are constants. In this paper, A = 0.17, B = 0.4348 [26].
The influence of imbibition, capillary pressure, and complex transport mechanisms on shale gas flow are considered through Equations (12), (14) and (17).

Field Data Matching
The gas production and flowback rate of the Woodford shale formation are used for history data matching. The field data and initial input parameters of the model are from Lu [16]. The input parameters are summarized in Table 4. The simulation setup is shown in Figure 8. The history data matching results for the well are shown in Figures 9 and 10.
where K B is the Boltzmann constant. λ c is the collision diameter of gas molecules. ζ is the dimensionless rarefaction coefficient of gas, which is also the function of the Knudsen number [26]: where A, B, and ζ 0 are constants. In this paper, A = 0.17, B = 0.4348 [26]. The influence of imbibition, capillary pressure, and complex transport mechanisms on shale gas flow are considered through Equations (12), (14) and (17).

Field Data Matching
The gas production and flowback rate of the Woodford shale formation are used for history data matching. The field data and initial input parameters of the model are from Lu [16]. The input parameters are summarized in Table 4. The simulation setup is shown in Figure 8. The history data matching results for the well are shown in Figures 9 and 10. The geometry of the simulated reservoir is shown in Figure 8. The hydraulic fracture is assumed as a rectangle and it fully penetrates the height of the reservoir. If the gravity effect is neglected, the reservoir is simplified to the 2D case. The geometry of the simulated reservoir is shown in Figure 8. The hydraulic fracture is assumed as a rectangle and it fully penetrates the height of the reservoir. If the gravity effect is neglected, the reservoir is simplified to the 2D case. In the injection process, two boundary conditions are given at the hydraulic fracture: where w Q is the injection rate of the fracturing fluid. f H and f L are, respectively, the height and the length of the hydraulic fracture. F Γ stands for the hydraulic fracture.
During the shut-in time, no fluids flow out from the formation. The bottom hole pressure changes as the time goes on. The "no flow" boundary allows no fluids to enter into the formation from the fracture surface. The pressure change in the fracture surface can be simulated when the "no flow" boundary is applied to the fracture surface. The following boundary condition is applied to the hydraulic fracture: During the production process, the following boundary condition is applied to the hydraulic fracture: The other outer boundaries are still set as nonflow boundaries, i.e., Equation (21). Before carrying out the data matching, there are four steps. First, the geometry of the stimulated reservoir is established in Figure 8, and then it should be meshed. Quadrilateral grids are adopted in this paper. Second, the basic parameters listed in Table. 4, including the relative permeability functions (Equations (11)-(12)), capillary pressure function (Equation (13)), and the imbibition function (Equation (14)) are inputted into the program. Third, three analysis steps are set in the program-i.e., injecting In the injection process, two boundary conditions are given at the hydraulic fracture: where Q w is the injection rate of the fracturing fluid. H f and L f are, respectively, the height and the length of the hydraulic fracture. Γ F stands for the hydraulic fracture. During the shut-in time, no fluids flow out from the formation. The bottom hole pressure changes as the time goes on. The "no flow" boundary allows no fluids to enter into the formation from the fracture surface. The pressure change in the fracture surface can be simulated when the "no flow" boundary is applied to the fracture surface. The following boundary condition is applied to the hydraulic fracture: During the production process, the following boundary condition is applied to the hydraulic fracture: The other outer boundaries are still set as nonflow boundaries, i.e., Equation (21). Before carrying out the data matching, there are four steps. First, the geometry of the stimulated reservoir is established in Figure 8, and then it should be meshed. Quadrilateral grids are adopted in this paper. Second, the basic parameters listed in Table 4, including the relative permeability functions (Equations (11)-(12)), capillary pressure function (Equation (13)), and the imbibition function (Equation (14)) are inputted into the program. Third, three analysis steps are set in the program-i.e., injecting fluids first, then shutting the well, and finally production-and the conditions are given in Equations (19)- (22). Finally, the timestep is set to 1 h before starting the solution program.
As can be seen from Figures 9 and 10, the simulation results have a good agreement with the field data. In Figure 9, the sudden drop in field gas rate is also predicted by the simulation results.
Energies 2020, 13, x FOR PEER REVIEW 10 of 18 fluids first, then shutting the well, and finally production-and the conditions are given in Equations (19)- (22). Finally, the timestep is set to 1 h before starting the solution program. As can be seen from Figures 9 and 10, the simulation results have a good agreement with the field data. In Figure 9, the sudden drop in field gas rate is also predicted by the simulation results. Figure 9. A comparison between the field data and the simulated gas production. Figure 10. A comparison between the field data and the simulated water production.

Results and Discussion
In this section, the simulation result was discussed as the effect of imbibition on the production and shut-in time optimization.

Imbibition Effect on Gas Production
Based on the experimental data from Woodford shale samples [17], it indicated that the initial water saturation was irreducible, which meant that the original saturated region was a single gas flow region. After fracturing, the near hydraulic fracture region creates a gas-water flow. The gas-water flow region lowers the gas production [18]. Figure 11 shows the influence of imbibition on gas production. The shut-in time is 1 day in Figures 11-13. It can be seen from Figure 11 that the imbibition effect improves the gas production and shortens the time of obtaining the gas flow. When we consider the influence of imbibition, the accumulated gas production is 2.53 × 10 8 m 3 , which is 12.68 percent higher than that neglecting the influence of imbibition. The scenario without considering the imbibition effect is calculated using the published research of Wei et al. [18]. fluids first, then shutting the well, and finally production-and the conditions are given in Equations (19)- (22). Finally, the timestep is set to 1 h before starting the solution program.
As can be seen from Figures 9 and 10, the simulation results have a good agreement with the field data. In Figure 9, the sudden drop in field gas rate is also predicted by the simulation results. Figure 9. A comparison between the field data and the simulated gas production. Figure 10. A comparison between the field data and the simulated water production.

Results and Discussion
In this section, the simulation result was discussed as the effect of imbibition on the production and shut-in time optimization.

Imbibition Effect on Gas Production
Based on the experimental data from Woodford shale samples [17], it indicated that the initial water saturation was irreducible, which meant that the original saturated region was a single gas flow region. After fracturing, the near hydraulic fracture region creates a gas-water flow. The gas-water flow region lowers the gas production [18]. Figure 11 shows the influence of imbibition on gas production. The shut-in time is 1 day in Figures 11-13. It can be seen from Figure 11 that the imbibition effect improves the gas production and shortens the time of obtaining the gas flow. When we consider the influence of imbibition, the accumulated gas production is 2.53 × 10 8 m 3 , which is 12.68 percent higher than that neglecting the influence of imbibition. The scenario without considering the imbibition effect is calculated using the published research of Wei et al. [18].

Results and Discussion
In this section, the simulation result was discussed as the effect of imbibition on the production and shut-in time optimization.

Imbibition Effect on Gas Production
Based on the experimental data from Woodford shale samples [17], it indicated that the initial water saturation was irreducible, which meant that the original saturated region was a single gas flow region. After fracturing, the near hydraulic fracture region creates a gas-water flow. The gas-water flow region lowers the gas production [18]. Figure 11 shows the influence of imbibition on gas production. The shut-in time is 1 day in Figures 11-13. It can be seen from Figure 11 that the imbibition effect improves the gas production and shortens the time of obtaining the gas flow. When we consider the influence of imbibition, the accumulated gas production is 2.53 × 10 8 m 3 , which is 12.68 percent higher than that neglecting the influence of imbibition. The scenario without considering the imbibition effect is calculated using the published research of Wei et al. [18]. shale core. Figure 13 shows the distribution of the dimensionless gas permeability of shale-i.e., 0 / rg kk k , where imbibition is neglected. Figure 12 indicates that the dimensionless gas permeability in the vicinity of hydraulic fracture is up to 8 when the imbibition happens. The enhanced gas permeability near hydraulic fractures is the reason for improving the gas production. However, when neglecting the imbibition effect, the dimensionless gas permeability near the hydraulic fracture is as low as 0.0025. The damaged gas permeability in the vicinity of the hydraulic fractures is responsible for the decline in gas production [18]. Comparing the region size of the enhanced gas permeability in Figure 12 and damaged gas permeability in Figure 13, we know that the effective distance of the capillary pressure-driven imbibition is much shorter than the pressure difference-driven water flow. . The shut-in time is 1 day. Figure 11. The influence of imbibition on gas production; the shut-in time is 1 day. Figure 11. The influence of imbibition on gas production; the shut-in time is 1 day. Figure 12 shows the distributions of the dimensionless permeability of shale considering the imbibition of the fracturing fluid-i.e., 0 / rg kk k , where 0 k is the initial matrix permeability of the dry shale core. Figure 13 shows the distribution of the dimensionless gas permeability of shale-i.e., 0 / rg kk k , where imbibition is neglected. Figure 12 indicates that the dimensionless gas permeability in the vicinity of hydraulic fracture is up to 8 when the imbibition happens. The enhanced gas permeability near hydraulic fractures is the reason for improving the gas production. However, when neglecting the imbibition effect, the dimensionless gas permeability near the hydraulic fracture is as low as 0.0025. The damaged gas permeability in the vicinity of the hydraulic fractures is responsible for the decline in gas production [18]. Comparing the region size of the enhanced gas permeability in Figure 12 and damaged gas permeability in Figure 13, we know that the effective distance of the capillary pressure-driven imbibition is much shorter than the pressure difference-driven water flow.   Figures 14 and 15 show the distributions of pore pressure (water phase pressure) and water saturation after shutting the well for various days. It can be seen from Figure 14 that the longer the shutin time, the larger the pressurized region is but the lower the BHP (Bottom Hole Pressure) is. Figure 15 points out that if the shut-in time is longer, the imbibed region becomes larger. When the well is only shut in for 1 day, the fracturing fluid is concentrated near the hydraulic fracture, and the water intrudes into the shale matrix 8.4 m after being shut in for 1 day. For the scenarios shut-in 3 days and shut-in 5 days, the water intrudes into the shale matrix 23.5 and 32 m, respectively. Thus, the flowback rate of the fracturing fluid is fast. After 3 days of shutting well, more fracturing fluid is imbibed into the shale matrix. This could cause more inner gas to be displaced by the fracturing fluid to obtain the initial gas production. Figure 13. The distribution of dimensionless gas permeability without considering the imbibition effect-i.e., k·k rg /k 0 . The shut-in time is 1 day. Figure 12 shows the distributions of the dimensionless permeability of shale considering the imbibition of the fracturing fluid-i.e., kk rg /k 0 , where k 0 is the initial matrix permeability of the dry shale core. Figure 13 shows the distribution of the dimensionless gas permeability of shale-i.e., kk rg /k 0 , where imbibition is neglected. Figure 12 indicates that the dimensionless gas permeability in the vicinity of hydraulic fracture is up to 8 when the imbibition happens. The enhanced gas permeability near hydraulic fractures is the reason for improving the gas production. However, when neglecting the imbibition effect, the dimensionless gas permeability near the hydraulic fracture is as low as 0.0025. The damaged gas permeability in the vicinity of the hydraulic fractures is responsible for the decline in gas production [18]. Comparing the region size of the enhanced gas permeability in Figure 12 and damaged gas permeability in Figure 13, we know that the effective distance of the capillary pressure-driven imbibition is much shorter than the pressure difference-driven water flow. Figures 14 and 15 show the distributions of pore pressure (water phase pressure) and water saturation after shutting the well for various days. It can be seen from Figure 14 that the longer the shut-in time, the larger the pressurized region is but the lower the BHP (Bottom Hole Pressure) is. Figure 15 points out that if the shut-in time is longer, the imbibed region becomes larger. When the well is only shut in for 1 day, the fracturing fluid is concentrated near the hydraulic fracture, and the water intrudes into the shale matrix 8.4 m after being shut in for 1 day. For the scenarios shut-in 3 days and shut-in 5 days, the water intrudes into the shale matrix 23.5 and 32 m, respectively. Thus, the flowback rate of the fracturing fluid is fast. After 3 days of shutting well, more fracturing fluid is imbibed into the shale matrix. This could cause more inner gas to be displaced by the fracturing fluid to obtain the initial gas production.  Figures 14 and 15 show the distributions of pore pressure (water phase pressure) and water saturation after shutting the well for various days. It can be seen from Figure 14 that the longer the shutin time, the larger the pressurized region is but the lower the BHP (Bottom Hole Pressure) is. Figure 15 points out that if the shut-in time is longer, the imbibed region becomes larger. When the well is only shut in for 1 day, the fracturing fluid is concentrated near the hydraulic fracture, and the water intrudes into the shale matrix 8.4 m after being shut in for 1 day. For the scenarios shut-in 3 days and shut-in 5 days, the water intrudes into the shale matrix 23.5 and 32 m, respectively. Thus, the flowback rate of the fracturing fluid is fast. After 3 days of shutting well, more fracturing fluid is imbibed into the shale matrix. This could cause more inner gas to be displaced by the fracturing fluid to obtain the initial gas production.    Figures 14 and 15 show the distributions of pore pressure (water phase pressure) and water saturation after shutting the well for various days. It can be seen from Figure 14 that the longer the shutin time, the larger the pressurized region is but the lower the BHP (Bottom Hole Pressure) is. Figure 15 points out that if the shut-in time is longer, the imbibed region becomes larger. When the well is only shut in for 1 day, the fracturing fluid is concentrated near the hydraulic fracture, and the water intrudes into the shale matrix 8.4 m after being shut in for 1 day. For the scenarios shut-in 3 days and shut-in 5 days, the water intrudes into the shale matrix 23.5 and 32 m, respectively. Thus, the flowback rate of the fracturing fluid is fast. After 3 days of shutting well, more fracturing fluid is imbibed into the shale matrix. This could cause more inner gas to be displaced by the fracturing fluid to obtain the initial gas production.

Shut-In Time Optimization
The fracturing fluid is pumped into the reservoir through hydraulic fractures and then imbibed into the shale matrix. When the well is shut after hydraulic fracturing, two effects of the fracturing fluid exist: (1) the fracturing fluid enters the matrix pores, thus the gas relative permeability declines; (2) the microfractures were clearly observed after the imbibition experiment in shale rocks, which had no observed fractures before the experiment. In addition, the permeability testing indicated that there was a big improvement in permeability after imbibition. In the testing, the permeability of shale after imbibition can be up to 153.57 times higher than the original permeability [17]. The gas production and flowback volume in the first 100 days are analyzed in the flowback optimization. The accumulated gas production after 100 days is chosen as the selection index of optimal shut-in time. Five scenarios of shut-in time are analyzed in this section; for more obvious contrast, only three scenarios are shown in Figures 16 and 17. In Figure 16, it is indicated that the accumulated gas production of the shut-in 3 days scenario is the largest, which means that the optimal shut-in time is 3 days. When the well is shut for 1 day, the gas flow is obtained from the reservoir after 7 days flowback of the fracturing fluid. When the well is shut for 3 days and 5 days, the initial gas productions are, respectively, 2.955 × 10 6 and 2.792 × 10 6 m 3 .
The fracturing fluid is pumped into the reservoir through hydraulic fractures and then imbibed into the shale matrix. When the well is shut after hydraulic fracturing, two effects of the fracturing fluid exist: (1) the fracturing fluid enters the matrix pores, thus the gas relative permeability declines; (2) the microfractures were clearly observed after the imbibition experiment in shale rocks, which had no observed fractures before the experiment. In addition, the permeability testing indicated that there was a big improvement in permeability after imbibition. In the testing, the permeability of shale after imbibition can be up to 153.57 times higher than the original permeability [17].
The gas production and flowback volume in the first 100 days are analyzed in the flowback optimization. The accumulated gas production after 100 days is chosen as the selection index of optimal shut-in time. Five scenarios of shut-in time are analyzed in this section; for more obvious contrast, only three scenarios are shown in Figures 16 and 17. In Figure 16, it is indicated that the accumulated gas production of the shut-in 3 days scenario is the largest, which means that the optimal shut-in time is 3 days. When the well is shut for 1 day, the gas flow is obtained from the reservoir after 7 days flowback of the fracturing fluid. When the well is shut for 3 days and 5 days, the initial gas productions are, respectively, 2.955 × 10 6 and 2.792 × 10 6 m 3 .   into the shale matrix. When the well is shut after hydraulic fracturing, two effects of the fracturing fluid exist: (1) the fracturing fluid enters the matrix pores, thus the gas relative permeability declines; (2) the microfractures were clearly observed after the imbibition experiment in shale rocks, which had no observed fractures before the experiment. In addition, the permeability testing indicated that there was a big improvement in permeability after imbibition. In the testing, the permeability of shale after imbibition can be up to 153.57 times higher than the original permeability [17].
The gas production and flowback volume in the first 100 days are analyzed in the flowback optimization. The accumulated gas production after 100 days is chosen as the selection index of optimal shut-in time. Five scenarios of shut-in time are analyzed in this section; for more obvious contrast, only three scenarios are shown in Figures 16 and 17. In Figure 16, it is indicated that the accumulated gas production of the shut-in 3 days scenario is the largest, which means that the optimal shut-in time is 3 days. When the well is shut for 1 day, the gas flow is obtained from the reservoir after 7 days flowback of the fracturing fluid. When the well is shut for 3 days and 5 days, the initial gas productions are, respectively, 2.955 × 10 6 and 2.792 × 10 6 m 3 .   It can be seen from Figure 17 that the shut-in time has little influence on the flowback water volume. Therefore, we will focus on the relationship between the shut-in time and gas production.
The shut-in time and corresponding accumulated gas production after 100 days are plotted in Figure 18. Due to the two opposite influences of fracturing fluid on the shale matrix (i.e., with the increase in water saturation, imbibition increases the absolute permeability of the shale matrix; however, the gas relative permeability declines), the accumulated gas production is not monotonically related to the shut-in time. The largest accumulated gas production can be achieved after shutting the well for 3 days.
The shut-in time and corresponding accumulated gas production after 100 days are plotted in Figure 18. Due to the two opposite influences of fracturing fluid on the shale matrix (i.e., with the increase in water saturation, imbibition increases the absolute permeability of the shale matrix; however, the gas relative permeability declines), the accumulated gas production is not monotonically related to the shut-in time. The largest accumulated gas production can be achieved after shutting the well for 3 days. Figure 18. The relation between 100 days accumulated gas production and the shut-in time.

Conclusions
This paper developed a numerical model coupling the imbibition of fracturing fluid with the gaswater flow model. This model is used to investigate the effect of imbibition on Woodford shale. The experimental and field data of Woodford shale formation were matched by the model to study the effect of imbibition on gas production and shut-in time optimization. The numerical model enables us to efficiently optimize the shut-in time to achieve a high engineering profit. Taking the Woodford Shale as an example, we contrasted the effect of imbibition on gas production and optimized the best shut-in time. The conclusions are the following.
Imbibition is shown to cause a matrix permeability decrease, while the unstable natural fracture permeability increases in the experiment. When the numerical model was developed, the relationship between the imbibed fluid volume and permeability change was considered and matched. Hence, this gas-water coupled imbibition model can be used to simulate the effect of imbibition on permeability, production, and water recovery. For Woodford Shale, the gas phase permeability is 8 times higher than the initial matrix permeability after imbibition, which shortens the time of obtaining the gas flow.
When applying for the experimental and field data of the Woodford shale formation as an example, the simulating result indicates that the model that considers the imbibition effect could provide a more accurate result of gas production. When considering the imbibition effect, the accumulated gas production after 100 days production is 12.68 percent higher than that neglecting the influence of imbibition.
According to the simulation, when the well is shut in during hydraulic fracturing, imbibition affects the permeability and gas production. During the early days, the imbibed fluid makes natural fractures unstable to increase the formation permeability. That increase can bring a higher gas production. Meanwhile, when larger volumes of fluid are imbibed in the formation, the gas production becomes worse because the gas flow in the porous medium is restrained by the imbibed fluid so that the relative gas permeability is low. Therefore, the shut-in time of hydraulic fracturing should be optimized to maximize gas production. Based on the simulation results, the optimal shut-in time of Woodford Shale is 3 days. Figure 18. The relation between 100 days accumulated gas production and the shut-in time.

Conclusions
This paper developed a numerical model coupling the imbibition of fracturing fluid with the gas-water flow model. This model is used to investigate the effect of imbibition on Woodford shale. The experimental and field data of Woodford shale formation were matched by the model to study the effect of imbibition on gas production and shut-in time optimization. The numerical model enables us to efficiently optimize the shut-in time to achieve a high engineering profit. Taking the Woodford Shale as an example, we contrasted the effect of imbibition on gas production and optimized the best shut-in time. The conclusions are the following.
Imbibition is shown to cause a matrix permeability decrease, while the unstable natural fracture permeability increases in the experiment. When the numerical model was developed, the relationship between the imbibed fluid volume and permeability change was considered and matched. Hence, this gas-water coupled imbibition model can be used to simulate the effect of imbibition on permeability, production, and water recovery. For Woodford Shale, the gas phase permeability is 8 times higher than the initial matrix permeability after imbibition, which shortens the time of obtaining the gas flow.
When applying for the experimental and field data of the Woodford shale formation as an example, the simulating result indicates that the model that considers the imbibition effect could provide a more accurate result of gas production. When considering the imbibition effect, the accumulated gas production after 100 days production is 12.68 percent higher than that neglecting the influence of imbibition.
According to the simulation, when the well is shut in during hydraulic fracturing, imbibition affects the permeability and gas production. During the early days, the imbibed fluid makes natural fractures unstable to increase the formation permeability. That increase can bring a higher gas production. Meanwhile, when larger volumes of fluid are imbibed in the formation, the gas production becomes worse because the gas flow in the porous medium is restrained by the imbibed fluid so that the relative gas permeability is low. Therefore, the shut-in time of hydraulic fracturing should be optimized to maximize gas production. Based on the simulation results, the optimal shut-in time of Woodford Shale is 3 days.
The flow back volume is also simulated in the model to explain the reason for the low water recovery after hydraulic fracturing in shale formations.
The contribution of this work is to close the gap between the laboratory and the field. It is helpful to understand the mechanism affecting the flow back and production after hydraulic fracturing in shale formations.