Development and Calibration of a Semianalytic Model for Shale Wells with Nonuniform Distribution of Induced Fractures Based on ES-MDA Method

Abstract

Given reliable parameters, a newly developed semianalytic model could offer an efficient option to predict the performance of the multi-fractured horizontal wells (MFHWs) in unconventional gas reservoirs. However, two major challenges come from the accurate description and significant parameters uncertainty of stimulated reservoir volume (SRV). The objective of this work is to develop and calibrate a semianalytic model using the ensemble smoother with multiple data assimilation (ES-MDA) method for the uncertainty reduction in the description and forecasting of MFHWs with nonuniform distribution of induced fractures. The fractal dimensions of induced-fracture spacing (d_fs) and aperture (d_fa) and tortuosity index of induced-fracture system (θ) are included based on fractal theory to describe the properties of SRV region. Additionally, for shale gas reservoirs, gas transport mechanisms, e.g., viscous flow with slippage, Knudsen diffusion, and surface diffusion, among multi-media including porous kerogen, inorganic matter, and fracture system are taken into account and the model is verified. Then, the effects of the fractal dimensions and tortuosity index of induced fractures on MFHWs performances are analyzed. What follows is employing the ES-MDA method with the presented model to reduce uncertainty in the forecasting of gas production rate for MFHWs in unconventional gas reservoirs using a synthetic case for the tight gas reservoir and a real field case for the shale gas reservoir. The results show that when the fractal dimensions of induced-fracture spacing and aperture is smaller than 2.0 or the tortuosity index of induced-fracture system is larger than 0, the permeability of induced-fracture system decreases with the increase of the distance from hydraulic fractures (HFs) in SRV region. The large d_fs or small θ causes the small average permeability of the induced-fracture system, which results in large dimensionless pseudo-pressure and small dimensionless production rate. The matching results indicate that the proposed method could enrich the application of the semianalytic model in the practical field.

1. Introduction

Unconventional reservoirs, especially the tight sand and shale reservoirs, have become more and more significant with the development of the horizontal wells and hydraulic fracturing technologies. The complex fracture network surrounding the multistage fractured horizontal wells (MFHWs) called stimulated reservoir volume (SRV), has obvious influences on the production performance [1,2,3]. Currently, a lot of effort has been done to model the pressure transient and rate transient of MFHWs aiming to forecast the production accurately. Through numerical methods, analytical methods, and empirical methods [4,5,6,7], production data analysis in unconventional reservoirs, such as Barnett, Bakken, and Eagle, shows that the linear flow in formation, especially in hydraulic fractures (HFs) and SRV, dominates the production of MFHWs. Meanwhile, the advantage of analytic models is their simplicity and less parameters compared with numerical and empirical models, although they cannot characterize the over-complex heterogeneity properties and fracture system. Thus, considering the main flow regimes in the life of MFHWs in unconventional gas reservoirs, various linear flow models based on linear-flow assumptions have been developed. Ozkan et al. [8] and Brown et al. [9] reported a trilinear model to study the performance of MFHWs in unconventional reservoirs, which divides the formation into three main regions. In order to describe the SRV size between HFs more reasonably, Stalgorova and Matter [10] improved the trilinear model to a five-linear model by simplifying SRV in a region with limited width. On this basis, many scholars studied the production rate or pressure transient performances of MFHWs in unconventional reservoirs [11,12]. Although the models mentioned above can simulate the main linear-flow regimes in both SRV and unstimulated reservoir volume (USRV), the heterogeneous distribution of complex induced fractures are not considered in SRV region, which is the major distinct properties between SRV and USRV as shown in Figure 1.

The distribution of induced fractures system in SRV affected by HFs and the pre-existing closed natural fractures (NFs) with self-similar characterization can be described by fractal theory [13,14,15]. Chang and Yortsos [16] quantified the heterogeneous permeability and porosity of NFs using fractal theory. Then, Cossio et al. [17] introduced the fractal permeability/porosity relation [18] into the fractal diffusivity equation (FDE), which proved to be more fundamental than the classical diffusivity equation. Wang et al. [19] combined FDE with dual-porosity media model in the trilinear model to describe oil flow in SRV region of MFHWs in tight reservoirs, and their model provided a suitable way to study the flow behaviors in heterogeneous induced fractures systems. Sheng et al. [20] extended the fractal trilinear model to shale gas reservoirs coupling multi-porosity media (porous kerogen, inorganic matter, and induced-fracture) with multiple gas transport mechanisms (viscous flow with slippage, Knudsen diffusion, and surface diffusion). However, they did not consider the unstimulated reservoir volume (USRV) between HFs, which cannot accurately describe the microseismic results. In addition, the fracture porosity and permeability distribution affected by the distribution of induced-fracture spacing and aperture, but there are limited reports about the calculation methods for the non-uniform distribution of induced fractures spacing and aperture using fractal theory [21].

History matching to estimate the critical parameters of formation and hydraulic fracturing and rate prediction of MFHWs plays important roles in unconventional reservoirs based on the numerical and analytic models. Samandarli et al. [22] presented a semianalytic method to estimate reservoir parameters by history matching the production data of hydraulically fractured shale gas wells using a least absolute value regression. Correction for gas properties changing with time and desorbed gas at low pressure was considered in the model. According to the algebraic equations, bilinear flow regime and linear flow regime were used to match both synthetic data and field data by doing regression on matrix permeability, fracture permeability, and length. Hategan and Hawkes [23] used a simple analytical model based on the solution of the pseudo-steady state equation to analyze the production of MFHWs in shale gas reservoirs, and thought most shale gas reservoirs fit different type curves by normalizing to reservoir pressure, permeability, and number of fracture stages. Considering the effect of pressure-dependent natural-fracture permeability on production from shale gas wells, Cho et al. [24] modified the trilinear flow model to match the production data of Haynesville and Barnett shale gas wells. The reservoir properties and the correlation coefficients used for the pressure-dependent permeability were matched ignoring the non-uniqueness issues caused by the large number of parameters. Ogunyomi et al. [25] developed an approximate analytical solution to the double-porosity model for MFHWs in unconventional oil reservoirs. Due to the solution obtained in real-time space, the model was used to match the field data with the least-square method mainly by changing the parameters related to time, and forecast the production rates. Clarkson et al. [26] combined analytical, semi-analytical, and empirical methods for history matching and production forecasting for MFHWs in tight and shale gas reservoirs. Linear-to-boundary (LTB) model [27], composite model [9,10], and semi-analytical model (SAM) [26] were used to match the field data by fitting the parameters of SRV and unstimulated region. Chen et al. [28] proposed a comprehensive semi-analytical model for MFHWs considering the shale gas flow mechanisms. The unknowns including the well length, fracture number, conductivity, and length were determined by history matching. Yao et al. [7] developed an analytical multi-linear model based on linear flow and radial flow assumption. Permeability and length of different regions, and fracture length were estimated by history matching with several field data. Therefore, the semianalytic models coupling with history matching methods can be applied for not only obtaining some critical parameters for MFHWs but also predicting the production rate.

The objective of this paper is development and calibration of a semianalytic model for shale wells with nonuniform distribution of induced fractures based on ensemble smoother with multiple data assimilation (ES-MDA). The fractal dimensions of induced-fracture spacing and aperture and tortuosity index of induced-fracture system are included based on fractal theory to describe the properties of SRV region. For shale gas reservoirs, gas transport mechanisms (viscous flow with slippage, Knudsen diffusion, and surface diffusion) among multi-media (porous kerogen, inorganic matter, and fracture system) are taken into account and the model is verified. The effects of the fractal dimensions and tortuosity index of induced fractures on MFHWs performances are analyzed. Then, we employ the ES-MDA method with the presented model to reduce uncertainty in the forecasting of gas production rate for MFHWs in unconventional gas reservoirs using a synthetic case for the tight gas reservoir and a real filed case for the shale gas reservoir. The matching results indicate that the proposed approach could enrich the application of the semianalytic model in the practical field. Finally, some conclusions are drawn and presented based on the results and analysis.

2. Fractal Semianalytic Model for Mfhws in Unconventional Gas Reservoirs

As shown in Figure 2, the basic workflow of a new semianalytic model construction in this study consisted of two main parts: firstly, for SRV, which dominates the productivity of MFHWs in unconventional gas reservoirs, the fractal fracture spacing distribution (FFSD) was taken into account to describe its heterogeneous properties instead of classical idealized dual-media model; secondly, based on the linear flow assumption, combining with the linear models for unstimulated reservoir volume (USRV) (original formation), a coupled fractal multi-linear flow model (FMFM) for MFHWs was developed. The detailed derivations of the FMFM for MFHWs in unconventional gas reservoirs are presented in following sections.

2.1. Fractal Fracture Network Permeability and Porosity in SRV

Hydraulic fracturing creates the fracture network surround the horizontal well-bore and hydraulic fracture. Due to differences in situ stress and fracture perforation position, the distribution of induced-fracture properties including induced-fracture spacing and aperture are heterogeneous [21], which affects the distribution of permeability and porosity of fracture network directly. Therefore, the fractal dimensions for both induced-fracture spacing and aperture (d_fs and d_fa) are first discussed in this section. Then the fracture network permeability and porosity are proposed (details are shown in Appendix A). As shown in Figure 2c, using the fractal theory; thus, the fracture permeability can be obtained by

$$k_f\left(y\right)=k_{i-f}\left(y\right)n\left(y\right)\frac{A_p}{A_r}\frac{\mathrm{d}x}{\mathrm{d}L}=k_w{\left(\frac{y}{y_w}\right)}^{3\left(d_{fa}-2\right)+\left(d_{fs}-2\right)-\theta },$$

(1)

where k_f is the fracture permeability, m²; μ is the fluid viscosity, Pa·s; A_r is the sectional-area of a REV, m²; A_p is the sectional-area of induced-fracture, m²; dL is the flow length in the induced-fracture, m; k_w is the fracture permeability at the reference point, m²; y_w is the position of reference point, m; n is the number of fracture sites per unit thickness, m⁻¹; d_fs is the fractal dimension of induced-fracture spacing; d_fa is the fractal dimension of induced-fracture aperture; and θ is the tortuosity index of induced fracture [30].

Moreover, the fracture porosity is given by [21]

$${\varphi }_f\left(y\right)={\varphi }_{i-f}\frac{b_f\left(y\right)}{b_f\left(y\right)+s_f\left(y\right)}={\varphi }_w{\left(\frac{y}{y_w}\right)}^{d_{fs}+d_{fa}-4},$$

(2)

where ${\varphi }_f$ is the fracture porosity; ${\varphi }_{i-f}$ is the single induced-fracture porosity; ${\varphi }_w$ is the fracture porosity at the reference point; s_f is induced-fracture spacing, m; and b_f is induced-fracture aperture, m.

Then, Equations (1) and (2) are the basic forms of fracture permeability and porosity considering the fractal distribution in SRV of MFHWs in unconventional reservoirs. When y_w = 0.1 m and k_w = 10 × 10⁻¹⁵ m², the fracture permeability distribution with different fractal dimensions of induced-fracture spacing and aperture (d_fs and d_fa) and tortuosity index of induced-fracture (θ) are shown in Figure 3. The fracture permeability k_f decreases if d_fs < 2 and increases if d_fs > 2 with the increase of the distance from the reference point y_w as shown in Figure 3a. Meanwhile, the fracture permeability k_f decreases if d_fa < 2 and increases if d_fa > 2 with the increase of the distance from the reference point y_w as shown in Figure 3b. Figure 3c shows that the tortuosity index of induced-fracture θ decreases when θ < 0 and increases when θ > 0 with the increase of the distance from the reference point y_w [31]. The reasons can be illustrated as follows: If d_fs < 2 or d_fa < 2, as the distance from the reference point increases, the induced-fracture number or aperture decreases and the induced-fracture spacing increases; if d_fs = 2 or d_fa = 2, as the distance from the reference point increases, the induced-fracture number or aperture and the induced-fracture spacing are constants; if d_fs > 2 or d_fa > 2, as the distance from the reference point increases, the induced-fracture number or aperture increases and the induced-fracture spacing decreases; if θ > 0, the tortuosity of pores is more complex than that when θ < 0. Therefore, for SRV, the fractal dimensions of nature fracture are usually reported from 1.3 to 1.7 [32].

2.2. Transient Diffusivity Equations in SRV And USRV

Based on the trilinear model, a newly coupled FMFM was developed for MFHWs in unconventional gas reservoirs as shown in Figure 2c, which followed the basic linear flow assumptions and some simplifying assumptions: (1) the reservoir formation consisted of five regions including the hydraulic fracture, the effective SRV and three USRV regions; (2) the fractal triple-media model (FTMM) was applied to describe the heterogeneous characteristics of induced fractures in SRV, which consisted of three media (fracture network, kerogen, and inorganic matter); (3) idealized dual-media model (IDMM) was used to characterize the homogeneous properties of formation in USRV, which included porous kerogen and inorganic matter; (4) porous kerogen and inorganic matter are uniform and continuous media in both SRV and USRV regions; (5) for shale gas reservoirs, viscous flow, Knudsen diffusion, surface diffusion, and slip factor were taken into account; (6) single-phase fluid flowed into the horizontal wellbore from reservoir only by the hydraulic fractures, and frictional loss within the horizontal well and effects of gravity and capillary forces were ignored in the reservoir; (7) the hydraulic fractures had identical properties, and were evenly distributed along the horizontal well; and (8) there were the no-flow boundaries parallel to the hydraulic fractures at the middle of two fractures in closed reservoir.

2.2.1. Diffusivity Equations of FTMM in SRV (Region II)

Since the induced-fracture system, porous kerogen and inorganic matter coexist in SRV, there are two types of cross-flow flux happening: one is from porous kerogen to inorganic matter, and another one is from inorganic matter to induced-fracture system, both of which can be described by Warren–Root pseudo-steady model. Thus, combining Equations (1) and (2), the diffusivity equation of fluid flow in the induced-fracture system can be expressed as

$$\begin{array}{ll}\frac{\partial \left[{\varphi }_w{\left(y/y_w\right)}^{d_{fs}+d_{fa}-4}\left(p_{f}M/ZRT\right)\right]}{\partial t}& =\frac{\partial }{\partial y}\left[\frac{p_{f}M}{ZRT{\mu }_g}{k}_w{\left(\frac{y}{y_w}\right)}^{3\left(d_{fa}-2\right)+\left(d_{fs}-2\right)-\theta }\frac{\partial p_f}{\partial y}\right]\\ & +\frac{1}{x_f}\frac{p_{f}M}{ZRT{\mu }_g}{k}_m{\frac{\partial p_{3m}}{\partial x}|}_{x=x_f}+q_{m-f}^{\prime }\end{array},$$

(3)

where p_f is the induced-fracture pressure in SRV, Pa; Z is the gas compressibility factor; R is the universal gas constant, 8.314 J/(K·mol); μ_g is the gas viscosity, Pa·s; k_m is the apparent permeability of the inorganic matter, m²; $p_m$ is the inorganic matter pressure, Pa; $x_f$ is the hydraulic fracture half-length, m; subscript, 3 represents region III; and $q_{m-f}^{\prime }$ is the mass transfer from inorganic matter to the induced-fracture system, kg/m³·s.

For porous kerogen, if the Knudsen diffusion, surface diffusion, viscous flow and slippage effect are included, the apparent permeability can be given by (Sheng et al., 2019b)

$$k_k=\frac{{\varphi }_k}{{\tau }_k}\frac{2r_k}{3}\sqrt{\frac{8ZRT}{\pi M}}{C}_g{\mu }_g+{\epsilon }_{ks}\left(1-{\varphi }_k\right)D_s\frac{c_{\mu s}ZRT}{{\left(p_L+p_k\right)}^2}{\mu }_g+\frac{{\varphi }_k}{{\tau }_k}\frac{r_k^2}{8}\left[1+\sqrt{\frac{8ZRT}{M}}\frac{{\mu }_g}{p_k{r}_k}\left(\frac{2}{f}-1\right)\right],$$

(4)

where $k_k$ is the apparent permeability of the porous kerogen, m²; ${\varphi }_k$ is the porous kerogen porosity; ${\tau }_k$ is the porous kerogen tortuosity; $r_k$ is the pore size in porous kerogen, m; $C_g$ is the gas compressibility, Pa⁻¹; ${\epsilon }_{ks}$ is proportion of solid kerogen volume in total interconnected matrix; $D_s$ is the surface diffusion coefficient, m²/s; $c_{\mu s}$ is the Langmuir volume on the porous kerogen, mol/m³; $p_L$ is the Langmuir pressure, Pa; $p_k$ is the porous kerogen pressure, Pa; and f is fraction of molecules striking pore wall which are diffusely reflected.

Here, the Warren–Root pseudo-steady-state model is used to describe the cross-flow between porous kerogen and inorganic matter. Therefore, the diffusivity equation for the porous kerogen can be obtained by

$$\frac{\partial \left[{\varphi }_k\left(p_{2k}M/ZRT\right)+{\epsilon }_{ks}\left(1-{\varphi }_k\right)c_s\right]}{\partial t}={\sigma }_k\frac{p_{2k}M{k}_k}{ZRT{\mu }_g}\left(p_{2m}-p_{2k}\right),$$

(5)

where $c_{\mu s}$ is the adsorbed gas volume on the porous kerogen, mol/m³, which follows Langmuir isotherm equation; ${\sigma }_k$ is the pseudo-steady-state shape factor for the porous kerogen, 1/m²; and subscript, 2 represents region II.

Similarly, if the Warren–Root pseudo-steady-state model is also used to describe the cross-flow between inorganic matter and induced-fracture system, the diffusivity equation for the porous kerogen can be expressed as

$$\frac{\partial \left({\varphi }_m{p}_{2m}M/ZRT\right)}{\partial t}={\sigma }_k\frac{p_{2k}M{k}_k}{ZRT{\mu }_g}\left(p_{2k}-p_{2m}\right)-{\sigma }_m\frac{p_{2m}M{k}_m}{ZRT{\mu }_g}\left(p_f-p_{2m}\right),$$

(6)

where ${\varphi }_m$ is the inorganic matter porosity; ${\sigma }_m$ is the pseudo-steady-state shape factor for the inorganic matter, 1/m²; $k_m$ is the apparent permeability of the inorganic matter, m², which includes the Knudsen diffusion, viscous flow, and the slippage effect are included, and can be given by

$$k_m=\frac{{\varphi }_m}{{\tau }_m}\frac{2r_m}{3}\sqrt{\frac{8ZRT}{\pi M}}{C}_g{\mu }_g+\frac{{\varphi }_m}{{\tau }_m}\frac{r_m^2}{8}\left[1+\sqrt{\frac{8ZRT}{M}}\frac{{\mu }_g}{p_m{r}_m}\left(\frac{2}{f}-1\right)\right],$$

(7)

where ${\tau }_m$ is the inorganic matter tortuosity; $r_m$ is the pore size in the inorganic matter, m.

Outer boundary conditions: flux continuity between SRV (region II) and USRV (region IV)

$$\frac{k_w}{{\mu }_g}{\left(\frac{y}{y_w}\right)}^{3\left(d_{fa}-2\right)+\left(d_{fs}-2\right)-\theta }{\frac{p_{f}M}{ZRT}\frac{\partial p_f}{\partial y}|}_{y=y_f}=\frac{k_m}{{\mu }_g}{\frac{p_{4m}M}{ZRT}\frac{\partial p_{4m}}{\partial y}|}_{y=y_f}.$$

(8a)

Inner boundary conditions: pressure continuity between SRV (region II) and HF (region I)

$${p_f|}_{y=y_w}={p_{HF}|}_{y=y_w}.$$

(8b)

In Equation (8a,b), $y_w$ is assumed as hydraulic fracture half-width, m; $y_f$ is the half-width of SRV in y-direction, m; and subscript, 4 represents region IV and HF represents hydraulic fracture.

2.2.2. Diffusivity Equations of FTMM in USRV (Region III, Region IV and Region V)

The fluid flow in both region III and region V can be described as linear flow in x-direction. Thus, according to Equations (5) and (6), the diffusivity equations of porous kerogen and inorganic matter in these two regions can be expressed as, respectively,

$$\{\begin{array}{l}\frac{\partial \left[{\varphi }_k\left(p_{3k}M/ZRT\right)+M{\epsilon }_{ks}\left(1-{\varphi }_k\right)c_s\right]}{\partial t}={\sigma }_k\frac{p_{3k}M{k}_k}{ZRT{\mu }_g}\left(p_{3m}-p_{3k}\right)\hfill \\ \frac{\partial \left({\varphi }_m{p}_{3m}M/ZRT\right)}{\partial t}=\frac{\partial }{\partial x}\left(\frac{p_{3m}M{k}_m}{ZRT{\mu }_g}\frac{\partial p_{3m}}{\partial x}\right)+{\sigma }_k\frac{p_{3k}M{k}_k}{ZRT{\mu }_g}\left(p_{3k}-p_{3m}\right)\hfill \end{array},\mathrm{for}\mathrm{region}\mathrm{III},$$

(9)

$$\{\begin{array}{l}\frac{\partial \left[{\varphi }_k\left(p_{5k}M/ZRT\right)+M{\epsilon }_{ks}\left(1-{\varphi }_k\right)c_s\right]}{\partial t}={\sigma }_k\frac{p_{5k}M{k}_k}{ZRT{\mu }_g}\left(p_{5m}-p_{5k}\right)\hfill \\ \frac{\partial \left({\varphi }_m{p}_{5m}M/ZRT\right)}{\partial t}=\frac{\partial }{\partial x}\left(\frac{p_{5m}M{k}_m}{ZRT{\mu }_g}\frac{\partial p_{5m}}{\partial x}\right)+{\sigma }_k\frac{p_{5k}M{k}_k}{ZRT{\mu }_g}\left(p_{5k}-p_{5m}\right)\hfill \end{array},\mathrm{for}\mathrm{region}\mathrm{V},$$

(10)

where subscript 5 represents region V.

Outer boundary conditions: no-flow conditions for region III and region V

$${\frac{\partial p_{3m}}{\partial x}|}_{x=x_e}=0,$$

(11a)

$${\frac{\partial p_{5m}}{\partial x}|}_{x=x_e}=0.$$

(11b)

Inner boundary conditions: pressure continuity for region III—SRV (region II), and region V—region IV

$${p_{3m}|}_{x=x_f}={p_f|}_{x=x_f},$$

(11c)

$${p_{5m}|}_{x=x_f}={p_{4m}|}_{x=x_f}.$$

(11d)

In Equation (11a–d), $x_e$ is the reservoirs half-size in x-direction, m; subscript 4 represents region IV.

Whereas, the fluid flow in y-direction in region IV, and the diffusivity equation can be given by

$$\{\begin{array}{l}\frac{\partial \left[{\varphi }_k\left(p_{4k}M/ZRT\right)+M{\epsilon }_{ks}\left(1-{\varphi }_k\right)c_s\right]}{\partial t}={\sigma }_k\frac{p_{4k}M{k}_k}{ZRT{\mu }_g}\left(p_{4m}-p_{4k}\right)\hfill \\ \frac{\partial \left({\varphi }_m{p}_{4m}M/ZRT\right)}{\partial t}=\frac{\partial }{\partial y}\left(\frac{p_{4m}M{k}_m}{ZRT{\mu }_g}\frac{\partial p_{4m}}{\partial y}\right)+{\sigma }_k\frac{p_{4k}M{k}_k}{ZRT{\mu }_g}\left(p_{4k}-p_{4m}\right)+\frac{1}{x_f}\frac{p_{4m}M}{ZRT{\mu }_g}{k}_m{\frac{\partial p_{5m}}{\partial x}|}_{x=x_f}\hfill \end{array}$$

(12)

Outer boundary conditions: no-flow conditions for region IV

$${\frac{\partial p_{4m}}{\partial y}|}_{y=y_e}=0.$$

(13a)

Inner boundary conditions: pressure continuity for region IV—SRV (region II)

$${p_{4m}|}_{y=y_f}={p_f|}_{y=y_f}.$$

(13b)

In Equation (13a,b), $y_e$ is the half-size of HF spacing in y-direction, m.

2.2.3. Diffusivity Equations of FTMM in HF (Region I)

Linear flow occurs in the hydraulic fractures with closed tip and constant production in x-direction, and the diffusivity equation can be expressed as

$$\frac{\partial \left({\varphi }_F{p}_{F}M/ZRT\right)}{\partial t}=\frac{\partial }{\partial x}\left(\frac{p_{F}M/ZRT}{{\mu }_g}{k}_F\frac{\partial p_F}{\partial x}\right)+\frac{1}{y_w}\frac{p_{F}M}{ZRT{\mu }_g}{k}_w{\left(\frac{y}{y_w}\right)}^{3\left(d_{fa}-2\right)+\left(d_{fs}-2\right)-\theta }{\frac{\partial p_f}{\partial y}|}_{y=y_w},$$

(14)

where ${\varphi }_F$ is the HF porosity; $p_F$ is the HF pressure, Pa; and $k_F$ is the HF permeability, m².

Outer boundary conditions: no-flow conditions for region I

$${\frac{\partial p_F}{\partial x}|}_{x=x_f}=0.$$

(15a)

Inner boundary conditions: constant production without wellbore storage and skin effect

$${\frac{\partial p_F}{\partial x}|}_{x=0}=-\frac{q_{f}ZRT{\mu }_g}{2M{k}_{F}h{y}_w{p_F|}_{x=0}},$$

(15b)

where $q_f$ is constant well production from HF, kg/s.

For solving the FMFM as shown in Appendix B, we can obtain the pseudo-pressure in well bottom-hole in Laplace domain as

$$\overline{{\phi }_{fD}}={\overline{{\phi }_{FD}}|}_{x_D=0}=\frac{\pi k_w{x}_f}{Z{k}_F{y}_{w}s·\tan \left(\sqrt{\frac{k_f}{k_F}\left({\omega }_{F}s-\frac{x_f}{y_w}{a}_1\right)}\right)\sqrt{\frac{k_f}{k_F}\left({\omega }_{F}s-\frac{x_f}{y_w}{a}_1\right)}},$$

(16)

where, $a_1=\sqrt{a_2}{y}_{wD}^{\frac{1+m_1-2n_1}{2}}\left(\frac{I_{\frac{-1-m_1}{2+m_1-n_1}}\left(\frac{2\sqrt{a_2}}{2+m_1-n_1}{y}_{wD}^{\frac{2+m_1-2n_1}{2}}\right)-\frac{h_1}{h_2}{K}_{\frac{-1-m_1}{2+m_1-n_1}}\left(\frac{2\sqrt{a_2}}{2+m_1-n_1}{y}_{wD}^{\frac{2+m_1-2n_1}{2}}\right)}{y_{wD}^{\frac{1-n_1}{2}}\left(I_{\frac{1-n_1}{2+m_1-n_1}}\left(\frac{2\sqrt{a_2}}{2+m_1-n_1}{y}_{wD}^{\frac{2+m_1-2n_1}{2}}\right)+\frac{h_1}{h_2}{K}_{\frac{1-n_1}{2+m_1-n_1}}\left(\frac{2\sqrt{a_2}}{2+m_1-n_1}{y}_{wD}^{\frac{2+m_1-2n_1}{2}}\right)\right)}\right)$ and other variables are shown in Appendix B.

Based on Equation (16), if the wellbore storage and skin factor are not considered, the dimensionless flow rate can be derived as [33,34]

$$\overline{q_D}=\frac{Z{k}_F{y}_w·\tan \left(\sqrt{\frac{k_f}{k_F}\left({\omega }_{F}s-\frac{x_f}{y_w}{a}_1\right)}\right)\sqrt{\frac{k_f}{k_F}\left({\omega }_{F}s-\frac{x_f}{y_w}{a}_1\right)}}{\pi s{k}_w{x}_f}$$

(17)

If we hope to obtain the solution in the real-time domain, here, the Stehfest algorithm is applied, and Equation (17) becomes

$$q_D=\frac{\ln 2}{t_D}{\displaystyle \sum }_{i=1}^N{\left(-1\right)}^{\frac{N}{2}+1}\times \left\{\begin{array}{c}\left[{\displaystyle {\displaystyle \sum }_{k=\frac{i+1}{2}}^{\min \left(i,\frac{N}{2}\right)}}\frac{k^{\frac{N}{2}}\left(2k\right)!}{\left(\frac{N}{2}-k\right)! k! \left(k-1\right)! \left(i-k\right)! \left(2k-i\right)!}\right]\\ \times \frac{Z{k}_F{y}_w·\tan \left(\sqrt{\frac{k_f}{k_F}\left({\omega }_{F}s-\frac{x_f}{y_w}{a}_1\right)}\right)\sqrt{\frac{k_f}{k_F}\left({\omega }_{F}s-\frac{x_f}{y_w}{a}_1\right)}}{\pi s{k}_w{x}_f}\times \left(\frac{\ln 2}{t_D}i\right)\end{array}\right\}.$$

(18)

3. Model Verification and Discussion

3.1. Model Verification

The SRV size was assumed to be the same as with the half-length of HF spacing (y_f = y_e), the presented FMFM model can be simplified as the model proposed by Sheng et al. [20]. The basic parameters used for FMFM and the Sheng’s model are listed in

Table 1. Parameters used in multi-linear flow model (FMFM) and Sheng’s model [20].

Parameters	Value	Parameters	Value
Fractal dimension of induced-fracture spacing, dfs	1.95	Pore size in porous kerogen, rk (m)	10 × 10−9
Fractal dimension of induced-fracture aperture, dfa	2.0	Portion of kerogen volume, εks	0.5
Tortuosity index of induced-fracture, θ	−0.05	Porous kerogen porosity, ϕk	0.1
Porosity of induced-fracture, ${\varphi }_f$	10−4	Porous kerogen tortuosity, τk	5
HF half-spacing, ye (m)	100	Langmuir pressure, pL (MPa)	13.78
HF half-length, xf (m)	50	Gas viscosity, μg (mPa·s)	0.0184
HF permeability, kF (m2)	103 × 10−15	Gas compressibility, Cg (MPa−1)	5 × 10−2
HF half-width, yw (m)	0.01	Surface diffusion coefficient, Ds (m2/s)	1 × 10−5
HF porosity, ϕF	10−3	Molecular mass of shale gas, M (kg/mol)	0.016
Total compressibility of the inorganic matter, Ctm (MPa−1)	7.5 × 10−3	Langmuir volume on kerogen surface, cμs (mol/m3)	700
Total compressibility of induced-fracture, Ctf (MPa−1)	4 × 10−3	Fraction of molecules striking pore wall which are diffusely reflected, f	0.8
Inorganic matter porosity, ϕm	0.1	Total compressibility of the porous kerogen, Ctk (MPa−1)	7.5 × 10−3
Inorganic matter tortuosity, τm	5	Reservoir thickness, h (m)	19
Induced-fracture permeability at yw, kw (m2)	2 × 10−15	Initial pressure, pi (MPa)	17
Pore size in inorganic matter, rm (m)	20 × 10−9	Well bottom pressure, pwf (MPa)	12
Formation temperature, T (K)	338	Reservoir half-width, xe (m)	200

. The comparison of two models are shown as Figure 4. The results indicate that the production rate and cumulative production calculated by the FMFM presented in this study can fit well with those of Sheng’s model. The classical trilinear flow model can be simplified by FMFM when the SRV size between HFs and the fractal distribution of induced-fracture are not considered. In addition, if the shale gas transport mechanisms are ignored in FMFM, the model can be used to predict the production of MFHWs in tight gas and even in tight oil reservoirs fast.

3.2. Sensitivity Analysis of Properties of SRV

Figure 5 and Figure 6 show that the dimensionless pseudo-pressure and dimensionless production rate with different tortuosity index and fractal dimensions of induced-fracture spacing by FMFM, which reflects the various properties of SRV for MFHWs in unconventional gas reservoirs. Figure 5 shows that the dimensionless pseudo-pressure is larger and the dimensionless production rate is smaller when the tortuosity index of induced-fracture θ increases. The larger tortuosity index of induced-fracture θ means the longer flow path of gas transport in SRV. When the tortuosity index of induced-fracture θ is larger than 0, the properties of induced-fracture will decrease far from HF. When the tortuosity index of induced-fracture θ is equal to 0, the properties of induced-fracture is homogeneous in SRV.

Figure 6 indicates that the dimensionless pseudo-pressure is smaller and the dimensionless production rate is larger when the fractal dimension of induced-fracture d_fs increases. The larger fractal dimension of induced-fracture d_fs leads to a larger equivalent permeability in SRV. In addition, the linear low in SRV can make the WHFWs produce more gas from the formation, which suggests that the induced-fracture system can improve the development effect for unconventional reservoirs.

4. Application of FMFM with ES-MDA Method

4.1. ES-MDA Method

In order to avoid the model update and parameter-state inconsistency, which makes the ensemble Kalman filter (EnKF) complicated, and achieve a better data matching with the ensemble smoother (ES), Emerick and Reynolds [35] reported ES-MDA by assimilating the same observed data multiple times with the covariance matrix of the measurement errors simultaneously. Numerical experiments and examples have shown that the ES-MDA method can provide a better data matching and reduce the uncertainty of model description than the ES method and the EnKF method. Based on the synthetic case and field case data, the application of ES-MDA history matching technology and FMFM for MFHWs is discussed in this section.

The reservoir simulation models are typically stable functions of the rock property fields, which is different with the oceanic and atmospheric models. Based on the common assumption in history-matching problems that the model uncertainty is ignored, we just need to take the parameter-estimation problem into account for ES. Thus, the analyzed vector of model parameters can be given by

$${\mathit{m}}_j^a={\mathit{m}}_j+{\mathit{C}}_{MD}{\left({\mathit{C}}_{DD}+{\mathit{C}}_D\right)}^{-1}\left({\mathit{d}}_{uc,j}-{\mathit{d}}_j\right), for j=1,\cdots ,N_e,$$

(19)

where ${\mathit{C}}_{MD}$ is the cross-covariance matrix between the prior vector of model parameters ${\mathit{m}}_j$ and the vector of predicted data $\mathit{d}$; ${\mathit{C}}_{DD}$ is the N_d—dimension auto-covariance matrix of predicted data; N_d is the total number of measurements assimilated; ${\mathit{d}}_{uc}~\mathcal{N}\left({\mathit{d}}_{obs},{\mathit{C}}_D\right)$; ${\mathit{d}}_{obs}$ is the N_d—dimensional vector of observed data; ${\mathit{C}}_D$ is the N_d—dimension covariance matrix of observed data measurement errors; and N_e is the number of ensemble member.

Emerick and Reynolds [35] have mentioned that ES-MDA can be used in the nonlinear case because ES is equivalent to a single Gauss–Newton iteration with a full step and an average sensitivity estimated from the prior ensemble, in which case the MDA can be interpreted as an “iterative” ensemble smoother. The ES-MDA algorithm is presented as follows:

(1)

Estimate the number of data assimilations N_a, and the coefficients ${\alpha }_i$ (${\displaystyle \sum }_{i=1}^{N_a}\frac{1}{{\alpha }_i}=1$), for $i=1,\cdots ,N_a$.

(2)

For $i=1 to N_a$

• Run the ensemble from time zero for obtaining the vector of predicted data

  $${\mathit{d}}_j=\mathit{g}\left({\mathit{m}}_j\right), for j=1,\cdots ,N_e,$$

  (20)

  where $\mathit{g}\left(·\right)$ is the nonlinear forward model, ${\mathit{d}}_j$ is assignment the model parameters to vector ${\mathit{m}}_j$ at time zero.
• For each ensemble member, perturb the observation vector by

  $${\mathit{d}}_{uc,j}={\mathit{d}}_{obs}+\sqrt{{\alpha }_i}{\mathit{C}}_D^{1/2}{z}_{d,j}, for j=1,\cdots ,N_e,$$

  (21)

  where $z_{d,j}~\mathcal{N}\left(0,{\mathit{{\rm I}}}_{N_d}\right)$.
• Update the ensemble

  $${\mathit{m}}_j^a={\mathit{m}}_j+{\mathit{C}}_{MD}{\left({\mathit{C}}_{DD}+{\alpha }_i{\mathit{C}}_D\right)}^{-1}\left({\mathit{d}}_{uc,j}-{\mathit{d}}_j\right), for j=1,\cdots ,N_e,$$

  (22)

  where $\mathit{C}={\mathit{C}}_{DD}+{\alpha }_i{\mathit{C}}_D\equiv \left[\begin{array}{cc}\begin{array}{cc}{\mathit{C}}_{DD}+{\alpha }_1{\mathit{C}}_D& {\mathit{C}}_{DD}\\ {\mathit{C}}_{DD}& {\mathit{C}}_{DD}+{\alpha }_2{\mathit{C}}_D\end{array}& \begin{array}{cc}\cdots & {\mathit{C}}_{DD}\\ \cdots & {\mathit{C}}_{DD}\end{array}\\ \begin{array}{cc}⋮& \\ {\mathit{C}}_{DD}& \cdots \end{array}& \begin{array}{cc}\ddots & ⋮\\ & {\mathit{C}}_{DD}+{\alpha }_{N_a}{\mathit{C}}_D\end{array}\end{array}\right]$; ${\mathit{C}}_{MD}=\frac{1}{N_e-1}{\displaystyle \sum }_{j=1}^{N_e}\left({\mathit{m}}_j-\overline{\mathit{m}}\right){\left({\mathit{d}}_j-\overline{\mathit{d}}\right)}^T$; ${\mathit{C}}_{DD}=\frac{1}{N_e-1}{\displaystyle \sum }_{j=1}^{N_e}\left({\mathit{d}}_j-\overline{\mathit{d}}\right){\left({\mathit{d}}_j-\overline{\mathit{d}}\right)}^T$; $\overline{\mathit{m}}$ and $\overline{\mathit{d}}$ represent the average values of model parameters and prediction parameter ensembles respectively.

In the procedure of MDA algorithm mentioned above, the data are assimilated N_a times continuously, and the ensemble needs to be rerun before each data assimilation. At the same time, in the ES-MDA algorithm, every time the data is assimilated repeatedly, the disturbance observation vector is resampled. The procedure will reduce the sampling problems due to the matching outliers that may occur when the observation data is perturbed.

4.2. Synthetic Case

Based on the simulator Eclipse, the numerical model of a MFHW in the unconventional gas reservoir is established. The numbers of grids are X × Y × Z = 80 × 260 × 3 and the grid steps are X × Y × Z = 5 m × 5 m × 5 m, respectively. A dual-porosity media model is used to describe the SRV (DUALPORO keyword), and a single-porosity medium model is used to describe the USRV. Hydraulic fractures are described by LGR keyword as shown in Figure 7. If only the gas viscous flow with slippage effect is considered, then the FMFM can be simplified as a semianalytic model for a MFHW with the homogeneous SRV in the tight gas reservoir. We set the same modeling parameters of FMFM with those of numerical model and then discuss the applicability of FMFM with ES-MDA.

Firstly, according to the dimensionless parameters listed in

Table A1. Basic parameters used in FMFM.

Parameters	Expression
Pseudo-pressure for gas phase	$\phi \left(p\right)=2{{\displaystyle \int }}_{p_i}^p\left(p/{\mu }_{g}Z\right)\mathrm{d}p$
Dimensionless pseudo-pressure	${\phi }_D=2\pi k_{w}h{c}_i\left({\phi }_i-\phi \right)/\left(p_i{q}_f\right)$
Dimensionless time	$t_D=k_{w}t/\left({\left(\varphi C_t\right)}_{f-m-k}{\mu }_g{x}_f^2\right)$
Dimensionless length	$x_D=x/x_f$, $y_D=y/y_w$, $x_{eD}=x_e/x_f$, $y_{eD}=y_e/x_f$, $y_{fD}=y_f/x_f$, $y_{wD}=y_w/x_f$
Total compressibility coefficient of porous kerogen	$C_{tk}=C_g+\left({\epsilon }_{ks}{c}_{\mu s}ZRT\right)/\left({\varphi }_k{\left(p_L+p_k\right)}^2\right)$
Transfer coefficient from porous kerogen to inorganic matter	${\lambda }_k={\sigma }_k{k}_k{x}_f^2/k_w$
Transfer coefficient from inorganic matter to induced-fracture	${\lambda }_m={\sigma }_m{k}_m{x}_f^2/k_w$
Storage coefficient of porous kerogen	${\omega }_k=\left({\varphi }_k{C}_{tk}\right)/{\left(\varphi C_t\right)}_{f-m-k}$
Storage coefficient of inorganic matter	${\omega }_m=\left({\varphi }_m{C}_g\right)/{\left(\varphi C_t\right)}_{f-m-k}$
Storage coefficient of HFs	${\omega }_F=\left({\varphi }_F{C}_g\right)/{\left(\varphi C_t\right)}_{f-m-k}$

Subscript i represents the initial conditions; f, m, k represents fracture, inorganic matter, and porous kerogen respectively.

, the production rate of the numerical model can be non-dimensionalized as $q_D=p{q}_f/\left(2\pi k_{w}h{c}_i\left(\phi -{\phi }_{wf}\right)\right)$.

Then, the N_d—dimensional vector of observed data can be chosen as:

$${\mathit{d}}_{obs}={\left[{\mathit{d}}_s^T,{\mathbf{q}}_o^T\right]}_{N_d\times 1}^T,$$

(23)

where d_s donates the N_d-dimensional vector of static parameters of the reservoir and MFHW; q_o donates the dimensionless production rate of numerical model.

The N_d × N_d covariance matrix of observed data measurement errors is given by

$${\mathit{C}}_D={\left[\begin{array}{ccc}{\sigma }_{d_1}^2& & \\ & \ddots & \\ & & {\sigma }_{d_{N_d}}^2\end{array}\right]}_{N_d\times N_d},$$

(24)

where σ_d is equal to 0.3 for real production rate.

The matching results of dimensionless production rate by ES-MDA after 4 iterations and error analysis are shown in Figure 8. The values of some variables are listed as N_p = 23, N_d = 115, N_e = 100, and $\sqrt{{\alpha }_i}$ = 2. Meanwhile, the root-mean-square error and the average objective function in Figure 7 can be expressed as, respectively,

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{j=1}^{Ne}{\left({\mathit{d}}_{sobs}-{\mathit{d}}_{sj}\right)}^2}{N_p}},$$

(25)

$$O\left(\mathit{m}\right)=\frac{1}{N_d}{\left[\mathit{g}\left(\mathit{m}\right)-{\mathit{d}}_{obs}\right]}^T{C}_D^{-1}\left[\mathit{g}\left(\mathit{m}\right)-{\mathit{d}}_{obs}\right].$$

(26)

Figure 8 shows that the error between the ensemble mean model and that of true model is large with one iteration, and the distribution range of ensemble members is large. All models of the ensemble (light blue curves) represent the results calculated by FMFM, the ensemble mean (blue curve) represents the calculated results matching with unperturbed observation data (average value of ensemble members), the true value (red curve) is the true dimensionless production rate obtained by FMFM, and the numerical results (red point) is the observed data in Figure 8. We can see that with two iterations, the distribution range of ensemble members becomes small, and ensemble mean model basically coincides with the true model. With three and four iteration, the distribution range of ensemble members is further reduced around the true model, and the ensemble mean model and the observed data are well-matching, which is almost the same with the true model. In addition, the RMSE and objective function values of the prior model are large, indicating that the model errors and data errors are large. After one iteration, the error decreases. The RMSE and objective function values becomes small and basically the same after three and four iterations, which indicates that the matching results are good as shown in

Table 2. Comparison of matching results of FMFM and numerical model parameters by ES-MDA.

Parameters	Eclipse	FMFM	Parameters	Eclipse	FMFM
Total compressibility coefficient, MPa−1	5.5 × 10−3	5.2714 × 10−3	USRV matrix permeability, m2	1 × 10−19	0.887 × 10−18
HF permeability, m2	15 × 10−13	7.82 × 10−13	Induced-fracture permeability, m2	1 × 10−16	1.33 × 10−16
SRV matrix permeability, m2	1 × 10−18	0.77 × 10−18	Other parameters	Equal values

. Therefore, the discussion mentioned above suggests that parameters of FMFM can be obtained by automatically matching production data of the numerical model by ES-MDA method. The technique can also obtain reservoir physical properties and fracturing parameters by matching the actual production data of the oilfield.

4.3. Field Case

The presented FMFM model with ES-MDA history matching method was applied based the shale gas production data of a MFHW in western China [36]. Compared with the presented models based on different simulating conditions, the matching results are shown in Figure 9. Noting that the actual shale gas production rate of a MFHW in western China is represented by red points; when d_fs = d_fa = 2, θ = 0, y_fD = y_eD, the FMFM can be assumed as a typical trilinear model with homogenous SRV (black line); the results of Sheng’s FMPM model [20] with fractal SRV are drawn by the light green line, and they calculated the d_fs= 1.90 by the box-counting method and the fractal random-fracture-network algorithm and θ = −0.05 by the random walk method; based on the presented FMFM with ES-MDA, the fractal dimensions, tortuosity index, and SRV size can be obtained as shown by the blue line. Figure 9 shows that the fractal dimension and tortuosity index of induced-fracture system matched by FMFM model based on ES-MDA approximate the results of Sheng’s FMPM model. The unmatched early-time data was caused by the early-time flow back process. In addition, the results dimensionless production rate calculated by FMFM were smaller but matched better with actual data than Sheng’s FMPM model when the SRV size was taken into account. This section mainly provides an application case of our presented approach.

5. Conclusions

In this paper, a semianalytic fractal multi-linear flow model (FMFM) for MFHWs in unconventional gas reservoirs with consideration of the heterogeneous distribution of the properties of induced-fracture system is proposed. Fractal dimensions of induced-fracture spacing and aperture and tortuosity index of induced-fracture system are included based on fractal theory to describe the properties of SRV region. For shale gas reservoirs, gas transport mechanisms (viscous flow with slippage, Knudsen diffusion, and surface diffusion) among various medium including porous kerogen, inorganic matter, and fracture system are take into account in FMFM. Then, combining with ES-MDA, the FMFM can be applied not only for prediction of gas production rate for MFHWs in unconventional gas reservoirs but also for main uncertainty parameters matching. Some findings can be drawn as follows:

(1)

The permeability and porosity of the induced-fracture system are affected by the heterogeneous distribution of induced-fractures spacing and aperture. When the fractal dimensions of induced-fracture spacing (d_fs) and aperture (d_fa) are smaller than 2.0, the permeability of the induced-fracture system decreases with the increase of the distance from HFs in SRV region, which will become increase if d_fs > 2.0 or d_fa > 2.0. When the tortuosity index of the induced-fracture system θ is larger than 0, the permeability of induced-fracture system also decreases with the increase of the distance from HFs in SRV region.

(2)

The FMFM divides the formation into three types of porous media in shale reservoirs (porous kerogen, inorganic matter, and fracture system). Triple-porosity media and dual-porosity media are used to describe the fractal SRV region and USRV region, respectively. Multiple gas transport mechanisms such as viscous flow with slippage, Knudsen diffusion, and surface diffusion are considered, and gas flow behaviors in different regions are coupled by pseudo-steady cross-flow among various media. The FMFM is verified by other presented model and the results show that the large d_fs or small θ causes the small average permeability of the induced-fracture system, which results in large dimensionless pseudo-pressure and small dimensionless production rate.

(3)

Combining the FMFM with ES-MDA history-matching method, the synthetic case for the tight gas reservoir and field case for the shale gas reservoir are discussed. Various main parameters inversions of HFs and NFs are conducted in SRV region. Thus, the presented method can be applied for gas production predicting and history-matching in unconventional gas reservoirs.