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

: Given reliable parameters, a newly developed semianalytic model could o ﬀ er an e ﬃ cient 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 signiﬁcant 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 ﬂow with slippage, Knudsen di ﬀ usion, and surface di ﬀ usion, among multi-media including porous kerogen, inorganic matter, and fracture system are taken into account and the model is veriﬁed. Then, the e ﬀ ects 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 ﬁeld 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 ﬁeld.


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

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 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].
Energies 2020, 13,3718 3 of 20 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.

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

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 (dfs and dfa) 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 where kf is the fracture permeability, m 2 ; μ is the fluid viscosity, Pa·s; Ar is the sectional-area of a REV, m 2 ; Ap is the sectional-area of induced-fracture, m 2 ; dL is the flow length in the induced-fracture, m; kw is the fracture permeability at the reference point, m 2 ; yw is the position of reference point, m; n is the number of fracture sites per unit thickness, m −1 ; dfs is the fractal dimension of induced-fracture spacing; dfa 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] ( ) = ( ) ( ) + ( ) = ,

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 where k f is the fracture permeability, m 2 ; µ is the fluid viscosity, Pa·s; A r is the sectional-area of a REV, m 2 ; A p is the sectional-area of induced-fracture, m 2 ; dL is the flow length in the induced-fracture, m; k w is the fracture permeability at the reference point, m 2 ; y w is the position of reference point, m; n is the number of fracture sites per unit thickness, m −1 ; 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] where φ f is the fracture porosity; φ i− f is the single induced-fracture porosity; φ 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 −15 m 2 , 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 Energies 2020, 13, 3718 5 of 20 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].
porosity at the reference point; sf is induced-fracture spacing, m; and bf 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 yw = 0.1 m and kw = 10×10 −15 m 2 , the fracture permeability distribution with different fractal dimensions of inducedfracture spacing and aperture (dfs and dfa) and tortuosity index of induced-fracture (θ) are shown in Figure 3. The fracture permeability kf decreases if dfs < 2 and increases if dfs > 2 with the increase of the distance from the reference point yw as shown in Figure 3a. Meanwhile, the fracture permeability kf decreases if dfa < 2 and increases if dfa > 2 with the increase of the distance from the reference point yw 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 yw [31]. The reasons can be illustrated as follows: If dfs < 2 or dfa < 2, as the distance from the reference point increases, the induced-fracture number or aperture decreases and the induced-fracture spacing increases; if dfs = 2 or dfa = 2, as the distance from the reference point increases, the induced-fracture number or aperture and the induced-fracture spacing are constants; if dfs > 2 or dfa > 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].

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 Energies 2020, 13, 3718 6 of 20 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.

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 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 2 ; 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 is the mass transfer from inorganic matter to the induced-fracture system, kg/m 3 ·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) where k k is the apparent permeability of the porous kerogen, m 2 ; φ k is the porous kerogen porosity; τ k is the porous kerogen tortuosity; r k is the pore size in porous kerogen, m; C g is the gas compressibility, Pa −1 ; ε ks is proportion of solid kerogen volume in total interconnected matrix; D s is the surface diffusion coefficient, m 2 /s; c µs is the Langmuir volume on the porous kerogen, mol/m 3 ; 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 where c µs is the adsorbed gas volume on the porous kerogen, mol/m 3 , which follows Langmuir isotherm equation; σ k is the pseudo-steady-state shape factor for the porous kerogen, 1/m 2 ; 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 where φ m is the inorganic matter porosity; σ m is the pseudo-steady-state shape factor for the inorganic matter, 1/m 2 ; k m is the apparent permeability of the inorganic matter, m 2 , which includes the Knudsen diffusion, viscous flow, and the slippage effect are included, and can be given by where τ 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) Inner boundary conditions: pressure continuity between SRV (region II) and HF (region I) 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.

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, , for region III, where subscript 5 represents region V. Outer boundary conditions: no-flow conditions for region III and region V Inner boundary conditions: pressure continuity for region III-SRV (region II), and region V-region IV 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 In Equation (13a,b), y e is the half-size of HF spacing in y-direction, m.

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 where φ F is the HF porosity; p F is the HF pressure, Pa; and k F is the HF permeability, m 2 .
Outer boundary conditions: no-flow conditions for region I Inner boundary conditions: constant production without wellbore storage and skin effect 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 where, 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] Energies 2020, 13, 3718 9 of 20 If we hope to obtain the solution in the real-time domain, here, the Stehfest algorithm is applied, and Equation (17) becomes

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. 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. 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, d fs  Figures 5 and 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 dfs increases. The larger fractal dimension of induced-fracture dfs 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.  Figures 5 and 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.  Figures 5 and 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 dfs increases. The larger fractal dimension of induced-fracture dfs 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.  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.

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 historymatching problems that the model uncertainty is ignored, we just need to take the parameterestimation problem into account for ES. Thus, the analyzed vector of model parameters can be given by where is the cross-covariance matrix between the prior vector of model parameters and the vector of predicted data ; is the Nd-dimension auto-covariance matrix of predicted data; Nd is the total number of measurements assimilated; ~( , ); is the Nd-dimensional vector of observed data; is the Nd-dimension covariance matrix of observed data measurement errors; and Ne 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 Na, and the coefficients (∑ = 1), for = 1, ⋯ , .
(2) For = 1 a. Run the ensemble from time zero for obtaining the vector of predicted data where (•) is the nonlinear forward model, is assignment the model parameters to vector at time zero.
b. For each ensemble member, perturb the observation vector by

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 where C MD is the cross-covariance matrix between the prior vector of model parameters m j and the vector of predicted data d; C DD is the N d -dimension auto-covariance matrix of predicted data; N d is the total number of measurements assimilated; d uc ∼ N (d obs , C D ); d obs is the N d -dimensional vector of observed data; 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 α i ( (2) For i = 1 to N a a.
Run the ensemble from time zero for obtaining the vector of predicted data where g(·) is the nonlinear forward model, d j is assignment the model parameters to vector m j at time zero. b.
For each ensemble member, perturb the observation vector by where z d,j ∼ N 0, I N d .
c. Update the ensemble where . . .
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.

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.
where , ~ 0, . c. Update the ensemble and represent the average values of model parameters and prediction parameter ensembles respectively. In the procedure of MDA algorithm mentioned above, the data are assimilated Na 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.

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, the production rate of the numerical model can be non-dimensionalized as = 2 ℎ − ⁄ . Then, the Nd-dimensional vector of observed data can be chosen as: where ds donates the Nd-dimensional vector of static parameters of the reservoir and MFHW; qo donates the dimensionless production rate of numerical model. The Nd × Nd covariance matrix of observed data measurement errors is given by where σd is equal to 0.3 for real production rate. Firstly, according to the dimensionless parameters listed in Table A1, the production rate of the numerical model can be non-dimensionalized as q D = pq f / 2πk w hc i ϕ − ϕ w f .
Then, the N d -dimensional vector of observed data can be chosen as: 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 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 √ α i = 2. Meanwhile, the root-mean-square error and the average objective function in Figure 7 can be expressed as, respectively, 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 Np = 23, Nd = 115, Ne = 100, and =2. Meanwhile, the root-mean-square error and the average objective function in Figure 7 can be expressed as, respectively, 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. 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.   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. 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.

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.

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 dfs = dfa = 2, θ = 0, yfD = yeD, 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 dfs= 1.90 by the box-counting method and the fractal random-fracturenetwork 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.

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 Figure 9. The matching results of actual shale gas production rate and different models.

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.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Derivation of Fractal Permeability and Porosity of Induced-Fracture
As shown in Figure 2c, using the fractal theory, the total number of induced-fracture, which is relevant to the fractal dimension of induced-fracture spacing d fs , can be expressed as [37] L x,y 0 n L x,y dL x,y ∝ L d f s x,y , where L x,y is the position of system, m. Then, Chang and Yortsos [16] presented the expression of the distribution of fracture site as n L x,y = αL where D is the Euclidean dimension, which is equal to 2 in Cartesian coordinate system; α is a constant related to fracture porosity. For the fracture site, we can obtain s f L x,y n L x,y + b f L x,y n L x,y = A SRV , where A SRV is the hydraulic fracture half-length, m. Substituting Equation (A2) into Equation (A3), Equation (A3) can be rewritten as In the actual unconventional reservoirs, the value of induced-fracture aperture is commonly far less than that of induced-fracture spacing (b f << s f ), and thus Equation (A4) becomes If the induced-fracture spacing at the plane of HF is known, which can be chosen as a reference point, the FFSD in Cartesian coordinate system can be obtained by Equation (A5) s f (y) = s f w y y w where s fw is the induced-fracture spacing at reference point, m. For the distribution of induced-fracture aperture (FFAD), the basic form can also be obtained according to Equation (A6) [21] b f (y) = b f w y y w where b fw is the induced-fracture aperture at reference point, m. Bai and Elsworth [38] reported the expression of a fracture with a slab shape as In a representative elementary volume (REV) of fracture network, the fluid flux in the fracture system is given by v = A r k f (y) µ dp dx = n(y)A p k i− f (y) µ dp dL , where v is flux of induced-fracture system in a REV, m 3 /s; µ is the fluid viscosity, Pa·s; and dp is the pressure-difference of a REV in y direction, Pa. Thus, the fracture permeability and porosity can be obtained as Equations (1) where a 4 = a 5 − k w k m y 2 wD √ a 5 tan h (y lD − y eD ) √ a 5 .
(2) Dimensionless forms and solutions for SRV.
(3) Dimensionless forms and solutions for HFs.
According to Equation (14), the dimensionless diffusivity equations in Laplace domain for HFs can be given by Solving Equation (A16) with boundary conditions, we can obtain the pseudo-pressure in well bottom-hole in Laplace domain as Equation (16).