Numerical Investigation of the E ﬀ ect of Partially Propped Fracture Closure on Gas Production in Fractured Shale Reservoirs

: Nonuniform proppant distribution is fairly common in hydraulic fractures, and di ﬀ erent closure behaviors of the propped and unpropped fractures have been observed in lots of physical experiments. However, the modeling of partially propped fracture closure is rarely performed, and its e ﬀ ect on gas production is not well understood as a result of previous studies. In this paper, a fully coupled ﬂuid ﬂow and geomechanics model is developed to simulate partially propped fracture closure, and to examine its e ﬀ ect on gas production in fractured shale reservoirs. Speciﬁcally, an e ﬃ cient hybrid model, which consists of a single porosity model, a multiple porosity model and the embedded discrete fracture model (EDFM), is adopted to model the hydro-mechanical coupling process in fractured shale reservoirs. In ﬂow equations, the Klinkenberg e ﬀ ect is considered in gas apparent permeability, and adsorption / desorption is treated as an additional source term. In the geomechanical domain, the closure behaviors of propped and unpropped fractures are described through two di ﬀ erent constitutive models. Then, a stabilized extended ﬁnite element method (XFEM) iterative formulation, which is based on the polynomial pressure projection (PPP) technique, is developed to simulate a partially propped fracture closure with the consideration of displacement discontinuity at the fracture interfaces. After that, the sequential implicit method is applied to solve the coupled problem, in which the ﬁnite volume method (FVM) and stabilized XFEM are applied to discretize the ﬂow and geomechanics equations, respectively. Finally, the proposed method is validated through some numerical examples, and then it is further used to study the e ﬀ ect of partially propped fracture closures on gas production in 3D fractured shale reservoir simulation models. This work will contribute to a better understanding of the dynamic behaviors of fractured shale reservoirs during gas production, and will provide more realistic production forecasts.


Introduction
In the last decade, the shale gas industry has received worldwide attention [1][2][3]. Horizontal drilling [4] and hydraulic fracturing [5][6][7][8][9][10] are two key techniques for economically improving the production of shale reservoirs. Compared with conventional fracturing techniques, the liquid nitrogen technique [8-10] has its unique advantages for fracturing hard and high-temperature formations. During the hydraulic fracturing treatment, large amounts of fluid and proppants are injected into the shale formation to create propped fractures which act as conductive paths for the fluid flowing from the

Fluid-Governing Equations
In this study, we assume that gas is stored as a free phase in the fractures, while it is stored as both free and adsorbed phases in the shale matrix. The governing equation for isothermal singlephase gas flow is obtained from the mass conservation law [38], and its general integral form is: where n is the unit normal vector of boundary Г, and A, F and q are the mass accumulation, mass flux and source terms on domain Ω, respectively. The mass accumulation is expressed as: where ϕ is the Lagrange porosity, and its definition can be found in references [3,28]. ρg is the gas density, and m is the mass of gas adsorbed per solid volume (only exists for matrix), and it is determined according to the Langmuir's isotherm [39,40] as a function of gas pressure: where p is the gas pressure, M is the gas molar mass, Z is the gas compression factor, R is the universal gas constant, T is the reservoir temperature, ρr is the rock bulk density, ρgstd is the gas density at standard condition, and VL and PL are the Langmuir volume and Langmuir pressure, respectively. The mass flux term in Equation (1) is given by: (4) where k is the absolute permeability, µ is the gas viscosity, g is the gravity acceleration, D is the depth, and b is the Klinkenberg factor, which incorporates the gas-slippage effect in gas-flow models by modifying the absolute permeability as a function of gas pressure [38,41] k k 1 1 4 , , 2 2.81708 where c is the gas compressibility, and Dk is the Knudsen diffusion coefficient for real gas [3].

Fluid-Governing Equations
In this study, we assume that gas is stored as a free phase in the fractures, while it is stored as both free and adsorbed phases in the shale matrix. The governing equation for isothermal single-phase gas flow is obtained from the mass conservation law [38], and its general integral form is: where n is the unit normal vector of boundary Γ, and A, F and q are the mass accumulation, mass flux and source terms on domain Ω, respectively. The mass accumulation is expressed as: where φ is the Lagrange porosity, and its definition can be found in references [3,28]. ρ g is the gas density, and m is the mass of gas adsorbed per solid volume (only exists for matrix), and it is determined according to the Langmuir's isotherm [39,40] as a function of gas pressure: where p is the gas pressure, M is the gas molar mass, Z is the gas compression factor, R is the universal gas constant, T is the reservoir temperature, ρ r is the rock bulk density, ρ gstd is the gas density at standard condition, and V L and P L are the Langmuir volume and Langmuir pressure, respectively. The mass flux term in Equation (1) is given by: where k is the absolute permeability, µ is the gas viscosity, g is the gravity acceleration, D is the depth, and b is the Klinkenberg factor, which incorporates the gas-slippage effect in gas-flow models by modifying the absolute permeability as a function of gas pressure [38,41] 2.81708 k/φ πRT 2M (5) where c is the gas compressibility, and D k is the Knudsen diffusion coefficient for real gas [3].
Energies 2020, 13, 5339 4 of 24 The fluid flow boundary conditions are v·n q = q on Γ q , p = p on Γ p (6) where v is the flow rate, and q and p are the prescribed flow rate and pressure on boundaries Γ q and Γ p , respectively, as shown in Figure 2. n q is the unit normal vector of Γ q .
where ∇ indicates the gradient operator, and the superscript T indicates transpose.
where t and u are the prescribed traction and displacement on boundaries Гt and Гu, respectively, nt is the unit normal vector of Г,; pHF is the gas pressure in the fracture, which is exerted on the inner fracture faces ГHF, HF n is the unit normal vector of ГHF pointing to Ω + with HF HF HF − + = =− n n n , as shown in Figure 2, and ps is the effective stress exerted on the inner fracture faces, which is caused by the compression of proppants (propped) or the contact of fracture faces (unpropped), and its formula will be given in Equations (12) and (13).

Constitutive Relations
In this study, we assume that a hydraulic fracture is represented by two parallel plates, and the closure of the propped hydraulic fractures mainly depends on proppant deformation, thus its general constitutive model is defined as

Geomechanics-Governing Equations
The governing equation for solid deformation under the assumption of quasi-static infinitesimal deformation is given as [42] (the sign convention of continuum mechanics is adopted here): where σ is the total stress tensor, and ρ b is the bulk density. In the hybrid model, the total stress tensor is defined as [3] σ = Cε − αpI, outside the SRV C up ε + l K dr b l p l I, in the SRV (8) where C and α are the elasticity tensor and Biot coefficient of the shale matrix, I is the identity tensor, C up and K dr are the upscaled elasticity tensor and drained bulk modulus at the level of gridblock, b is a coupling coefficient, and the subscript l indicates the l-th sub-gridblock within each matrix gridblock.
Detailed definitions of C up , K dr and b l can be found in reference [3]. The strain tensor ε is given as the symmetric gradient of displacement vector u from infinitesimal deformation: where ∇ indicates the gradient operator, and the superscript T indicates transpose. The geomechanics boundary conditions are σ·n HF = −(p HF + p s )·n HF on Γ HF (11) where t and u are the prescribed traction and displacement on boundaries Γ t and Γ u , respectively, n t is the unit normal vector of Γ; p HF is the gas pressure in the fracture, which is exerted on the inner fracture faces Γ HF , n HF is the unit normal vector of Γ HF pointing to Ω + with n HF = n − HF = −n + HF , as shown in Figure 2, and p s is the effective stress exerted on the inner fracture faces, which is caused by the compression of proppants (propped) or the contact of fracture faces (unpropped), and its formula will be given in Equations (12) and (13).

Constitutive Relations
In this study, we assume that a hydraulic fracture is represented by two parallel plates, and the closure of the propped hydraulic fractures mainly depends on proppant deformation, thus its general constitutive model is defined as where f s indicates the nonlinear stress-strain relationship of proppant deformation, and it can be obtained from mechanical experiments, ε s = − 1 〚 u 2 〛 ·n HF /d HF0 is the normal proppant strain, 〚 u 〛 = u + − u − represents the displacement difference between fracture faces, and d HF0 is the initial fracture aperture.
For the closure of unpropped hydraulic fractures, a constitutive model modified from the self-contact condition and penalty method [43] is written as where d HF = d HF0 + 〚 u 〛 ·n HF is the residual fracture aperture, and E n 1 is a normal penalty parameter. This constitutive model allows for a small interpenetration of the fracture faces because of the finite value of E n [43]. Based on the Kozeny-Carman equation [44,45], the dynamic permeability of matrix k m can be written as where k m and φ m are the matrix permeability and porosity, respectively. The subscript 0 represents the initial state. The dynamic permeability of natural fractures is [3] where K f is the drained bulk modulus of natural fractures, and ε v is the volumetric strain of a gridblock. According to the parallel plate assumption and cubic law [46][47][48], the dynamic permeability of unpropped hydraulic fractures can be written as A number of experimental results have reported that the dynamic permeability of propped hydraulic fractures is related to effective stress [49,50] k HF = f k (p s ) (17) where f k indicates the relationship between fracture permeability and effective stress, and it can be obtained from conductivity experiments [51][52][53]. In this study, only the normal deformation of hydraulic fractures is considered, and the shear slip and related dilation will be included in future publications.

Numerical Schemes
The mixed space discretization is applied in this study. The space discretization for fluid flow is achieved by using the FVM, in which pressure is located at the elements' centers. On the other hand, the stabilized XFEM is employed for geomechanics, where the displacement vectors are located at grid vertices. The geometry discretization is achieved by using orthogonal grids, and hydraulic fractures Energies 2020, 13, 5339 6 of 24 are discretized according to the intersections of hydraulic fractures and orthogonal grids. Besides, the nested grids based on the improved MINC method are generated inside each matrix block in the SRV. More details about the geometry discretization can be found in reference [3].

Flow Equations Discretization
The FVM is applied to discretize the gas flow-governing equation (i.e., Equation (1)), and the standard first-order approximation is used for time discretization, while accumulation and mass flux terms are evaluated fully implicitly. The residual of Equation (1) can be written as: where n is the time level, i and j are element numbers, ∆t is the time step, V is the element volume, G i is the adjacent element set of element I, ij + 1/2 represents the upstream weight value on the interface between elements i and j, λ = 1/µ is the gas mobility, and T ij is the transmissibility between elements i and j, defined as where A ij is the interface area, d i and d j are vertical distances from cell center to the interface, and k ij+1/2 is the harmonic-weighted permeability. To include the Klinkenberg effect on gas flow, the permeability in Equation (19) should be evaluated as (k(1 + b/p)) ij+1/2 . More details about the calculation of T ij in different domains can be found in reference [3]. Finally, Equation (18) can be solved by using the Newton-Raphson method, which iterates until convergence where i denotes the iterative step within one time step.

Geomechanics Equations Discretization
The governing equation for solid deformation is discretized by using the XFEM. The weak form of Equation (7) after incorporating Equations (8) and (11) is given as [44] Ω δε : C a ε dΩ = Ω δε : (21) where C a indicates C and C up for the regions outside SRV and in SRV, respectively, and p a = p in the region outside SRV, while we define p a of a gridblock in SRV as The equal-order linear interpolation of displacement does not satisfy the discrete LBB condition [30], and this leads to the unstable approximation (i.e., displacement oscillation) for Equation (21). Therefore, the PPP technique [30] is used to establish a stabilized XFEM formulation by adding a stabilizing term to Equation (21) Ω δε : C a εdΩ = Ω δε : p a I dΩ + Γ HF Stabilizing term (23) where τ > 0 is a stabilization parameter [3], which typically has a dimensionless value in the order of 1.0, and M is bulk modulus, equal to K dr in this study. Π denotes a projection operator, and its definition can be found in references [3,30].
To represent the displacement discontinuity between hydraulic fracture faces, the displacement field is approximated by adding enrichment functions to standard finite element space [54]. By substituting the enriched displacement field and test function of the displacement field into Equation (23), and satisfying the necessity that weak form should hold for all kinematically admissible test functions, Equation (23) can be written as: with Stabilizing term (26) where u, a, b, and c indicate the standard and enriched DOFs, respectively, f int and f ext indicate internal force vector and external force vector, respectively, m = [1,1,0] T (2D) and m = [1,1,1,0,0,0] T (3D). Because the effective stress p s is nonlinear, Equation (24) would be solved by using the Newton-Raphson method, and its discretized form can be written as where superscript k is the current iteration step, and k-1 is the previous iteration step.

Sequential Implicit Method
To solve the coupled flow and geomechanical model, the sequential implicit method is used because of its stability and flexibility [32,33]. Specifically, the flow model is solved first, then the geomechanics model is solved using the intermediate solution information (gas pressure) obtained from the flow model. This sequential procedure is iterated at each time step until the solution converges to within an acceptable tolerance, and more details can be seen in reference [3]. A flow chart showing the framework for solving the coupled flow and geomechanics model is given in Figure 3.
Energies 2020, 13 Figure 3. Schematic of the fixed stress sequential implicit method.

Model Verifications
To show the accuracy of the proposed method, two numerical examples are presented in this section. Firstly, we test the stabilized XFEM iterative formulation for partially propped fracture deformation. After that, the simulation of gas production in a 3D deformable shale reservoir is validated. The reference results of these two examples are both from the standard finite element model (SFEM), in which hydraulic fractures are explicitly modeled by using refined grids, and the fully coupled fully implicit schemes are used (implemented by COMSOL Multiphysics [55]). In addition, more verifications of the proposed method can be found in our previous studies [3,28].

Partially Propped Fracture Deformation Problem
The model used to test the stabilized XFEM (S-XFEM) formulation for partially propped fracture deformation is shown in Figure 4a. This medium is homogenous (Young's modulus: 1 GPa, Poisson's ratio: 0.2) and contains one vertical fracture (1.68 × 1.68 × 0.005 m). Only one-half of the fracture (below the red line) is filled with proppant, and the nonlinear stress-strain curve of proppant is shown in Figure 4b. A uniform vertical downward constant load (100 MPa) is applied on the top boundary,

Model Verifications
To show the accuracy of the proposed method, two numerical examples are presented in this section. Firstly, we test the stabilized XFEM iterative formulation for partially propped fracture deformation. After that, the simulation of gas production in a 3D deformable shale reservoir is validated. The reference results of these two examples are both from the standard finite element model (SFEM), in which hydraulic fractures are explicitly modeled by using refined grids, and the fully coupled fully implicit schemes are used (implemented by COMSOL Multiphysics [55]). In addition, more verifications of the proposed method can be found in our previous studies [3,28]. The model used to test the stabilized XFEM (S-XFEM) formulation for partially propped fracture deformation is shown in Figure 4a. This medium is homogenous (Young's modulus: 1 GPa, Poisson's ratio: 0.2) and contains one vertical fracture (1.68 × 1.68 × 0.005 m). Only one-half of the fracture (below the red line) is filled with proppant, and the nonlinear stress-strain curve of proppant is shown in Figure 4b. A uniform vertical downward constant load (100 MPa) is applied on the top boundary, while other boundaries are fixed in their normal directions. A displacement observation line (black dotted line) is located at the middle of the fracture. In this model, we simulate pure mechanical deformation, and compare the results of different methods.

1.68m
(a) Schematic of the fractured medium (b) The nonlinear stress-strain curve Computational grids of the SFEM are refined near the fracture, and its number is 56,190. On the other hand, uniform grids (15 × 25 × 25) are used for the S-XFEM. The x-displacements obtained from the SFEM and S-XFEM are shown in Figure 5, and it can be seen that they are qualitatively close. This is because the S-XFEM can accurately simulate a partially propped fracture closure without requiring highly refined grids. Besides, it also shows that the deformation of the upper half of the fracture (unpropped) is greater than that of the lower half (propped). The computational times for simulating this model used by the SFEM and the proposed method are 15 s and 10 s, respectively. The residual fracture apertures of the whole fracture, obtained from the SFEM (black dots) and S-XFEM (surface), are shown in Figure 6, and the excellent agreement between these results can be seen. In addition, Figure 6 also gives the residual fracture apertures along the displacement observation line, which are obtained from different methods (i.e., SFEM, S-XFEM and conventional XFEM). We can see that the result of S-XFEM is in good agreement with that of SFEM, while there is some oscillation in the result from conventional XFEM, especially at the fracture boundary. This is caused by the deficiency in displacement approximation (i.e., it does not satisfy the discrete LBB condition) of the conventional XFEM.  Computational grids of the SFEM are refined near the fracture, and its number is 56,190. On the other hand, uniform grids (15 × 25 × 25) are used for the S-XFEM. The x-displacements obtained from the SFEM and S-XFEM are shown in Figure 5, and it can be seen that they are qualitatively close. This is because the S-XFEM can accurately simulate a partially propped fracture closure without requiring highly refined grids. Besides, it also shows that the deformation of the upper half of the fracture (un-propped) is greater than that of the lower half (propped). The computational times for simulating this model used by the SFEM and the proposed method are 15 s and 10 s, respectively.

1.68m
(a) Schematic of the fractured medium (b) The nonlinear stress-strain curve Computational grids of the SFEM are refined near the fracture, and its number is 56,190. On the other hand, uniform grids (15 × 25 × 25) are used for the S-XFEM. The x-displacements obtained from the SFEM and S-XFEM are shown in Figure 5, and it can be seen that they are qualitatively close. This is because the S-XFEM can accurately simulate a partially propped fracture closure without requiring highly refined grids. Besides, it also shows that the deformation of the upper half of the fracture (unpropped) is greater than that of the lower half (propped). The computational times for simulating this model used by the SFEM and the proposed method are 15 s and 10 s, respectively. The residual fracture apertures of the whole fracture, obtained from the SFEM (black dots) and S-XFEM (surface), are shown in Figure 6, and the excellent agreement between these results can be seen. In addition, Figure 6 also gives the residual fracture apertures along the displacement observation line, which are obtained from different methods (i.e., SFEM, S-XFEM and conventional XFEM). We can see that the result of S-XFEM is in good agreement with that of SFEM, while there is some oscillation in the result from conventional XFEM, especially at the fracture boundary. This is caused by the deficiency in displacement approximation (i.e., it does not satisfy the discrete LBB condition) of the conventional XFEM. The residual fracture apertures of the whole fracture, obtained from the SFEM (black dots) and S-XFEM (surface), are shown in Figure 6, and the excellent agreement between these results can be seen. In addition, Figure 6 also gives the residual fracture apertures along the displacement observation line, which are obtained from different methods (i.e., SFEM, S-XFEM and conventional XFEM). We can see that the result of S-XFEM is in good agreement with that of SFEM, while there is some oscillation in the result from conventional XFEM, especially at the fracture boundary. This is caused by the deficiency in displacement approximation (i.e., it does not satisfy the discrete LBB condition) of the conventional XFEM.

Gas Production and Stress Evolution in a 3D Fractured Shale Reservoir
As shown in Figure 7a, a 3D model is designed to validate the proposed method for predicting gas production and stress evolution in fractured shale reservoirs. A hydraulic fracture filled with proppant is at the center of this reservoir. For simplicity, this reservoir is assumed to be a dualporosity and single-permeability medium, and the proppant deformation is assumed to be linear elastic. The production well is located at the center of the hydraulic fracture. For flow, we have no flow at the boundaries. For geomechanics, the uniform constant force (25 MPa) is applied on the top, right and back boundaries, and the other boundaries are fixed in their normal direction, respectively. The stress-dependent normalized fracture conductivity is shown in Figure 8, and other model properties are listed in Table 1. Computational grids for the SFEM (169,082) and the proposed method (6417) are shown in Figure 7b,c, respectively.  Four different methods (i.e., SFEM, the proposed method, and the simplified methods 1 and 2) are applied to simulate the depletion production. The difference between the proposed method and simplified method 1 is that the displacement discontinuity at the hydraulic fracture was neglected in simplified method 1, while in simplified method 2, only the gas flow is simulated, and the geomechanical effects are simplified by correlating the fracture conductivity to pressure, which is based on the assumption that the total stress acting on the hydraulic fracture remains constant during depletion production. To make a better comparison of these methods, the permeabilities of the matrix

Gas Production and Stress Evolution in a 3D Fractured Shale Reservoir
As shown in Figure 7a, a 3D model is designed to validate the proposed method for predicting gas production and stress evolution in fractured shale reservoirs. A hydraulic fracture filled with proppant is at the center of this reservoir. For simplicity, this reservoir is assumed to be a dual-porosity and single-permeability medium, and the proppant deformation is assumed to be linear elastic. The production well is located at the center of the hydraulic fracture. For flow, we have no flow at the boundaries. For geomechanics, the uniform constant force (25 MPa) is applied on the top, right and back boundaries, and the other boundaries are fixed in their normal direction, respectively. The stress-dependent normalized fracture conductivity is shown in Figure 8, and other model properties are listed in Table 1. Computational grids for the SFEM (169,082) and the proposed method (6417) are shown in Figure 7b,c, respectively.

Gas Production and Stress Evolution in a 3D Fractured Shale Reservoir
As shown in Figure 7a, a 3D model is designed to validate the proposed method for predicting gas production and stress evolution in fractured shale reservoirs. A hydraulic fracture filled with proppant is at the center of this reservoir. For simplicity, this reservoir is assumed to be a dualporosity and single-permeability medium, and the proppant deformation is assumed to be linear elastic. The production well is located at the center of the hydraulic fracture. For flow, we have no flow at the boundaries. For geomechanics, the uniform constant force (25 MPa) is applied on the top, right and back boundaries, and the other boundaries are fixed in their normal direction, respectively. The stress-dependent normalized fracture conductivity is shown in Figure 8, and other model properties are listed in Table 1. Computational grids for the SFEM (169,082) and the proposed method (6417) are shown in Figure 7b,c, respectively.  Four different methods (i.e., SFEM, the proposed method, and the simplified methods 1 and 2) are applied to simulate the depletion production. The difference between the proposed method and simplified method 1 is that the displacement discontinuity at the hydraulic fracture was neglected in simplified method 1, while in simplified method 2, only the gas flow is simulated, and the geomechanical effects are simplified by correlating the fracture conductivity to pressure, which is based on the assumption that the total stress acting on the hydraulic fracture remains constant during depletion production. To make a better comparison of these methods, the permeabilities of the matrix

Name
Value Initial absolute permeability of matrix, m 2 2 · 10 −20 Initial absolute permeability of natural fracture, m 2 1 · 10 −17 Initial absolute permeability of hydraulic fracture, m 2 1 · 10 −17 Natural fracture spacing, m 1.0 Initial aperture of natural fracture and hydraulic fracture, m 5 · 10 −6 , 0.01 Initial porosity of matrix, natural fracture and hydraulic fracture 0.05, 1.0, 0. The pressure, x-stress (σxx) and x-displacement (ux) fields after 1000 days, obtained from different methods, are shown in Figure 9, and it can be seen that the results of the proposed method and SFEM are qualitatively close, while there are some differences between the results of simplified method 1 and SFEM, especially near the hydraulic fracture in the x-stress field, and this is caused by the neglect of displacement discontinuity. The computational times used by the SFEM and the proposed method are 220 s and 95 s, respectively. It can also be seen that the computational performance of the proposed method is better than that of the SFEM.

Name Value
Initial absolute permeability of matrix, m 2 2 × 10 −20 Initial absolute permeability of natural fracture, m 2 1 × 10 −17 Initial absolute permeability of hydraulic fracture, m 2 1 Four different methods (i.e., SFEM, the proposed method, and the simplified methods 1 and 2) are applied to simulate the depletion production. The difference between the proposed method and simplified method 1 is that the displacement discontinuity at the hydraulic fracture was neglected in simplified method 1, while in simplified method 2, only the gas flow is simulated, and the geomechanical effects are simplified by correlating the fracture conductivity to pressure, which is based on the assumption that the total stress acting on the hydraulic fracture remains constant during depletion production. To make a better comparison of these methods, the permeabilities of the matrix and natural fractures are assumed to be constant.
The pressure, x-stress (σ xx ) and x-displacement (u x ) fields after 1000 days, obtained from different methods, are shown in Figure 9, and it can be seen that the results of the proposed method and SFEM are qualitatively close, while there are some differences between the results of simplified method 1 and SFEM, especially near the hydraulic fracture in the x-stress field, and this is caused by the neglect of displacement discontinuity. The computational times used by the SFEM and the proposed method are 220 s and 95 s, respectively. It can also be seen that the computational performance of the proposed method is better than that of the SFEM.  Comparisons of the cumulative gas and the evolution of σxx at the observation point (100, 215, 30) located on the hydraulic fracture are plotted in Figure 10. They show that the cumulative gas of simplified method 2 is significantly lower than that of the reference. This is because simplified method 2 overestimates the total stress acting on the hydraulic fracture, as shown in Figure 10b. The results of simplified method 1 show certain improvements, but there are still some differences compared with the references because of its neglect of displacement discontinuity. On the other hand, the results of the proposed method are almost the same as those of the references (the average errors of cumulative gas and the evolution of σxx are 2.26% and 0.5%, respectively). Therefore, the Comparisons of the cumulative gas and the evolution of σ xx at the observation point (100, 215, 30) located on the hydraulic fracture are plotted in Figure 10. They show that the cumulative gas of simplified method 2 is significantly lower than that of the reference. This is because simplified method 2 overestimates the total stress acting on the hydraulic fracture, as shown in Figure 10b. The results of simplified method 1 show certain improvements, but there are still some differences compared with the references because of its neglect of displacement discontinuity. On the other hand, the results of the proposed method are almost the same as those of the references (the average errors of cumulative gas and the evolution of σ xx are 2.26% and 0.5%, respectively). Therefore, the geomechanical and displacement discontinuity have significant impacts on gas production in shale reservoirs, and the proposed method can accurately simulate their effects.
Energies 2020, 13, x 13 of 23 geomechanical and displacement discontinuity have significant impacts on gas production in shale reservoirs, and the proposed method can accurately simulate their effects.
(a) Cumulative gas (b) Observation point total stress

Applications and Discussions
In this section, a series of case studies is conducted to study the effect of partially propped fracture closure on gas production in shale reservoirs. In order to reduce the simulation time, only one single stage (120 × 250 × 30 m) of a multi-stage hydraulic fractured shale reservoir is modeled, as shown in Figure 11a Table 2. The stress-dependent normalized fracture conductivity is shown in Figure 8, and the nonlinear stress-strain curve of proppant deformation is shown in Figure 12. Unless otherwise noted, these parameters and stress-strain curve are used for the following simulations. The dynamic permeabilities of the matrix, natural fracture and hydraulic fracture are calculated by using Equations (14)- (17). The computational grids used for this model are shown in Figure 11b. Please note that the surrounding rock (impermeable) is simulated for the better consideration of geomechanical effects in this model. For convenience, only the reservoir part is plotted in the following analysis.

Applications and Discussions
In this section, a series of case studies is conducted to study the effect of partially propped fracture closure on gas production in shale reservoirs. In order to reduce the simulation time, only one single stage (120 × 250 × 30 m) of a multi-stage hydraulic fractured shale reservoir is modeled, as shown in Figure 11a Table 2. The stress-dependent normalized fracture conductivity is shown in Figure 8, and the nonlinear stress-strain curve of proppant deformation is shown in Figure 12. Unless otherwise noted, these parameters and stress-strain curve are used for the following simulations. The dynamic permeabilities of the matrix, natural fracture and hydraulic fracture are calculated by using Equations (14)- (17). The computational grids used for this model are shown in Figure 11b. Please note that the surrounding rock (impermeable) is simulated for the better consideration of geomechanical effects in this model. For convenience, only the reservoir part is plotted in the following analysis.
Energies 2020, 13, x 13 of 23 geomechanical and displacement discontinuity have significant impacts on gas production in shale reservoirs, and the proposed method can accurately simulate their effects.
(a) Cumulative gas (b) Observation point total stress

Applications and Discussions
In this section, a series of case studies is conducted to study the effect of partially propped fracture closure on gas production in shale reservoirs. In order to reduce the simulation time, only one single stage (120 × 250 × 30 m) of a multi-stage hydraulic fractured shale reservoir is modeled, as shown in Figure 11a Table 2. The stress-dependent normalized fracture conductivity is shown in Figure 8, and the nonlinear stress-strain curve of proppant deformation is shown in Figure 12. Unless otherwise noted, these parameters and stress-strain curve are used for the following simulations. The dynamic permeabilities of the matrix, natural fracture and hydraulic fracture are calculated by using Equations (14)- (17). The computational grids used for this model are shown in Figure 11b. Please note that the surrounding rock (impermeable) is simulated for the better consideration of geomechanical effects in this model. For convenience, only the reservoir part is plotted in the following analysis.

Effect of initial Aperture Distribution
We first study the effect of initial aperture distribution on gas production through the following two different cases, as shown in Figure 13: (a) Case 1 with a more realistic aperture distribution obtained from geomechanical calculation, in which the initial pressure inside the hydraulic fracture is set as 42 MPa to mimic the fracturing process; (b) Case 2 with a uniform aperture of 0.005 m, which means Case 2 has the same cumulative gas as Case 1 when the geomechanical effect is not considered. Figure 14 illustrates the comparison of σxx distribution between Case 1 and Case 2 after 10 years. The results of cumulative gas and bottomhole closure stress are provided in Figure 15. We can see that the cumulative gas of Case 2 is slightly lower than that of Case 1. This is because that the aperture distribution has some effects on the total stress around the hydraulic fracture, and this results in a higher closure stress (acting on the hydraulic fracture) in Case 2, as shown in Figure 15b. Because the initial aperture distribution in Case 1 is more realistic than that in Case 2, it will be used in the following simulations.

Effect of initial Aperture Distribution
We first study the effect of initial aperture distribution on gas production through the following two different cases, as shown in Figure 13: (a) Case 1 with a more realistic aperture distribution obtained from geomechanical calculation, in which the initial pressure inside the hydraulic fracture is set as 42 MPa to mimic the fracturing process; (b) Case 2 with a uniform aperture of 0.005 m, which means Case 2 has the same cumulative gas as Case 1 when the geomechanical effect is not considered. Figure 14 illustrates the comparison of σ xx distribution between Case 1 and Case 2 after 10 years. The results of cumulative gas and bottomhole closure stress are provided in Figure 15. We can see that the cumulative gas of Case 2 is slightly lower than that of Case 1. This is because that the aperture distribution has some effects on the total stress around the hydraulic fracture, and this results in a higher closure stress (acting on the hydraulic fracture) in Case 2, as shown in Figure 15b. Because the initial aperture distribution in Case 1 is more realistic than that in Case 2, it will be used in the following simulations.

Effect of Propped Fracture Length
As shown in Figure 16, three cases, with propped fracture lengths of 140 m (Case 3), 100m (Case 4) and 60 m (Case 5), along with Case 1 (180 m), are used to investigate the effect of propped fracture length on gas production. Figure 17 illustrates the comparison of σxx distribution among the four cases after 10 years of depletion. We can see that the σxx around the propped fracture is slightly higher than that around the unpropped fracture. This is because the propped fracture part is stiffer than the unpropped part, which must suffer higher stress to support the boundary force. Figure 18 shows the comparisons of cumulative gas and production rate among the different cases. It can be seen that gas production correlates positively with the propped fracture length because of the increase in the contact area between the high-conductivity propped fracture and the matrix. Therefore, the

Effect of Propped Fracture Length
As shown in Figure 16, three cases, with propped fracture lengths of 140 m (Case 3), 100m (Case 4) and 60 m (Case 5), along with Case 1 (180 m), are used to investigate the effect of propped fracture length on gas production. Figure 17 illustrates the comparison of σxx distribution among the four cases after 10 years of depletion. We can see that the σxx around the propped fracture is slightly higher than that around the unpropped fracture. This is because the propped fracture part is stiffer than the unpropped part, which must suffer higher stress to support the boundary force. Figure 18 shows the comparisons of cumulative gas and production rate among the different cases. It can be seen that gas production correlates positively with the propped fracture length because of the increase in the contact area between the high-conductivity propped fracture and the matrix. Therefore, the

Effect of Propped Fracture Length
As shown in Figure 16, three cases, with propped fracture lengths of 140 m (Case 3), 100m (Case 4) and 60 m (Case 5), along with Case 1 (180 m), are used to investigate the effect of propped fracture length on gas production. Figure 17 illustrates the comparison of σxx distribution among the four cases after 10 years of depletion. We can see that the σxx around the propped fracture is slightly higher than that around the unpropped fracture. This is because the propped fracture part is stiffer than the unpropped part, which must suffer higher stress to support the boundary force. Figure 18 shows the comparisons of cumulative gas and production rate among the different cases. It can be seen that gas production correlates positively with the propped fracture length because of the increase in the contact area between the high-conductivity propped fracture and the matrix. Therefore, the

Effect of Propped Fracture Length
As shown in Figure 16, three cases, with propped fracture lengths of 140 m (Case 3), 100m (Case 4) and 60 m (Case 5), along with Case 1 (180 m), are used to investigate the effect of propped fracture length on gas production. Figure 17 illustrates the comparison of σ xx distribution among the four cases after 10 years of depletion. We can see that the σ xx around the propped fracture is slightly higher than that around the unpropped fracture. This is because the propped fracture part is stiffer than the unpropped part, which must suffer higher stress to support the boundary force. Figure 18 shows the comparisons of cumulative gas and production rate among the different cases. It can be seen that gas production correlates positively with the propped fracture length because of the increase in the contact area between the high-conductivity propped fracture and the matrix. Therefore, the completion design should strive to increase the length of the propped fracture so as to produce more gas.
Energies 2020, 13, x 16 of 23 completion design should strive to increase the length of the propped fracture so as to produce more gas. (a) Cumulative gas (b) Production rate Figure 18. Comparisons of cumulative gas and production rate obtained from different cases.

Effect of Propped Fracture Height
To analyze the effect of propped fracture height on gas production, four cases with different propped fracture heights (Case 6: 15 m, Case 7: 10.5 m, Case 8: 9.5 m and Case 9: 5 m) are designed, as shown in Figure 19. It should be noted that the wellbore is placed on the propped fracture in Case 6 and Case 7, but on the unpropped fracture in Case 8 and Case 9.
The comparisons of σxx distribution and conductivity distribution among different cases after 10 years are illustrated in Figures 20 and 21, respectively. The results of cumulative gas and production rate for different cases are provided in Figure 22. We can see that in Case 6, Case 7 and Case 9, gas production correlates positively with the propped fracture height. However, there is little difference between the gas production of Case 6 and Case 7, because the wells of both cases are located on the propped part of the hydraulic fracture, which leads to a large well production index. However, the completion design should strive to increase the length of the propped fracture so as to produce more gas. (a) Cumulative gas (b) Production rate Figure 18. Comparisons of cumulative gas and production rate obtained from different cases.

Effect of Propped Fracture Height
To analyze the effect of propped fracture height on gas production, four cases with different propped fracture heights (Case 6: 15 m, Case 7: 10.5 m, Case 8: 9.5 m and Case 9: 5 m) are designed, as shown in Figure 19. It should be noted that the wellbore is placed on the propped fracture in Case 6 and Case 7, but on the unpropped fracture in Case 8 and Case 9.
The comparisons of σxx distribution and conductivity distribution among different cases after 10 years are illustrated in Figures 20 and 21, respectively. The results of cumulative gas and production rate for different cases are provided in Figure 22. We can see that in Case 6, Case 7 and Case 9, gas production correlates positively with the propped fracture height. However, there is little difference between the gas production of Case 6 and Case 7, because the wells of both cases are located on the propped part of the hydraulic fracture, which leads to a large well production index. However, the completion design should strive to increase the length of the propped fracture so as to produce more gas. (a) Cumulative gas (b) Production rate Figure 18. Comparisons of cumulative gas and production rate obtained from different cases.

Effect of Propped Fracture Height
To analyze the effect of propped fracture height on gas production, four cases with different propped fracture heights (Case 6: 15 m, Case 7: 10.5 m, Case 8: 9.5 m and Case 9: 5 m) are designed, as shown in Figure 19. It should be noted that the wellbore is placed on the propped fracture in Case 6 and Case 7, but on the unpropped fracture in Case 8 and Case 9.
The comparisons of σxx distribution and conductivity distribution among different cases after 10 years are illustrated in Figures 20 and 21, respectively. The results of cumulative gas and production rate for different cases are provided in Figure 22. We can see that in Case 6, Case 7 and Case 9, gas production correlates positively with the propped fracture height. However, there is little difference between the gas production of Case 6 and Case 7, because the wells of both cases are located on the propped part of the hydraulic fracture, which leads to a large well production index. However, the

Effect of Propped Fracture Height
To analyze the effect of propped fracture height on gas production, four cases with different propped fracture heights (Case 6: 15 m, Case 7: 10.5 m, Case 8: 9.5 m and Case 9: 5 m) are designed, as shown in Figure 19. It should be noted that the wellbore is placed on the propped fracture in Case 6 and Case 7, but on the unpropped fracture in Case 8 and Case 9.
conductivity of the unpropped fracture where the wellbore is located. Another interesting observation is that although the propped fracture height of Case 8 is less than that of Case 6 and Case 7, the gas production of Case 8 is higher. This is because the wellbore of Case 8 is placed on a high conductivity arch, which is usually formed at the top of the proppant bank, as shown in Figure 21. The fracture at this area may remain open throughout the production, which leads to an extremely large well production index.  The comparisons of σ xx distribution and conductivity distribution among different cases after 10 years are illustrated in Figures 20 and 21, respectively. The results of cumulative gas and production rate for different cases are provided in Figure 22. We can see that in Case 6, Case 7 and Case 9, gas production correlates positively with the propped fracture height. However, there is little difference between the gas production of Case 6 and Case 7, because the wells of both cases are located on the propped part of the hydraulic fracture, which leads to a large well production index. However, the gas production of Case 9 is much lower than that of Case 6 and Case 7 because of the limited conductivity of the unpropped fracture where the wellbore is located. Another interesting observation is that although the propped fracture height of Case 8 is less than that of Case 6 and Case 7, the gas production of Case 8 is higher. This is because the wellbore of Case 8 is placed on a high conductivity arch, which is usually formed at the top of the proppant bank, as shown in Figure 21. The fracture at this area may remain open throughout the production, which leads to an extremely large well production index.
Energies 2020, 13, x 17 of 23 gas production of Case 9 is much lower than that of Case 6 and Case 7 because of the limited conductivity of the unpropped fracture where the wellbore is located. Another interesting observation is that although the propped fracture height of Case 8 is less than that of Case 6 and Case 7, the gas production of Case 8 is higher. This is because the wellbore of Case 8 is placed on a high conductivity arch, which is usually formed at the top of the proppant bank, as shown in Figure 21.
The fracture at this area may remain open throughout the production, which leads to an extremely large well production index.  Energies 2020, 13, x 17 of 23 gas production of Case 9 is much lower than that of Case 6 and Case 7 because of the limited conductivity of the unpropped fracture where the wellbore is located. Another interesting observation is that although the propped fracture height of Case 8 is less than that of Case 6 and Case 7, the gas production of Case 8 is higher. This is because the wellbore of Case 8 is placed on a high conductivity arch, which is usually formed at the top of the proppant bank, as shown in Figure 21.
The fracture at this area may remain open throughout the production, which leads to an extremely large well production index.

Effect of Proppant Distribution
Recently, the alternating (i.e., pulsed) injection of proppant-free fluids and proppant-containing fluids has been used to form heterogeneous proppant distributions (i.e., proppant pillars) [56,57]. The proppant pillars may withstand the closure stress to prevent the unpropped fracture between them being closed, and highly conductive flow channels can be generated. As shown in Figure 23, four cases (Cases 10-13) with different proppant distribution patterns are constructed to analyze the effects of proppant distribution on gas production.

Effect of Proppant Distribution
Recently, the alternating (i.e., pulsed) injection of proppant-free fluids and proppant-containing fluids has been used to form heterogeneous proppant distributions (i.e., proppant pillars) [56,57]. The proppant pillars may withstand the closure stress to prevent the unpropped fracture between them being closed, and highly conductive flow channels can be generated. As shown in Figure 23, four cases (Cases 10-13) with different proppant distribution patterns are constructed to analyze the effects of proppant distribution on gas production.
The comparisons of σxx distribution and conductivity distribution among different cases after 10 years are illustrated in Figures 24 and 25, respectively. We can see that in the first two cases, the highly conductive flow channels between proppant pillars are not connected, while they are connected in the other two cases. The results of cumulative gas and production rate for different cases are provided in Figure 26. It can be seen from the results of Case 10 and Case 11 that cumulative gas decreases as the number of flow channels increases. This is because gas flow in the fracture is mainly controlled by proppant pillars when the flow channels are not connected. The proppant pillars will withstand more closure stress as the flow channels increase, leading to the lower conductivity of the proppant pillars. Compared with Case 10 and Case 11, the cumulative gas is much higher in Case 12 and Case 13. However, the gas production curves are nearly overlapped for Case 12 and Case 13. This is because when flow channels are connected, the gas flow in the fracture is mainly controlled by the flow channels, and the fracture can be regarded as infinitely conductive in this situation. Therefore, improving the connectivity between flow channels is more important than increasing the number of flow channels.   The comparisons of σ xx distribution and conductivity distribution among different cases after 10 years are illustrated in Figures 24 and 25, respectively. We can see that in the first two cases, the highly conductive flow channels between proppant pillars are not connected, while they are connected in the other two cases. The results of cumulative gas and production rate for different cases are provided in Figure 26. It can be seen from the results of Case 10 and Case 11 that cumulative gas decreases as the number of flow channels increases. This is because gas flow in the fracture is mainly controlled by proppant pillars when the flow channels are not connected. The proppant pillars will withstand more closure stress as the flow channels increase, leading to the lower conductivity of the proppant pillars. Compared with Case 10 and Case 11, the cumulative gas is much higher in Case 12 and Case 13. However, the gas production curves are nearly overlapped for Case 12 and Case 13. This is because when flow channels are connected, the gas flow in the fracture is mainly controlled by the flow channels, and the fracture can be regarded as infinitely conductive in this situation. Therefore, improving the connectivity between flow channels is more important than increasing the number of flow channels.    (a) Cumulative gas (b) Production rate Figure 26. Comparisons of cumulative gas and production rate obtained from different cases.

Conclusions
The main contribution of this study is that we developed an efficient hybrid model, which includes the effect of each hydraulic fracture explicitly, with non-matching grids, and represents the effect of natural fractures implicitly, to simulate the closure of partially propped fractures and to examine the effect of this on gas production in fractured shale reservoirs. Besides, a stabilized XFEM iterative formulation based on the PPP technique is developed to accurately simulate a partially propped fracture closure with consideration of displacement discontinuity at fracture interfaces, and the modified fixed-stress sequential implicit method is applied to solve the coupled problem. We have validated our proposed method with the SFEM through some numerical examples, and the following conclusions are obtained from the case studies: 1.
The displacement discontinuity at hydraulic fractures and the initial aperture distribution of the hydraulic fractures have significant impacts on the evolution of stress near hydraulic fractures, thus a good performance can be predicted more accurately when these effects are considered; 2.
Gas production is positively correlated with the propped fracture length and height, and the connectivity between the well bore and the high conductivity arch plays an important role in production; 3.
Proppant distribution can significantly affect well performance. A fracture with alternating

Conclusions
The main contribution of this study is that we developed an efficient hybrid model, which includes the effect of each hydraulic fracture explicitly, with non-matching grids, and represents the effect of natural fractures implicitly, to simulate the closure of partially propped fractures and to examine the effect of this on gas production in fractured shale reservoirs. Besides, a stabilized XFEM iterative formulation based on the PPP technique is developed to accurately simulate a partially propped fracture closure with consideration of displacement discontinuity at fracture interfaces, and the modified fixed-stress sequential implicit method is applied to solve the coupled problem. We have validated our proposed method with the SFEM through some numerical examples, and the following conclusions are obtained from the case studies:

1.
The displacement discontinuity at hydraulic fractures and the initial aperture distribution of the hydraulic fractures have significant impacts on the evolution of stress near hydraulic fractures, thus a good performance can be predicted more accurately when these effects are considered; 2.
Gas production is positively correlated with the propped fracture length and height, and the connectivity between the well bore and the high conductivity arch plays an important role in production; 3.
Proppant distribution can significantly affect well performance. A fracture with alternating propped-unpropped-propped sections is favored because the highly conductive flow channels can be generated. In this fracture, improving the connectivity between flow channels is more important than increasing the number of flow channels.

Conflicts of Interest:
The authors declare no conflict of interest.