Prediction of Production Performance of Refractured Shale Gas Well considering Coupled Multiscale Gas Flow and Geomechanics

Production simulation is an important method to evaluate the stimulation effect of refracturing. Therefore, a production simulation model based on coupled fluid flow and geomechanics in triple continuum including kerogen, an inorganic matrix, and a fracture network is proposed considering the multiscale flow characteristics of shale gas, the induced stress of fracture opening, and the pore elastic effect. The complex transport mechanisms due to multiple physics, including gas adsorption/desorption, slip flow, Knudsen diffusion, surface diffusion, stress sensitivity, and adsorption layer are fully considered in this model. The apparent permeability is used to describe the multiple physics occurring in the matrix. The model is validated using actual production data of a horizontal shale gas well and applied to predict the production and production increase percentage (PIP) after refracturing. A sensitivity analysis is performed to study the effects of the refracturing pattern, fracture conductivity, width of stimulated reservoir volume (SRV), SRV length of new and initial fractures, and refracturing time on production and the PIP. In addition, the effects of multiple physics on the matrix permeability and production, and the geomechanical effects of matrix and fracture on production are also studied. The research shows that the refracturing design parameters have an important influence on the PIP. The geomechanical effect is an important cause of production loss, while slippage and diffusion effects in matrix can offset the production loss.


Introduction
After initial fracturing of the shale gas well, monitoring results of production logging show that the local fracturing stage produces no gas and about one third of perforation cluster produces no gas or contributes little to gas production [1,2]. A limited region of initial fracturing results in insufficient gas supply for propped hydraulic fractures, and fracture conductivity loss resulted from the geomechanical effect during the production process causes a rapid production decline [3][4][5][6][7]. To solve the problem of rapid production decline of shale gas wells, a large number of refracturing technologies and field experiments have been carried out in North America and achieved a well stimulation effect [8][9][10][11][12].
Refracturing is an important technology for further exploiting the stimulation effect of shale reservoirs between fractured intervals where primary fracturing is not fully utilized. Productivity prediction of a refractured shale gas well is an important method to effectively evaluate the stimulation effect. Some scholars have carried out productivity performance simulation of refractured shale gas wells using a numerical simulation method. Based on a gas-water, twophase, dual-permeability model, Tavassoli et al. [13] studied the production performance and the optimal time of refracturing in shale gas wells using numerical simulation and parameter sensitivity analysis methods. The results show that the final production of a shale gas well can be increased by about 30%. Huang et al. [14] established a numerical model to predict the production performance of a shale gas well considering different refracturing scenarios and analyzed the effects of the fracture length and conductivity on gas production. Guo et al. [15] proposed a production prediction model of refractured horizontal well in tight reservoirs and investigated the optimal refracturing scenario and time. Jayakumar et al. [16] adopted a synthetic model to evaluate the effects of fracture spacing, matrix permeability, fracture conductivity, and orientation on gas production and economic benefits of refractured wells considering different fracturing scenarios. Considering the effect of geomechanical effect on production, a multistep productivity model is presented, and the most important factors affecting the productivity of refractured well are analyzed in detail [17]. Urban et al. [18] put forward a shale gas production model with consideration of multiple pore characteristic to calculate the optimal refracturing time, production, and recovery and compared the benefits of refracturing and infill drilling. However, the multiscale flow characteristic of shale gas and the complex transport mechanism in different scale media (kerogen, inorganic matrix, and fracture) and coupling fluid flow and geomechanics are not fully considered in the abovementioned mathematical models.
Only by clarifying the mechanisms of shale gas storage and flow can we accurately evaluate the production performance of refracturing in shale gas wells. It is confirmed that shale gas is mainly stored in the shale reservoir with adsorbed gas and free gas, and the shale gas transport behaves complex multiscale characteristics. Many mathematical models have been proposed to describe the gas transport and simulate the initial fracturing production of shale gas. The current models concerning shale gas production evaluation are established based on multiple-porosity characteristics of shale reservoir and complex flow behavior in multiscale media [19][20][21][22], including Knudsen diffusion, surface diffusion, adsorption/desorption and viscous flow, and the effects of complex gas flow mechanisms and parameters related to SRV are investigated using single-, dual-, triple-and fourporosity models. Zhao et al. [23] presented a single-porosity medium model to analyze the contribution of multiple transport mechanisms on gas production. Azom and Javadpour [24] established a dual-porosity/dual-permeability continuum model to simulate shale gas production, in which Knudsen diffusion and slip flow effects in matrix pores were incorporated. Hu et al. [25] applied the dual-continuum model considering the complex gas flow behaviors and stress sensitivity to predict the shale gas production performance. However, some researchers considered that the conventional single-porosity medium and dual-continuum cannot accurately describe the multiscale gas transport process. Hence, two types of triple-continuum models are proposed. One is the triple-continuum model based on kerogen, an inorganic matrix, and fracture [26][27][28]; the other is the triple-porosity model composed of macrofracture, microfracture, and matrix [29][30][31][32]. Other researchers have presented a quadrupleporosity medium model [33][34][35] by dividing the fractures into natural microfractures and hydraulic macrofractures. Hinkley et al. [36] developed a quadruple-porosity model incorporating macrofractures, microfractures, kerogen, and inorganic matrix to simulate unconventional reservoir production. Li et al. [37] put forward a quadruple-porosity medium model composed of organic kerogen, an inorganic matrix, natural fracture, and a hydraulic fracture network to predict the shale gas production.
The treatment of a complex fracture network in a shale reservoir mainly includes a continuum model and a discrete fracture model [38]. In the continuum model, the complex fracture system is treated as continuum field. Therefore, it is not suitable for simulating large-scale fractures with strong local heterogeneity [39]. However, discrete fracture simulation based on unstructured grids can be employed to clearly describe the fluid flow in each fracture [40][41][42][43] because the complex fractures are explicitly characterized using the discrete fracture model (DFM). Therefore, it is necessary to represent the specific orientation of each fracture in this simulation method. By approximating the complex fracture network as vertical and orthogonal discrete fracture, Cipolla et al. [44,45] applied the numerical simulator to study the influence of fracture parameters on the gas production. However, it is neither practical nor advantageous to simulate a large number of natural fractures with the DFM. Thus, Moinfar et al. [46] recommended adopting the coupled continuum model and the DFM to model the unconventional reservoir production. In their article, the DFM is presented to describe the large-scale complex fractures, while the dual-continuum model is applied to simulate a large number of natural fractures. Xu et al. [47] first described the simulation method of an embedded discrete fracture model (EDFM) and then applied the EDFM to mimic the fluid flow in complex hydraulic fractures. The mixed-continuum method and DFM have been widely used to simulate shale gas production. The large-scale hydraulic fractures are handled as discrete fractures, and the dual-continuum and single-porosity media are adopted to describe the stimulated and unstimulated regions in the composite models, respectively [48][49][50][51].
In the process of shale gas production, the closure stress of main and secondary fractures increases, and the fracture network conductivity decreases gradually as the reservoir pore pressure decreases, resulting in the rapid decline of shale gas production. Aybar et al. [4,52] established a trilinear flow model and a numerical simulation model to quantitatively analyze the effect of stress-sensitive natural fracture and hydraulic fracture conductivity loss in the SRV region on long-term shale gas production, respectively. Based on reservoir simulation software, Yu and Sepehrnoori [53,54] calculated the impact of fracture conductivity reduction on shale gas production. However, the coupling mechanism of fluid flow and solid deformation is not considered in these models. Some scholars put forward the coupling fluid flow and geomechanical models (including single-porosity, double-porosity, and triple-porosity models) to simulate the coupling process of shale gas seepage and reservoir stress variation according to the theory of effective stress and pore elasticity of multiple-porosity media. Fan et al. [55] developed a fully coupled fluid flow and solid deformation model incorporating matrix and discrete hydraulic fracture systems to analyze the effect of geomechanics on matrix permeability and gas production. Gao et al. (2016) proposed a fluid-solid coupling model considering gas transport in kerogen and inorganic matrix to simulate shale gas production and established an apparent kerogen permeability model and an inorganic 2 Geofluids matrix permeability model with consideration of the effective stress variation, respectively. Peng et al. [56] presented a coupled fluid flow and geomechanical deformation model considering multiscale shale gas flow. In this mathematical model, continuity equations of kerogen, inorganic matrix, and fracture system are established, respectively, and the correlations between the permeability, porosity, and volume strain are derived, respectively. Based on the single-porosity medium approach, Kim et al. [57] developed a coupled gas seepage and stress model, in which two porosity and permeability correlations by calculating the effective stress are treated as the coupling terms. Sang et al. [58] put forward a coupled gas flow and geomechanical model based on triplecontinuum approach and derived the correlations between the porosity, permeability of kerogen, inorganic matrix, natural fracture, and volume strain. Based on the theory of pore elasticity, Peng et al. [59] proposed an apparent permeability model considering the adsorption strain and presented a coupled gas flow and geomechanical model to study the evolution process of apparent shale gas permeability under different stress boundary conditions. It is still a great challenge to accurately evaluate shale gas production due to multiple physics, coupled gas flow in multiple-porosity media, and effective stress change during shale gas production. However, at present, the factors considered in shale gas production simulation are not comprehensively enough. In this paper, a coupled flow and geomechanical model of shale gas in a kerogen-inorganic matrix-fracture system is established to evaluate shale gas production performance. The effects of stress-dependent permeability of matrix and fracture network, desorption, surface diffusion, and slippage flow on gas production are taken into account in the model. The actual production data of a shale gas well are employed to validate the model. The production performance and stimulation effect of a horizontal shale gas well are predicted considering different refracturing scenarios and fracturing design parameters based on the field production performance data of initial fracturing. In addition, the effects of multiple physics on the matrix permeability and production, and the geomechanical effects of matrix and fracture on gas production are also studied. This study and research results provide vital theoretical and engineering significance for understanding shale gas complex flow mechanisms, refracturing design, and production evaluation.

Mathematical Model
A shale reservoir should be considered a dual-porosity medium system composed of fracture and matrix after reservoir stimulation. According to the theory of dual-porosity medium, matrix pore is the main fluid storage space, while a fracture system is the main channel of fluid seepage. However, a shale matrix can be divided into organic kerogen and inorganic matrix, and organic kerogen is dispersed in an inorganic matrix system. Therefore, according to the theory of triple continuum [26,60], the migration paths of shale gas are as follows [27] during the production process, as shown in Figure 1: (1) desorption of adsorption gas in kerogen and migration through kerogen pore itself together with free gas, (2) transscale transport between kerogen and inorganic matrix, (3) internal gas transport in inorganic matrix, (4) transscale flow between inorganic matrix and fracture system, (5) and viscous flow in fracture network system. According to the abovementioned shale gas migration process, a mathematical model is established based on triplecontinuum considering multiple transport mechanisms.

Gas Transport and Mass Conservation Equation in
Kerogen Matrix 2.1.1. Apparent Matrix Permeability. As shown in Figure 2, gas transport in kerogen matrix includes multiple physics, incorporating slippage effect, Knudsen diffusion, and surface diffusion. During the shale gas production, the stresssensitive effect results in the decrease of pore radius, while the thinning adsorption layer due to gas desorption can enlarge the pore radius [61]. Because of multiple physical effects, the gas flow behavior becomes more complex. The apparent gas permeability is an important parameter to characterize the above physics. Therefore, considering slip flow, Knudsen diffusion, and surface diffusion in kerogen matrix, the total gas mass flow in kerogen is proposed by introducing a contribution coefficient to describe the contribution of slip flow and Knudsen diffusion, respectively [20].
where J t is the total mass flow (kg/m 2 /s), Kn is the Knudsen number (dimensionless), J vs , J kk , and J ks are the mass flux of slippage flow, Knudsen diffusion, and surface diffusion (kg/m 2 /s), which are expressed with Equations (2), (3), and (4), respectively [37]: Kerogen matrix Inorganic matrix Fracture system Gas flow Kerogen matrix Inorganic matrix Fracture system Gas flow

Geofluids
The apparent gas permeability in the kerogen is derived by expressing Equation (1) as the form of Darcy flow mass flowrate, as described in where K kapp is the gas apparent permeability in the porous kerogen system (μm 2 ), τ is the pore tortuosity (dimensionless), B k is the gas slip coefficient in kerogen (dimensionless), r eff is the effective kerogen pore radius (m), μ g is the gas viscosity (mPa·s), ρ g is the gas density (kg/m 3 ), M g is the molar molecular mass (kg/mol), R is the gas constant (8.314 J/mol/K), T is the reservoir temperature (K), ε kp is the percentage of kerogen pore volume in total matrix pore volume (dimensionless), φ m is the total matrix porosity (dimensionless), φ f is the natural fracture porosity (dimensionless), D kk is the Knudsen diffusion coefficient in kerogen (m 2 /s), D s is the surface diffusion coefficient (m 2 /s), V L is the Langmuir volume (m 3 /kg), P L is the Langmuir pressure (MPa), ρ s is the shale core density (kg/m 3 ), and V std is the molar gas volume at the standard condition (m 3 /mol), and P k is the kerogen system pressure (MPa). The adsorption layer and stress-dependent matrix permeability result in the decrease of nanopore radius that will also affect gas microflow behavior in matrix pores [20]. The effective kerogen pore radius is mathematically expressed by the following equation: where r ki is the initial nanopore radius in porous kerogen (m), σ me is the mean effective stress (MPa), C φ is the material constant (dimensionless), and d c is the molecular diameter of methane (4 × 10 −10 m). The slip factor, B k , is defined as [62]: where α is the tangential momentum adjustment coefficient (dimensionless) and the value of α is set to 0.8 in this study.

Mass Conservation Equation in
Kerogen Matrix. A part of shale gas is stored as free gas and adsorbed gas in kerogen bulk with a well-developed nanopore network, and gas transport mechanisms in kerogen incorporates adsorption and desorption, Knudsen diffusion, surface diffusion, and slippage flow. Therefore, with consideration of above transport mechanisms and quasisteady gas flow from kerogen to inorganic matrix, the continuity equation of gas migration in the kerogen system is written as follows: where C k is the moles of free gas per kerogen pore volume (mol/m 3 ), ε ks denotes the proportion of kerogen grain volume in total shale grain volume (dimensionless), W km is the mass transfer term between kerogen and inorganic matrix (mol/m 3 /s), and C μ denotes the adsorbed gas concentration per kerogen volume (mol/m 3 ). It is assumed that the adsorbed gas in kerogen is adsorbed on the pore wall in the form of monolayer, and the adsorbed phase occupies the pore space of kerogen. The Langmuir isotherm adsorption equation is used to describe the absolute adsorption gas volume in kerogen [26,63], which is represented by The apparent permeability of kerogen differs slightly from that of inorganic matrix, and the steady-state flow will be achieved in a relatively short time. Hence, the Warren-Root, dual-porosity model based on pseudosteady-state flow is applied to describe the gas transport between kerogen and inorganic matrix. The steady-state mass transfer term between kerogen and inorganic matrix is [26] where σ km is the shape factor of porous kerogen system (1/m 2 ) and P m is the inorganic matrix system pressure (MPa).

Gas Continuity Equation in
Inorganic Matrix. Free gas is stored in the inorganic matrix system, and the gas transport mechanism in the inorganic matrix system includes slippage effect and Knudsen diffusion considering matrix stress sensitivity. In addition, there are cross-scale mass transfer between inorganic matrix and kerogen and fracture network, because kerogen provides gas source to inorganic matrix, and inorganic matrix supplies gas to fracture system. Consequently, based on the principle of mass conservation, the gas continuity equation in shale inorganic matrix is written as where D km is the Knudsen diffusion coefficient in inorganic matrix (m 2 /s), B m is the gas slip coefficient in inorganic matrix (dimensionless), W mf is the mass transfer term between inorganic matrix and fracture network system (mol/m 3 /s), C m is the mole of free gas per inorganic matrix pore volume (mol/m 3 ), and K mi represents intrinsic permeability of inorganic matrix (μm 2 ). It takes a long time for fluid flow from the inorganic matrix to the fracture network system to get into a steady state since the permeability of inorganic matrix and hydraulic fracture network is quite different. The Warren-Root, unsteady dual-porosity model is adopted to describe the gas migration between the inorganic matrix and the fracture system. The transient mass transfer term between inorganic matter and fracture system is where P f is the pore pressure in the natural fracture system (MPa), K mapp is the gas apparent permeability in inorganic matrix pore (μm 2 ), σ mf is the transient shape factor between the inorganic matrix and the fracture network system (1/m 2 ). According to the definition [24,26], the unsteady shape factor, σ mf , is defined as where L fx and L fy represent the spacing of fractures in the x-and y-directions (m), respectively, and P i is the initial reservoir pressure (MPa). Similarly, the apparent permeability of the inorganic matrix system considering slippage effect, Knudsen diffusion, and matrix stress-sensitive effect is given by [27]: where C g is gas compressibility (MPa -1 ). The inorganic matrix pore radius, Knudsen diffusion coefficient, and slip factor in inorganic matrix are expressed by Equations (15), (16), and (17), respectively: where r mi is the initial nanopore radius in inorganic matrix (m) and r m is the effective inorganic matrix pore radius (m).

Gas Continuity Equation in
Fracture Network System. The natural fractures are activated, and then, a connected complex fracture network incorporating main and secondary fractures is created near the wellbore in the process of hydraulic fracturing in the shale reservoir. In order to facilitate the simulation calculation, we simplify the complex fracture network into a vertical orthogonal fracture network [44,45,64,65] and adopt the dual-continuum approach to characterize the fluid seepage in the fracture network area [49,51]. Considering the gas supply and production well, the gas continuity equation in the fracture network system is written as follows: where K f is the fracture permeability (μm 2 ), C f is the moles of free gas per pore volume in the fracture network system (mol/m 3 ), and Q well is the well production term (mol/m 3 /s), which is defined as where W f is the hydraulic fracture width (m), P wf is the bottom-hole flowing pressure (MPa), r e is the effective well radius (m), r w is the well radius (m), and V b is the well grid volume (m 3 ).

Geomechanical Model.
Induced stress field will be formed around the propped hydraulic fractures [66,67] and the in situ stress field will change due to the poroelastic effect and opening of fractures [68][69][70]. The reservoir permeability will be affected by both of the additional stresses. Thus, it is necessary to establish a geomechanical model considering the influence of poroelastic effect and induced stress [71]. Shale gas well production is a process of coupling gas seepage and reservoir deformation. The governing equation for multiscale shale deformation is proposed based on 5 Geofluids the multimedium pore elasticity theory. Shale deformation is mainly determined by effective stress and gas desorption. The relationship between stress and strain is presented as [56,58]: where G is the shear modulus (MPa); λ is Merla constant (MPa); ε v is the volume strain (dimensionless); σ ij represents the total stress (MPa); δ ij represents the Kroneker symbol (when i = j, then δ ij = 1, else δ ij = 0); ε s is the Langmuir volume strain (dimensionless); α k , α m , and α f are the effective stress coefficients of kerogen system, inorganic matrix system, and fracture system (dimensionless), respectively; and ε ij is the solid strain (dimensionless), which is defined as where u represents the displacement (m). The total stress expressed in Equation (20) satisfies the stress balance equation where f is the body force (MPa). Equations (20) and (21) are substituted into Equation (22), the displacement equation is written as To calculate the induced stress around the hydraulic fracture, a mathematical model of induced stress field for hydraulic fracture is presented based on the 2D displacement discontinuity theory (DDM) [72]. The hydraulic fracture is divided into N elements, and the shear and normal displacement discontinuities of any element are calculated by introducing the boundary condition Equation (25) into Equation (24). Finally, the displacement discontinuity is introduced into Equation (24) to calculate the induced stress of any element in the x-y plane.
The boundary condition of any element j is presented as follows where σ i xx and σ i yy are the induced normal stresses in the x-and y-directions (MPa), respectively; σ i xy is the induced shear stress (MPa); , and A i,j sy represent the elastic coefficients (MPa/m) [72]; D j s and D j n represent the shear and normal displacement discontinuities on the fracture element j (m), respectively; N represents the number of fracture element (dimensionless); and ΔP f represents the net pressure of opening fracture (MPa).
The increasing closure pressure during shale gas production will result in a gradual reduction in hydraulic fracture conductivity. The relationship between the hydraulic fracture conductivity and closure stress has been obtained based on experimental studies [4,5]. According to relevant experiments and previous research results [73,74], the permeability of the shale reservoir fracture system owing to increasing effective stress can be calculated by where K f s is the fracture permeability under the current effective stress (mD), K f0 is the initial fracture permeability (mD), d f is the stress sensitivity coefficient of the fracture system (MPa -1 ), and σ me represents the mean effective stress considering the impact of induced stress and pore elastic stress (MPa).

Initial and Boundary Conditions.
The initial condition is written as It is assumed that the outer boundary of the model is closed and the inner boundary is a constant bottom-hole pressure. Thus, the inner boundary condition is The outer boundary condition is where Γ I and Γ o represent the inner boundary and outer boundary, respectively.

Model Solution
Due to the high nonlinearity of the overall mathematical models, the differential Equations (8), (11), (18), and (23) 6 Geofluids are discretized using the block-center difference scheme and solved numerically using the finite difference method. The specific solution procedures are as follows: Step 1. Divide the simulation area into rectangular grids and input reservoir and initial fracture parameters Step 2. Calculate the pore pressure distribution of the fracture system, inorganic matrix, and kerogen system by the iterative method in each time step. Then, the pore pressure is substituted into Equation (23) to calculate the displacement and strain. The gas transport Equations (8), (11), and (18) are coupled to solid deformation Equation (23) through volume strain and permeability. More details can be found in Zhang et al. [35] and Charoenwongsa et al. [75] Step 3. Calculate the gas production at this time step Step 4. Determine whether the simulation time is equal to the refracturing time, and if the condition is satisfied, refracturing and reservoir parameters are updated Step 5. Continue to the next time step until the simulation time is over

Comparison and Verification of Model Calculation
Results. Based on the fracture and reservoir parameters of a typical shale gas reservoir [76], the calculation results of this model are compared with those of Grieser et al. [77], Yu et al. [76], and field data. It is worth noting that the simulation results of Yu et al. [76] are calculated from the commercial reservoir simulator CMG-IMEX. The simulation results are shown in Figure 3. It can be seen that the calculation results of our model are relatively consistent with the simulation results of Grieser et al. [77], the reservoir simulator results of Yu et al. [76], and the actual shale gas production data. It also can be seen that the daily production without considering the stress-sensitive effect is quite different from the actual field production data. Hence, the stress-sensitive effect is an important factor affecting shale gas production, which should be considered in the simulation of shale gas production performance.

Production Matching and Prediction.
The model was calibrated with monthly average production data of a horizontal shale gas well within 3 years, in southwest China [37]. The initial fracturing and reservoir parameters are listed in Table 1, and the matching and prediction results are provided in Figure 4. It can be seen that the predicted results using the model are in good agreement with the actual production performance. Next, we will perform a sensitivity analysis to study the effects of refracturing scenarios, fracture conductivity, SRV width, SRV length, and the refracturing time on the gas production and the PIP based on the initial fracturing production data. In order to improve the calculation efficiency, a single stage with several hydraulic fractures is investigated to calculate daily production and cumulative production. The total production of the horizontal well is the product of the total number of fracturing stage and the gas production of a single fracturing stage.
The length, width, and height of the representative segment are 400 m × 100 m × 40 m, respectively. The length, width, and height of each grid block are set at 5 m × 2 m × 4 m, with a local grid refinement of 1 × 7 × 1 for each hydraulic fracture. In this simulation work, the hydraulic fracture and SRV region are characterized by modifying the grid permeability of the fracture system [25,78].

Influence of Matrix Multiphysics on Permeability and
Production. To study the influence of different physics in the matrix on permeability and production, six cases are designed as shown in Table 2. When the pore pressure decreases from 30 MPa to 10 MPa, it is calculated based on Figure 5(a) that the permeability of Case#2 is lower than that of Case#1 about 50% when the stress-sensitive effect is considered. However, considering the slippage effect in Case#3 can compensate the permeability loss caused by the stresssensitive effect, and the apparent permeability increases further when the Knudsen diffusion is taken into consideration. For Case#5 and Case#6, the pore radius and intrinsic permeability decrease on account of the existence of the adsorption layer. However, the adsorption layer becomes thinner owing to gas desorption, and the apparent permeability increases gradually with the decreasing pore pressure. According to Figure 5(b), it can be observed that when the pressure is reduced from 10 MPa to 2 MPa, the apparent permeability of Case#5 and Case#6 considering the Knudsen diffusion and surface diffusion increases rapidly, especially the effect of surface diffusion. Figure 6 illustrates the contributions of slip flow, Knudsen diffusion, and surface diffusion to apparent permeability in the organic matrix as pore pressure decreases. The contribution of slip flow to apparent permeability decreases, while the contributions of surface diffusion and Knudsen diffusion increase with the decreasing pore pressure, respectively. When the pore pressure drops to approximately 16 MPa,  [76], and field production data.

Case#1
Without any physics Case#2 Consider the stress sensitivity Case#3 Consider the stress sensitivity and slip flow Case#4 Consider the stress sensitivity, slip flow, and Knudsen diffusion Case#5 Consider the stress sensitivity, slip flow, Knudsen diffusion, and adsorption/desorption Case#6 Consider the stress sensitivity, slip flow, Knudsen diffusion, adsorption/desorption, and surface diffusion 8 Geofluids the contribution of surface diffusion is greater than that of slip flow, and the contribution of Knudsen diffusion will be greater than that of slip flow when the pore pressure is further reduced to about 6 MPa. At the condition of low pore pressure, surface diffusion is the main gas transport phenomenon. As shown in Figure 7, the loss of apparent permeability considering only stress sensitivity is approximate in the range of 30% to 25% when the pore pressure declines from 30 MPa to 2 MPa. When only the adsorption layer is taken into consideration, the apparent permeability reduction declines from 45% to 12%. If both the stress sensitivity and adsorption layer are considered, the loss of apparent permeability decreases from 65% to 30%, approximately. With the decrease of pore pressure, the pore radius decreases due to stress sensitivity. However, the thinning of the adsorption layer results in an increase in the pore radius. Thus, the loss of permeability decreases as the pore pressure declines.

Geofluids
The cumulative production of shale gas considering different physical phenomena are displayed in Figure 8. By calculating the difference between the cumulative production, it is discovered that only considering stress sensitivity in Case#2 results in 17.33% less cumulative gas production than the case without any physics. However, when the slippage effect and Knudsen diffusion are considered, the cumulative production reductions corresponding to Case#3 and Case#4 are reduced to 5.3% and 1.9%, respectively. The results show that slippage and Knudsen diffusion can compensate for the production loss caused by stress sensitivity. The cumulative production increases about 31.0% when the gas adsorption/desorption effects are considered. This is because that the adsorbed gas is an important source of shale gas production. When the surface diffusion effect in Case#6 is taken into account, the cumulative production is increased 42.3% approximately compared with Case#1, which indicates that surface diffusion has a certain contribution to shale gas production.
The effects of matrix and fracture stress sensitivity on cumulative production are illustrated in Figure 9. By calculating the cumulative production, it can be discovered that the production loss caused by considering both matrix and fracture stress sensitivity is about 23.9% compared with that without any geomechanics considered. The consideration of only the effect of fracture or matrix stress sensitivity results in 17.5% and 8.4% loss in the cumulative production, respectively. The effect of geomechanics is an important cause of shale gas well production loss, and the effect of fracture geomechanics is much more serious than that of matrix stress sensitivity.

Refracturing Scenario.
In order to study the effect of the refracturing pattern on stimulation effect, the following four types of refracturing scenarios are evaluated, as shown in Figure 10. Case 1 is refracturing of the initial fractures; Case 2 indicates that new fractures are created between the initial fractures; Case 3 represents refracturing of initial fractures and propagation of the initial fractures; and Case 4 is refracturing of the initial and new fractures.
The gas flow rate and cumulative gas production of the horizontal well based on different refracturing scenarios are shown in Figures 11(a) and 11(b), respectively. In addition, Figure 12 shows the distribution of reservoir gas content after 10 years, corresponding to the four different refracturing scenarios. Production increase percentages for the four different refracturing scenarios are listed in Tables 3(a) and 3(b), respectively. It should be noted that percentages of production increase in Table 3 are calculated according to the daily production and cumulative production after 1 and 7 years of refracturing, respectively. It is observed that the PIP of daily production for the four refracturing scenarios after 1 year is 33.5%, 112.4%, 130.7%, and 117.1%, respectively, and the PIP of cumulative production is 12.8%, 37.5%, 30.5%, and 41.0%, respectively. After 7 years of refracturing production, the PIP of daily production is 1.9%, 36.0%, 56.6%, and 32.3%, and the percentage of cumulative production increment is 11.2%, 48.5%, 52.3%, and 49.5%, respectively, in comparison to the case when no refracturing is considered. It is found that the PIP of Case 1 is the smallest, and the percentage differences between Case 2, Case 3, and Case 4 are small, but much larger than that of Case 1. According to Figure 11(b) and Tables 3(a) and 3(b), it can be discovered that Case 4 is the largest as for the PIP in a certain production period, but Case 3 will become the largest until the late production stage. In general, the final differences of PIP in cumulative production between Case 2, Case 3, and Case 4 are small. However, the treatment cost of Case 2 is the smallest, because only new fractures need to be created, while it is necessary to refracture the initial fractures except for the new ones in both Case 3 and Case 4.

Influence of Initial Fracture Conductivity on
Refracturing Production. In order to study the influence of initial fracture conductivity on refracturing production. We calculated the daily production and cumulative production of shale gas wells in the two scenarios of refracturing and no refracturing when the initial fracture conductivity is 0.025 D·cm and 0.2 D·cm, respectively, as illustrated in Figures 13(a) and 13(b). For both scenarios, the refracturing fracture conductivity is 0.2 D·cm. It is observed that when the initial fracture conductivity is 0.025 D·cm, the initial daily production and cumulative production increase significantly after refracturing, and the cumulative production increases about 1 time. However, when the initial fracture conductivity is 0.2 D·cm, the daily production increase after refracturing is small, and the cumulative production hardly increases.

SRV Length of Refracturing Initial
Fracture. It is also important to investigate the effect of the refractured SRV length of initial fracture on daily production and cumulative production, and the model scenario is illustrated in Case 3 of Figure 10. The daily production, cumulative production, and the percentages of production increase according to different fracture lengths of initial fracture refracturing are provided in Figures 14(a) and 14(b) and Table 4. When the SRV length of initial fracture refracturing is 100 m, 125 m, and  Figure 9: Effects of matrix and fracture stress sensitivity on cumulative production. 26.8%, and 56.6%, and the PIP of cumulative production is 11.2%, 32.7%, and 52.3%, respectively, after 7 years of refracturing. It should be noted that when the fracture length is 100 m, this case indicates that only the initial fracture is fractured without creating new fractures. The PIP is very small when only the initial fracture is refractured without increasing the initial SRV length. However, increasing the SRV length of initial fracture refracturing can achieve a better stimulation effect. Fracture. Figures 15(a) and 15(b) illustrate the effect of SRV bandwidth of initial fracture on daily and cumulative gas production. Table 5 shows the PIP of daily production, and cumulative production after 10 years of production based on different initial SRV widths. According to the data in Table 5, it is observed that when the SRV width of initial fracture increases from 6 m to 18 m, the PIP of daily production declines from 36.0% to 17.1%, and the PIP of cumulative production dwindles from 48.5% to 20.3% after 10 years of production. It is also found that shale

Geofluids
gas well production increases with the increasing initial SRV width, but the percentage of production increases after refracturing decreases with the increasing initial SRV width.

Fracture Conductivity of Refractured New Fracture.
Based on the refracturing scenario of Case 2 shown in Figure 10, we investigated different refractured new fracture conductivities on daily production and cumulative production of the horizontal shale gas well, and the results are provided in Figures 16(a) and 16(b). The percentages of production increase after 7 years of refracturing according to different refractured new fracture conductivities are provided in Table 6. It can be seen that the daily production increment is remarkable at the initial stage of refracturing, but it will decrease with the production time. It is also observed that when the fracture conductivity increases from 0.025 D·cm to 0.2 D·cm, the percentage of daily production increment grows from 19.1% to 36.0%, and the percentage of cumulative       13 Geofluids production increment enhances from 21.2% to 48.5% after 7 years of refracturing. Increasing the fracture conductivity of refractured new fracture can markedly improve the stimulation effect of the refractured shale gas well.

SRV Length of Refractured New
Fracture. It is critical to investigate various SRV lengths of refracturing scenario based on Case 2, as displayed in Figure 10. For this purpose, three different SRV lengths are considered, and the results are given in Figures 17(a) and 17(b), and the percentage of production increase after 7 years of refracturing corresponding to different refractured new SRV lengths are listed in Table 7.
It is observed that the PIP of daily production is 8.7%, 36%, and 74.2%, and the PIP of cumulative production is 23.9%, 48.5%, and 71.4% corresponding to the SRV length of 50 m, 100 m, and 150 m, respectively, when compared with no refracturing. It is also discovered that shale gas production increases greatly in the initial stage after refracturing, as illustrated in Figure 17(a). However, the decrease of fracture conductivity and reservoir pressure with the increase of production time will result in the decrease of PIP. Generally speaking, increasing the length of SRV can significantly improve the stimulation effect of refracturing in shale gas wells. 5.3.6. SRV Width of Refractured New Fracture. The effects of SRV width of refractured new fracture on daily and cumulative productions are illustrated in Figures 18(a) and 18(b), respectively, and the incremental percentages of daily and cumulative production after 10 years of production corresponding to different SRV widths of new fracture of refracturing are displayed in Table 8. It is observed that the cumulative gas production increases with an increase in the SRV width, but the incremental extent decreases. According to the data in Table 8, it is also observed that the percentage of incremental daily production is 34.7%, 35.5%, and 36.0%, respectively, and the percentage of incremental cumulative production is 37.4%, 48.5% and 54.9% corresponding to the SRV width of 0 m, 6 m, and 12 m, respectively, in comparison to the case that no refracturing is considered. Hence, creating SRV during refracturing and increasing the SRV width of    14 Geofluids new fracture can increase the gas production of the refractured shale gas well.
5.3.7. Refracturing Time. Figures 19(a) and 19(b) show the effects of different refracturing times on the daily and cumulative productions, and Table 9 provides the percentages of incremental daily production, cumulative production after 10 years of production. Compared with no refracturing, refracturing in the third, fourth, and fifth year resulted in 36.0%, 41.7%, and 51.7% increase in daily production, and 48.5%, 43.8%, and 41.0% increase in cumulative production, respectively. According to the production data in Table 9, it can be seen that the refracturing time has an impact on the PIP of cumulative production of the horizontal shale gas well. The percentage of incremental cumulative production of refracturing in the third year is larger than that of refracturing in the fifth year. With the increase of production time, the difference of incremental cumulative production percentage corresponding to different refracturing times will become smaller.

Conclusions
In this paper, a production simulation model based on coupled fluid flow/geomechanics in triple-continuum including kerogen, inorganic matrix, and fracture network is proposed to predict the production performance of refracturing shale gas well. The mathematical model is solved numerically and validated using actual production data of a horizontal shale gas well. On this basis, the model is applied to predict daily production, cumulative production, and percentage of incremental production based on different refracturing parameters. The following conclusions are drawn: (1) Refracturing of new fracture and propagation of initial fracture can significantly increase the gas production compared with refracturing of the initial fractures. The incremental percentage of initial fracture refracturing is approximately 10%~13%, while the incremental percentage of refracturing new fracture and propagation of initial fracture is about          16 Geofluids 30%~50%. However, when the conductivity of the initial fracture is small, refracturing the initial fracture can still get a better stimulation effect (2) Increasing fracture conductivity of refractured new fracture can dramatically increase shale gas well production. The PIP after 7 years of refracturing is approximately 20%~50% when fracture conductivity is enhanced from 0.025 D·cm to 0.2 D·cm (3) The increase in shale gas production is also obvious by increasing the SRV length of refractured new fracture. The percentage of incremental cumulative production after 7 years of refracturing is about 24%~70% when the SRV length of refractured new fracture is increased from 50 m to 150 m (4) The increase in shale gas production is substantial by refracturing the initial fracture and propagating the initial fracture length. The percentage of incremental cumulative production after 7 years of refracturing is about 11.2%~52.3% as the initial fracture length is extended from 100 m to 150 m (5) The production of the shale gas well increases with an increase in the SRV width of the refractured new fracture. The incremental production percentage after 7 years is approximately 37.4%~54.9% when the SRV    17 Geofluids width of the refractured new fracture increases from 0 m to 12 m (6) Shale gas well production increases with the increasing initial SRV width, but the percentage of incremental cumulative production after refracturing declines with the increase of initial SRV width (7) Shale gas production is influenced by the refracturing time. The earlier the refracturing time is, the higher the initial production and the cumulative production will be. However, the percentage difference of cumulative production increment after refracturing will become smaller as production time increases (8) The contribution of slip flow decreases, while the contributions of surface diffusion and Knudsen diffusion become more remarkable as the pore pressure declines, respectively. The decrease in permeability loss becomes smaller as the pore pressure decreases, because the apparent permeability will increase greatly owing to the significant enhancement of diffusion and slippage effects (9) The stress-sensitive effect is an important reason for the gas production loss, but slippage and diffusion effects can offset the production loss caused by matrix stress sensitivity to a certain extent. The contribution of gas desorption and surface diffusion is significant to shale gas production

Data Availability
The data used to support the findings of this study are included within the article and are available from the corresponding authors.