A Combination Approach of the Numerical Simulation and Data-Driven Analysis for the Impacts of Refracturing Layout and Time on Shale Gas Production

Abstract

Refracturing can alleviate the rapid decline of shale gas production with a low drilling cost, but an appropriate fracture layout and optimal refracturing time have been unclear without a heavy computation load. This paper proposes a combination approach with a numerical simulation and data-driven analysis to quickly evaluate the impacts of the refracturing layout and refracturing time on shale gas production. Firstly, a multiphysical coupling model with the creep of natural fractures is established for the numerical simulation on shale gas production. Secondly, the effects of the refracturing layout and refracturing time on the shale gas production are investigated through a single factor sensitivity analysis, but this analysis cannot identify the fracture interaction. Thirdly, the influence of fractures interaction on shale gas production is explored through a combination of a global sensitivity analysis (GSA) and an artificial neural network (ANN). The GSA results observed that the adjacent fractures have more salient interferences, which means that a denser fracture network will not significantly increase the total gas production, or will reduce the contribution from each fracture, resulting in higher fracturing costs. The new fractures that are far from existing fractures have greater contributions to cumulative gas production. In addition, the optimal refracturing time varies with the refracturing layout and is optimally implemented within 2–3 years. A suitable refracturing scale and time should be selected, based on the remaining gas reserve. These results can provide reasonable insights for the refracturing design on the refracturing layout and optimal time. This ANN-GSA approach provides a fast evaluation for the optimization of the refracturing layout and time without enormous numerical simulations.

1. Introduction

A shale gas reservoir usually has an extremely low porosity and permeability, and hydraulic fracturing is one of the most successful technologies to stimulate shale gas reservoirs for the permeability enhancement [1,2,3]. However, a key problem is that shale gas production declines rapidly after hydraulic fracturing. In order to solve this quick decline problem and to avoid the uneconomical drilling of new horizontal wells, refracturing these fractured shale gas reservoirs is an effective option for the gas production enhancement. Refracturing refers to the injection of fracturing fluid into the same perforation or new hole to expand the geometry of the fracture network for the reservoir conductivity enhancement. According to the estimation by Jacobs [4], the refracturing cost is about one-third of the cost of drilling a new horizontal well. Therefore, refracturing is an economic option for the shale gas production enhancement.

The US has the fourth-largest shale gas reserves in the world and has the most advanced shale gas exploitation technology, such as the refracturing technique, which has been applied to wells with less productive and poor reservoir qualities [5]. Different refracturing treatments have been numerically and analytically investigated. For example, Tavassoli et al. [6] compared the change of gas production before and after refracturing through a sensitivity analysis on the permeability, porosity and conductivity. They concluded that the enhancement of permeability, porosity and conductivity can increase the flow rate and cumulative production of shale gas to some extent. They further introduced a long-term refracturing efficiency index to evaluate the effect of the refracturing stimulation. This index is the ratio of the cumulative natural gas production in a 30-year production cycle after and before refracturing. Comparing the effects of refracturing at different initial fractures shows that the cumulative gas production from wells with a big initial fracture spacing is similar to that from wells with a small initial fracture spacing if the final fracture spacing is the same. Dahl et al. [7] analyzed the influence of both the reservoir and fracture parameters (fracture length, width and number of fractures) on the effects of a refractured well in the Eagle Ford Shale. They found that the shale gas production after 10 years still increases with new fractures. Therefore, refracturing is an effective option in both technical and economic considerations, but the refracturing layout and refracturing time should be optimized under the reservoir pressure and stress conditions. The new fracture geometry should be determined in a spatial-temporal domain. China contains the current world’s largest shale gas reserves [8]. However, due to the diversity and complexity of the geology in China, the exploitation costs and environmental costs of China are obviously higher than those of the US. The refracturing technology in China is still in the engineering experiment stage. Research on the optimization of the refracturing design for shale gas is conducive to ensuring China’s energy security and reducing production costs. Shi et al. [9]. studied the influence of the fracturing strategy on the gas production in the Fuling shale gas field. The results show that the combination of time-phased staged fracturing and refracturing treatments can fracture each stage more effectively. Huang et al. [10] numerically simulated the gas production of the shale gas wells by considering the reservoir volume and multi-scale gas transport. Their results show that the flow rate from the existing perforation cluster is higher than that from a new perforation zone in a short term, but refracturing a new perforation layer obtains more cumulative gas production than refracturing the original fractures in the long run.

Refracturing may be affected by the mechanical properties of a fracture. However, the above studies on refracturing are based on the assumption of the elastic shale reservoir and does not consider the effect of fracture creep. Yin et al. [11] measured the creep properties of shale by indentation tests under triaxial loadings. They found that the creep behavior of a shale sample is dominated by the organics and clay constituents and follows a logarithmic law. Following the hydraulic fracturing, fractures may experience creep deformation, especially normal creep. This deformation may change the shale gas production. Rassouli et al. [12] explored the short-term and long-term creep phenomenon of a tight reservoir through experiments and investigated the influence of reservoir creep on permeability. They found that permeability decreases not only with the increase of effective stress, but also with time. Refracturing is a key measure to prevent the production from decline. Incorporating the stress creep mechanism helps provide more accurate permeability results for the field development [13]. In order to describe the change of reservoir permeability with time, this study introduces a normal creep of a natural fracture closure to investigate the relationship between the permeability decline and the optimal refracturing time.

An optimal refracturing scheme should fully consider both the reservoir properties and the gas production history. The geometric parameters (number, spacing, and length) of the fractures may significantly affect the shale gas production. However, the current optimization analysis on refracturing was carried out by only a local sensitivity analysis where the impact of a single parameter on the simulation outputs is analyzed by keeping other parameters fixed. This local sensitivity analysis cannot characterize the interactions among all input parameters, thus being of a significant drawback. The interactions are sometimes even much more significant than the individual effects [14]. For example, there is an interactive effect between fractures, and the shale gas production will be not enhanced endlessly with the increase of the fracture geometric parameters. Compared with the local sensitivity analysis, a global sensitivity analysis (GSA) can quantify the impact of the interactions among all parameters on the simulation results. One of the most widely used GSA methods is the Sobol method [15]. This is a variance-based sensitivity analysis approach which uses sensitivity indices to characterize the importance of uncertain parameters and their interactions. We will apply the Sobol GSA method in this paper to quantify the impact of the fractures and their interactions on shale gas production.

The Sobol method requires a big amount of gas production data as samples. These data on reservoir formation and gas production may be obtained from field observations or numerical simulations. Directly performing large-scale calculations in traditional numerical simulations is a heavy, even an overburdened, job and can easily lead to insufficient memory and system crash. Thus, a proxy model with a high accuracy and a high computation speed is necessary. For oil and gas productions, a variety of machine learning models have been established [16,17], including multiple linear regression (MLR), support vector machine (SVM), extreme gradient boosting (XGBoost), random forest (RF) and artificial neural network (ANN). These machine learning models are fast and accurate in the prediction of shale gas production and shale reservoir characteristics. They are often used as proxy models for shale gas production, providing a large quantity of samples for statistical analysis. Li et al. [18] applied ANN to the Sobol global sensitivity analysis and random sampling high dimensional model representations (RS-HDMRs). They used the ANN model as a proxy model of combustion kinetics and output lots of prediction results for the Sobol method and the RS-HDMR method. These accelerate the operation of the two sensitivity analyses. Wang and Chen [19] used the ANN and Sobol global sensitivity analysis to study the impact of eighteen parameters on the cumulative production of shale gas. Their work demonstrated that a well-trained ANN model can remarkably reduce the computation cost for a global sensitivity analysis. However, these artificial intelligence approaches have not been applied to any refracturing scheme design in shale gas reservoirs. ANN has the advantages in simple structure, short time consuming and strong fitting ability. It can quickly generate a large amount of sample data for the Sobol global sensitivity analysis. Hence, this study chooses the ANN as the proxy model to predict the shale gas production in the post-refracturing stage.

This paper proposes a combination approach of the numerical simulation and data-driven analysis to investigate the influence of the refracturing layout and refracturing time on shale gas production. Figure 1 presents the workflow in this study. First, a fluid-solid coupling model is formulated for a refractured shale gas reservoir. This model considers the normal creep of the natural fractures and the gas flow in a fractured three-zone reservoir. Then, the influences of the fracture layout on the gas production are numerically simulated by a single factor sensitivity analysis. Third, 500 sets of production curves under different refracturing schemes are generated for the ANN training. The ANN model is used as a proxy model for the shale gas production prediction. Finally, the shale gas production data obtained from the ANN are used for the global sensitivity analysis to explore the influence of the fractures and their interactions on the shale gas production. This study does not consider the fracture propagation and the interactions among the fractures during hydraulic fracturing. It only considers the effects of the fracture length variation on the gas production in the post refracturing stage.

2. A Numerical Simulation Model

2.1. A Three-Zone Model in a Fractured Reservoir

A fractured reservoir after hydraulic fracturing can be divided into three zones: a matrix zone, a fractured zone, and hydraulic fractures [20]. The fractured zone contains natural and artificial micro-fractures. In shale gas production, shale gas flows from the matrix to the natural fractures, then to hydraulic fractures, and finally to the well [21]. Fractures as the main channels of shale gas flow play an important role in the shale gas extraction. A three-zone model for the horizontal well plus hydraulic fracturing is shown in Figure 2, where the arrows represent the direction of the gas flow. Zone A is the rock matrix formation area far from the horizontal well. Hydraulic fracturing does not affect this zone, thus this zone is in a natural state or a single-porosity medium zone with an extremely low porosity and permeability. Zone B is a near-horizontal well zone, where hydraulic fracturing generates many micro-fractures except for natural fractures. In this zone, a dual-porosity medium model is used to describe the gas flow in the shale matrix and micro-fractures. Both adsorbed gas and free gas are stored in the matrix but only free gas is in the micro-fractures. Zone C has hydraulic fractures which have only free gas. This zone connects to Zone B through both sides and the horizontal well through its lower part. In this zone, hydro-fracturing generates a lot of visible fractures. The shale gas in zones A and B flows to zone C.

2.2. Deformation of a Shale Gas Reservoir

The closure of an unsupported natural fracture is the main deformation of reservoirs. The natural fractures are more prone to creep during shale gas production [22]. In order to simplify numerical simulations, we assume that the deformation is viscoelastic for natural fractures and elastic for shale matrix. For single-porosity medium, the constitutive law for deformation is

$${\epsilon }_{ij}=\left(\frac{1+v}{E}{\sigma }_{ij}-\frac{v}{E}{\sigma }_{kk}{\delta }_{ij}\right)-\frac{\alpha }{3K}{p}_k{\delta }_{ij}-\frac{{\epsilon }_s}{3}{\delta }_{ij}$$

(1)

where the subscript $i$ and $j$ represent any two orthogonal directions of the $x$, $y$ and $z$ directions. The compressive strain and stress are positive. ${\epsilon }_{ij}$ is the strain tensor, $v$ is Poisson’s ratio, $E$ is Young’s modulus, and ${\sigma }_{ij}$ is the stress tensor. $p_k$ is the pore pressure in the matrix. ${\sigma }_{kk}={\sigma }_{11}+{\sigma }_{22}+{\sigma }_{33}$ is the sum of normal stresses. ${\delta }_{ij}$ is the Kronecker function. $\alpha$ is Biot’s coefficient of the matrix. $K$ is the bulk modulus of the shale matrix. ${\epsilon }_s$ is the adsorption volumetric strain.

The Navier equation for the deformation of the rock mass under body force, pore pressure, and adsorption/desorption swelling is obtained as

$$G{u}_{i,jj}+\frac{G}{1-2v}{u}_{j,ji}+\alpha p_{k,j}+K{\epsilon }_{s,i}+f_i=0$$

(2)

where $G$ is the shear modulus of the shale matrix.

In Zone B, the matchstick model is used to represent the dual-porosity medium of the fractures and matrix (see Figure 3). Fractures and matrix have different mechanical properties: The shale matrix blocks are in elastic deformation and are connected by viscoelastic fractures. Their mechanical connection is shown in Figure 4 and the virtual box is considered as a representative elementary volume (REV). This dual-porosity medium has the following constitutive equation:

$${\epsilon }_{ij}=\left(\frac{1+v}{E_t}{\sigma }_{ij}-\frac{v}{E_t}{\sigma }_{kk}{\delta }_{ij}\right)-\frac{\alpha }{3K_t}{p}_m{\delta }_{ij}-\frac{\beta }{3K_t}{p}_f{\delta }_{ij}-\frac{{\epsilon }_s}{3}{\delta }_{ij}$$

(3)

where $E_t$ is the elastic modulus of the dual-porosity medium. $\beta$ is Biot’s coefficient in fracture, and $p_f$ is the gas pressure in the fracture.

The fractures have a great influence on the mechanical properties of shale reservoir. The elastic moduli of the matrix block and REV have the following relationship:

$$\frac{1}{E_{\mathrm{t}}}=\frac{1}{E}+\frac{1}{K_{\mathrm{n}}a}$$

(4)

where $K_{\mathrm{n}}$ is the normal stiffness of the fracture and of the time-dependence.

The stress ${\sigma }_{ij}$ is expressed in terms of strain ${\epsilon }_{ij}$ as

$${\sigma }_{ij}=2G_t{\epsilon }_{ij}-{\epsilon }_s{K}_t{\delta }_{ij}+\left(\alpha p_m+\beta p_f\right){\delta }_{ij}+\lambda {\epsilon }_{kk}{\delta }_{ij}$$

(5)

where $E_t$ is the shear modulus of the dual-porosity medium, and $\lambda$ is the lame constant.

Therefore, the Navier equation in the dual-porosity medium is

$$G_t{u}_{i,jj}+(\lambda +G_t)u_{j,ji}+\alpha p_{m,i}+\beta p_{f,i}-K_t{\epsilon }_{s,i}+f_i=0$$

(6)

2.3. Gas Flow Equation in the Matrix

The mass conservation equation of shale gas in a single-porosity medium is

$$\frac{\partial m_k}{\partial t}+\nabla \cdot \left(-{\rho }_g\frac{k_k}{\mu }\nabla p_k\right)=Q_{\mathrm{s}}$$

(7)

where ${\rho }_g$ is the density of shale gas, $k_k$ is the matrix permeability, $\mu$ is the dynamic viscosity coefficient of the reservoir gas, and $m_k$ is the gas mass including both free gas and adsorbed gas as

$$m_k={\rho }_g{\varphi }_k+{\rho }_{ga}{\rho }_s\frac{V_L{p}_k}{P_L+p_k}$$

(8)

where ${\varphi }_k$ is the porosity of the single-porosity medium zone.

Substituting Equation (8) into Equation (7) gets the gas flow equation in shale matrix as

$${\varphi }_k\frac{\partial P_k}{\partial t}+P_k\frac{\partial {\varphi }_k}{\partial t}+P_a{\rho }_s\frac{V_L{P}_L}{{\left(P_L+p_k\right)}^2}\frac{\partial p_k}{\partial t}-\nabla \left(p_k\left(\frac{k_k}{\mu }\right)\nabla p_k\right)=0$$

(9)

The change of the matrix porosity is described as [23]

$$\frac{d{\varphi }_k}{{\varphi }_k}=\left(\frac{1}{K}-\frac{1}{K_p}\right)\left(d\overline{\sigma }-d{p}_k\right)$$

(10)

where $K$ is the bulk modulus of shale and $K_p$ is the bulk modulus of pores.

The bulk modulus of shale is much larger (that is $1/K\ll 1/K_p$) and the compression coefficient of the matrix is $c_m=1/K_p$. Integrating Equation (10) obtains

$${\varphi }_k={\varphi }_{k0}\exp \left\{-c_m\left[\left(\overline{\sigma }-{\overline{\sigma }}_0\right)-\left(p_k-p_0\right)\right]\right\}$$

(11)

where the subscript 0 represents the initial state.

For a plane strain problem (${\epsilon }_{zz}=0$), ${\sigma }_{zz}$ is expressed as

$${\sigma }_{zz}=v\left({\sigma }_{xx}+{\sigma }_{yy}\right)+\alpha (1-2v)p_k+\frac{E{\epsilon }_s}{3}$$

(12)

The mean stress $\overline{\sigma }$ is calculated by

$$\overline{\sigma }=\frac{1}{3}{\sigma }_{kk}=\frac{1}{3}({\sigma }_{xx}+{\sigma }_{yy}+{\sigma }_{zz})$$

(13)

Thus,

$$\overline{\sigma }-{\overline{\sigma }}_0=\frac{1-2v}{3}\alpha \left(p-p_0\right)+\frac{E{\epsilon }_s}{9}$$

(14)

Substituting Equation (14) into Equation (11) obtains the porosity ratio of the single-porosity medium as

$$\frac{{\varphi }_k}{{\varphi }_{k0}}=\exp \{\left(-c_m\right)[\left(\frac{1-2v}{3}\alpha -1\right)\left(p_k-p_0\right)+\frac{E{\epsilon }_L}{9}(\frac{p_k}{p_k+P_L}-\frac{p_0}{p_0+P_L})]\}$$

(15)

The shale gas flow in the matrix has more complicated mechanisms. In this study, three flow mechanisms are included: viscous flow, Knudsen diffusion and surface diffusion. The permeability of a single-porosity medium is expressed as

$$k_k=\frac{4\mu J_t}{\pi d_p^2{\rho }_g}\cdot \frac{1}{\nabla p_k}$$

(16)

where $J_t$ is the total flow flux of the matrix pores including the viscous flow, the Knudsen diffusion and the surface diffusion [20].

The continuity equation of the gas flow in the matrix of the dual-porosity medium is

$$\frac{\partial m_m}{\partial t}-\nabla \left({\rho }_g\frac{k_m}{\mu }\nabla p_m\right)=-\frac{{\pi }^2{\rho }_g{k}_m}{\mu }\left(\frac{1}{l_x^2}+\frac{1}{l_y^2}\right)\left(p_m-p_f\right)$$

(17)

where $m_m$ is the gas mass including the adsorbed gas and free gas:

$$m_m={\rho }_g{\varphi }_m+{\rho }_{ga}{\rho }_s\frac{V_L{p}_m}{p_L+p_m}$$

(18)

The governing equation of the gas flow in the matrix is then obtained as

$$\left({\varphi }_m+\frac{p_a{\rho }_s{V}_L{P}_L}{{\left(P_L+p_m\right)}^2}\right)\frac{\partial p_m}{\partial t}+p_m\frac{\partial {\varphi }_m}{\partial t}-\nabla \left(\frac{k_m}{\mu }{p}_m\nabla p_m\right)=\frac{{\pi }^2{p}_m{k}_m}{\mu }\left(\frac{1}{l_x^2}+\frac{1}{l_y^2}\right)\left(p_f-p_m\right)$$

(19)

where ${\phi }_m$ and $k_m$ are the porosity and permeability of the matrix, respectively.

The porosity of the matrix in a dual-porosity medium zone has the same form as that in a single-porosity medium zone. It is expressed as

$$\frac{{\varphi }_m}{{\varphi }_{m0}}=\exp \left\{\left(-c_m\right)\left[\begin{array}{l}\frac{1-2v}{3}\left[\alpha \left(p_m-p_0\right)+\beta \left(p_f-p_0\right)\right]\\ -\left(p_m-p_0\right)+\frac{E_t{\epsilon }_L}{9}\left(\frac{p_m}{p_m+P_L}-\frac{p_0}{p_0+P_L}\right)\end{array}\right]\right\}$$

(20)

The permeability ratio is obtained through the cubic law as

$$\frac{k_m}{k_{m0}}=\exp \left\{\left(-3c_m\right)\left[\begin{array}{l}\frac{1-2v}{3}\left[\alpha \left(p_m-p_0\right)+\beta \left(p_f-p_0\right)\right]\\ -\left(p_m-p_0\right)+\frac{E_t{\epsilon }_L}{9}\left(\frac{p_m}{p_m+P_L}-\frac{p_0}{p_0+P_L}\right)\end{array}\right]\right\}$$

(21)

2.4. Gal Flow in the Fractures and Creep of the Natural Fractures

The governing equation of the gas flow in the fractures is

$${\varphi }_f\frac{\partial p_f}{\partial t}+p_f\frac{\partial {\varphi }_f}{\partial t}-\nabla \left(\frac{k_f}{\mu }{p}_f\nabla p_f\right)=\frac{{\pi }^2{p}_m{k}_m}{\mu }\left(\frac{1}{l_x^2}+\frac{1}{l_y^2}\right)\left(p_m-p_f\right)$$

(22)

where ${\phi }_f$ and $k_f$ are the porosity and permeability of the fracture network, respectively.

This equation ignores the desorption of shale gas on the fractures. The deformation of the REV along the horizontal direction (x or y direction) is the summation of the deformation of the fractures and matrix in that direction:

$$\Delta u_{\mathrm{t}}=\Delta b+\Delta a$$

(23)

where $b$ is the fracture width, $a$ is the matrix block width, and $\Delta u_{\mathrm{t}}$ is the total deformation of the REV in the horizontal direction.

Under b << a, Equation (23) is rewritten in the form of strain:

$$\Delta b=\Delta u_{\mathrm{t}}-\Delta a=(a+b)\Delta {\epsilon }_{\mathrm{t}}-a\Delta {\epsilon }_m\approx a\Delta {\epsilon }_{\mathrm{t}}-a\Delta {\epsilon }_m$$

(24)

where $\Delta {\epsilon }_{\mathrm{t}}$ and $\Delta {\epsilon }_m$ are the linear strains of the REV and the matrix block along the horizontal direction, respectively.

The linear deformation of the REV and the shale matrix block in the x direction is

$$\Delta u_{\mathrm{t}}\approx a\Delta {\epsilon }_t=\frac{a}{E_{\mathrm{t}}}\left[\Delta {\sigma }_{\mathrm{ex}}-\mu \left(\Delta {\sigma }_{\mathrm{e}y}+\Delta {\sigma }_{\mathrm{e}z}\right)\right]-\Delta {\epsilon }_{\mathrm{s}}$$

(25)

$$\Delta a\approx a\Delta {\epsilon }_m=\frac{a}{E}\left[\Delta {\sigma }_{\mathrm{ex}}-\mu \left(\Delta {\sigma }_{\mathrm{e}y}+\Delta {\sigma }_{\mathrm{e}z}\right)\right]-\Delta {\epsilon }_{\mathrm{s}}$$

(26)

where, $\Delta {\sigma }_{\mathrm{e}y}$, and $\Delta {\sigma }_{ez}$ are the effective stress increments in the x, y and z directions, respectively.

Substituting Equations (25) and (26) into Equation (24) gets the change of fracture width as

$$\Delta b=\left(\frac{a}{E_{\mathrm{t}}}-\frac{a}{E}\right)\left[\Delta {\sigma }_{\mathrm{ex}}-\mu \left(\Delta {\sigma }_{\mathrm{e}y}+\Delta {\sigma }_{ez}\right)\right]$$

(27)

Substituting Equation (4) into Equation (27) has

$$\Delta b=\left(\frac{1}{K_{\mathrm{n}}}\right)\left[\Delta {\sigma }_{\mathrm{ex}}-\mu \left(\Delta {\sigma }_{\mathrm{e}y}+\Delta {\sigma }_{ez}\right)\right]$$

(28)

The effective stress is

$$\Delta {\sigma }_e=\Delta \sigma -\beta \cdot \Delta P_f$$

(29)

For a plane strain problem (${\epsilon }_{zz}=0$), the ${\sigma }_{zz}$ in the dual-porosity medium is

$${\sigma }_{zz}=v\left({\sigma }_{xx}+{\sigma }_{yy}\right)+\alpha (1-2v)p_m+\beta (1-2v)p_f+\frac{E{\epsilon }_s}{3}$$

(30)

Therefore, the total stress increment is

$$\Delta \sigma =\left[\begin{array}{c}\Delta {\sigma }_{xx}\\ \Delta {\sigma }_{yy}\\ \Delta {\sigma }_{zz}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \alpha (1-2v)(p_m-p_0)+\beta (1-2v)(p_f-p_0)+\frac{E({\epsilon }_s-{\epsilon }_{s0})}{3}\end{array}\right]$$

(31)

The effective stress is obtained as

$$\Delta {\sigma }_e=\Delta \sigma -\beta \cdot \Delta P_f=\left[\begin{array}{c}-\beta \cdot \Delta P_f\\ -\beta \cdot \Delta P_f\\ \Delta {\sigma }_{zz}-\beta \cdot \Delta P_f\end{array}\right]$$

(32)

The fracture normal deformation $\Delta b_x$ and $\Delta b_y$ can be obtained by substituting Equation (32) into Equation (28). The deformation $\Delta a_x$ and $\Delta a_y$ of the matrix block can be obtained by substituting Equation (32) into Equation (26). For the matchstick model, the porosity is

$$\phi =1-\frac{a^2}{{(a+b)}^2}\approx \frac{2b}{a}$$

(33)

The horizontal permeability is

$$k_x=k_y=\frac{b^3}{12a}$$

(34)

The natural fracture is a significant factor affecting the reservoir deformation. Wang et al. [24] described the general constitutive relation for a fracture as

$$\Delta {\delta }_I=D_{IJ}\Delta {\tau }_J$$

(35)

$$\left[\begin{array}{c}\Delta {\delta }_n\\ \Delta {\delta }_s\\ \Delta {\delta }_t\end{array}\right]=\left[\begin{array}{lll}{D}_{nn}\hfill & D_{ns}\hfill & D_{nt}\hfill \\ D_{sn}\hfill & D_{ss}\hfill & D_{st}\hfill \\ D_{tn}\hfill & D_{ts}\hfill & D_{tt}\hfill \end{array}\right]\left\{\begin{array}{c}\Delta {\sigma }_n\\ \Delta {\tau }_s\\ \Delta {\tau }_t\end{array}\right\}$$

(36)

where $\Delta {\sigma }_n$ is the normal stress; $\Delta {\tau }_s$ and $\Delta {\tau }_t$, are the shear stresses; $\Delta {\delta }_n$ is the closure displacement of the joint; $\Delta {\delta }_s$ and $\Delta {\delta }_t$ are the shear displacements. The tensor $D_{IJ}$ ($I,J=n,s,t$) represents nine components of the joint flexibility matrix.

In the matchstick model, the fracture closure is mainly caused by the normal deformation and the influence of tangential creep and the dilatancy effect is ignored. Equation (36) is rewritten as

$$\Delta {\delta }_n=D_{nn}\Delta {\tau }_n$$

(37)

where $D_{nn}=1/K_n$ is the normal flexibility.

The deep complex environment has an impact on the mechanical properties of shale and induces a nonlinear deformation. During the production of shale gas, the fractures will be normally closed, due to the increase of closure stress. When a normal closure stress is applied, the fractures first have an instantaneous elastic deformation, and then a viscoelastic closure for a long time. The fracture closure process is described by the generalized Kelvin model [25] in Figure 5 as

$$\epsilon =\frac{{\sigma }_0}{K_1}+\frac{{\sigma }_0}{K_2}\left(1-\exp (-K_{2}t/{\eta }_1)\right)$$

(38)

where $K_1$ and $K_2$ are the stiffness coefficients of two springs, and ${\eta }_1$ is the viscous coefficient.

Therefore, the normal stiffness is also a function of time as

$$\frac{1}{K_n(t)}=\frac{1}{G_{1n}}+\frac{1}{G_{2n}}\left[1-\exp \left(-G_{2n}t/{\eta }_n\right)\right]$$

(39)

$$K_n(t)=\frac{1}{\frac{1}{G_{1n}}+\frac{1}{G_{2n}}\left[1-\exp \left(-G_{2n}t/{\eta }_n\right)\right]}$$

(40)

where $G_{1n}$, $G_{2n}$ and ${\eta }_n$ are the normal elastic modulus, viscoelastic modulus and viscous coefficient of fracture, respectively.

The above analysis indicates that the pressure in the fractures and pores decreases and the effective stress of the reservoir increases in the shale gas production process. The viscoelastic deformation of the fractures and the elastic deformation of the matrix are the main causes for the reduction of the reservoir permeability. Figure 6 presents the cross-coupling processes for the shale gas production model. In the process of shale gas production, gas flows from the matrix to the fracture network and finally to the horizontal well. The gas pressure in the fracture network and the shale matrix gradually decreases. Due to the increase of the effective stress, the viscoelastic closure of the fractures and the elastic deformation of the matrix in the reservoir are more obvious, which leads to the further decrease of the permeability in the reservoir. Therefore, the shale gas production rate gradually decreases after the initial fracturing. The refracturing is conducive to restoring the reservoir permeability and increasing the gas production.

2.5. Gas Flow Equations in the Hydraulic Fractures

The aperture of the hydraulic fractures is usually large, but the gas flow in the hydraulic fractures can be still described by Darcy’s law as

$$v_{hf}=-\frac{k_{hf}}{\mu }{b}_{hf}{\nabla }_T{p}_{hf}$$

(41)

where $v_{hf}$ and $p_{hf}$ are the gas velocity and gas pressure in the hydraulic fractures.$b_{hf}$ and $k_{hf}$ are the width and permeability of the hydraulic fractures.

The continuity equation of the gas flow in the hydraulic fractures is

$$b_{hf}\frac{\partial \left({\rho }_g{\varphi }_{hf}\right)}{\partial t}+{\nabla }_T\left({\rho }_g\cdot v_{hf}\right)=0$$

(42)

where ${\rho }_g$ is the density of shale gas.

A hydraulic fracture filled with proppant is regarded as a porous medium. Its permeability is evolving in an exponential relationship:

$$k_{hf}=k_{hf0}{e}^{-3c_{hf}\left({\sigma }_e-{\sigma }_{e0}\right)}$$

(43)

where $c_{hf}$ is the compression coefficient of the hydraulic fracture filled with proppant, which is usually obtained through experiment. In this study, the field fracturing data in the Southeastern Sichuan Basin [26] are used to fit the permeability-pressure curve of the hydraulic fracture filled with high-strength ceramsite of 40/70 mesh. The compression coefficient is obtained as 0.0094.

2.6. Model Validation

The above governing equations are numerically solved by the finite element method within the platform of COMSOL Multiphysics. The equation for the shale reservoir deformation is solved by the solid mechanics module and the equations for gas flow are solved by the PDE module. This study uses field data from a horizontal shale gas well in Sichuan Basin, Southwest China for the model verification [8]. The horizontal wellbore is initially stimulated with six stages. The average half-length of the hydraulic fractures is 70 m, and the fracture spacing is about 70 m. Following 4 years of extraction, the gas production rate dropped sharply, and the refracturing job for this well was carried out to improve the output of shale gas. According to the microseismic monitoring data, the third to sixth hydraulic fractures from the original stimulation treatment were refractured. The computational parameters are listed in

Table 1. List of parameters of the shale reservoirs in the Sichuan Basin.

Parameters	Value
Elastic modulus of shale, $E$ (Pa)	2.2 × 1010
Elastic modulus of shale particles, $E_s$ (Pa)	5.5 × 1010
Methane viscosity, $\mu$ (Pa s)	2 × 10−5
Poisson ratio of shale, $v$	0.2
Reservoir temperature, $T$ (K)	358
$\mathrm{Langmuir}\mathrm{volume}\mathrm{constant},V_L$ (m3/kg)	0.003
Langmuir pressure constant, $P_L$ (MPa)	6 × 106
Langmuir volumetric strain constant, ${\epsilon }_L$	0.00185
Shale reservoir density, ${\rho }_z$ (kg/m3)	2580
Matrix pore flexion, $\tau$	11
Avogadro constant, $N_A$ (1/mol)	6.02 × 1023
Methane molar mass, $M_g$ (kg/mol)	0.016
Fractal Coefficient, $D_f$	2.6
Natural fracture spacing, $l_x=l_y$ (m)	1
Initial matrix porosity, ${\varphi }_{m0}$	0.005
Initial fracture porosity, ${\varphi }_{f0}$	0.001
Porosity of hydraulic fractures, ${\varphi }_{hf}$	1.78 × 10−6
Initial fracture permeability, $K_f$ (m2)	1 × 10−17
Initial matrix permeability, $K_m$ (m2)	1 × 10−19
Initial hydraulic fracture permeability, $K_{hf}$ (m2)	7.5 × 10−9
Initial pressure of shale reservoir, $p_0$ (MPa)	20
Compressibility coefficient of fractures, $c_f$ (1/Pa)	1 × 10−7
Compressibility coefficient of matrix, $c_m$ (1/Pa)	1 × 10−9
Compressibility coefficient of hydraulic fractures, $c_{hf}$ (1/Pa)	1 × 10−8
$\mathrm{Surface}\mathrm{diffusion}\mathrm{coefficient},D_{s0}$ (m2/s)	1 × 10−7
Initial opening degree of hydraulic fractures, $b_{hf}$ (m)	0.0003
Bottomhole flowing pressure, $p_{wf}$ (MPa)	15
Normal elastic modulus of natural fracture surface, $G_{1\mathrm{n}}$(GPa/m)	80
Viscoelastic modulus of natural fracture surface, $G_{2\mathrm{n}}$ (GPa/m)	18
Viscous coefficient of natural fracture surface, ${\eta }_n$ (GPa d/m)	9000

.

The production data in 10 years are calculated. The comparison between the field data and the simulated results is presented in Figure 7a for the gas production rate and in Figure 7b for the cumulative gas production. The gas production rate after refracturing is significantly higher than that without refracturing. The simulated results are generally in good agreement with the field data. This verified the effectiveness of this computation model.

3. ANN Based Global Sensitivity Analysis (ANN-GSA) Algorithm

The global sensitivity analysis can estimate the influence of the interaction among the parameters on the results, but a large number of samples are required. An artificial neural network (ANN) has the advantages of a high accuracy and fast running speed. This study introduces the ANN as a proxy model to formulate an ANN-GSA for the shale gas production prediction. This ANN-GSA works as follows: (1) Randomly generate multiple input parameters (fracture length), and use numerical simulation to calculate the shale gas production under different refracturing designs. (2) Determine the appropriate ANN structure and train the ANN model with the samples generated by the numerical simulation. (3) Use the well-trained ANN model as a proxy model to predict shale gas production and generate the sample data for the GSA estimation of the sensitivity index. The workflow of the ANN-GSA is shown in Figure 8.

3.1. Artificial Neural Network (ANN)

ANN has been a hot topic of artificial intelligence since its birth. Its development is inspired by human biology and can simulate some simple performances of the human brain. Prior to making predictions, the network should be trained with actual data. The training process has two stages: the signal forward propagation to obtain the output value, and the error back propagation to correct the weight matrix. The output value is continuously approaching to the actual value when the model is being trained [27,28]. In recent years, ANN has made great progress in oil and gas engineering and is widely used in the prediction of production and reservoir characteristics, the development of techno-economic framework for shale gas exploitation, and the simulation of gas flow in porous media [29].

The multilayer perceptron (or multilayer feed-forward network) is the most basic and important neural network which is connected by several groups of neurons (or nodes). As shown in Figure 9, the basic structure of a neural network is represented by a single-layer perception algorithm and includes the weighted summation and function processing [30]. The input-output relationship of a single neuron is generally expressed as

$$y=f\left({\displaystyle \sum _{i=1}^n{x}_i{\omega }_i}+b\right)$$

(44)

where $x$ is the input and $y$ is the output. ${\omega }_i$ is the weight of the connection between cells to reflect the strength of the information transmission between neurons.

Figure 10 presents a complete ANN with three layers: input layer, output layer and one hidden layer. The input layer receives data from the outside. The output layer realizes the output of the calculation results. The hidden layer is located between the input layer and the output layer and cannot be directly observed from outside of the system. In the forward propagation process, the input layer receives the data signal and passes to the hidden layer with the prescribed operation, and then passes from the hidden layer to the output layer to produce the result.

The neural network training is to minimize the error of the output by continuously correcting the weights. The error of the neural network is measured by the mean square error (MSE):

$$\mathrm{MSE}=\frac{1}{n}{\displaystyle \sum _{m=1}^n{\left(y_m-{\stackrel{⌢}{y}}_m\right)}^2}$$

(45)

where ${\stackrel{⌢}{y}}_m$ is the prediction obtained from the neural network and $y_m$ is the actual value.

The back-propagation (BP) algorithm is the most commonly used algorithm for the multilayer perceptron training. The weight changes in the direction of the maximum error gradient until the error reaches the minimum.

3.2. Global Sensitivity Analysis (GSA)

A sensitivity analysis is a common approach to evaluate the importance of each input model parameter to the output. A global sensitivity analysis does not only provide the contribution of each parameter to the model output, but also evaluates the influence of the interaction among the parameters on the model output. In this paper, the Sobol method is used to analyze the effects of nine fracture parameters and their interactions on the cumulative gas production. In the Sobol method, the variance of the model output is decomposed into fractions caused by the input and its combination. We introduce a d-dimensional model $Y=f(x)$, where $x$ represents the set of input variables $\left\{x_1,\dots ,x_d\right\}$. This function is decomposed into [31]

$$f(x)=f_0+{\displaystyle \sum _{i=1}^d{f}_i}\left(x_i\right)+{\displaystyle \sum _{i<j}^d{f}_{ij}}\left(x_i,x_j\right)+\cdots +f_{1,2,\dots ,d}\left(x_1,x_2,\dots ,x_d\right)$$

(46)

The variance of each term in Equation (31) is expressed as

$$\mathrm{Var}[Y]={\displaystyle \sum _{i=1}^d{V}_i}+{\displaystyle \sum _{i<j}^d{V}_{ij}}+\cdots +V_{1,2,\dots ,d}$$

(47)

where $V_i={\mathrm{Var}}_{xi}\left[E_{x~i}\left(Y|x_i\right)\right]$, and $V_{ij}={\mathrm{Var}}_{xij}\left[E_{x~ij}\left(Y|x_i,x_j\right)\right]-V_i-V_j$. The subscript $x_{\sim i}$ denotes the set of all independent variables except $x_i$. The effect of the ith input parameter to the output variance can be expressed by the first order index of a single input variable as

$$S_i=\frac{V_i}{\mathrm{Var}[Y]}$$

(48)

The second order index $S_{ij}$ denotes the contribution of the interaction between the ith and jth input parameters to the output variance.

$$S_{ij}=\frac{V_{ij}}{\mathrm{Var}[Y]}$$

(49)

Finally, the total order index of the ith parameter, which means the contribution of $x_i$ and its interactions to the output variance, is defined as

$$S_{Ti}=\frac{E_{x_{~i}}\left({\mathrm{Var}}_{x_i}\left[Y|x_{~i}\right]\right)}{\mathrm{Var}[Y]}$$

(50)

where ${\displaystyle \sum {}_{i=1}^d}{S}_i+{\displaystyle \sum {}_{i<j}^d{S}_{ij}}+\cdots +S_{1,2,\dots ,d}=1$ and ${\displaystyle \sum {}_{i=1}^d}{S}_{Ti}\ge 1$.

4. Sensitivity Analysis for the Refracturing Scheme

4.1. Computational Model and the Refracturing Scheme

The reservoir parameters are listed in

Table 2. Main parameters used in the simulations.

Parameter	Value
Elastic modulus of shale, $E$ (Pa)	2 × 1010
Elastic modulus of shale particles, $E_s$ (Pa)	5.5 × 1010
Methane viscosity, $\mu$ (Pa s)	2 × 10−5
Poisson ratio of shale, $v$	0.33
Reservoir temperature, $T$ (K)	350
$\mathrm{Langmuir}\mathrm{volume}\mathrm{constant},V_L$ (m3/kg)	0.0025
Langmuir pressure constant, $P_L$ (MPa)	3 × 106
Langmuir volumetric strain constant, ${\epsilon }_L$	0.00185
Shale reservoir density, ${\rho }_z$ (kg/m3)	2580
Matrix pore flexion, $\tau$	11
Avogadro constant, $N_A$ (1/mol)	6.02 × 1023
Methane molar mass, $M_g$ (kg/mol)	0.0195
Fractal Coefficient, $D_f$	2.6
Natural fracture spacing, $l_x=l_y$ (m)	1
Initial matrix porosity, ${\varphi }_{m0}$	0.065
Initial fracture porosity, ${\varphi }_{f0}$	0.012
Porosity of hydraulic fractures, ${\varphi }_{hf}$	1.78 × 10−6
Initial fracture permeability, $K_f$ (m2)	1 × 10−17
Initial matrix permeability, $K_m$ (m2)	1 × 10−19
Initial hydraulic fracture permeability, $K_{hf}$ (m2)	7.5 × 10−9
Initial pressure of shale reservoir, $p_0$ (MPa)	20
Compressibility coefficient of fractures, $c_f$ (1/Pa)	1 × 10−7
Compressibility coefficient of matrix, $c_m$ (1/Pa)	1 × 10−9
Compressibility coefficient of hydraulic fractures, $c_{hf}$ (1/Pa)	1 × 10−8
$\mathrm{Surface}\mathrm{diffusion}\mathrm{coefficient},D_{s0}$ (m2/s)	1 × 10−7
Initial opening degree of hydraulic fractures, $b_{hf}$ (m)	0.0003
Normal elastic modulus of natural fracture surface, $G_{1\mathrm{n}}$ (GPa/m)	80
Viscoelastic modulus of natural fracture surface, $G_{2\mathrm{n}}$ (GPa/m)	18
Viscous coefficient of natural fracture surface, ${\eta }_n$ (GPa d/m)	9000

[20]. The geometric model is 1100 m long, 500 m wide and 50 m thick. The model includes three hydraulic fractures with the spacing of 300 m, the length of 120 m. Because the shale reservoir and fracture distribution are symmetrical, only one-half of the production domain is simulated. The left and lower sides are constrained for deformation. The right and upper sides are free for deformation but a compressive stress of 40 MPa is used to simulate the in-situ stress condition. Zero flux on both sides is setting for the flow boundary. The initial pore pressure in the shale reservoir is 20 MPa. The following assumptions are made:

(1)

The closure of natural fractures is viscoelastic and the deformation of the matrix and hydraulic fractures is elastic.

(2)

Refracturing activates the opening of the existing natural and hydraulic fractures to the initial values.

The creep of the natural fractures heavily affects the reservoir production attenuation and the refracturing scheme. In order to illustrate the effect of the fracture creep on shale gas production, the cumulative gas production from elastic and viscoelastic reservoirs in 10 years was calculated without refracturing. Figure 11 compares their gas productions. In the early stage of production, their gas production curves are relatively close. With the continuous production of shale gas, the viscoelastic reservoir shows a more obvious production decline, and the production curve approaches to flat in the later stage. The refracturing technology can reopen the closed fractures, increase the complexity of the fracture network and the reservoir permeability. The key issue is to reasonably set the refracturing layout and time. According to the distribution of the hydraulic fractures and remaining shale gas, two refracturing schemes are proposed: refracturing in existing perforations and refracturing new fractures. Thus, three cases for refracturing are generally considered: (1) Fracturing new fractures between existing fractures. (2) Fracturing new fractures on both sides of existing fractures. (3) Refracturing in existing perforations [32,33]. The three refracturing cases are shown in Figure 12. The fracturing length is 120 m and the refracturing time is three years after the initial fracturing.

4.2. Single Factor Sensitivity Analysis for the Geometry of the Fractures on Gas Production

The gas flow rate and cumulative gas production from the three cases are compared in Figure 13. Case 3 has the best performance and case 1 has the worst performance. This is because the existing fractures extend into more unfractured reservoirs with abundant shale gas reserves and are less affected by other fractures. Due to the small remaining shale gas reserve and the low gas pressure in the middle of the existing fractures, the new fractures between the existing fractures do not significantly enhance gas production. In contrast, the new fractures on both sides of the existing fractures perform better. However, from the perspective of a fracture propagation mechanism, refracturing in existing perforations may be not the best option. On the one hand, when shale gas is exploited, the gas pressure and total stress near the existing fractures will decrease, so that the fluid in the fractures will encounter less resistance in the depletion zone. It is more likely to generate wider and shorter fractures when refracturing the existing fractures. On the other hand, when new fractures are generated around the existing fractures, the reservoir depletion will also offset the impact of the stress shadow to a certain extent and make it more difficult to generate new fractures than refracturing existing fractures [33]. The new fractures in a refracturing scheme should be evaluated according to the reservoir pressure, historical production data and fracture distribution.

Figure 14 shows the impact of the fracturing length on shale gas production from these three cases. Increasing the half-length of the hydraulic fracture can obviously enhance the cumulative gas production until the fracture half-length reaches a certain length. The effect of increasing the fracture length will be reduced due to the limited total shale gas content in the reservoir. Moreover, a too long hydraulic fracture length will cause serious production decline, which is not conducive to the long-term shale gas extraction. Figure 15 presents the effect of the number of new fractures on the cumulative production. As the number of new fractures increases, the cumulative production of shale gas also increases. The production decline might be more dramatic in cases with a bigger number of fractures. This is due to the limited total shale gas reserve and the interference between the fractures.

4.3. Global Sensitivity Analysis (ANN-GSA) for the Optimization of the Fracture Geometry

4.3.1. Refracturing Cases and Data Processing

This section assumes an ideal refracturing layout to analyze the effects of the fractures at different locations and their interactions on the cumulative production. As shown in Figure 16, fractures are numbered as fracture 1–fracture 9 from left to right. In this refracturing layout, three initial fractures (fracture 3, fracture 5 and fracture 7) will be continuously fractured, and six new fractures will be generated. The range of fracture length is determined, according to the literature [33,34]. This study assumes that the maximum half-length of the fractures is 250 m. The refracturing extended the length of these old fractures in a range from 0 to 130 m. The new fractures have a length between 50 m and 250 m. This length is used as a parameter of the fracture extension in the refracturing layout. The refracturing time ranges from one year to five years after the initial stimulation. The random algorithm with a uniform distribution is used to generate 500 random refracturing schemes. A refracturing scheme includes the refracturing layout and the refracturing time. Then the cumulative productions of these refracturing schemes in 10 years are calculated by the numerical model. The data before the neural network training are normalized as

$$x^{\prime }=\frac{x-\mathrm{mean}(x)}{\max (x)-\min (x)}$$

(51)

4.3.2. Design of the ANN Structure

A neural network with a simple structure can achieve good results. The ANN with a single hidden layer can fit the arbitrary nonlinear relations with arbitrary precision [35]. The structure of the neural network needs to be adjusted before training to achieve an optimal model performance. The input layer has nine neurons and the output layer has one neuron of gas production. Hence, this study optimized the ANN structure by adjusting the number of hidden layers and the hidden layer neurons. In order to evaluate the performance of the models with different structures, the error of the 5-fold cross validation is calculated and the result is shown in

Table 3. MSE value of the ANN models with different structures (×10⁻⁵).

Number of Hidden Layers	5 Hidden Neurons	10 Hidden Neurons	15 Hidden Neurons	20 Hidden Neurons	25 Hidden Neurons	30 Hidden Neurons	35 Hidden Neurons
1	33.4	27.8	20.2	10.1	10.1	9.9	9.7
2	30.6	22.7	14.5	10.5	10.0	9.2	9.3
3	27.4	21.3	13.8	9.7	9.7	9.1	8.7
4	25.6	19.7	13.6	9.6	9.7	8.7	8.3
5	25.8	17.2	13.1	9.7	8.9	8.2	8.5
6	24.9	15.1	12.5	9.4	8.4	7.8	8.1
7	22.5	12.9	10.9	9.2	7.7	7.6	7.6
8	21.2	9.4	9.7	8.6	8.2	7.7	7.2
9	19.6	12.4	10.3	8.7	7.6	7.6	7.2

. The K-fold cross validation is a method for the model performance evaluation. In this approach, all data sets are divided into K groups. The K-1 group are used as the training set, and the remaining ones are used as the testing set. Upon cycling through the K groups, the average error of the model on the testing set is taken as the generalization error. With the increase of the number of hidden layers and hidden neurons in the neural network, the prediction performance of the ANN model will be improved. In addition, the activation function provides the neural network with the ability of nonlinear mapping, and the appropriate activation function is beneficial to improve the prediction performance of the model [36]. Therefore, this study also compared the test errors of the ANN models with different activation functions. As shown in Figure 17, the MSE value of the ANN model with the activation function of ‘relu’ is the lowest. Therefore, in this study, the ‘relu’ is used as the activation function.

The ANN model is trained by 80% of the data and evaluated by 20% of the data. The MSE of the ANN model is 0.00007 on the training set and 0.0001 on the test set. Figure 18 depicts the comparison between the simulations and the model predictions at a 10-year cumulative production for the training set and the testing set. Their difference is presented in Figure 19 for 500 data points. The differences are mostly within ±1 × 10⁶ m³. The ANN model performs very well on both training and test sets. Its accuracy is within an acceptable range. In computation time, each numerical simulation takes 18s and the total time of the numerical simulations for 500 samples is about 1.5 h. The training time of the ANN model is about 3 min. Once well trained, the ANN model can complete a prediction in only a few milliseconds. The ANN model can provide a large number of samples for the global sensitivity analysis of the refracturing optimization and save a lot of computation time. Therefore, the ANN model has the advantages in short time-consuming and a high accuracy and is a good tool for the optimization of the subsequent refracturing design.

4.3.3. Sobol Global Sensitivity Analysis

The Sobol sequence generation algorithm is used to generate enough sets of fracture parameters for the proposed ANN model training. These fracture parameters and the outputs calculated by the ANN model are used as the samples for the Sobol global sensitivity analysis. In our experiment, we found that the calculated sensitivity index is unreliable when the number of samples is less than 10,000. The calculated sensitivity indices tend to be stable when the number of samples is greater than 30,000. In this study, the number of samples is 50,000.

Table 4. Some sensitivity indices varying with sample sizes.

Sensitivity Indices	3000 Samples	5000 Samples	10,000 Samples	30,000 Samples	50,000 Samples	150,000 Samples	500,000 Samples
S1	0.3441	0.3466	0.3500	0.3496	0.3496	0.3495	0.3496
S2	0.0384	0.0384	0.0386	0.0388	0.0387	0.0388	0.0388
S3	0.0797	0.0765	0.0755	0.0770	0.0768	0.0767	0.0767
S4	0.0017	0.0028	0.0024	0.0025	0.0025	0.0025	0.0024
S5	0.0263	0.0270	0.0273	0.0274	0.0272	0.0273	0.0274
S6	0.0038	0.0029	0.0033	0.0034	0.0033	0.0034	0.0033
S7	0.0677	0.0679	0.0715	0.0716	0.0708	0.0714	0.0712
S8	0.0458	0.0439	0.0457	0.0457	0.0451	0.0453	0.0453
S9	0.3513	0.3503	0.3508	0.3515	0.3513	0.3514	0.3513
Sum (Si)	0.9592	0.9567	0.9654	0.9675	0.9655	0.9665	0.9664
Sum (Sij)	0.0844	0.0598	0.0346	0.0334	0.0331	0.0320	0.0318

Where S_i represents the first-order sensitivity index of fracture i; Sum (S_i) is the sum of all first-order sensitivity indices; S_ij represents the second-order sensitivity index of fracture i and fracture j; Sum (S_ij) is the sum of all second-order sensitivity indices.

lists the nine first-order sensitivity indices and the sum of all first and second order sensitivity indices. They vary with sample sizes.

Figure 20 shows the first-order sensitivity indices and the total-order sensitivity indices of the fractures to the cumulative gas production, respectively. They indicate the main impact and the total impact of the fractures on the cumulative gas production. The new fractures far from the center of the fracture network (fracture 1 and fracture 9) contribute more to the shale gas production. The contribution of fracture 4 and fracture 6 to the cumulative gas production can be ignored. This is because the shale gas in the center of the reservoir is gradually depleted during the initial production, and the effect of fracturing new fractures at this location is limited. When the initial fracture spacing is small, fracturing new fractures between the existing fractures is not a good choice. The refracturing effect of the existing fractures (fracture 3, fracture 5 and fracture 7) is better than that of the surrounding new fractures, which is consistent with the results of the numerical simulation analysis. In addition, the total-order sensitivity indices of the fractures are higher than the first-order sensitivity indices, which indicate that the interaction between fractures will affect the production of shale gas.

Following the refracturing operations, the interaction between fractures may heavily affect the shale gas production, thus being not ignorable.

Table 5. Second-order index of the Sobol method.

	Fracture 1	Fracture 2	Fracture 3	Fracture 4	Fracture 5	Fracture 6	Fracture 7	Fracture 8	Fracture 9
Fracture 1		8.40	2.06	0.17	0.35	0.17	0.32	0.10	−0.03
Fracture 2	8.40		3.06	0.22	0.14	0.11	0.10	0.14	0.06
Fracture 3	2.06	3.06		1.03	1.68	0.19	0.60	0.56	0.45
Fracture 4	0.17	0.22	1.03		0.56	0.00	−0.02	−0.03	−0.12
Fracture 5	0.35	0.14	1.68	0.56		0.78	1.00	0.55	0.33
Fracture 6	0.17	0.11	0.19	0.00	0.78		0.70	0.10	0.10
Fracture 7	0.32	0.10	0.60	−0.02	1.00	0.70		3.78	1.29
Fracture 8	0.10	0.14	0.56	−0.03	0.55	0.10	3.78		8.46
Fracture 9	−0.03	0.06	0.45	−0.12	0.33	0.10	1.29	8.46

shows the sensitivity analysis results for the interaction between fractures on the cumulative production of shale gas, which are expressed as second-order indices. In addition, since shale gas reserves are limited, increasing the length of a fracture will reduce the production of its surrounding fractures. This phenomenon is called competitive relationship. The second-order indices also indirectly reflect the competitive relationship. The change of one fracture length has a great influence on the gas production of the other fracture when the second-order indices between two fractures is large. Figure 21 shows the 10 largest second-order indices. For new fractures, the S₁₂, S₂₃, S₇₈ and S₈₉ have the highest values. The interaction between fractures at the edge are more obvious. The remaining amount of shale gas in the fringe part of the fractured reservoir is more than that in the middle of the reservoir, so the fractures at the fringe can obtain more shale gas production, and have a greater impact on the production of the surrounding fractures. For existing fractures, the S₃₅ and S₅₇ have also relatively high values, indicating that the existing fractures have a great influence on each other when the existing fractures are refractured. In addition, Figure 20 shows that the first-order index of fracture 5 is lower than that of fracture 3 and fracture 7. The contribution of fracture 5 is small. Therefore, it is better to control the fracture spacing when refracturing the existing fractures, such as only refracturing fracture 3 and fracture 7. Such a finding can guide us to make reasonable arrangements for the refracturing layout: The fractures should extend to more undeveloped reservoirs for higher gas production. The competitive relationship between fractures can be reduced by increasing the fracture spacing and reducing the number of fractures. An overly dense, small-scale fracture network should be avoided.

4.3.4. Selection of the Optimal Refracturing Time

In the long-term production of shale gas, due to the increasing closure stress and normal creep, the natural fracture network will be obviously closed after the initial stimulation. The selection of the appropriate refracturing time for fracturing new hydraulic fractures and restimulating these closed natural fractures is key to improve the reproduction effect of the shale reservoirs. In order to analyze the effect of the refracturing time on the refracturing performance, the cumulative production of shale gas in 10 years is calculated under the same fracture layout but with a different refracturing time. Figure 22 shows the cumulative gas production after a different refracturing time. It is found that refracturing in the second year has the best performance. If refracturing is carried out six months after the initial fracturing, the production will decline more significantly in the later period. This is because in the early stage of production, the amount of gas in the reservoir is large, and the closure of fractures is small. Refracturing does not significantly enhance its permeability. On the contrary, the creep of the fractures will be more obvious in the subsequent long-term exploitation. Refracturing in the fifth year has the most obvious effect on the restimulation of the reservoir fractures. Although the reservoir permeability has increased significantly, a low shale gas volume will lead to poor shale gas production. As shown in Figure 23, refracturing too early or too late will adversely affect the production of shale gas. An optimal refracturing time should exist to achieve the highest cumulative gas production in a production cycle. For this fracture layout, refracturing in the second to third year is more appropriate.

According to the assumption (2) of Section 4.1, at the same time for refracturing, the effect of different fracture layouts on restoring the permeability of the existing fractures is the same, and the refracturing effect of the different fracture layouts is mainly related to the layout of the hydraulic fractures. The contribution of refracturing to shale gas production is mainly reflected in restoring the permeability of the existing natural and hydraulic fractures when the refracturing effect is small. On the contrary, when the refracturing effect is large, the contribution of refracturing to shale gas production is mainly reflected in the generation of new hydraulic fractures. Figure 23 also shows that refracturing in later stage has a better ability in restoring the permeability of existing fractures. Therefore, it is reasonable to assume that the optimal fracturing time of small-scale refracturing which is mainly used to restore the permeability of existing fractures may be slower than that of large-scale refracturing. In other words, the optimal time for refracturing the same well varies with the fracture network. When the large-scale refracturing is adopted, the refracturing time should be appropriately advanced. In order to explain this point in detail, we quantified the contribution of the refracturing layout to the cumulative gas production to study its relationship with the optimal refracturing time. It is unreasonable to directly define the reproduction efficiency by length and number of the fractures. Therefore, we defined the average reproduction efficiency (ARE) to evaluate the refracturing effects of the different fracture layouts as

$$ARE=\frac{{\displaystyle \underset{t_0}{\overset{t_1}{\int }}y(t,l_1,l_2,\dots ,l_9)dt}-(t_1-t_0)p_u}{(t_1-t_0)p_u}$$

(52)

where $p_u$ is the cumulative gas production in 10 years without refracturing, $y(t,l_1,l_2,\dots ,l_9)$ is the cumulative gas production in 10 years predicted by the ANN model. ARE refers to the average value of the reproduction efficiency at the same fracture layout within the time from $t_0$ to $t_1$. It is an indicator only related to the fracture distribution.

We generated 200 fracture network parameters with a uniform distribution, and calculated their AREs and optimal refracturing time (see Figure 24). The results show that the optimal refracturing time is in the second to fourth year. Different fracture layouts often have different optimal fracturing times. In the early production, shale gas reserve is large and the closure of fractures is not obvious. Large-scale refracturing with a high ARE value can enhance gas production more obviously. If the ARE value of the refracturing layout is too small at this time, the shale gas recovery efficiency will not increase significantly, and the permeability will continue to rapidly decline in the next long-term production, thus the shale gas production performance is poor. In the middle and late stages of the shale gas production, the fractures are significantly closed and the shale gas reserve is relatively small. Although refracturing with a high ARE value can effectively enhance the reservoir permeability, the reproduction performance will be not as good as the early fracturing. In the middle and late stages of production, large-scale refracturing is more costly and less effective than small-scale refracturing due to a low residual reserve. At this time, small-scale refracturing is more suitable.

The above investigations show that the optimal refracturing time is not fixed and should be adjusted, based on the production history and refracturing layout to achieve the maximum cumulative production. This discovery can provide a beneficial reference in predicting the shale gas production performance and refracturing operation. In actual production, an ideal hydraulic fracture layout is difficultly obtained. There are often some alternative refracturing layouts. These refracturing layouts (perforation position, fracture length and so on) are usually determined, based on the realistic engineering geological environment, in-situ stress distribution and the influence from the existing fracture network. Then the combination approach of the numerical simulation and data-driven analysis can be used to further identify the suitable refracturing layout and its optimal refracturing time.

However, the study also has some limitations: Firstly, this study quantifies the scale of fracturing by the length and number of fractures, and the real generated hydraulic fracture network contains more complex features, such as fracture width and tortuosity varying with length. The effect of these complex features on shale gas production after refracturing needs to be further studied. Secondly, how to determine the optimal layout of hydraulic fracture network by a data-driven method is still a problem. The fracturing location and geometry of the hydraulic fractures are often limited by complex in-situ stress conditions and rock physical properties. The results of this study are suitable for selecting the best solution among various feasible schemes, that is, selecting the best fracturing scheme after determining the perforation positions. Third, this proposed numerical model is greatly simplified and only considered the influence of natural fracture creep for the optimization of the refracturing layout and the optimal refracturing time. In future research, the influence of the reservoir viscoelastic-plasticity and temperature change can be comprehensively included for a more perfect shale gas production model.

5. Conclusions

This study investigated the impacts of the refracturing layout and refracturing time on shale gas production through a combination approach with numerical simulations and an artificial neural network. A multiphysical coupling model was established for the numerical simulations on the shale gas production. The influences of the fracture geometry parameters (fracture length and number) on the production were studied through a single factor sensitivity analysis. An ANN model as a proxy sampling model was further developed for the Sobol global sensitivity analysis. The influences of the fracture interactions on shale gas production were explored by a global sensitivity analysis. Finally, the effect of the refracturing time on the cumulative production was analyzed. The following conclusions can be drawn, based on these investigations.

(1)

The ANN model can quickly and accurately predict the cumulative shale gas production. The combination approach of the numerical simulations and the artificial neural network can largely reduce the computation loads in the optimization of the refracturing layout and time. Our examples show that the computation load can be reduced by more than 90%, compared to conventional optimization methods.

(2)

Refracturing the existing fractures and fracturing new fractures from both sides of the existing fractures have a better performance than fracturing new fractures between the existing fractures. Increasing the number and length of hydraulic fractures can significantly enhance the refracturing efficiency, but its effect becomes less obvious with the continuous increase of the number and length of hydraulic fractures.

(3)

The Sobol global sensitivity analysis can identify the impact of the fractures interaction (especially adjacent fractures) on shale gas production. The full extension of the refractures to the undeveloped reservoir zone with reasonable fracture spacing can largely enhance the cumulative shale gas production.

(4)

This study found that an optimal refracturing time for a maximum cumulative gas production in ten years is within 2–4 years after the initial fracturing. The optimal refracturing time depends on the refracturing layout. The fracture layout with a higher average reproduction efficiency generally has an earlier optimal refracturing time.

(5)

The normal creep of natural fractures has a vital impact on the reservoir permeability reduction and refracturing time. Compared with the elastic reservoir model, the gas production and permeability from the viscoelastic reservoir decline more rapidly, but the viscoelastic reservoir productivity increases more obviously after refracturing.