A Fully Three Dimensional Semianalytical Model for Shale Gas Reservoirs with Hydraulic Fractures

: Two challenges exist for modeling gas transport in shale. One is the existence of complex gas transport mechanisms, and the other is the impact of hydraulic fracture networks. In this study, a truly three dimensional semianalytical model was developed for shale gas reservoirs with hydraulic fractures of various shapes. Using the instantaneous point source solution, the pressure are solved for a bounded reservoir with fully 3D, partially penetrated hydraulic fractures of different strike angles and dip angles. The fractures could have various shapes such as rectangles, disks and ellipses. The shale gas diffusion equations considers complex transport mechanism such as gas slippage and gas diffusion. This semianalytical model is veriﬁed with a commercial software and an analytical method for single fully penetrated rectangle fracture, and the production results of shale gas are consistent. The impacts of fracture height and strike angles are investigated by ﬁve systematically constructed models. The comparison shows that the production increases proportionally with the fracture height, and decreases with the increase of strike angles. The method proposed in this study could also be applied in well testing to analyze the reservoir properties and used to forecast the production for tight oil and conventional resources.


Introduction
Shale formations are characterized by low permeability, nanopores, high total organic carbon, and natural fractures. In shale, the gas content is mainly composed of free gas and adsorbed gas. Because of its low permeability, shale gas did not gain much success until the application of a horizontal well drilling and hydraulic fracturing technique. Currently, shale gas accounts for about half of the natural gas produced in the United States of America (U.S.A.) and that ratio will increase to 70% in 2040 according to a prediction by the U.S. Energy Information Agency [1,2] To simulate the gas production process in shales, many complex gas transport mechanisms such as gas slippage and desorption and the presence of natural fractures need to be considered, which makes the simulation of shale gas production a challenging work.
There are mainly three gas coexisting transport mechanisms during the gas production process, including gas diffusion, gas slippage and gas desorption [3]. The complexity of these mechanisms makes the classical free gas diffusivity equations inadequate to model the gas transport in shales. Despite the availability of some existing models for simulation of gas transport in shales [4,5], most of them do not fully capture the complex mechanisms, and a comprehensive gas transport equation is still needed.
Another challenge presented when dealing with shale gas transport simulation is the modeling of complex fracture patterns. The dual-porosity or dual permeability model [6,7] is widely used in numerical simulations of the fractures, but the model is usually over-simplified compared to the actual field complexities, and the resolution of this approach is not enough to capture the details of fractures. In that respect, discrete-fracture models (DFM), using finite-volume or finite-difference methods, have been developed [8]. To be compatible with the complex geometries of fractures, such as nonplanar fractures and fractures with variable aperture, reservoir simulators using unstructured grids have also been designed [9][10][11][12] in order to explicitly model the fractures. However, the use of unstructured grids leads to high computational cost so that it is only used in limited applications in real field studies. As a solution, a technique that directly incorporates complex fractures in a conventional structured grid was developed [13,14], and Moinfar et al. [15] named it embedded discrete fracture model (EDFM). The EDFM was further developed for 3D reservoir simulation by Moinfar et al. [16]. Recently, the EDFM has been widely used for naturally and hydraulically fractured reservoirs [17][18][19][20]. Despite the accurate capture of the fracture geometries, the EDFM method needs extra computation time to calculate the extra transmissibility caused by the non-neighborhood connections between fractures and matrix cells.
Assuming the rock properties are homogeneous, there are many analytical and semianalytical models derived to efficiently simulate the gas transport in shales. Gringarten et al. [21] used Green's functions to solve unsteady flow problems for a wide variety of reservoir flow problems. Cinco-Lay and Samaniego [22] presented a new technique for analyzing pressure transient data for wells intercepted by a finite conductivity vertical fracture. Wan and Aziz [23] described a new semianalytical solution for horizontal wells with multiple fractures with different strike angles and fully/partially penetrating the formation in the vertical direction. With that solution, the authors calculated the well index combined with a numerically computed gridblock pressure. Lin and Zhu [24] studied the well performance for fractured horizontal wells in an infinite slab reservoir, by using the instantaneous solutions of plane sources. Zhou et al. [25] designed a semianalytical model to simulate the gas production in Barnett shales with fully penetrated planar fractures. Yu et al. [26] extended Zhou's model to simulate well production for reservoirs with fully penetrated non-planar fractures, which have varying fracture width and fracture permeability. Their model was verified with a numerical reservoir simulator for single fracture case before being applied to several case studies based on Marcellus shale properties. Yang et al. [27] developed a robust and comprehensive model for real gas transport in shales with complex non-planar fracture network, with gas diffusion and gas slippage effect considered. He et al. [28] proposed a semianalytical method to diagnose the locations of underperforming fully penetrated hydraulic fractures through pressure transient analysis in tight gas reservoirs. Wang et al. [29] constructed a semianalytical model for simulating shale gas reservoirs with complex fully penetrated fractures.
In most of the above studies, due to the complexity of deriving the explicit pressure solutions for partially penetrated fractured in a bounded reservoir, the fractures are assumed to be fully penetrated rectangle shapes in a slab reservoir. Because coupling pressure and flow rate could bring an extra difficulty, most of the previous work concentrate on the pressure profile without being concerned about the production. In real situations, the fractures are most likely partially penetrated and have various shapes, such as disk and ellipse [30][31][32]. In this study, a novel semianalytical model is designed for bounded reservoir with multiple fully or partially penetrated fractures with various shapes such as rectangles, disks, and ellipses. First, the pressure solutions for fractures with dip angles or strike angles, and various shapes of fractures are derived, including rectangle, ellipse and disk shapes. Based on the pressure solutions and the mathematical equation developed for shale gas reservoirs with proper assumptions, the semianalytical model is introduced and described. The new model is verified against a commercial reservoir simulator and an analytical model for a bounded reservoir with a horizontal well and three fully penetrated fractures. In order to investigate the impact of fracture heights and strike angles on gas production, this model is then used to model the shale gas production for five bounded reservoirs with different hydraulic fracture patterns, varying from fully penetrated to partially penetrated, from fractures without strike angles to fractures with different strike angles. To the best of the authors' knowledge, this study is the first to design a semianalytical model for shale gas reservoirs with partially penetrated fractures of various shapes. The methodology proposed in this study could be expanded and implemented in existing reservoir simulation tools to investigate various fracture patterns.

Model Development
For gas transport from shale matrix into fracture segment, the diffusivity equation is revised for conventional gas reservoirs by considering the important gas transport mechanisms such as gas slippage, gas diffusion, and gas desorption as follows: where ρ g is gas density, φ is rock porosity, k is reservoir permeability, C g is the isothermal gas compressibility factor, p is pressure, µ g is gas viscosity, V a is pore volume fraction of adsorbed gas, K n is Knudsen number, α is a constant and close to 1, 8αK n is used to describe the gas slippage effect [33], D g is the Fickian diffusivity of gas component through the pore and used to describe the gas diffusion process, which might be the combination of three distinct mechanisms acting individually or simultaneously: bulk diffusion, Knudsen diffusion and surface diffusion [34]. δ is a dimensionless constrictivity factor (≤1), which accounts for variation of the pore size along its length caused by small pores and can be calculated by the empirical equation developed by Satterfield et al. [35]: where d gas is gas molecule diameter (methane molecular diameter is 0.38 nm) and d pore is the pore diameter, τ is a dimensionless tortuosity factor (≥1), which can be estimated from porosity and gas saturation by [36]: where S g is gas saturation, K a is the differential equilibrium partitioning coefficient of gas at a constant temperature, which is used to describe the gas desorption mechanism and defined as [37]: where ρ a is the adsorbed gas mass per unit shale volume (kilograms of adsorbed gas per cubic meter of solid). A general form of Brunauer-Emmett-Teller (BET) isotherm is used to model gas desorption effect [38], which is given below: where ν(p) is the gas volume of adsorption at pressure p, p o is the saturation pressure of the gas, ν m is the maximum adsorption gas volume when the entire adsorbent surface is being covered with a complete monomolecular layer, C is a constant related to the net heat of adsorption, n is number of the adsorption layers. Note that when n = 1, Equation (5) will reduce to the standard Langmuir isotherm [39]. For the general form of BET isotherm, the differential equilibrium partitioning coefficient of gas can be expressed as: For gas flow from fracture to wellbore, the non-Darcy flow effect is considered by using the Forchheimer [40] modification to Darcy's law below: where v is velocity, β is the non-Darcy Forchheimer coefficient, which can be determined using the correlation proposed by Evans and Civan [41] as given below: β = 1.485 × 10 9 k f 1.021 (10) where k f is the fracture permeability.

3D Semianalytical Methods
Using the method of images and considering the geometry, we can generate the solution for a unit-strength, instantaneous point source in an infinite-slab reservoir with impermeable boundaries at z = 0 and h, which is given below [42]: where η x , η y , η z is the hydraulic conductivity defined as η i = k i φµc t (i = x, y, z), and M is the interest point and M is the source point.
If we assume the hydraulic conductivity is isotropic, then Equation (11) could be rewritten as follows: Using the superposition theory, we can get the pressure solution for rectangles with no strike angles ( Figure 1a) as follows: dz dy dτ Rewriting Equation (13) we get: where x w , y w , z w is the fracture center. Calculating the integration with respect to Z gives: Similarly: x y z is the fracture center. Calculating the integration with respect to Z gives: Similarly: Thus Equation (14) could be simplified to double integral expression, as follows: Thus Equation (14) could be simplified to double integral expression, as follows:

Slanted Rectangles Fracture with Strike Angle
Sometimes the fracture plane is not orthogonal to the x-y plane, but rather has a strike angle θ as plotted in Figure 1b. In this case, the pressure solution corresponding to Equation (14) will be derived as follows: It is necessary to point out that "csc θ" term in Equation (18) is obtained when the plane integration is performed for each fracture plane. Applying the same technique as in Equations (15) and (16), we can simplify this to the following expression: The corresponding bounded domain formula will become: Equation (20) will be the primary equation to be used hereafter for rectangle fractures.

Slanted Rectangles Fracture with Dip Angle
Next we study the cases in which the fracture plane is not necessarily orthogonal to the x-y plane, and instead the angle between the fracture plane and the x-y plane is set to be θ.
In this case, the pressure drop solution corresponding to Equation (14) will be as follows: dz dy dτ which could be written as: dz dy dτ If θ = 90, then Equation (21) will be simplified to Equation (20). The corresponding bounded domain formula will then become:

Pressure Solutions for Rectangle, Ellipse, and Disk Fractures
To show the versatility of our method, and for the purpose of faster computation, the pressure solutions for three other shapes of fractures are constructed in a slab reservoir. Following the same logic as for the rectangle fracture in a bounded reservoir, the pressure solution of other shapes could also be derived. The pressure solution for a disk fracture in a slab reservoir is given as below: Similarly, the pressure solution for an ellipse fracture in a slab reservoir is given as: Three cases are designed with three different shapes of fractures. In all three cases, the depth of the slab reservoir is 30 m, and the center of all three fractures is at (0, 15). The dimensionless pressure results after 1000 days are plotted in Figure 2. The dimensionless pressure is defined as the actual minimum and maximum pressure being translated to [0, 1] range by the linear translation. Figure 2a is a fracture with rectangular shape, where the length of the rectangle equals 20 m and the height equals 10 m. Figure 2b has a disk fracture, with radius equal to 10 m. Figure 2c has an elliptical fracture, with a long radius equal to 10 m and short radius equal to 5 m. Using Equation (17) and Equations (23) and (24), the pressure solution could be calculated. To concentrate on different pressure distribution contours caused by different shapes of the fractures, only the dimensionless pressure values are plotted. The dimensionless pressure solutions are shown in Figure 2.
Comparing Figure 2a-c, there are two observations to make. First, it shows that our pressure solutions capture the pressure contours accurately, with the boundary effect clearly seen from the top and bottom of the reservoir. Second, the fracture shape impact the pressure distribution thus the shale gas production significantly, so that it is worthwhile obtaining better measurement about the fracture shapes. Comparing Figures 2a-c, there are two observations to make. First, it shows that our pressure solutions capture the pressure contours accurately, with the boundary effect clearly seen from the top and bottom of the reservoir. Second, the fracture shape impact the pressure distribution thus the shale gas production significantly, so that it is worthwhile obtaining better measurement about the fracture shapes.

Semianalytical Model for 3D Fractures in Shale Gas Reservoir
For the non-planar fracture geometry with varying fracture width and fracture permeability along fracture length, a semi-analytical model is developed to efficiently simulate the shale gas production. The semi-analytical model discretizes the non-planar fracture geometry into a number of small fracture segments in order to approximately describe the non-planar fracture geometry. Each fracture segment can have different fracture properties such as fracture length, width, and permeability to approximately represent the non-planar fracture geometry. The semi-analytical model mainly consists of two parts for simulating shale gas production: gas flow from shale matrix into fracture and gas flow from fracture to wellbore. First, we used an analytical solution to approximately solve the revised diffusivity equation for gas transport from shale matrix into fracture segments. Second, gas flow in each fracture segment is described by the non-Darcy flow equation. For the first part, an analytical solution is used to solve the gas diffusivity equation with the following assumptions: (1) The reservoir is bounded; (2) The reservoir is isotropic and homogeneous.

Semianalytical Model for 3D Fractures in Shale Gas Reservoir
For the non-planar fracture geometry with varying fracture width and fracture permeability along fracture length, a semi-analytical model is developed to efficiently simulate the shale gas production. The semi-analytical model discretizes the non-planar fracture geometry into a number of small fracture segments in order to approximately describe the non-planar fracture geometry. Each fracture segment can have different fracture properties such as fracture length, width, and permeability to approximately represent the non-planar fracture geometry. The semi-analytical model mainly consists of two parts for simulating shale gas production: gas flow from shale matrix into fracture and gas flow from fracture to wellbore. First, we used an analytical solution to approximately solve the revised diffusivity equation for gas transport from shale matrix into fracture segments. Second, gas flow in each fracture segment is described by the non-Darcy flow equation. For the first part, an analytical solution is used to solve the gas diffusivity equation with the following assumptions: (1) The reservoir is bounded; (2) The reservoir is isotropic and homogeneous; (3) It only allows gas flow from fractures to wellbore; (4) There is no gravity effect.
According to the superposition principle, one can calculate the pressure change at any location in the system by adding the contributions from all fracture segments. The pressure solution p jc at the center of the j-th fracture segment x jc , y jc , z jc for fluid flow in shale with the following expressions: where p j x jc , y jc , z jc , t is defined as Equation (20) and: In all cases of this study, the constant pressure well condition is assumed. In this case, the semianalytical model will discretize a fracture into N f fracture segments with N v nodes (N ν = N f + 1). All unknown that need to be solved include: (1) N f gas flux at fracture segments, q f i = i = 1, 2, · · · , N f .
(3) N ν pressure at all nodes p i (i = 1, 2. · · · N ν ), but the pressure at the well node is known, so that there are N f pressure unknowns at nodes.
So there are 2N f + N ν unknown variables. The corresponding equations needed to solve include: (1) N ν mass balance equations at all nodes: (2) N f non-Darcy flow equations at all fracture segments, which are described by Equation (9).
Totally there are 2N f + N ν equations, equal to the number of unknowns. Thus the model is well-defined. The Newton-Raphson iterative method was used to solve the final nonlinear system of governing equations, caused by the non-Darcy effect of gas flow. In the simulation, it is assumed that there is only gas flow under the condition of residual water saturation and the reservoir is homogeneous with a constant height, porosity, and permeability. It should be noted that the current semi-analytical model did not capture the effect of pressure-dependent matrix permeability with reservoir depletion. In order to analyze this effect, numerical reservoir simulation is recommended.

Model Validation
Before the newly derived 3D semianalytical model is used to investigate the effects of fracture height and fracture strike angles, a commercial simulator and an analytical model are used to verify the accuracy of this model. The verification case has three bi-wing hydraulic fractures and each fracture is discretized equally into 10 segments. The model dimension is 260 m × 600 m × 50 m, which corresponds to the reservoir length, width and thickness, respectively. The bottomhole pressure is held at 7.5 MPa for all simulations. Other basic reservoir, fluid and fracture parameters are listed in Table 1. The simulation was run for 10,000 days, for the newly developed 3D semianalytical model, an analytical model (Kappa Topaz) and a commercial reservoir simulator (CMG). The comparison results of the gas flow rates are plotted in Figure 3. The three curves in Figure 3 are so close to each other, with relative error less than 3%, illustrating that the newly developed 3D semianalytical model can accurately calculate the shale gas production for reservoirs with fully penetrated fractures. For the detailed description of the comparison, readers could refer to the authors' another work [29]. The simulation was run for 10,000 days, for the newly developed 3D semianalytical model, an analytical model (Kappa Topaz) and a commercial reservoir simulator (CMG). The comparison results of the gas flow rates are plotted in Figure 3. The three curves in Figure 3 are so close to each other, with relative error less than 3%, illustrating that the newly developed 3D semianalytical model can accurately calculate the shale gas production for reservoirs with fully penetrated fractures. For the detailed description of the comparison, readers could refer to the authors' another work [29].

Model Applications
After the accuracy and effectiveness of the 3D semianalytical model are verified, five cases of shale gas reservoirs with different fracture patterns are constructed, in order to study the impact of fracture height and strike angles on shale gas production. In general, the gas production is controlled by many properties of the fracture and matrix, such as permeability, fracture conductivity, but to focus on the major innovations of this study, which is the construction of a truly three dimensional hydraulic fracture model, with different strike/dip angles, we only analyze the impacts of fracture heights and strike angles. The five cases are summarized as follows:

Model Applications
After the accuracy and effectiveness of the 3D semianalytical model are verified, five cases of shale gas reservoirs with different fracture patterns are constructed, in order to study the impact of fracture height and strike angles on shale gas production. In general, the gas production is controlled by many properties of the fracture and matrix, such as permeability, fracture conductivity, but to focus on the major innovations of this study, which is the construction of a truly three dimensional hydraulic The fracture patterns of these five cases and the corresponding reservoir models are plotted in Figure 4. The green planes are the fractures, with different strike angles and fracture heights. The blue lines are the horizontal wells. The reservoir size and well properties are the same for all five cases, but the fracture heights and strike angles are different. The fracture patterns of these five cases and the corresponding reservoir models are plotted in Figure 4. The green planes are the fractures, with different strike angles and fracture heights. The blue lines are the horizontal wells. The reservoir size and well properties are the same for all five cases, but the fracture heights and strike angles are different.    Table 1. Using the semianalytical model we just developed and verified, the gas production of five cases are computed.
First, by just analyzing Case 2, the impact of fracture height could be investigated. Figure 5 shows the production data for each fracture in Case 2. To show the difference, the semi-log plot on y axis is used. The blue line is the daily production rate of the first and third fractures with 45 m height, and the red line is for the second fracture with 40 m height. It is clear that the higher fractures have more shale gas produced than the lower one. In fact, at the early stage (from day 1 to day 50), the increase of production for higher fracture is about 12.5%. But in late stage (after 4000 days), that increment could become 400%. The production data of single fractures for other cases have similar comparison results (15% and 500%). The differences between Case 2 and Case 3 are the strike angles. Case 2 has three fractures with 0° strike angles, but Case 3 has tree fractures with strike angles −10°, 0°and 10°. The same strike angle difference exists for Case 4 (Figure 4d) and Case 5 (Figure 4e).
The model dimensions (length, width and thickness) for the five cases are the same: 260 m × 600 m × 50 m. The bottomhole pressure is held at 7.5 MPa for all simulations. Other basic reservoir, fluid and fracture parameters are still the same with Table 1. Using the semianalytical model we just developed and verified, the gas production of five cases are computed.
First, by just analyzing Case 2, the impact of fracture height could be investigated. Figure 5 shows the production data for each fracture in Case 2. To show the difference, the semi-log plot on y axis is used. The blue line is the daily production rate of the first and third fractures with 45 m height, and the red line is for the second fracture with 40 m height. It is clear that the higher fractures have more shale gas produced than the lower one. In fact, at the early stage (from day 1 to day 50), the increase of production for higher fracture is about 12.5%. But in late stage (after 4000 days), that increment could become 400%. The production data of single fractures for other cases have similar comparison results (15% and 500%).   To investigate the impact of fracture height and strike angles, we compare the cumulative production data for these five cases, which are plotted in Figure 6, which shows that the cumulative production data for these five cases are decreasing as the fracture heights decrease and the strike angle variates. In fact, after 2500 days, the decreasing production caused by fracture height between Case 1 and Case 2 is around 10% to 15%. The corresponding ratio from Case 2 to Case 4 is about 10-20%. The decrease ratio caused by strike angles change is approximately 6-8% between Case 2 and Case 3, while the corresponding decrease ratio between Case 4 and Case 5 is approximately 3% to 8%. To investigate the impact of fracture height and strike angles, we compare the cumulative production data for these five cases, which are plotted in Figure 6, which shows that the cumulative production data for these five cases are decreasing as the fracture heights decrease and the strike angle variates. In fact, after 2500 days, the decreasing production caused by fracture height between Case 1 and Case 2 is around 10% to 15%. The corresponding ratio from Case 2 to Case 4 is about 10-20%. The decrease ratio caused by strike angles change is approximately 6-8% between Case 2 and Case 3, while the corresponding decrease ratio between Case 4 and Case 5 is approximately 3% to 8%. Figure 6. Comparisons of shale gas cumulative production for five cases. Figure 6 shows clearly that the impacts of fracture areas for the shale gas productions. With decreasing fracture heights from Case 1 to Case 2 and Case 4, while their lengths are the same, the fracture areas are correspondingly decreasing. It is seen from Figure 6 that the cumulative productions of Case 1, Case 2 and Case 4 are decreasing too, with the similar ratio as the fracture heights. For example, the total fracture heights of Case 1 and Case 2 are 150 m and 130 m, respectively, and the ratio of the total heights are 87%, which is very close to the cumulative production ratio of Case 1 and Case 2. The same comparison result holds true for Case 4.
The well interference effect has also been widely studied recently [43][44][45]. In Figure 6, comparing Case 2 and Case 3, the only difference is that Case 3 has nonzero strike angles. The difference between their cumulative productions is very clear. The cumulative production of Case 3 is about 7% less than Case 4. This quantitatively verifies the impact of strike angles on shale gas production, which provides measurement criteria for fracturing design and evaluations.

Discussion
Shale gas production has become an important part of the energy supply of the world, especially for the United States. It will become even more dominant in the next 20 years, but modelling shale gas reservoirs production has proved challenging due to the complex fluid flow mechanisms and the existence of various fracture patterns. Although many techniques have been proposed, such as numerical methods (Finite Volume Method, Finite Element Method) based on structured or unstructured grids, analytical methods, and decline curve analysis, there is still a need to design a method that honors the complex fracture geometry while having faster computation speed. In this study, a new truly 3D semianalytical model for bounded shale gas reservoirs with various fracture patterns is designed and validated with a commercial simulator and an analytical model. Then this model is applied to investigate the impacts of fracture heights and strike angles on daily production rate and cumulative production.  Figure 6 shows clearly that the impacts of fracture areas for the shale gas productions. With decreasing fracture heights from Case 1 to Case 2 and Case 4, while their lengths are the same, the fracture areas are correspondingly decreasing. It is seen from Figure 6 that the cumulative productions of Case 1, Case 2 and Case 4 are decreasing too, with the similar ratio as the fracture heights. For example, the total fracture heights of Case 1 and Case 2 are 150 m and 130 m, respectively, and the ratio of the total heights are 87%, which is very close to the cumulative production ratio of Case 1 and Case 2. The same comparison result holds true for Case 4.
The well interference effect has also been widely studied recently [43][44][45]. In Figure 6, comparing Case 2 and Case 3, the only difference is that Case 3 has nonzero strike angles. The difference between their cumulative productions is very clear. The cumulative production of Case 3 is about 7% less than Case 4. This quantitatively verifies the impact of strike angles on shale gas production, which provides measurement criteria for fracturing design and evaluations.

Discussion
Shale gas production has become an important part of the energy supply of the world, especially for the United States. It will become even more dominant in the next 20 years, but modelling shale gas reservoirs production has proved challenging due to the complex fluid flow mechanisms and the existence of various fracture patterns. Although many techniques have been proposed, such as numerical methods (Finite Volume Method, Finite Element Method) based on structured or unstructured grids, analytical methods, and decline curve analysis, there is still a need to design a method that honors the complex fracture geometry while having faster computation speed. In this study, a new truly 3D semianalytical model for bounded shale gas reservoirs with various fracture patterns is designed and validated with a commercial simulator and an analytical model. Then this model is applied to investigate the impacts of fracture heights and strike angles on daily production rate and cumulative production.
To the best of the authors' knowledge, this study is the first to rigorously derive the pressure solution for three dimensional fracture plane in a bounded shale gas reservoir. The fractures could be fully penetrated or partially penetrated. They can also have any strike angle or dip angle. The shapes are not constrained either. The pressure solutions for fractures with shapes of rectangles, disks, and ellipses have been derived explicitly. The same derivation technique could be adopted to derive pressure solutions for any shape of fractures with explicit expressions of the planes.
For complicate shapes of fractures such as disks and ellipses, the pressure solutions have the form of triple integrals. The computation of these triple integrals will be time-consuming because of the singularity existing at the end points of integral interval. In order to efficiently model shale gas reservoirs with complex fracture shapes and patterns, a better numerical integration algorithm needs to be developed and implemented into the current semianalytical model.
This study concentrates on single phase, shale gas reservoirs, but the technique proposed here could be easily expanded to simulate multiphase reservoirs. For multiphase reservoirs with fully penetrated, rectangle shape fractures, several researchers have made great effort to develop semianalytical models to model the pressure profile and production [46]. The pressure solution derived in this study could be directly implemented into multiphase semianalytical model to model more complex fluid physics.
Finally, the fractures in this study are all synthetic with uniform fracture widths, permeability, and fracture conductivity, in the purpose of concentrating on investigating the impact of fracture heights and angles systematically and quantitatively. In the future, more practical fracture patterns will be implemented into the semianalytical model by integrating seismic data or well logging data. The fracture patterns could be generated based on microseismic data and fracture outcrop figures [47,48], which will make the semianalytical model more realistic. The fracture widths, fracture conductivity could also be variables as what the authors have done before for fully penetrated nonplanar fractures [26]. The classical gas transport equations used in this study could also be replaced with fractional gas diffusion equations to include the historical dependence effect of the shales [49][50][51][52][53].

Conclusions
With the rapidly growing interest in unconventional resources such as shale gas, studying the impact of hydraulic fractures on shale gas production has become an urgent task for petroleum engineers. In this study, based on the gas transport equation considering multiple fluid mechanisms such as gas desorption, gas diffusion and gas slippage, a semianalytical model was constructed and solved for shale gas reservoirs with fully 3D hydraulic fractures. This method could model hydraulic fractures with different strike/dip angles and different shapes. This semianalytical model was verified with a commercial software and an analytical method. The shale gas production of reservoirs with different fracture heights and strike angles are simulated and compared. Based on our study, the following conclusions could be drawn:

•
The pressure for partially penetrated fractures with different strike angles, dip angles and shapes could be derived by integrating the instantaneous point source solution on the fracture plane; • Based on comprehensive shale gas transport equations considering complex gas transport mechanisms and complex fracture patterns, a novel semianalytical model to simulate the gas transport and production is developed by discretizing each fracture to multiple segments; • This semianalytical mode has been verified with a commercial simulator and an analytical model for a shale gas reservoir with three fully penetrated fractures, and the relative error is within 3%; • Five shale gas reservoir cases with different fracture heights and strike angles have been investigated by using this semianalytical model. It is found that the decreasing of fracture height will lead to decreasing production rate and cumulative shale gas production, by a ratio close to the fracture areas ratio; and strike angles could cause well interference and decrease the daily production rate and cumulative shale gas production.

•
The technique proposed in this study and the partial results shed a light on the importance of fracture shapes and heights. In most of the shale gas models, the hydraulic fractures are always assumed to be rectangle and fully penetrated in the pay zone. But from the results of this study, the fracture shapes and heights play essential roles in influencing the shale gas production. On the other hand, these explicit pressure results derived in this study for fracture planes with different shapes and angles could also be used in well testing to investigate the fracture and reservoir properties. Thirdly, although the cases studied here are all synthetic and simple, the methods could still be used for complex fracture patterns with the help of more accurate geophysical measurement. In conclusion, the method proposed in this work will be widely used by petroleum engineers dealing with unconventional reservoirs with complex fractures. Total length of the reservoir, x-direction L y Total length of the reservoir, y-direction L z Total length of the reservoir, z-direction