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

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.


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.

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.

Conceptual Model 1
The matrix is divided into two sections.Some further assumptions are made to derive the following formulas: 1.
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);

2.
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: Based on the assumptions mentioned above, the pore volume ratio χ is introduced to describe this characterization [24], which satisfies: The matrix is artificially divided into two sections.Some further assumptions are made to derive the following formulas: 1.
Sub-matrix m1 and sub-matrix m2 are strictly separated in the vertical direction, and the whole matrix is divided into two layers.2.
Both of them have the same width L f , length L F and porosity (see Figure 2b), which is given as: Sub-matrix m1 and sub-matrix m2 have different thicknesses denoted as H 1 and H 2 , and the total thickness of the matrix is denoted as H, which is the sum of those two parts: Based on the assumptions, the thickness ratio γ is introduced to describe this character, and the definition of this parameter is derived as:

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), and then the final dimensionless governing equations are derived.

Dimensionless Variables Definitions
Dimensionless pseudo pressure Dimensionless time Dimensionless storability ratio Dimensionless production data For Model 1, the corresponding equation governing fluid flowing is presented as follows: 1.
The pressure equation governing fluid flow in hydraulic fracture is given as: The pressure equation governing fluid flow in natural fracture is given as: 3.
The pressure equation governing fluid flow in matrix 1 is given as: The pressure equation governing fluid flow in matrix 2 is given as:

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

1.
The pressure equation governing fluid flow in hydraulic fracture is given as: The pressure equation governing fluid flow in natural fracture is given as: 3.
The pressure equation governing fluid flow in matrix 2 is given as: In the lower layer unit, the mathematical model is described by the following differential equations.

4.
The pressure equation governing fluid flow in hydraulic fracture is given as: The pressure equation governing fluid flow in natural fracture is given as: The pressure equation governing fluid flow in matrix 1 is given as: In addition, the corresponding conditions are list in Table 2.
Table 2. Qualitative description of initial and boundary conditions.

Condition Model 1 Model 2
Initial Outer boundary

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 For Model 2, the dimensionless production rate is given as 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].

Results and Discussion
In this model, some characterized parameters includes: ω F , ω f , ω m , χ, γ.

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.Comparison of our model with alternative analytical approach in the dimensionless coordinate.

Validation with Numerical Approach
Furthermore, a numerical simulation is conducted to compare with the models (the input data are listed in Table 3).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.

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 γ.

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.

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): 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: 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 0 can be calculated as: and the standard deviation can be computed from:

Procedure
Based on above analysis, the procedure would be divided into three steps: 1.
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; 2.
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); 3.
Based on Equation ( 20), we can calculate the production rate and cumulative rate at a certain production time T 0 , and then compare those results with the numerical simulation by calculating standard deviation using Equation ( 21).

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.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 1 = 23.06,L 2 = 19.88.Furthermore, based on Equation (19), the weights for those two models are as follows: w 1 = 0.54, w 2 = 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 4 m 3 /day, SD[G t = 1000] = 12.4 × 10 4 m 3 .Therefore, this new aggregation model is a practical tool in terms of predicting shale gas and tight gas rate.

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.

1.
For hydraulic fracture: diffusivity equation that control HF-horizontal well communication is described as For Natural fracture: diffusivity equation that control NF-HF communication is described as 3.
For sub-matrix 1: diffusivity equation that control matrix-HF communication is described as For sub-matrix 2: diffusivity equation that control matrix-NF communication is described as In Upper Layer, the pressure governing equations in hydraulic fracture (HF), natural fracture (NF) and sub matrix 2 are given as: In Lower Layer, the pressure governing equations in hydraulic fracture (HF), natural fracture (NF) and sub matrix 1 are given as:

. Initial and Boundary Condition
In addition, initial conditions in HF, NF and matrix are given as:

•
In model 1, inner boundary conditions are given as:

•
In model 2, inner boundary condition in HF, NF and matrix are given as:  In addition, the inner and outer boundary conditions are satisfied as: In the unit of lower layer, the pressure governing equations in hydraulic fracture, natural fracture and matrix 1 are given as: = ω F s ψ DF1 , in hydraulic fracture (A24) In terms of conceptual model 2, the total production rate is the sum of upper layer and upper layer, therefore, the final solution of production rate combing Equations (A36) and (A39):

Figure 1 .
Figure 1.(a) The illustration of multi-fractured horizontal well; (b) The two segments of matrix in stimulated reservoir volume (SRV) zone.

Figure 3 .
Figure 3.Comparison of our model with alternative analytical approach in the dimensionless coordinate.

Figure 4 .
Figure 4. (a) Depleted shale gas reservoir model established using Eclipse 2010; (b) Comparison between numerical model and those two new models proposed in this paper.

Figure 5 .
Figure 5. Type curve of flow regime: (a) Cartesian plot for model 1; (b) log-log plot for model 1; (c) Cartesian plot for model 2; (d) log-log plot for model 2.

Figure 6 .
Figure 6.Natural seismic event for this well, where the spots with different color stand for the individual stages.

Figure 8 .
Figure 8.Comparison between aggregation model and numerical simulation: (a) gas rate; (b) cumulative rate.

Table 1 .
Definition of dimensionless variables.

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

Table 4 .
Input parameters used in the field example.