Trilinear fuzzy seepage model of shale gas reservoirs with multiple fuzzy parameters

Abstract

This paper focuses on establishing the trilinear fuzzy seepage model with multiple fuzzy parameters for shale gas reservoirs. Different from the conventional seepage models of shale gas reservoirs, the multiple fuzzy parameters seepage model uses fuzzy numbers to describe some parameters with uncertainty. Firstly, the multiple fuzzy parameters seepage model is constructed based on fuzzy concepts. The fuzzy structure element method and the centroid method are used to solve the fuzzy seepage model and defuzzifier, respectively. Secondly, the advantages of the development fuzzy model over the conventional seepage model are discussed and illustrated through numerical examples and simulations. Finally, to further study the seepage laws inside shale gas reservoirs, this paper explores the sensitivity of relevant main control parameters to gas production based on the development model.

Keywords

Shale gas reservoirs fuzzy parameter fuzzy differential equation fuzzy structural element fuzzy modeling

1 Introduction

Since modern times, the industry has developed rapidly, and the demand for energy is increasing. As a result, conventional energy sources such as petroleum and coalbed methane are gradually depleted. The shale gas reservoir is expected to become a new energy pillar because of its large storage capacity and long exploitation period [1 –6]. Most shale gas producing areas have a wide distribution range, large thickness, and generally contain gas, which enables shale gas wells to create gas at a stable rate for an extended period. Compared with conventional natural gas, shale gas has poor physical properties and has extremely low permeability. These physical properties of shale gas reservoirs directly lead to high exploitation costs and development difficulties. It is necessary to have a deep understanding of the seepage conditions inside the reservoirs to realize the economic exploitation of shale gas reservoirs. The innovation of horizontal well fracturing technology has played a significant role in the rapid growth of the shale gas industry. Research on the hydraulic fracturing and seepage model of shale gas reservoirs has made significant progress.

The production stage of shale gas reservoirs is dominated by linear flow, so many linear flow models have emerged to study the seepage law of shale gas reservoirs. Sheng-Tai Lee [7] first developed a new analytical solution based on an approximate trilinear flow model to study the transient behavior of wells intercepted by vertical fractures with limited conductivity. This solution considers the effects of skin, wellbore storage, and fracture storage, as well as constant pressure and constant speed conditions. Erdal Ozkan [8] studied the performance of horizontally fractured wells in millidarcy permeability (conventional) and micro-nano-permeability (unconventional) reservoirs. The trilinear flow model proposed in this work and the information obtained from it should be helpful for the design and performance prediction of multiple fractured horizontal wells in shale gas reservoirs. Yu-long Zhao [9] proposed a “three-hole” mathematical model to describe fluid flow from shale gas formation to multi-fractured horizontal wells. Based on the new model, the effects of horizontal well length, fracture number, and Langmuir volume on the type curve are analyzed.

However, these conventional seepage models of shale gas reservoirs generally have two problems.

The traditional seepage models took into account the complexity of the model in the modelling process and idealized uncertain parameters to a specific value.

The parameters in the process of nodeling were determined by the laboratory or by expert experience.

Therefore, to solve these problems is the main motivation for the current study. Fuzzy set theory is a powerful tool for solving uncertainty problems [10 –20]. Since 1965, the fuzzy set theory has been in development for more than 50 years, showing substantial advantages in many fields, such as automatic control, computer, and information processing [21 –23]. Based on the IVHFHSWA(interval-valued hesitant fuzzy Hamacher synergetic weighted average) operator, Guoliang Li proposed a new method for shale gas region selection with interval value hesitant fuzzy information. According to this, more choices can be provided to decision-makers. The fuzzy seepage theory [24 –28] was successfully introduced into the modeling analysis of the stope gas [29]. In 2019, Zhang Duo, Qiu Dong, et al. first proposed the concept of fuzzy permeability and introduced the fuzzy set theory into the study of the fuzzy seepage model [30]. Based on this, the fuzzy seepage model was established, which took into account the uncertainties in the mining process. However, these models only consider one fuzzy parameter, and there are many other parameters with uncertainty in the modeling process.

This paper mainly studies the trilinear fuzzy seepage models of shale gas reservoirs with multiple fuzzy parameters.With the above in mind, we summarize the novelties of this work as below:

Introduce the fuzzy concepts(fuzzy formation initial pressure, fuzzy reservoir thickness, fuzzy permeability, fuzzy fractures number, and fuzzy main fracture width) into the seepage model of shale gas reservoirs and treat the difficult-to-measure parameters as fuzzy numbers;

Based on the proposed fuzzy concepts, the fuzzy seepage model for shale gas reservoirs is constructed. The fuzzy structure element method and the centroid method are used to solve fuzzy differential equations and defuzzification, respectively. A numerical experiment is given to show the effectiveness and superiority of the improved fuzzy seepage model;

Sensitivity analysis is performed for specific parameters that have a significant impact on shale gas reservoir production. In this way, the study of seepage theory in shale gas reservoirs is promoted.

The remainder of this note is organized as follows. In Section 2, we recall some relevant concepts. In Section 3, the fuzzy concepts are considered, fuzzy seepage models of shale gas reservoir are built, and the fuzzy structural element method is used to solve those models. For the practical application, the centroid method of defuzzification should be further used to obtain a single numerical solution corresponding to the fuzzy solution set. In Section 4, Sensitivity analyses are performed for certain parameters that have a significant impact on shale gas reservoir production. In Section 5, some concluding remarks are drawn.

2 Basic concepts

2.1 Fuzzy structure element

Definition 2.1. [30] Let E be the fuzzy set on the real number field R, the membership function is E (x), x ∈ R, if E (x) satisfies the following properties:

E (0) =1, E (1 +0) = E (-1 - 0) =0;

E (x) is a single-increasing right continuous function on the interval [-1, 0), and E (x) is a single-drop left continuous function in the interval (0, 1];

When - ∞ < x < -1 or 1< x < + ∞, E (x) =0.

Then the fuzzy set E is the fuzzy structural element on R.

Definition 2.2. [30] If the membership function of the fuzzy structure element is E (x) =0 on the interval (-1, 1), and it is a continuous strictly single increment function within (-1, 0), it is a continuous strictly one - down function within the interval (0, 1), then E is a regular fuzzy structure element.

Definition 2.3. [30] Let E have membership function $E (x) = {\begin{matrix} 1 + x & x \in [- 1, 0], \\ 1 - x & x \in (0, 1], \\ 0 & other . \end{matrix}$ (1)E is called triangular structural element. Trigonometric structuring element is a regular fuzzy structural element.

Theorem 2.4. [30] Let E be a fuzzy structural element on R and have a membership function E (x). The set of function f (x) is monotonically bounded on the interval [-1, 1], then the fuzzy transformation f (E) of fuzzy structural element E is a finite fuzzy number on R and the membership function of f (E) is E (f^-1 (x)). Here f^-1 (x) is f (x) rotation symmetric function of x, y. If f (x) is strictly monotonous on the interval [-1, 1], then f^-1 (x) is the inverse of f (x).

2.2 Derivative of fuzzy value function based on structural element

Fuzzy numbers and fuzzy value functions are the most basic concepts in the fuzzy analysis.

Definition 2.5. [30] If for ∀λ ∈ (0, 1], the λ-level sets f_λ (x) of the fuzzy value function $\tilde{f} (x)$ is derivable on D. Then the fuzzy value function $\tilde{f} (x)$ is derivable on D, and ${\tilde{f}}^{'} (x) = ⋃_{λ \in [0, 1]} λ * f_{λ}^{'} (x),$ (2) is the derivative of $\tilde{f} (x)$ on D.

Definition 2.6. [30] Let $\tilde{f} (x) = g (x, E)$ be the fuzzy value function generated by the fuzzy structure element E, where g (x, y) is a monotone bounded function with respect to y on [-1, 1], and constitute an ordered cone or a locally ordered cone on X. If the function g (x, y) is derivable with respect to x on D ⊆ X, then $\tilde{f} (x)$ is said to be derivable on D, and its derivative is defined as $\frac{d \tilde{f} (x)}{dx} = \frac{\partial g (x, y)}{\partial x} |_{y = E} .$ (3)

The following theorem proves that Definition 2.5 and Definition 2.6 are equivalent.

Theorem 2.7. [30] The necessary and sufficient condition for the fuzzy value function $\tilde{f} (x) = g (x, E)$ to be derivable under the condition of Definition 2.6 is that ∀λ ∈ (0, 1], the λ-level sets f_λ (x) of $\tilde{f} (x)$ is derivable on D, and the derivative function $\tilde{f^{'}} (x)$ is still a fuzzy value function on D.

2.3 Definition of fuzzy finite differential equation

With the development of society, the research system has become more and more complex. Due to various reasons, cannot be fine to measure some parameters of the system, and can only use fuzzy numbers or fuzzy value function estimate. Therefore, the solution of the differential equation produces a certain fuzziness. According to the background of the actual problem, put forward the concept of the fuzzy differential equation.

Definition 2.8. [30] The differential equation is the conditional equation between the unknown function and its derivative and the known function. By solving the equation, the unknowns can be determined. Differential equations including ordinary differential equations and partial differential equations, unified to remember as follows $F (t; x; D^{m} x; G (t)) = 0,$ (4) where t = t (t₁, t₂, . . . , t_n) is an independent variable; x = x (t) is an unknown function; D^mx is a set of unknown functions until the derivative of the m order, and $Dx = (\frac{\partial x}{\partial t_{1}}, \frac{\partial x}{\partial t_{2}}, . . ., \frac{\partial x}{\partial t_{n}})$ , $D^{k} x = \frac{\partial^{k} x}{\partial {t_{1}}^{k_{1}} \partial {t_{2}}^{k_{2}} . . . \partial {t_{n}}^{k_{n}}}$ , k₁, k₂, . . . k_n are nonnegative integers, and meet $\sum_{i = 1}^{n} k_{i} = k$ , k = 1, 2, . . . , m. G (t) is a set of known functions and constants.

Definition 2.9. [30] The fuzzy differential equation refers to the conditional equation between the unknown fuzzy value function and its fuzzy value derivative function and the known fuzzy value function (or the known fuzzy number). The differential equation of fuzzy function is recorded as $F (t; \tilde{x}; D^{m} \tilde{x}; \tilde{G_{1}} (t)) = \tilde{G_{2}} (t) .$ (5) In which $\tilde{G_{1}} (t)$ , $\tilde{G_{2}} (t)$ represent fuzzy parameter terms, $D^{m} \tilde{x}$ is the same order (reverse order) fuzzy derivative function of the unknown fuzzy value function $\tilde{x}$ .

Definition 2.10. [30] Let $\tilde{g} (t)$ for the given fuzzy valued function, and λ ∈ (0, 1], and the λ cut set of $\tilde{g} (t)$ is the interval value function $g_{λ} (t) = [g_{λ}^{-} (t), g_{λ}^{+} (t)]$ . The function $g_{λ}^{-} (t)$ and $g_{λ}^{+} (t)$ are the λ qualified functions of $\tilde{g} (t)$ . The total of all λ qualified functions of the fuzzy value function $\tilde{g} (t)$ is C_λ.

Definition 2.11. [30] Let $\tilde{G} (t)$ be a known function of fuzzy value, and g (t) is a known common function, with $μ_{\tilde{G} (t)} (g (t)) = λ$ . For a given differential equation F (t ; x ; D^mx ; g (t)) = 0 with a parameter of g (t), the existence of the solution of the differential equation under a given initial condition D₀. The solution is x (t, g (t)), and $F (t; x; {\tilde{D}}^{m} x; \tilde{G} (t)) = 0,$ (6) is a fuzzy finite differential equation. Its solution is $\tilde{x} (t) = ⋃_{g (t) \in C_{λ}} λ * x (t, g (t)) .$ (7)

Definition 2.12. [30] If the fuzzy parameter term in the fuzzy finite differential equations is replaced by the explicit parameter, the ordinary differential equation is called the characterization equation of the fuzzy finite differential equations, and the explicit parameter of replacement is the characterization parameter.

Theorem 2.13. Let $F (t; x; {\tilde{D}}^{m} x; \tilde{G} (t)) = 0,$ (8) be a single parameter fuzzy finite differential equation, and the corresponding portrayal equation is $F (t; x; D^{m} x; g (t)) = 0 .$ (9) The existence of solutions for (9) under a given initial condition D₀, and the solution is f (t) = f_s (g ; D₀), the solution of the equation (8) is $\tilde{x} (t; D_{0}) = f_{s} (g; D_{0}) |_{g = \tilde{G} (t)} = f_{s} (\tilde{G} (t); D_{0}) .$ (10)

Definition 2.14. [30] Set the fuzzy differential differential equation $F (t; x; {\tilde{D}}^{m} x; \tilde{G} (t)) = 0$ , in which $\tilde{G} (t) (t \in Ω)$ is a bounded closed fuzzy value function. If the solution $\tilde{x} (t; D_{0})$ exists under the given initial condition D₀ and is a fuzzy value function on the Ω, the definite solution $\tilde{x} (t; D_{0})$ of the fuzzy finite differential equation can be expressed.

Theorem 2.15. [30] For a fuzzy finite differential equation with single fuzzy parameter, If the portrayed parameters and solutions of the corresponding portrayal equations are monotonous, the solution of the fuzzy finite differential equation can be expressed.

2.4 Multiple fuzzy parameters differential equation

Without loss of generality, this section uses double fuzzy parameters as an example to illustrate.

Let $G (t, x, D^{m} x, \tilde{L} (t)) = 0$ be a fuzzy differential equation with two fuzzy parameter terms ${\tilde{Q}}_{1}$ and ${\tilde{Q}}_{2}$ . For ∀λ ∈ (0, 1], the λ-level sets of ${\tilde{Q}}_{1}$ and ${\tilde{Q}}_{2}$ are interval value functions, which are recorded as ${\tilde{Q}}_{λ}^{(1)} (t) = [γ_{1 (λ)} (t), γ_{2 (λ)} (t)]$ , ${\tilde{Q}}_{λ}^{(2)} (t) = [δ_{1 (λ)} (t), δ_{2 (λ)} (t)]$ , respectively. Let us denote that $D_{Q λ} (t) = {\tilde{Q}}_{λ}^{(1)} (t) \times {\tilde{Q}}_{λ}^{(2)} (t)$ . For a given (γ (t) , δ (t)) ∈ D_Qλ (t), the interval function terms ${\tilde{Q}}_{λ}^{(1)} (t)$ and ${\tilde{Q}}_{λ}^{(2)} (t)$ in $G (t, x, D^{m} x, \tilde{L} (t)) = 0$ are replaced by γ (t) and δ (t), respectively, to obtain a general differential equation containing only ordinary function terms. The differential equation satisfies the existence of the solution of the condition D₀, and the solution is g (t) = Fs (γ (t) , δ (t) ; D₀). Then the solution of the interval differential equation with two interval function terms ${\tilde{Q}}_{λ}^{(1)} (t)$ and ${\tilde{Q}}_{λ}^{(2)} (t)$ is

$\begin{matrix} g (t; {\tilde{Q}}_{λ}^{1} (t), {\tilde{Q}}_{λ}^{2} (t); D_{0}) \\ = \cup_{(γ (t), δ (t)) \in D_{Q λ} (t)} Fs (γ (t), δ (t); D_{0}) . \end{matrix}$ (11)

Theorem 2.16. [30] The characterization parameters of the characterization equations of the interval differential equations with two interval function terms and the solutions of the equations are monotonic. Assume that both parameters and solutions are in the same order, and ∀ (γ (t) , δ (t)) ∈ D_Pλ (t), we have $\begin{matrix} Fs (γ_{1 (λ)} (t), δ_{1 (λ)} (t); D_{0}) \leq Fs (γ (t), δ (t); D_{0}) \\ \leq Fs (γ_{2 (λ)} (t), δ_{2 (λ)} (t); D_{0}) . \end{matrix}$ Then the solution of the interval equation (11) can be formally written as an interval value function $\begin{matrix} g (t; {\tilde{Q}}_{λ}^{(1)} (t), {\tilde{Q}}_{λ}^{(2)} (t); D_{0}) = [Fs (γ_{1 (λ)} (t), δ_{1 (λ)} (t); \\ D_{0}), Fs (γ_{2 (λ)} (t), δ_{2 (λ)} (t); D_{0})] . \end{matrix}$ From the performance theorem, the fuzzy solution of fuzzy differential equations with two fuzzy parameter terms is $\begin{matrix} \tilde{g} (t) = & \cup_{λ \in [0, 1]} λ * {\tilde{g}}_{λ} (t) \\ = & \cup_{λ \in [0, 1]} λ * g (t; {\tilde{P}}_{λ}^{(1)} (t), {\tilde{P}}_{λ}^{(2)} (t); D_{0}) \\ = & [Fs (γ_{1 (λ)} (t), δ_{1 (λ)} (t); D_{0})] \\ = & \cup_{λ \in [0, 1]} λ * [Fs (γ_{1 (λ)} (t), δ_{1 (λ)} (t); D_{0}), \\ Fs (γ_{2 (λ)} (t), δ_{2 (λ)} (t); D_{0})] . \end{matrix}$

3 Shale gas reservoir fuzzy seepage models

In order to efficiently and reasonably develop shale gas reservoirs, the trilinear seepage model of fracturing horizontal wells in shale gas reservoirs was proposed [36]. The fracturing of horizontal wells in shale gas reservoirs will form a fracture network with a hydraulic fracture. Shale gas reaches the fracture network through bedrock desorption, diffusion and seepage, and then flows into the well through hydraulic fractures, which forms a complete seepage process. However, there are obvious differences in the seepage of bedrock, natural fractures and hydraulic fractures. The bedrock is the storage space for shale gas, and its multi-level seepage mode is prominent; natural fractures increase the seepage channels, forming a bridge connecting the bedrock and hydraulic fracturing fractures. Hydraulic fracturing fractures are the main seepage channel and the contact bridge between the well and the reservoir, with high permeability. The seepage flow in the above three kinds of media is very complicated, so it is simplified according to the relationship. It is divided into three partial seepage modes as shown in the Fig. 1. Since the seepage velocity of natural gas in the three parts is relatively slow, it is assumed that all three parts follow linear seepage, which is defined as trilinear seepage. The relevant model assumptions are as follows:

The research model is proposed for the single-phase flow of a constant compressible fluid. The flow of liquid from the reservoir to a horizontal well depends solely on hydraulic fracturing. The output produced directly from the horizontal well surface is negligible;

It is assumed that hydraulic fractures have the same characteristics and are equally spaced along the horizontal well, d_F. Both the flow in the internal reservoir between hydraulic fracturing and the flow in the external reservoir outside the hydraulic fracturing tip are assumed to be linear.

Fig.1

Diagram of trilinear seepage flow model for multi-stage fractured horizontal wells.

Under the conditions assumed in paper [36], a trilinear seepage model under the polar coordinate system is derived below (parameters were defined as Table 1 shows).

Table 1

Parameters for numerical simulation

Parameter	Unit	Value	Parameter	Unit	Value
Hydraulic fracture half length, x_F	m	75	Fracture porosity,φ_f	%	0.45
Hydraulic fracture porosity,φ_F	%	0.6	Matrix compressibility,C_tm	1/MPa	0.00022
Fracture compressibility,C_tf	1/MPa	0.00035	Well depth,H	m	1600
Matrix permeability,K_m	md	0.0000001	Fracture permeability,K_f	md	0.0001
Hydraulic fracture permeability,K_F	md	300	Formation thickness, h	m	46
Horizontal section length, x_F	m	1006	Compression factor,Z	number 0.89
Reservoir size in x-direction, x_e	m	300	Matrix porosity,φ_m	%	4.5
Hydraulic fracture width, w_F	mm	2	Langmuir volume, V_L	m³/m³	2
Number of hydraulic fractures, n_F	-	33	Langmuir pressure, P_L	MPa	10
Hydraulic fracture compressibility,C_tF	1/MPa	0.00045	Initial reservoir pressure,P_i	MPa	21
Initial reservoir temperatuer,T	°C	65	Molar gas constant, R	Jmol ^-1 K ^-1	8.314
Reservoir size in y-direction, y_e	m	15.625	Gas viscosity, μ	mPas	0.027

Liner flow in matrix

The seepage equation of the matrix system in polar coordinates is:

$\begin{matrix} \frac{\partial}{\partial x} (\frac{P_{m} M k_{m}}{RTZ μ} \frac{\partial P_{m}}{\partial r}) = \frac{\partial}{\partial t} (\frac{φ_{m} P_{m} M}{RTZ} \\ + ρ_{gsc} \frac{V_{L} P_{m}}{P_{L} + P_{m}}) \end{matrix}$ (12) Define the following pseudo pressure functions and pressure transmiting coefficients: $\begin{matrix} m_{m} = 2 \int_{0}^{P_{m}} \frac{P}{μ Z} dP, \\ η_{m} = 2.637 \times 10^{- 4} \frac{k_{m}}{φ_{m} C_{tm} μ} \end{matrix}$ where $C_{tm} = \frac{1}{P_{m}} - \frac{1}{Z} \frac{\partial Z}{\partial P_{m}} + \frac{ρ_{gsc} V_{L} P_{L} RTZ}{φ_{m} M P_{m} (P_{L} + P_{m})^{2}}$ . So the seepage equation in the form of pseudo pressure can be obtained as follows: $\frac{\partial^{2} m_{m} (x, t)}{\partial x^{2}} = \frac{1}{η_{m}} \frac{\partial m_{m}}{\partial t}$ In order to facilitate the calculation, the following dimensionless parameters are further defined: $\begin{matrix} m_{mD} = \frac{k_{m} y_{e}}{1422 q_{F} T}, \\ t_{D} = 1.69 \times 10^{- 7} \frac{k_{f} t}{φ_{f} μ C_{tf} x_{e}^{2}} \end{matrix}$ $\begin{matrix} x_{D} = \frac{x}{x_{F}}, y_{D} = \frac{y}{x_{F}}, η_{mD} = \frac{η_{m}}{η_{f}}, \\ η_{f} = 1.82 \times 10^{- 6} \frac{k_{f}}{φ_{f} μ C_{tf}} \end{matrix}$ Then the dimensionless form of seepage equation is: $\frac{\partial^{2} m_{mD} (x_{D}, t_{D})}{\partial x_{D}^{2}} = \frac{1}{η_{mD}} \frac{\partial m_{mD}}{\partial t_{D}}$ By laplace transformation, the seepage equation can be converted into the following form: $\frac{\partial^{2} {\bar{m}}_{mD} (x_{D}, t_{D})}{\partial x_{D}^{2}} - \frac{s}{η_{mD}} {\bar{m}}_{mD} = 0$

Liner flow in fracture network.

The seepage equation in a fracture network is: $\frac{\partial}{\partial y} (\frac{P_{f} M k_{f}}{RTZ μ} \frac{\partial P_{f}}{\partial r}) + 2 q_{mf} = \frac{\partial}{\partial t} (\frac{φ_{f} P_{f} M}{RTZ})$ (13) where $q_{mf} = \frac{P_{m} M k_{m}}{RTZ μ} \frac{\partial P_{m}}{\partial x} |_{x = x_{F}}$ Similar to matrix systems, define pressure and dimensionless parameters: $\begin{matrix} m_{f} = 2 \int_{0}^{P_{f}} \frac{P}{μ Z} dP, m_{fD} = \frac{k_{f} x_{F}}{1422 q_{F} T} \\ [m_{f} (P_{i}) - m_{f} (P)], η_{fD} = \frac{η_{f}}{η_{f}} = 1 \end{matrix}$ By dimensionless analysis and Laplace transform, we can get $\frac{\partial^{2} {\bar{m}}_{fD} (y_{D}, s)}{\partial y_{D}^{2}} + 2 \frac{k_{m}}{k_{f}} \frac{\partial {\bar{m}}_{mD}}{\partial x_{D}} |_{x_{D} = 1} - s {\bar{m}}_{fD} = 0$ (14)

Liner flow in hydraulic fractures

By dimensionless analysis and Laplace transform, we can get

$\begin{matrix} \frac{\partial^{2} {\bar{m}}_{FD} (y_{D}, s)}{\partial x_{D}^{2}} + \frac{4}{C_{FD}} \frac{\partial {\bar{m}}_{fD}}{\partial y_{D}} |_{y_{D} = W_{D} / 2} \\ - \frac{s}{η_{FD}} {\bar{m}}_{FD} = 0 \end{matrix}$ (15) where $m_{F} = 2 \int_{0}^{P_{F}} \frac{P}{μ Z} dP, η_{FD} = \frac{η_{F}}{η_{f}}$ , $η_{F} = 1.82 \times 10^{- 6} \frac{k_{F}}{φ_{F} μ C_{tF}}$ , $C_{FD} = \frac{k_{F} W}{k_{f} x_{F}}$ .

So the trilinear seepage problem of shale gas is as follows:

{\begin{matrix} \frac{\partial^{2} {\bar{m}}_{mD} (x_{D}, s)}{\partial x_{D}^{2}} - \frac{s}{η_{mD}} {\bar{m}}_{mD} (x_{D}, s) = 0, \\ \frac{\partial^{2} {\bar{m}}_{fD} (y_{D}, s)}{\partial y_{D}^{2}} + 2 \frac{K_{m}}{K_{f}} \frac{\partial {\bar{m}}_{mD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 1} = s {\bar{m}}_{fD} (y_{D}, s), \\ \frac{\partial^{2} {\bar{m}}_{FD} (x_{D}, s)}{\partial x_{D}^{2}} + \frac{4}{C_{FD}} \frac{\partial {\bar{m}}_{fD} (y_{D}, s)}{\partial y_{D}} |_{y_{D} = \frac{w_{D}}{2}} = \frac{s}{η_{FD}} {\bar{m}}_{FD} (x_{D}, s), \\ \frac{\partial {\bar{m}}_{mD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = x_{eD}} = \frac{\partial {\bar{m}}_{fD} (y_{D}, s)}{\partial y_{D}} |_{y_{D} = y_{eD}} = 0 \\ \frac{\partial {\bar{m}}_{FD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 1} = 0, q = q_{F} n_{F} \\ {\bar{m}}_{mD} (x_{D}, s) |_{x_{D} = 1} = {\bar{m}}_{fD} (x_{D}, s) |_{x_{D} = 1} \\ {\bar{m}}_{fD} (y_{D}, s) |_{y_{D} = w_{D} / 2} = {\bar{m}}_{FD} (y_{D}, s) |_{y_{D} = w_{D} / 2}, \\ {\bar{m}}_{FD} (x_{D}, s) |_{x_{D} = 0} = {\bar{m}}_{wD}, \frac{\partial {\bar{m}}_{FD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 0} = - \frac{π}{C_{FD} s} . \end{matrix}

(16) So

$\begin{matrix} q = A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}) . \end{matrix}$ (17) where $A = \frac{6.89 \times 10^{6} x_{e} K_{m}^{2} t_{D}^{- \frac{1}{4}}}{π y_{eD} η_{m} T μ {zK}_{f} x_{F}}$ , B = 2η_f (x_eD - 1) + C_fDy_eDη_m, $C = \frac{y_{eD} η_{m}}{K_{f}^{2} x_{F}^{2}}$ .

3.1 Single fuzzy parameter seepage models

It is well known that the distribution of shale gas reservoirs is exceptionally uneven due to environmental factors such as temperature and topography, which leads to some parameters of the reservoir with uncertainty. Therefore, this section considers fuzzy numbers to describe some parameters with an uncertainty of shale gas reservoirs.

3.1.1 Fuzzy initial formation pressure model

The initial formation pressure, that is, the reservoir pressure measured when the reservoir is in a pressure equilibrium state before the reservoir is mined. Figure 2 shows the trend of gas production for different initial formation pressure, P_i. It can be seen from the figure that the shale gas reservoir is significantly affected by the initial pressure of the shale gas reservoir.

Fig.2

The trend of gas production for different formation initial pressure P_i.

Reservoir pressure generally increases with increasing burial depth and temperature. The study of reservoir pressure is of great significance. It has not only an essential impact on the gas content of the reservoir and the state of gas storage but also the energy of gas flowing from the fracture to the wellbore. Reservoir pressure is generally measured by well testing, and the commonly used method is the water injection/pressure drop method. Based on the fact that the burial depth and temperature at different locations in the shale reservoir are different, this paper plans to use fuzzy number ${\tilde{P}}_{i}$ to describe the initial formation permeability.

${\begin{matrix} \frac{\partial^{2} {\tilde{\bar{m}}}_{mD} (x_{D}, s)}{\partial x_{D}^{2}} - \frac{s}{η_{mD}} {\tilde{\bar{m}}}_{mD} (x_{D}, s) = 0, \\ \frac{\partial^{2} {\tilde{\bar{m}}}_{fD} (y_{D}, s)}{\partial y_{D}^{2}} + 2 \frac{K_{m}}{K_{f}} \frac{\partial {\tilde{\bar{m}}}_{mD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 1} = s {\tilde{\bar{m}}}_{fD} (y_{D}, s), \\ \frac{\partial^{2} {\bar{m}}_{FD} (x_{D}, s)}{\partial x_{D}^{2}} + \frac{4}{C_{FD}} \frac{\partial {\tilde{\bar{m}}}_{fD} (y_{D}, s)}{\partial y_{D}} |_{y_{D} = \frac{w_{D}}{2}} = \frac{s}{η_{FD}} {\tilde{\bar{m}}}_{FD} (x_{D}, s), \\ \frac{\partial {\tilde{\bar{m}}}_{mD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = x_{eD}} = \frac{\partial {\tilde{\bar{m}}}_{fD} (y_{D}, s)}{\partial y_{D}} |_{y_{D} = y_{eD}} = 0 \\ \frac{\partial {\tilde{\bar{m}}}_{FD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 1} = 0, q = q_{F} n_{F} \\ {\tilde{\bar{m}}}_{mD} (x_{D}, s) |_{x_{D} = 1} = {\tilde{\bar{m}}}_{fD} (x_{D}, s) |_{x_{D} = 1} \\ {\tilde{\bar{m}}}_{fD} (y_{D}, s) |_{y_{D} = w_{D} / 2} = {\tilde{\bar{m}}}_{FD} (y_{D}, s) |_{y_{D} = w_{D} / 2}, \\ {\tilde{\bar{m}}}_{FD} (x_{D}, s) |_{x_{D} = 0} = {\tilde{\bar{m}}}_{wD}, \frac{\partial {\tilde{\bar{m}}}_{FD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 0} = - \frac{π}{C_{FD} s} . \end{matrix}$ (18)

(16) is the characteristic equation corresponding to (18). Take P_i as the variable of (17), and the other parameters and variables are constants. The differential of the (17) about P_i is as follows $q_{1}^{'} = {AP}_{i} (y_{eD} - \frac{w}{2 x_{F}}) K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}) .$ (19) (19) is always greater than zero. So (17) is an P_i as a variable, an increasing function. Therefore, according to Difinition 2.13 (see [30]), the solution of the fuzzy defined differential equation (18) is piecewise representation, and the solution is $\begin{matrix} {\tilde{q}}_{1} = 2 A ({\tilde{P}}_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}) . \end{matrix}$ Let $\tilde{P_{i}} = a_{1} + b_{1} \cdot E$ , a₁ > 0, b₁ > 0, where E is a triangular structure element, and the membership function of E is $E (x) = {\begin{matrix} 1 + x & x \in [- 1, 0], \\ 1 - x & x \in (0, 1], \\ 0 & other . \end{matrix}$ (20) Let the function $\begin{matrix} f_{1} (x) = A (x^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}) . \end{matrix}$ then $\begin{matrix} f_{1} (\tilde{P_{i}}) = A ({\tilde{P}}_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}) . \end{matrix}$ Let $\begin{matrix} g_{1} (x) = A [(a_{1} + b_{1} \cdot x)^{2} - P_{ω f}^{2}] (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}) . \end{matrix}$ Where g₁ (x) in [-1, 1] monotone bounded, $f (\tilde{P_{i}}) = g_{1} (E)$ . According to the local mapping theorem, we can obtain $\begin{matrix} s_{1} = g_{1}^{- 1} (x) = - \frac{a_{1}}{b_{1}} + \\ \sqrt{\frac{x}{{Ab}_{1}^{2} (y_{eD} - \frac{w}{2 x_{F}}) K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2})} + \frac{P_{wf}^{2}}{b_{1}^{2}}} . \end{matrix}$ So $μ_{f_{1} ({\tilde{P}}_{i})} (x) = E (g_{1}^{- 1} (x)) = {\begin{matrix} 1 + g_{1}^{- 1} (x) & x \in [α_{1}, α_{2}], \\ 1 - g_{1}^{- 1} (x) & x \in [α_{2}, α_{3}], \\ 0 & other . \end{matrix}$ Where $\begin{matrix} α_{1} = A [(a_{1} - b_{1})^{2} - P_{ω f}^{2}] (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}), \\ α_{2} = A (a_{1}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}), \\ α_{3} = A [(a_{1} + b_{1})^{2} - P_{ω f}^{2}] (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} {hwn}_{F} (B + {CK}_{F}^{2} w^{2}) . \end{matrix}$ The centroid method is also called area centroid or barycenter method. The centroid method is the most popular and exciting method of all non-fuzzy techniques. The algebraic representation of the centroid method is as follows [30]: $z^{*} = \int μ_{\tilde{c}} (z) \cdot zdz / \int μ_{\tilde{c}} (z) dz .$ (21) so $q_{1}^{*} = \int μ_{f_{1} ({\tilde{P}}_{i})} (x) \cdot xdx / \int μ_{f_{1} ({\tilde{P}}_{i})} (x) dx .$ (22) According to the model assumption, initial pressure is a fuzzy number $\tilde{P_{i}} = a_{1} + b_{1} \cdot E$ . In conjunction with the data in Table 1 and the algebraic representation of the centroid method, we can find that $q_{1}^{*}$ is the single numerical solution corresponding to the fuzzy solution set. Figure 3 presents the comparison of gas production between the fuzzy initial formation pressure model and the conventional seepage model. It can be inferred from the figure that the fuzzy seepage model is consistent with the shape of the characteristic curve of the conventional seepage model, but the shale gas yield predicted by the fuzzy seepage model is slightly lower than the conventional seepage model. The reason for this phenomenon is that the conventional seepage model does not take into account the heterogeneity of the reservoir.

Fig.3

Comparison of the gas production q between the conventional seepage model and the fuzzy initial formation pressure seepage model.

3.1.2 Fuzzy formation thickness model

Reservoir thickness refers to the effective thickness of a reservoir whose top and bottom are restricted by non-permeable rock formations whose physical parameters exceed or equal to the lower limit of the reservoir. Figure 4 shows the trend of gas production for different reservoir thickness. It can be seen from the figure that as the thickness of the shale gas reservoir increases, the gas production of the shale gas reservoir also increases.

Fig.4

The trend of gas production for different reservoir thickness h.

As with the formation of conventional oil and gas, the formation of commercial shale oil and gas requires the effective thickness of organic black shale to reach a specific limit. When the thickness of the reservoir is 30–50 m, it is enough to produce commercial airflow. The larger the effective thickness, the larger the total organic matter, and the higher the shale gas enrichment degree. Given the irregular shape of shale gas reservoirs and the extreme heterogeneity of the reservoir itself, this paper uses fuzzy number $\tilde{h}$ to describe the formation thickness.

(16) is the characteristic equation corresponding to (23). Take h as the variable of (17), and the other parameters and variables are constants. Find the derivative of the (17) about h $q_{2}^{'} = A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) K_{F}^{2} {wn}_{F} (B + {CK}_{F}^{2} w^{2})$ (24) The value of (24) is positive. So (17) is an h as a variable, an increasing function. Therefore, the solution of the fuzzy defined differential equation (23) is piecewise representation, and the solution is ${\tilde{q}}_{2} = A (P_{i}^{2} - P_{ω f}^{2}) (y_{e D} - \frac{w}{2 x_{F}}) K_{F}^{2} \tilde{h} w n_{F} (B + C K_{F}^{2} w^{2}) .$ Let $\tilde{h} = a_{2} + b_{2} \cdot E$ , a₂ > 0, b₂ > 0. Let us denote the function as $f_{2} (h) = A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) K_{F}^{2} h {wn}_{F} (B + {CK}_{F}^{2} w^{2}),$ (25) then $\begin{matrix} g_{2} (x) = A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} (a_{2} + b_{2} \cdot x) {wn}_{F} (B + {CK}_{F}^{2} w^{2}) . \end{matrix}$ Where g₂ (x) is monotonous and bounded in [-1, 1], $f_{2} (\tilde{h}) = g_{2} (E)$ . According to the local mapping theorem (see [30]), we obtain $\begin{matrix} s_{2} = g_{2}^{- 1} (x) = - \frac{a_{2}}{b_{2}} + \\ \frac{x}{{Ab}_{2} (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) K_{F}^{2} {wn}_{F} (B + {CK}_{F}^{2} w^{2})}, \end{matrix}$ so $μ_{f_{2} (\tilde{h})} (x) = E (g_{2}^{- 1} (x)) = {\begin{matrix} 1 + g_{2}^{- 1} (x) & x \in [β_{1}, β_{2}], \\ 1 - g_{2}^{- 1} (x) & x \in [β_{2}, β_{3}], \\ 0 & other . \end{matrix}$ Where $\begin{matrix} β_{1} = A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} (a_{2} + b_{2} \cdot x) {wn}_{F} (B + {CK}_{F}^{2} w^{2}), \\ β_{2} = A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} a_{2} {wn}_{F} (B + {CK}_{F}^{2} w^{2}), \\ β_{3} = A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{w}{2 x_{F}}) \\ K_{F}^{2} (a_{2} - b_{2} \cdot x) {wn}_{F} (B + {CK}_{F}^{2} w^{2}) . \end{matrix}$ Using the centroid method to defuzzifier, we get [30]. $q_{2}^{*} = \int μ_{f_{2} (\tilde{h})} (x) \cdot xdx / \int μ_{f_{2} (\tilde{h})} (x) dx .$ (26) According to the above assumption, formation thickness is a fuzzy number $\tilde{h} = a_{2} + b_{2} \cdot E$ . In conjunction with the data in Table 1 and the algebraic representation of the centroid method, we can find that $q_{2}^{*}$ is the single numerical solution corresponding to the fuzzy solution set.

Figure 5 presents a comparison between the fuzzy formation thickness model and the conventional seepage model. It can be inferred that the fuzzy formation thickness seepage model is identical to the curve shape of the conventional seepage model. This phenomenon shows that the fuzzy formation thickness seepage model accords with the seepage law of shale gas reservoirs. Also, it can be seen from the figure that the gas production in the fuzzy formation thickness seepage model is about 5% lower than that of the conventional seepage model. The fuzzy seepage model solves the problem of high prediction data of conventional seepage models.

Fig.5

Comparison of the gas production q between the conventional seepage model and the fuzzy reservoir thickness seepage model.

3.1.3 Fuzzy formation permeability model

Shale gas reservoirs that have been fractured by horizontal wells are composed of the matrix, fractures, and fracturing fractures. The matrix is the primary storage site of shale gas reservoirs, and fractures and fracturing fractures are the main migration sites of shale gas reservoirs. Among them, the permeability of fracturing fractures is the largest, and the permeability of matrix and fractures is several orders of magnitude different from that of fracturing fractures. Figure 6 shows the trend of gas production for different reservoir permeability of fracturing fractures. It can be seen from the figure that as the permeability of the shale gas reservoir increases, the gas production of the shale gas reservoir also increases. This indicates that the reservoir permeability is a significant factor affecting the productivity of shale gas reservoirs. Therefore, this paper uses fuzzy numbers to describe the permeability of fracturing fractures. Given the irregular shape of shale gas reservoirs and the extreme heterogeneity of the reservoir itself, this paper uses fuzzy number ${\tilde{K}}_{F}$ to describe the reservoir permeability of fracturing fractures.

Fig.6

The trend of gas production for different permeability K₀ of large fractures.

Figure 7 presents a comparison between the fuzzy formation permeability model and the conventional seepage model. It can be inferred that the fuzzy formation permeability model is identical to the curve shape of the conventional seepage model. These curves show that the fuzzy formation permeability model accords with the seepage law of shale gas reservoirs. Also, it can be seen from Fig. 7 that the reservoir pressure in the conventional seepage model is about 10% higher than that of the fuzzy formation permeability model. The difference between the two models is mainly due to the fuzzy seepage model using fuzzy numbers to describe the formation permeability of fracturing fractures.

Fig.7

Comparison of the gas production between the fuzzy formation permeability model and the conventional seepage model.

3.1.4 Fuzzy fractures number model

Horizontal well drilling and staged fracturing are currently the core technologies used for shale gas production [31, 32]. Horizontal well fracturing technology [31 –35]: When the fracturing fluid enters the reservoir through the perforation hole at high speed, perforation friction will occur and increase with the increase of pump displacement, which will increase the bottom hole pressure. When the bottom hole pressure exceeds the fracture pressure of multiple fractured intervals, fractures are fractured at each interval. Figure 8 shows the trend of gas production for different fractures numbers. It can be seen from the figure that as the number of fractures increases, the gas production of shale gas reservoirs also increases significantly. This phenomenon indicates that fracturing technology is of vital importance to the exploitation of shale gas reservoirs.

Fig.8

The trend of gas production for different fractures number n_F.

Because the heterogeneity of the shale gas reservoir makes the number of fractured fractures uncertain, therefore, this paper uses fuzzy number ${\tilde{n}}_{F}$ to describe the number of fractures.

Figure 9 presents a comparison between the fuzzy fractures number model and the conventional seepage model. It can be inferred that the fuzzy seepage model is identical to the curve shape of the conventional seepage model. These curves show that the fuzzy gas viscosity seepage model accords with the seepage law of shale gas reservoirs. Also, it can be seen from the figure that the reservoir pressure in the conventional seepage model is about 30% higher than that of the fuzzy fractures number model.

Fig.9

Comparison of the gas production between the fuzzy fractures number model and the conventional seepage model.

3.1.5 Fuzzy main fracture width model

The main fracture width refers to the width of hydraulic fracturing fractures in shale gas reservoirs. Figure 10 shows the trend of gas production for different fracture width. It can be seen from the figure that a slight increase in the width of fractures has caused a significant increase in the production of shale gas reservoirs. For the extreme non-uniformity of shale reservoirs and the irregularity of reservoir shape, this paper uses the fuzzy number $\tilde{w}$ to describe the main fracture width of shale reservoirs. Based on this, the fuzzy main fracture width seepage model is established as follows.

Fig.10

The trend of gas production for different fracture width w.

${\begin{matrix} \frac{\partial^{2} {\bar{m}}_{mD} (x_{D}, s)}{\partial x_{D}^{2}} - \frac{s}{η_{mD}} {\bar{m}}_{mD} (x_{D}, s) = 0, \\ \frac{\partial^{2} {\bar{m}}_{fD} (y_{D}, s)}{\partial y_{D}^{2}} + 2 \frac{K_{m}}{K_{f}} \frac{\partial {\bar{m}}_{mD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 1} = s {\bar{m}}_{fD} (y_{D}, s), \\ \frac{\partial^{2} {\bar{m}}_{FD} (x_{D}, s)}{\partial x_{D}^{2}} + \frac{4}{{\tilde{C}}_{FD}} \frac{\partial {\bar{m}}_{fD} (y_{D}, s)}{\partial y_{D}} |_{y_{D} = \frac{{\tilde{w}}_{D}}{2}} = \frac{s}{η_{FD}} {\bar{m}}_{FD} (x_{D}, s), \\ \frac{\partial {\bar{m}}_{mD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = x_{eD}} = \frac{\partial {\bar{m}}_{fD} (y_{D}, s)}{\partial y_{D}} |_{y_{D} = y_{eD}} = 0 \\ \frac{\partial {\bar{m}}_{FD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 1} = 0, q = q_{F} n_{F} \\ {\bar{m}}_{mD} (x_{D}, s) |_{x_{D} = 1} = {\bar{m}}_{fD} (x_{D}, s) |_{x_{D} = 1} \\ {\bar{m}}_{fD} (y_{D}, s) |_{y_{D} = {\tilde{w}}_{D} / 2} = {\bar{m}}_{FD} (y_{D}, s) |_{y_{D} = {\tilde{w}}_{D} / 2}, \\ {\bar{m}}_{FD} (x_{D}, s) |_{x_{D} = 0} = {\bar{m}}_{wD}, \frac{\partial {\bar{m}}_{FD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 0} = - \frac{π}{{\tilde{C}}_{FD} s} . \end{matrix}$ (35) (16) is the characteristic equation corresponding to (35). Take w as the variable of the (17), and the other parameters and variables are constants. Find the derivative of the (17) about w $q_{5}^{'} = - \frac{2 A_{5} C_{5} w^{3}}{x_{F}} + 3 A_{5} C_{5} y_{e D} w^{2} - \frac{A_{5} B_{5} w}{x_{F}} + A_{5} B_{5} y_{e D}$ (36) where ${\begin{matrix} A_{5} = \frac{6.89 \times 10^{6} (P_{i}^{2} - P_{ω f}^{2}) K_{F}^{2} {hx}_{e} K_{m}^{2} n_{F} t_{D}^{- \frac{1}{4}}}{π y_{eD} η_{m} T μ {zK}_{f} x_{F}}, \\ B_{5} = 2 η_{f} (x_{eD} - 1) + C_{fD} y_{eD} η_{m}, \\ C_{5} = \frac{K_{F}^{2} y_{eD} η_{m}}{K_{f}^{2} x_{F}^{2}} . \end{matrix}$ The value of (36) is positive. So (17) is an w as a variable, an increasing function. Therefore, the solution of the fuzzily defined differential equation (35) is piecewise representation, and the solution is ${\tilde{q}}_{5} = A (P_{i}^{2} - P_{ω f}^{2}) (y_{e D} - \frac{\tilde{w}}{2 x_{F}}) K_{F}^{2} \tilde{w} h n_{F} (B + C K_{F}^{2} {\tilde{w}}^{2}) .$ Let $\tilde{w} = a_{5} + b_{5} \cdot E$ , a₅ > 0, b₅ > 0. If $f_{5} (\tilde{w}) = A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{\tilde{w}}{2 x_{F}}) K_{F}^{2} \tilde{w} h n_{F} (B + {CK}_{F}^{2} {\tilde{w}}^{2}),$ (37) then $\begin{matrix} g_{5} (x) = A (P_{i}^{2} - P_{ω f}^{2}) [y_{eD} - \frac{(a_{5} + b_{5} \cdot x)}{2 x_{F}}] \\ (a_{5} + b_{5} \cdot x) K_{F}^{2} h n_{F} [B + {CK}_{F}^{2} (a_{5} + b_{5} \cdot x)^{2}], \end{matrix}$ Where g₅ (x) is monotonous and bounded in [-1, 1], $f_{5} (\tilde{w}) = g_{5} (E)$ . According to the local mapping theorem (see [30]), we obtain $\begin{matrix} s_{5} = g_{5}^{- 1} (x) = \frac{H_{5}^{\frac{1}{2}}}{6 H_{4} \frac{1}{6}} - \frac{4 A_{5} C_{5} a_{5} b_{5}^{3} - 2 A_{5} C_{5} b_{5}^{3} x_{F} y_{ed}}{4 A_{5} C_{5} b_{5}^{4}} \\ + \frac{1}{6 H_{4}^{\frac{1}{6}} H_{5}^{\frac{1}{4}}} (12 H_{2} H_{5}^{\frac{1}{2}} - 9 H_{4}^{\frac{2}{3}} H_{5}^{\frac{1}{2}} \\ + 3 \sqrt{6} H_{7} H_{6}^{\frac{1}{2}} - H_{1}^{2} H_{5}^{\frac{1}{2}} - 12 H_{1} H_{4}^{\frac{1}{3}} H_{5}^{\frac{1}{2}})^{\frac{1}{2}} . \end{matrix}$ (38) Where $\begin{matrix} H_{1} = \frac{B_{5}}{C_{5} b_{5}^{2}} - \frac{3 x_{F}^{2} y_{eD}^{2}}{2 b_{5}^{2}}, \\ H_{2} = \frac{3 x_{F}^{4} y_{eD}^{4}}{16 b_{5}^{4}} - \frac{2 x x_{F}}{A_{5} C_{5} b_{5}^{4}} + \frac{3 B_{5} x_{F}^{2} y_{eD}^{2}}{4 C_{5} b_{5}^{4}}, \\ H_{3} = 16 H_{1}^{4} H_{2} + 4 (\frac{x_{F}^{3} y_{eD}^{3}}{b_{5}^{3}} + \frac{B_{5} x_{F} y_{eD}}{C_{5} b_{5}^{3}})^{2} \\ H_{1}^{3} + 27 (\frac{x_{F}^{3} y_{eD}^{3}}{b_{5}^{3}} + \frac{B_{5} x_{F} y_{eD}}{C_{5} b_{5}^{3}})^{4} + 128 H_{1}^{2} H_{2}^{2} + 256 H_{2}^{3} \\ + 144 (\frac{x_{F}^{3} y_{eD}^{3}}{b_{5}^{3}} + \frac{B_{5} x_{F} y_{eD}}{C_{5} b_{5}^{3}})^{2} H_{1} H_{2}, \\ H_{4} = \frac{1}{2} (\frac{x_{F}^{3} y_{eD}^{3}}{b_{5}^{3}} + \frac{B_{5} x_{F} y_{eD}}{C_{5} b_{5}^{3}})^{2} + \frac{H_{1}^{3}}{27} + \frac{4 H_{1} H_{2}}{3} + \frac{\sqrt{3 H_{3}}}{18}, \\ H_{5} = 9 H_{4}^{\frac{2}{3}} + H_{1}^{2} - 6 H_{1} H_{4}^{\frac{1}{3}} - \frac{9 x_{F}^{4} y_{eD}^{4}}{4 b_{5}^{4}} + \frac{24 x x_{F}}{A_{5} C_{5} b_{5}^{4}} - \frac{9 B_{5} x_{F}^{2} y_{eD}^{2}}{C_{5} b_{5}^{4}}, \\ H_{6} = 27 (\frac{x_{F}^{3} y_{eD}^{3}}{b_{5}^{3}} + \frac{B_{5} x_{F} y_{eD}}{C_{5} b_{5}^{3}})^{2} + 2 H_{1}^{3} + 72 H_{1} H_{2} + 3 \sqrt{3 H_{3}}, \\ H_{7} = \frac{x_{F}^{3} y_{eD}^{3}}{b_{5}^{3}} + \frac{B_{5} x_{F} y_{eD}}{C_{5} b_{5}^{3}} . \end{matrix}$ so $μ_{f_{5} (\tilde{w})} (x) = E (g_{5}^{- 1} (x)) = {\begin{matrix} 1 + g_{5}^{- 1} (x), & x \in [ɛ_{1}, ɛ_{2}] \\ 1 - g_{5}^{- 1} (x), & x \in [ɛ_{2}, ɛ_{3}] \\ 0 & other . \end{matrix}$ $\begin{matrix} ɛ_{1} & = & A (P_{i}^{2} - P_{ω f}^{2}) [y_{eD} - \frac{(a_{5} - b_{5})}{2 x_{F}}] K_{F}^{2} (a_{5} - b_{5}) h n_{F} \\ [B + C K_{F}^{2} (a_{5} - b_{5})^{2}], \\ ɛ_{2} & = & A (P_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{a_{5}}{2 x_{F}}) K_{F}^{2} a_{5} h n_{F} (B + C K_{F}^{2} a_{5}^{2}), \\ ɛ_{3} & = & A (P_{i}^{2} - P_{ω f}^{2}) [y_{eD} - \frac{(a_{5} + b_{5})}{2 x_{F}}] K_{F}^{2} (a_{5} + b_{5}) h n_{F} \\ [B + C K_{F}^{2} (a_{5} + b_{5})^{2}] . \end{matrix}$ And $q_{5}^{*} = \int μ_{f_{5} (\tilde{w})} (x) \cdot xdx / \int μ_{f_{5} (\tilde{w})} (x) dx .$ (39) According to the above assumption, the main fracture width is a fuzzy number $\tilde{w} = a_{5} + b_{5} \cdot E$ . In conjunction with the data in Table 1 and the algebraic representation of the centroid method, we can find that $q_{5}^{*}$ is the single numerical solution corresponding to the fuzzy solution set. Figure 11 presents the comparison between the fuzzy main fracture width model and the conventional seepage model. It can be inferred that the fuzzy seepage model is identical to the curve shape of the conventional seepage model. These curves show that the fuzzy main fracture width model accords with the seepage law of shale gas reservoirs. Also, it can be seen from the figure that the reservoir pressure in the conventional seepage model is about 16% higher than that of the fuzzy main fracture width model.

Fig.11

Comparison of the gas production between the fuzzy main fracture width model and the conventional seepage model.

3.2 Multiple fuzzy parameters trilinear seepage model

The five parameters (fuzzy initial formation pressure, fuzzy reservoir thickness, fuzzy formation permeability, fuzzy fractures number, and fuzzy main fracture width) mentioned in Section 3.1 seepage model exert a significant impact on shale gas reservoir mining. Five single fuzzy parameter seepage models are obtained by describing the uncertainty parameters with fuzzy numbers. After in-depth research and comparative analysis, it is more appropriate to use fuzzy numbers to describe the five parameters. Therefore, based on the previous study, this paper establishes the multiple fuzzy parameters seepage model to further enchance the accuracy of the production prediction of a shale gas reservoir. ${\begin{matrix} \frac{\partial^{2} {\tilde{\bar{m}}}_{mD} (x_{D}, s)}{\partial x_{D}^{2}} - \frac{s}{η_{mD}} {\tilde{\bar{m}}}_{mD} (x_{D}, s) = 0, \\ \frac{\partial^{2} {\tilde{\bar{m}}}_{fD} (y_{D}, s)}{\partial y_{D}^{2}} + 2 \frac{K_{m}}{K_{f}} \frac{\partial {\tilde{\bar{m}}}_{mD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 1} = s {\tilde{\bar{m}}}_{fD} (y_{D}, s), \\ \frac{\partial^{2} {\bar{m}}_{FD} (x_{D}, s)}{\partial x_{D}^{2}} + \frac{4}{{\tilde{C}}_{FD}} \frac{\partial {\tilde{\bar{m}}}_{fD} (y_{D}, s)}{\partial y_{D}} |_{y_{D} = \frac{{\tilde{w}}_{D}}{2}} = \frac{s}{η_{FD}} {\tilde{\bar{m}}}_{FD} (x_{D}, s), \\ \frac{\partial {\tilde{\bar{m}}}_{mD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = x_{eD}} = \frac{\partial {\tilde{\bar{m}}}_{fD} (y_{D}, s)}{\partial y_{D}} |_{y_{D} = y_{eD}} = 0 \\ \frac{\partial {\tilde{\bar{m}}}_{FD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 1} = 0, q = q_{F} {\tilde{n}}_{F} \\ {\tilde{\bar{m}}}_{mD} (x_{D}, s) |_{x_{D} = 1} = {\tilde{\bar{m}}}_{fD} (x_{D}, s) |_{x_{D} = 1} \\ {\tilde{\bar{m}}}_{fD} (y_{D}, s) |_{y_{D} = {\tilde{w}}_{D} / 2} = {\tilde{\bar{m}}}_{FD} (y_{D}, s) |_{y_{D} = {\tilde{w}}_{D} / 2}, \\ {\tilde{\bar{m}}}_{FD} (x_{D}, s) |_{x_{D} = 0} = {\tilde{\bar{m}}}_{wD}, \frac{\partial {\tilde{\bar{m}}}_{FD} (x_{D}, s)}{\partial x_{D}} |_{x_{D} = 0} = - \frac{π}{{\tilde{C}}_{FD} s} . \end{matrix}$ (40)

In fact, from the definition of the characteristic equation, (16) is also the characteristic equation corresponding to (40). Moreover, the solution of the characteristic equation is (17). According to the above research, parameters P_i, h, k_F, n_F, and w satisfy the conditions of monotonic. Therefore, the solution of the seepage model of the multiple fuzzy parameters shale gas reservoir can be accurately expressed, and the solution is $\tilde{q} = A ({\tilde{P}}_{i}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{\tilde{w}}{2 x_{F}}) {\tilde{K}}_{F}^{2} \tilde{w} \tilde{h} {\tilde{n}}_{F} (B + C {\tilde{K}}_{F}^{2} {\tilde{w}}^{2}) .$ If $\begin{matrix} g (x) = A [(a_{1} + b_{1} \cdot x)^{2} - P_{ω f}^{2}] [y_{eD} - \frac{(a_{5} + b_{5} \cdot x)}{2 x_{F}}] \\ (a_{3} + b_{3} \cdot x)^{2} (a_{5} + b_{5} \cdot x) (a_{2} + b_{2} \cdot x) \\ (a_{4} + b_{4} \cdot x) [B + C (a_{3} + b_{3} \cdot x)^{2} (a_{5} + b_{5} \cdot x)^{2}] . \end{matrix}$ and g (x) in [-1, 1] monotone bounded, we obtain g^-1 (x) are the root. $\sum_{i = 0}^{12} Q_{i} z_{i} = 0$ (41) The specific form of Q_i is shown in Appendix A (where J = A · C, I = A · B and $Q_{i} = M_{i 1} x_{F} y_{eD} + M_{i 2} + M_{i 3} P_{wf}^{2} + M_{i 4} P_{wf}^{2} x_{F} y_{eD}$ ,i = 0, . . . , 12). So $μ_{g (E)} (x) = E (g^{- 1} (x)) = {\begin{array}{l} 1 + g^{- 1} (x) & x \in [ζ_{1}, ζ_{2}], \\ 1 - g^{- 1} (x) & x \in [ζ_{2}, ζ_{3}], \\ 0 & other . \end{array}$ Where $\begin{matrix} ζ_{1} = & A [(a_{1} - b_{1})^{2} - P_{ω f}^{2}] [y_{eD} - \frac{(a_{5} - b_{5})}{2 x_{F}}] (a_{3} - b_{3})^{2} \\ (a_{5} - b_{5}) (a_{2} - b_{2}) (a_{4} - b_{4}) \times [B + C (a_{3} - b_{3})^{2} (a_{5} - b_{5})^{2}], \\ ζ_{2} = & A (a_{1}^{2} - P_{ω f}^{2}) (y_{eD} - \frac{a_{5}}{2 x_{F}}) a_{3}^{2} a_{5} a_{2} a_{4} (B + C a_{3}^{2} a_{5}^{2}), \\ ζ_{3} = & A [(a_{1} + b_{1})^{2} - P_{ω f}^{2}] [y_{eD} - \frac{(a_{5} + b_{5})}{2 x_{F}}] (a_{3} + b_{3})^{2} \\ (a_{5} + b_{5}) (a_{2} + b_{2}) (a_{4} + b_{4}) \times [B + C (a_{3} + b_{3})^{2} (a_{5} + b_{5})^{2}] . \end{matrix}$ Similar to Section 3.1, the fuzzy solution set is defuzzified by the centroid method. $P^{*} = \int μ_{g (E)} (x) \cdot xdx / \int μ_{g (E)} (x) (x) dx .$ (42) Based on the previous assumptions, some reservoir parameters are fuzzy number, ${\tilde{P}}_{i} = a_{1} + b_{1} \cdot E$ , $\tilde{h} = a_{2} + b_{2} \cdot E$ , ${\tilde{k}}_{0} = a_{3} + b_{3} \cdot E$ , ${\tilde{n}}_{F} = a_{4} + b_{4} \cdot E$ , and $\tilde{w} = a_{5} + b_{5} \cdot E$ . In conjunction with the data in Table 1 and the algebraic representation of the centroid method, we can find that ${\tilde{q}}^{*}$ is the single numerical solution corresponding to the fuzzy solution set.

The comparison analysis of single parameter fuzzy seepage models, multiple parameter fuzzy seepage model, and conventional seepage model is shown in Fig. 12. It can be inferred that the fuzzy seepage model is identical to the curve shape of the conventional seepage model. In addition, in the early times, the choice of the model does not affect gas production. Figure 12 also reveals that the gas production in fuzzy seepage models is lower than that of the conventional seepage model, and the multiple parameters fuzzy seepage model is lower than that of the single parameter fuzzy seepage models. The results show that compared with the conventional seepage model and the single fuzzy parameter seepage model, the multiple fuzzy parameters seepage model is closer to the shale gas reservoir on-site production.

Fig.12

Comparison of the gas production among single fuzzy parameter seepage models, multiple fuzzy parameters seepage model and the conventional seepage model.

4 Sensitivity of pressure distribution

Figure 13 presents the effect of bottom hole pressure, P_wf, on gas production. Four cases are investigated, P_wf = 8MPa, P_wf = 9MPa, P_wf = 10MPa, and P_wf = 11MPa. The figure reveals that increasing the temperature results in greater formation pressure. Increasing the bottom hole pressure from P_wf = 8MPa to P_wf = 11MPa results in about 33% decrease in gas production.

Fig.13

The effect of bottomhole pressure on the gas well productivity.

Figure 14 presents the relationship between the time and the gas production under different temperature, T. It can be inferred that the temperature has little effect on gas production in early mining. However, as mining time increases, the effect of the temperature, T, on gas production becomes more pronounced. When the time reaches a certain value, the influence of the temperature on the gas production gradually becomes stable. Figure 14 also reveals that as the temperature increases, the gas production gradually decreases.

Fig.14

The effect of temperature on the gas well productivity.

5 Conclusions

This paper primarily employs the fuzzy set theory into the study on seepage characteristics of shale gas reservoirs. The following conclusions are obtained from work presented in this paper.

Introduce the fuzzy concepts (fuzzy initial pressure, fuzzy reservoir thickness, fuzzy permeability, fuzzy fractures number, and fuzzy main fracture width) into the seepage model of shale gas reservoirs, on the one hand, the heterogeneity of the reservoir and the complex variability in the mining process can be fully considered. On the other hand, errors due to inaccurate laboratory measurements can be avoided.

Based on the proposed fuzzy concepts, a seepage model of multiple fuzzy parameters shale gas reservoir is established. The seepage model of the multiple fuzzy parameters is essentially a second-order fuzzy differential problem. There are many methods for solving fuzzy differentials. This paper uses the fuzzy structure element method to solve the seepage model to handle ergodic operations. After obtaining the fuzzy solution set, the fuzzy solution set is further defuzzification. By comparing the new model with the other five single fuzzy parameter models, and the conventional seepage model, the accuracy and superiority of the new model can be obtained.

To investigate the impact of various parameters on the pressure distribution of shale gas reservoirs, we conduct a series of sensitivity studies.

To further improve the accuracy of shale gas reservoir models, a perfect oil and gas reservoir seepage model should be taken into account the process of gas dispersion, temperature propagation, and other models of calculation. These problems still exist with ambiguity, which has to be solved in future research. There are many main control parameters with uncertainty except the five parameters(fuzzy initial pressure, fuzzy reservoir thickness, fuzzy permeability, fuzzy fractures number, and fuzzy main fracture width). So the seepage model of multiple fuzzy parameters also needs our further research.

Footnotes

Acknowledgments

The authors would like to thank the anonymous referees for their valuable comments and suggestions. This work is supported by the National Natural Science Foundation of China (No. 11361050 and No.11671284).

A Statement of the Q i

See Tables 2 to 11.

References

Rogers

J.H.

, Natural gas from shale formation – The evolution, evidences and challenges of shale gas revolution in United States, Renewable and Sustainable Energy Reviews, 2014.

Curtis

J.B.

, Fractured shale-gas systems, AAPG Bulletin 86(11) (2002), 1921–1938.

Ross

D.J.K.

and Bustin

R.M.

, The importance of shale composition and pore structure upon gas storage potential of shale gas reservoirs, Marine and Petroleum Geology 26(6) (2009), 916–927.

Vengosh

, Jackson

R.B.

, Warner

, et al., A critical review of the risks to water resources from unconventional shale gas development and hydraulic fracturing in the United States, Environmental Science and Technology 48(15) (2014), 8334–8348.

Jackson

R.B.

, Vengosh

, Darrah

T.H.

, et al., Increased Stray Gas Abundance in a Subset of Drinking Water Wells near Marcellus Shale Gas Extraction, Proceedings of the National Academy of Sciences 110(28) (2013).

Barenblatt

G.I.

, Monteiro

P.J.M.

and Rycroft

C.H.

, On a boundary layer problem related to the gas flow in shales, Journal of Engineering Mathematics 84(1) (2014), 11–18.

Lee

S.T.

and Brockenbrough

J.R.

, A New Approximate Analytic Solution for Finite-Conductivity Vertical Fractures, SPE Formation Evaluation 1(01) (1986), 75–88.

Ozkan

, Brown

M.L.

, Raghavan

, et al., Comparison of Fractured-Horizontal-Well Performance in Tight Sand and Shale Reservoirs, SPE Reservoir Evaluation and Engineering 14(02) (2011), 248–259.

Zhao

Y.L.

, Zhang

L.H.

, Zhao

J.Z.

, et al., “Triple porosity” modeling of transient well test and rate decline analysis for multi-fractured horizontal well in shale gas reservoirs, Journal of Petroleum Science and Engineering 110 (2013), 253–262.

10.

Qiu

, Lu

, Zhang

, et al., Algebraic properties and topological properties of the quotient space of fuzzy numbers based on Mare? equivalence relation, Fuzzy Sets and Systems 245 (2014), 63–82.

11.

Avner, Vengosh, Robert, et al., A Critical Review of the Risks to Water Resources from Unconventional Shale Gas Development and Hydraulic Fracturing in the United States, Environmental Science and Technology: ES and T 48(15) (2014), 8334–8348.

12.

Bede

and Gal

S.G.

, Generalizations of the differentiability of fuzzy-number-valued functions with applications to fuzzy differential equations, Fuzzy Sets and Systems 151(3) (2005), 581–599.

13.

Hai

, Gong

and Li

, Generalized differentiability for n-dimensional fuzzy-number-valued functions and fuzzy optimization, Information Sciences 2016, 374.

14.

Allahviranloo

and Chehlabi

, Solving fuzzy differential equations based on the length function properties, Soft Computing 19(2) (2015), 307–320.

15.

Song

and Wu

, Existence and uniqueness of solutions to Cauchy problem of fuzzy differential equations, Fuzzy Sets and Systems 110(1) (2000), 55–67.

16.

Bede

and Stefanini

, Generalized differentiability of fuzzy-valued functions, Fuzzy Sets and Systems 230(Complete) (2013), 119–141.

17.

Chehlabi

and Allahviranloo

, Positive or negative solutions to first-order fully fuzzy linear differential equations under generalized differentiability, Applied Soft Computing 70 (2018), 359–370.

18.

, Song

and Lee

E.S.

, Approximate Solutions, Existence, and Uniqueness of the Cauchy Problem of Fuzzy Differential Equations, Journal of Mathematical Analysis and Applications 202(2) (1996), 629–644.

19.

Vorobiev

and Seikkala

, Towards the theory of fuzzy differential equations, Fuzzy Sets and Systems 125(2) (2002), 231–237.

20.

Chalco-Cano

and Román-Flores

, Comparation between some approaches to solve fuzzy differential equations, Fuzzy Sets and Systems 160(11) (2009), 1517–1527.

21.

Lozano

and Wang

P.P.

, Fuzzy Set Theory Applied to Accounting Sciences, New Mathematics and Natural Computation (NMNC), 2020, 16.

22.

Arya

and Kumar

, A new picture fuzzy information measure based on shannon entropy with applications in opinion polls using extended VIKOR– TODIM approach, Computational and Applied Mathematics 39(3) (2020), 1–24.

23.

Yan, Zhijia, Yang, et al., A Reputation-Enhanced Hybrid Approach for Supplier Selection with Intuitionistic Fuzzy Evaluation Information, Mathematics, 2019.

24.

Zhang

, Shu

and Li

, Fuzzy structural element method for solving fuzzy dual medium seepage model in reservoir, Soft Computing 2020(3).

25.

Zhang

, Yu

, Shu

, et al., Fuzzy Structural Element Method for Solving Nonlinear Fuzzy Seepage Model of Shale Gas Reservoir, 2019 IEEE 14th International Conference on Intelligent Systems and Knowledge Engineering (ISKE), IEEE, 2019.

26.

, Liang

and Guo

Y.H.

, The Fuzzy Seepage Theory Based on Fuzzy Structuring Element Method, International Conference on Fuzzy Information and Engineering, Liaoning Technical University, 123000, Fu-xin, Liao-ning, China, 2007.

27.

NOVAKOVI, Aleksandar, RANKOVI, et al., DEVELOPMENT OF NEURO-FUZZY MODEL FOR DAM SEEPAGE ANALYSIS, Annals of the Faculty of Engineering – Hunedoara International Journal of Engineering 12(2) (2014), 133–136.

28.

Chongshi

G.U.

, Zhang

, Cai

, et al., Application of entropy-based fuzzy matter-element analysis in seepage monitoring of RCC dam, Frontiers of Architecture and Civil Engineering in China 5(001) (2011), 105–111.

29.

L.G.

and Peng

D.H.

, Interval-Valued Hesitant Fuzzy Hamacher Synergetic Weighted Aggregation Operators and Their Application to Shale Gas Areas Selection, Mathematical Problems in Engineering 2014 (2014), 1–25.

30.

Zhang

, Qiu

and Shu

, Fuzzy Structure Element Method for Solving Fuzzy Trilinear Seepage Model of Shale Gas Reservoir, International Journal of Fuzzy Systems 21(5) (2019), 1597–1612.

31.

Xie

, Li

, Zhang

, et al., Two-phase pressure transient analysis for multi-stage fractured horizontal well in shale gas reservoirs, Journal of Natural Gas Science and Engineering 21 (2014), 691–699.

32.

Chen

, Jiang

, Chen

, et al., Pressure response and production performance of volumetric fracturing horizontal well in shale gas reservoir based on boundary element method, Engineering Analysis with Boundary Elements 87(feb.) (2018), 66–77.

33.

Wei-Yang

, Xiao-Ping

, Lie-Hui

, et al., Production Decline Analysis for Two-Phase Flow in Multifractured Horizontal Well in Shale Gas Reservoirs, Journal of Chemistry 2015, (2015-9-27), 2015 (2015), 1–10.

34.

Tavassoli

, Yu

, Javadpour

, et al., Well Screen and Optimal Time of Refracturing: A Barnett Shale Well, Journal of Petroleum Engineering 2013, (2013-4-18), 2013 (2013), 1–10.

35.

Olmstead

S.M.

, Muehlenbachs

L.A.

, Shih

J.S.

, et al., Shale gas development impacts on surface water quality in Pennsylvania, Proceedings of the National Academy of Sciences of the United States of America 110(13) (2013), 4962–4967.

36.

Sang

, Chen

, Yang

, et al., A new mathematical model considering adsorption and desorption process for productivity prediction of volume fractured horizontal wells in shale gas reservoirs, Journal of Natural Gas Science and Engineering 19 (2014), 228–236.