Hybrid Coupled Multifracture and Multicontinuum Models for Shale Gas Simulation by Use of Semi-Analytical Approach

Abstract

Combining the advantages of multicontinuum and multifracture representations provides an easy-to-use tool to adequately capture the characteristic of the multiscaled fracture system in shale gas reservoir. A hybrid model is established on the basis of simplified conceptual productivity assumption, where the matrix volume is divided into two sub-domains (triple-porosity model and dual-depletion flowing model) and the fracture volume is represented by discrete finite conductivity fracture. In addition, the mechanisms of instant desorption, viscous flow and dual-depletion in matrix are taken into account. The rate transient responses are then obtained by use of semi-analytical approach. Based on the model, type curves are plotted and verified by comparing with alternative reliable methods. Different flow regimes in shale gas reservoirs can be identified and detected. The Generalized Likelihood Uncertainty Estimation methodology, based on probabilistic aggregation theory, is employed to integrating those two productivity models together such that the production can be predicted more accurately. A field example is applied to validate the applicability of this new model. Finally, it is concluded that the proposed model can predict the rate and cumulative rate more easily and practically.

1. Introduction

Shale gas reservoir is typical unconventional reservoir due to its ultra-low permeability and porosity. Generally speaking, there is no natural productive capacity for such reservoirs. Multi-stage fracturing technology for horizontal wells can activate the pre-existing natural fractures and generate a complex fracture network called stimulated reservoir volume (SRV), which has become the most effective method to develop shale resources [1]. Generally, SRV is defined as the part of the drainage area impacted by hydraulic fractures and concomitant reactive fractures, where the productivity of shale gas well could be effectively improved [2,3,4,5,6].

To characterize the transient performance for shale gas reservoirs, many scholars have made some efforts to simulate the fluids flow and transport behaviors in such complexly-structured fracture network through analytical and semi-analytical solutions [7,8,9,10]. These solutions are mostly proposed in the terms of multi-porosity model and multi-linear model. For the multi-porosity model, Barenblatt [11] and Warren and Root [12] firstly proposed the dual-porosity model to simulate the fluid exchange between matrix and fractures, where the exchange is assumed in the pseudo steady state. Then, alternative transient models were proposed for simulating the fluid transfer between matrix and fractures [13,14,15], which relaxed the pseudo-steady assumption. In addition, various triple-porosity models were subsequently established to overcome the drawback of the dual-porosity model in terms of multi-scale porosity description, such as dual-fracture/triple-porosity models [16] and improved triple-porosity/permeability model [17,18,19]. Next, for multi-linear model, numerous field practices demonstrate that many fractured reservoirs exhibit linear flow, which could last for several years. El-Banbi [20] proposed a linear dual-porosity model for linear fractured reservoirs without considering the impact of desorption, diffusion. Hasan and Al-Ahmadi [21] proposed a triple-porosity linear flow model with consideration of the impact of shale gas desorption and diffusion. Zhao et al [22] proposed a triple-porosity spherical flow model for the fractured infinite shale gas reservoirs by considering the impact of diffusion and desorption. Obinna and Hassan [23] proposed the quadrilinear flow mode to simultaneously the depletion from the matrix into natural and hydraulic fractures.

Nevertheless, it is very difficult to use these multi-porosity/linear models mentioned above to investigate the actual performance for shale gas reservoirs. It is caused by the fact that the characterization of diffusion and adsorption is not taken into account. Therefore, it is necessary to comprehensively take the impact of various flowing mechanisms and multi-porosity characteristics into account to obtain critical dynamic parameters for shale gas reservoirs. In this paper, for the purpose of obtaining accurate results of performance analysis, the main work is twofold: proposing two improved models to simulate production performance, and presenting a reliable probabilistic algorithm to integrate these two models. It is noted that the fundamental model is established on the basis of the work presented by Obinna and Hassan [23]. Some improvements are made: (1) SRV zone is simulated as a triple-porosity cubic model by taking into account instant adsorption and viscous flow; (2) two kinds of simplified conceptual models are generated to divide the matrix volume into two sub-domains for the purpose of allowing simultaneous matrix-hydraulic fracture (HF) and matrix-natural fracture (NF) depletion. Subsequently, the Generalized Likelihood Uncertainty Estimation (GLUE) approach is built, based on probabilistic aggregation theory, in order to reduce the uncertainty of performance analysis by integrating those two models. Finally, the numerical simulation and field application are applied to validate this new simplified aggregation model. The improved method in this paper could provide a robust methodology approach to quantify the uncertainty of production prediction.

2. Productivity Model

A schematic for the improved model is presented in Figure 1. Figure 1a shows the model consists of SRV region, unstimulated reservoir volume (USRV) region, and fracture region. Figure 1b shows the simultaneous matrix-depletion into HF and NF. In addition, some assumptions are emphasized as follows:

• SRV region is simplified as a cubic triple-porosity model, containing natural fractures, hydraulic fractures and the matrix;
• Hydraulic fractures are perpendicular to the horizontal well and evenly distributed along the wellbore, and the natural fractures are perpendicular to the hydraulic fracture. Horizontal wellbore are equal to L, and the length of hydraulic fracture and the width of reservoir are equal to y_e;
• Hydraulic fracture is finite conductivity and assumed to be penetrated fully;
• Only the fluid flowing from hydraulic fractures to wellbore is considered;
• Simultaneous matrix-depletion into HF and NF is assumed pseudo-steady state, and the exchange between HF and NF is assumed unsteady state;
• The effect of gravity and capillary pressure is not taken into account;
• Gas is slightly compressible and the compressibility coefficient is constant;
• This paper considers the simultaneous depletion from the matrix into HF and NF. The matrix in SRV region is artificially divided into two distinct segments which are denoted as sub-matrix m1 and sub-matrix m2 respectively. The depletion process from the matrix to HF and NF in SRV region (seen in Figure 1b);
• Sub-matrix m1 feeds the HF via inter-porosity exchange;
• Sub-matrix m2 feeds the NF via inter-porosity exchange.

The conceptual model can get some further detailed divisions based on different methods. Hence, two different approaches would be developed to characterize dual-depletion in the following section.

2.1. Conceptual Model

2.1.1. Conceptual Model 1

The matrix is divided into two sections. Some further assumptions are made to derive the following formulas:

• Sub-matrix m1 and sub-matrix m2 mix up with each other evenly, which means that both of them have the same width L_f, length L_F and thickness H (see the Figure 2a);
• Sub-matrix m1 and sub-matrix m2 have different porosities denoted as Φ_m1 and Φ_m2, and the total porosity of the total matrix denoted as Φ_m is the sum, presented as:

  $${\Phi }_m={\Phi }_{m1}+{\Phi }_{m2}$$

  (1)

Based on the assumptions mentioned above, the pore volume ratio χ is introduced to describe this characterization [24], which satisfies:

$$\chi =\frac{{\Phi }_{m1}}{{\Phi }_{m1}+{\Phi }_{m2}}$$

(2)

2.1.2. Conceptual Model 2

The matrix is artificially divided into two sections. Some further assumptions are made to derive the following formulas:

• Sub-matrix m1 and sub-matrix m2 are strictly separated in the vertical direction, and the whole matrix is divided into two layers.
• Both of them have the same width L_f, length L_F and porosity (see Figure 2b), which is given as:

  $${\Phi }_m={\Phi }_{m1}={\Phi }_{m2}$$

  (3)
• Sub-matrix m1 and sub-matrix m2 have different thicknesses denoted as H₁ and H₂, and the total thickness of the matrix is denoted as H, which is the sum of those two parts:

  $$H=H_1+H_2$$

  (4)

Based on the assumptions, the thickness ratio γ is introduced to describe this character, and the definition of this parameter is derived as:

$$\gamma =\frac{H_1}{H_1+H_2}$$

(5)

2.2. Mathematical Model

Based on the mass balance principle, the governing equations considering the adsorption, desorption and viscous flowing can be established respectively (see Appendix A).To simplify these equations, dimensionless variables are introduced (see in

Table 1. Definition of dimensionless variables.

Dimensionless Variables	Definitions
Dimensionless pseudo pressure	${\psi }_D=\frac{{\psi }_i-{\psi }_j}{{\psi }_i-{\psi }_{wf}}$, where $j=F,f_1,f_2,m1,m2$
Dimensionless time	$t_D=\frac{3.6k_{Fi}t}{\mu {\Lambda }_1{A}_{cw}}$, ${\Lambda }_1={(\phi c_t)}_{F+f_1+m1+m2}+\frac{Z{P}_{sc}T{V}_L{P}_L}{T_{sc}\overline{P}{(P_L+\overline{P})}^2}$
Dimensionless space	$y_D=\frac{y}{A_{cw}}$, $A_{cw}=2x_{e}h$, $x_D=\frac{x}{L_F/2}$
Dimensionless inter-porosity index	${\lambda }_{F-f}=\frac{12}{L_F^2}\frac{k_f}{k_F}{A}_{cw}$, ${\lambda }_{f-m1}=\frac{12}{L_F^2}\frac{k_{m1}}{k_F}{A}_{cw}$, ${\lambda }_{f-m2}=\frac{12}{L_F^2}\frac{k_{m2}}{k_F}{A}_{cw}$
Dimensionless storability ratio	${\omega }_j=\frac{{(\phi c_t)}_j}{{\Lambda }_1}$, where $j=F,f_1,f_2,m1,m2$
Dimensionless conductivity ratio	$\eta =\frac{k_{Fi}/A_{cw}}{k_{f1}/{(L_F/2)}^2}$
Dimensionless production data	$q_D=\frac{8.589\times {10}^{-6}{q}_{sc}T}{({\psi }_i-{\psi }_{wf})k_{Fi}\sqrt{A_{cw}}}$
Dimensionless permeability ratio	$\tau =k_{Fi}/k_{f2}$

), and then the final dimensionless governing equations are derived.

2.2.1. Conceptual Model 1

For Model 1, the corresponding equation governing fluid flowing is presented as follows:

• The pressure equation governing fluid flow in hydraulic fracture is given as:

  $$\frac{{\partial }^2{\psi }_{DF}}{\partial y_D^2}-\frac{{\lambda }_{F-f}}{3}{\left(\frac{\partial {\psi }_{Df}}{\partial x_D}\right)|}_{x_D=1}-\frac{{\lambda }_{F-m1}}{3}{\left(\frac{\partial {\psi }_{Df}}{\partial x_D^{*}}\right)|}_{x_D^{*}=1}={\omega }_F\frac{\partial {\psi }_{DF}}{\partial t_D}$$

  (6)
• The pressure equation governing fluid flow in natural fracture is given as:

  $$\frac{{\partial }^2{\psi }_{Df}}{\partial y_D^2}-\frac{{\lambda }_{F-m2}}{{\lambda }_{F-f}}{\left(\frac{\partial {\psi }_{Dm2}}{\partial y_D^{*}}\right)|}_{y_D^{*}=1}=\frac{3{\omega }_F}{{\lambda }_{F-f}}\frac{\partial {\psi }_{Df}}{\partial t_D}$$

  (7)
• The pressure equation governing fluid flow in matrix 1 is given as:

  $$\frac{{\partial }^2{\psi }_{Dm1}}{\partial x_D^{*}{}^2}=\frac{3}{{\lambda }_{F-m}}\chi \left({\omega }_m+{\omega }_d\right)\frac{\partial {\psi }_{Dm2}}{\partial t_D}$$

  (8)
• The pressure equation governing fluid flow in matrix 2 is given as:

  $$\frac{{\partial }^2{\psi }_{Dm2}}{\partial y_D^{*}{}^2}=\frac{3}{{\lambda }_{f-m2}}\left[\left(1-\chi \right){\omega }_m+{\omega }_d\right]\frac{\partial {\psi }_{Dm2}}{\partial t_D}$$

  (9)

2.2.2. Conceptual Model 2

In the upper layer unit, the mathematical model is described by the following differential equations.

• The pressure equation governing fluid flow in hydraulic fracture is given as:

  $$\frac{{\partial }^2{\psi }_{DF2}}{\partial y_D^2}-\frac{{\lambda }_{F-f}}{3}{\left(\frac{\partial {\psi }_{Df}}{\partial x_D}\right)|}_{x_D=1}={\omega }_F\frac{\partial {\psi }_{DF2}}{\partial t_D}$$

  (10)
• The pressure equation governing fluid flow in natural fracture is given as:

  $$\frac{{\partial }^2{\psi }_{Df2}}{\partial x_D^2}-\frac{{\lambda }_{F-m2}}{{\lambda }_{F-f}}{\left(\frac{\partial {\psi }_{Dm2}}{\partial y_D^{*}}\right)|}_{y_D^{*}=1}=\frac{3{\omega }_F}{{\lambda }_{F-f}}\frac{\partial {\psi }_{Df2}}{\partial t_D}$$

  (11)
• The pressure equation governing fluid flow in matrix 2 is given as:

  $$\frac{{\partial }^2{\psi }_{Dm2}}{\partial y_D^{*}{}^2}=\frac{3}{{\lambda }_{f-m2}}\left({\omega }_m+{\omega }_d\right)\frac{\partial {\psi }_{Dm2}}{\partial t_D}$$

  (12)
  In the lower layer unit, the mathematical model is described by the following differential equations.
• The pressure equation governing fluid flow in hydraulic fracture is given as:

  $$\frac{{\partial }^2{\psi }_{DF1}}{\partial y_D^2}-\frac{{\lambda }_{F-f}}{3}{\left(\frac{\partial {\psi }_{Df1}}{\partial x_D^{*}}\right)|}_{x_D^{*}=1}-\frac{{\lambda }_{F-m1}}{3}{\left(\frac{\partial {\psi }_{Dm1}}{\partial x_D^{*}}\right)|}_{x_D^{*}=1}={\omega }_F\frac{\partial {\psi }_{DF1}}{\partial t_D}$$

  (13)
• The pressure equation governing fluid flow in natural fracture is given as:

  $$\frac{{\partial }^2{\psi }_{Df1}}{\partial x_D^2}=\frac{3{\omega }_f}{{\lambda }_{F-f}}\frac{\partial {\psi }_{Df1}}{\partial t_D}$$

  (14)
• The pressure equation governing fluid flow in matrix 1 is given as:

  $$\frac{{\partial }^2{\psi }_{Dm1}}{\partial x_D^{*}{}^2}=\frac{3}{{\lambda }_{f-m1}}\left({\omega }_m+{\omega }_d\right)\frac{\partial {\psi }_{Dm1}}{\partial t_D}$$

  (15)

In addition, the corresponding conditions are list in

Table 2. Qualitative description of initial and boundary conditions.

Condition		Model 1	Model 2
Initial		${\psi }_{DF}(y_D,0)={\psi }_{DF1}(y_D,0)={\psi }_{DF2}(y_D,0)=0$
		${\psi }_{Df}(x_D,0)={\psi }_{Df1}(x_D,0)={\psi }_{Df2}(x_D,0)=0$
		${\psi }_{Dm1}(x_D^{*},0)={\psi }_{Dm2}(y_D^{*},0)=0$
Inner boundary	HF	${\frac{\partial {\psi }_{DF}}{\partial y_D}|}_{y_D=0}=-2\mathsf{\pi }$	${\frac{\partial {\psi }_{DF1}}{\partial y_D}|}_{y_D=0}=-\frac{2\mathsf{\pi }}{\gamma }$ ${\frac{\partial {\psi }_{DF2}}{\partial y_D}|}_{y_D=0}=-\frac{2\mathsf{\pi }}{1-\gamma }$
	NF	${\frac{\partial {\psi }_{Df}}{\partial x_D}|}_{x_D=0}=0$, ${{\psi }_{Df}|}_{x_D=1}={\psi }_{DF}$	${\frac{\partial {\psi }_{Df1}}{\partial x_D}|}_{x_D=0}=0$, ${{\psi }_{Df1}|}_{x_D=1}={\psi }_{DF1}$ ${\frac{\partial {\psi }_{Df2}}{\partial x_D}|}_{x_D=0}=0$, ${{\psi }_{Df2}|}_{x_D=1}={\psi }_{DF1}$
	Matrix	${\frac{\partial {\psi }_{Dm1}}{\partial x_D^{*}}|}_{x_D^{*}=0}=0$, ${{\psi }_{Dm1}|}_{x_D^{*}=1}={\psi }_{DF}$ ${\frac{\partial {\psi }_{Dm2}}{\partial y_D^{*}}|}_{y_D^{*}=0}=0$, ${{\psi }_{Dm2}|}_{y_D^{*}=1}={\psi }_{Df}$	${\frac{\partial {\psi }_{Dm1}}{\partial x_D^{*}}|}_{x_D^{*}=0}=0$, ${{\psi }_{Dm1}|}_{x_D^{*}=1}={\psi }_{DF1}$ ${\frac{\partial {\psi }_{Dm2}}{\partial y_D^{*}}|}_{y_D^{*}=0}=0$, ${{\psi }_{Dm2}|}_{y_D^{*}=1}={\psi }_{Df2}$
Outer boundary		${\frac{\partial {\psi }_{DF}}{\partial y_D^{}}|}_{y_D^{}=y_e/\sqrt{A_{cw}}}=0$, ${\frac{\partial {\psi }_{DF1}}{\partial y_D^{}}|}_{y_D^{}=y_e/\sqrt{A_{cw}}}=0$, ${\frac{\partial {\psi }_{DF2}}{\partial y_D^{}}|}_{y_D^{}=y_e/\sqrt{A_{cw}}}=0$

.

2.3. Model Solution

The Laplace transformation is used to solve Equations (6)–(15), the details of derivation are presented in Appendix B. As shown from Equations (A37) and (A41), the final solutions of these equations in Laplace space are derived as:

For Model 1, the dimensionless production rate is given as

$${\tilde{q}}_{Dsc}=\frac{\sqrt{f_F(s)}}{2\mathsf{\pi }s}\left[\frac{1-\exp \left(-2\sqrt{f_F(s)}{y}_{De}\right)}{1+\exp \left(-2\sqrt{f_F(s)}{y}_{De}\right)}\right]$$

(16)

For Model 2, the dimensionless production rate is given as

$${\tilde{q}}_{Dsc}=\frac{\gamma \sqrt{f_F(s)}}{2\mathsf{\pi }s}\left[\frac{1-\exp \left(-2\sqrt{f_{F1}(s)}{y}_{De}\right)}{1+\exp \left(-2\sqrt{f_{F1}(s)}{y}_{De}\right)}+\frac{(1-\gamma )\sqrt{f_{F2}(s)}}{2\mathsf{\pi }s}\frac{1-\exp \left(-2\sqrt{f_{F2}(s)}{y}_{De}\right)}{1+\exp \left(-2\sqrt{f_{F2}(s)}{y}_{De}\right)}\right]$$

(17)

However, to analyze the impact of relative parameters and identify flow regimes, the rate transient responses must be inverted into real time space by using the Stehfest numerical inversion algorithm [25].

3. Results and Discussion

In this model, some characterized parameters includes: ω_F, ω_f, ω_m, χ, γ.

3.1. Model Validation

3.1.1. Validation with Analytical Approach

When these parameters are equal to a special value, the new model proposed in this paper can be varied to alternative models which have been verified reasonable. For example, when ω_f is equates to 0, the model can be simplified as dual porosity slab model [26]; when χ = 0 or γ = 0, it indicates that the gas only depletes from matrix to NF [21]. Put another way, the new model can be equivalent to the traditional triple-linear model [20]. In summary, this new model can be universally applied to different formations only if some characterized parameters are equal to the corresponding values. Figure 3 displays the comparison of our model with alternative approach in the case of linear flow. Our model has a significant agreement with El-Banbi model [20] during the whole production cycle.

3.1.2. Validation with Numerical Approach

Furthermore, a numerical simulation is conducted to compare with the models (the input data are listed in

Table 3. Input parameters used in the case of constant flowing pressure.

Parameters	Symbol	Unit	Value
Initial pressure	Pi	MPa	20
Downhole pressure	Pwf	MPa	15
Formation temperature	T	K	333
Horizontal length	L	m	1200
Formation thickness	H	m	100
Macro-fracture length	yf	m	200
Total compressibility of HF	CtF	MPa−1	5 × 10−4
Total compressibility of NF	CtF	MPa−1	5 × 10−4
Total compressibility of Matrix	CtF	MPa−1	5 × 10−4
Porosity of HF	ΦF	Dimensionless	0.0005
Porosity of NF	Φf	Dimensionless	0.005
Total porosity of matrix	Φm	Dimensionless	0.08
Permeability of HF	kF	D	2
Permeability of natural fracture	kf	D	10−6
Permeability of sub-matrix m1	km1	D	10−8
Permeability of sub-matrix m2	km2	D	10−8
Langmuir pressure	PL	MPa	5
Langmuir volume	VL	sm3/m3	5

). The commercial software Eclipse 2010 (Eclipse Foundation, Ottawa, ON, Canada) is selected to simulate the behavior of triple-porosity reservoirs and verify the semi-analytical solutions in this paper.

The numerical model considers the flow influx towards the horizontal wellbore. Considering the symmetrical structure of the multi-fractured horizontal well, one representative segment is selected to represent one quadrant of the reservoir volume around a hydraulic fracture. This segment contains twenty micro fractures orthogonal to the hydraulic fractures at 10 m fracture spacing. As shown in Figure 4a, the model is built to be a two dimension model with 100 grid cells in the x-direction, 100 grid cells in y-direction and only one grid cell in the z-direction. The multi-porosity method is applied to further divide the matrix, thus the transient flow in the matrix can be simulated. Additionally, the desorption model is instant. The local grid refinement (LGR) method is employed to create the hydraulic fractures. Only then is the multi-stage hydraulic fracturing simulated. From the comparison results of transient rate in Figure 4b, the semi-analytical results in this paper are consistent with the results from numerical model.

Based on the model verification, it is proved that our model is accurate and can be further used to calculate more complex case.

3.2. Sensitivity Analysis

In this subsection, the parameters used are put into Table 3. The transient rate behavior is achieved by using Stehfest algorithm [25]. Figure 5 demonstrates the type curves of transient rate behavior for new model, where these curves are respectively shown in log-log and regular plots.

Seen from the type curve in Figure 5b,d, some well-known flowing regimes can be identified, such as the bilinear flow (the slope of curve is identify by −0.25), the linear flow (the slope of curve is identified by −0.5) and the boundary dominated flow (the slope of curve is approximately identified by −0.25).

Due to the limitation of χ and γ, it is possible that some flow regimes could not be identified. Under the assumption that (1) χ = 0 and γ = 0 or (2) χ = 1 and γ = 1, the two models are equivalent, which indicates that the characterization of flow regimes caused by two models is identical.

For model 1 seen in Figure 5a,b, with the increase of the pore volume ratio χ, the gas rate firstly increases and then decreases. When matrix and NF simultaneously deplete into HF, desorption occurs earlier and the rate is higher. However, when gas depletes from the matrix to HF, the boundary dominated flow occurs earlier, so the gas rate become lower. In term of model 2 seen in Figure 5c,d, with the increase of the thickness ratio γ, gas rate decreases. Model 2 assumes that there is no connection between the upper and lower layers. Meanwhile, the matrix and NF have no connection in the lower layer, so the impact of boundary is more severe and the gas rate decreases.

It is noted that this phenomenon mentioned above is not the most interesting part in this paper. By comparing Figure 5a,c, it is easily found that the gas rate is distinct at the early production period (less than 2 years), dependent on the pore volume ratio χ and the thickness ratio γ. During late production period, the gas rate is approximately identical under different values of χ and γ. It is a significant discovery in the field of reservoir engineering when conducting production prediction and estimating recovery in the case that the available data are limited, especially for shale gas wells. If pore volume ratio χ or thickness ratio γ are presented in the reasonable range, model 1 or 2 can perfectly match with the actual production data at the early stage, and then the future performance of gas well can be predicted based on the matched model. It is the reason that the intermediate- and late-time production performance is almost independent on the values of pore volume ratio χ or thickness ratio γ.

4. Model Aggregation

The numerical and semi-analytical results have been presented in Figure 4b. It is emphasized that the semi-analytical transient rate behavior presented by model 1 is higher than model 2. However, the rate predicted from the numerical simulation is in the range of two models. It indicates that the simulated curve of model 1 overlies model 2 during the whole production period when the value of pore volume ratio χ or thickness ratio γ is reasonable. Therefore, it makes sense to integrate the two models together to develop a generalized framework model for eliminating the error when predicting future performance of gas well. Here, an approach of Generalized Likelihood Uncertainty Estimation (GLUE) methodology is employed [27]. The approach of GLUE is presented based on the probabilistic aggregation theory. Using this approach, these two models can be integrated according to the personal weight which reduces the uncertainties.

4.1. GLUE Method

The theory of the GLUE method will be discussed in this section. When the pore volume ratio χ and thickness ratio γ can be reasonably chosen, the resulting simulation based on new model could provide better matching of limited history data. However, there exists remarkable error on rate prediction.

This study proposed an approach to integrate these two models, which includes two steps. (1) Two models are used to match the data to obtain reasonable values of the pore volume ratio χ and the thickness ratio γ; (2) A weight is assigned to each model based on some goodness-of-fit statistics, and then weighted mean and standard deviation of desired performance measures (the rate and cumulative rate) could be calculated. Here, a simple likelihood measure in the GLUE method can be defined using the root-mean-square-error (RMSE):

$$L_j={(1/RMS{E}_j)}^2$$

(18)

where, RMSE_j represents the j-th model, j = 1, 2 in this paper. Normalizing the likelihoods, the sum is equal to one, which provides the GLUE weight for the two models:

$$w_j=\frac{P{R}_j\cdot L_j}{{\displaystyle \sum _{j=1}^2\left(P{R}_j\cdot L_j\right)}}=\frac{P{R}_j/RMS{E}_j^2}{{\displaystyle \sum _{j=1}^2\left(P{R}_j/RMS{E}_j^2\right)}}$$

(19)

where, L_j is the likelihood functions, PR_j is the prior weight of each model (i.e., PR_j = 1/2 if all are assumed to be equally likely prior weighting). From these weights, the integrated gas rate q and cumulative rate G on certain production time T₀ can be calculated as:

$$q_{t=T_0}={\displaystyle \sum _{j=1}^2{w}_j{(q_{t=T_0})|}_j}$$

(20)

$$G_{t=T_0}={\displaystyle \sum _{j=1}^2{w}_j{(G_{t=T_0})|}_j}$$

(21)

and the standard deviation can be computed from:

$$SD[q_{t=T_0}]=\sqrt{{\displaystyle \sum _{j=1}^2{w}_j{[q_{t=T_0}(\mathrm{avg})-{(q_{t=T_0})|}_j]}^2}}$$

(22)

$$SD[G_{t=T_0}]=\sqrt{{\displaystyle \sum _{j=1}^2{w}_j{[G_{t=T_0}(\mathrm{avg})-{(G_{t=T_0})|}_j]}^2}}$$

(23)

4.2. Application Illustration

4.2.1. Procedure

Based on above analysis, the procedure would be divided into three steps:

• Actual production data are matched based on model 1 and model 2 respectively. As a result, the most reliable value of pore volume ratio χ and thickness ratio γ can be obtained;
• Based on history and predicted curves, the likelihoods can be calculated respectively for model 1 and model 2 using Equation (18), and then the likelihoods can be normalized using Equation (19);
• Based on Equation (20), we can calculate the production rate and cumulative rate at a certain production time T₀, and then compare those results with the numerical simulation by calculating standard deviation using Equation (21).

4.2.2. Field Example

A multi-stage fractured horizontal well is chosen from certain shale gas zone in Sichuan Basin, China. The geometry and distribution of hydraulic fractures can be diagnosed with seismic events during hydraulic-fracturing treatment. Fracture locations are interpreted by using the event amplitude as an intensity measure, and fracture sizes are estimated by the relationship between the amplitude and the seismic moment. It can be seen from Figure 6, the direction of fracture is N570E for this well. The basic data from the well testing and the core analysis are listed in

Table 4. Input parameters used in the field example.

Parameters	Symbol	Unit	Value
Initial pressure	Pi	MPa	22.5
Downhole pressure	Pwf	MPa	15
Formation temperature	T	K	314
Horizontal length	L	m	1450
Formation thickness	H	m	48.3
Macro-fracture length	yf	m	214
Total compressibility of HF	CtF	MPa−1	5.2 × 10−4
Total compressibility of NF	CtF	MPa−1	5.2 × 10−4
Total compressibility of Matrix	CtF	MPa−1	5.2 × 10−4
Porosity of HF	ΦF	Dimensionless	4.6 × 10−4
Porosity of NF	Φf	Dimensionless	5.3 × 10−3
Total porosity of matrix	Φm	Dimensionless	0.078
Permeability of HF	kF	D	1.8
Permeability of natural fracture	kf	D	2.1 × 10−6
Permeability of sub-matrix m1	km1	D	1.2 × 10−8
Permeability of sub-matrix m2	km2	D	1.2 × 10−8
Langmuir pressure	PL	MPa	5.2
Langmuir volume	VL	sm3/m3	5.2

.

In Figure 7, the red pots indicates the actual production data (production period of this well is less than 300 days), and colorful solid lines represent the simulated results from two semi-analytical models. It is not difficult to observe that (1) when χ = 0.55, the predictive (cumulative) rate can effectively match with the history data from model 1; and (2) when γ = 0.15, the (cumulative) rate can match with history data from model 2.

Thus, the likelihood for those two models can be calculated based on the difference between the history data and predictive data (see Figure 7), where L₁ = 23.06, L₂ = 19.88. Furthermore, based on Equation (19), the weights for those two models are as follows: w₁ = 0.54, w₂ = 0.46. Finally, the aggregation model can be presented based on Equation (20). Although the history data (less than 300 days) is limited from the viewpoint of establishing the integrated model, this new model has the capacity to predict the rate and cumulative rate for this well (see Figure 8).

To validate this proposed integrated model, a numerical model is developed to simulate the performance within 1000 days. At the same time, the new model simultaneously predicts the rate and cumulative rate within 1000 production days. We can apparently find that the rate can perfectly match with each other, while the cumulative rate predicted from the aggregation model can approximately match with the results calculated from the numerical model. In addition, those two models can provide an approximate range for the cumulative rate. It has the potential of decreasing the uncertainty for the prediction of estimation ultimate recovery, especially for unconventional reservoirs where the transient flow regime lasts for extremely long time. Based on Equation (21), the standard deviations of the rate and cumulative rate are obtained as follows: SD[q_t = 1000] = 0.012 × 10⁴ m³/day, SD[G_t = 1000] = 12.4 × 10⁴ m³. Therefore, this new aggregation model is a practical tool in terms of predicting shale gas and tight gas rate.

5. Conclusions

This paper presents two new simple models for shale gas reservoirs with multi-stage fractured horizontal well, where the USRV zone is regarded as a dual-porosity & dual-depletion system. The resulting solutions are more universal for conducting type curve analysis in homogeneous and naturally fractured reservoirs. Numerical simulation model is correspondingly developed to validate the semi analytical solutions.

Two characteristic parameters, the pore volume ratio and thickness ratio are defined to characterize the dual-depletion mechanism. A set of type curves is generated under different cases of those two parameters. Several regular regimes can be identified, such as bilinear flow, linear flow, and outer boundary-dominated flow, etc.

The Generalized Likelihood Uncertainty Estimation methodology is introduced to aggregate those two models to reduce the uncertainties on future performance forecast. The numerical simulation and field application further validate the new simplified aggregation model. This paper also demonstrates that the proposed uncertainty assessment protocol using statistical model averaging concepts (GLUE) is a robust method to quantify the uncertainty of production prediction and cumulative prediction.