*3.1. Model Verification*

The SRV size was assumed to be the same as with the half-length of HF spacing (*yf* = *ye*), 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].


**Figure 4.** Comparison of production rate and cumulative production between FMFM and Sheng's model [20].

#### *3.2. Sensitivity Analysis of Properties of 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 5.** Comparison of production rate and pseudo-pressure with different tortuosity index of induced-fracture θ by FMFM (*dfs* = 2 and *dfa* = 2).

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.** Comparison of production rate and pseudo-pressure with different fractal dimension of induced-fracture spacing *dfs* by FMFM (*dfa* = 2 and θ = −0.05).

#### **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

$$m\_j^a = m\_j + \mathcal{C}\_{\rm MD} (\mathcal{C}\_{\rm DD} + \mathcal{C}\_D)^{-1} \{d\_{\rm nc,j} - d\_j\}, \text{ for } j = 1, \dots, N\_\varepsilon,\tag{19}$$

where *CMD* is the cross-covariance matrix between the prior vector of model parameters *<sup>m</sup>j* and the vector of predicted data *d*; *CDD* is the *Nd*—dimension auto-covariance matrix of predicted data; *Nd* is the total number of measurements assimilated; *duc* ∼ N(*dobs*,*CD*); *dobs* is the *Nd*—dimensional vector of observed data; *CD* 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:

$$\text{(1)}\quad \text{Estimate the number of data assumptions } \mathbf{N\_{a}}, \text{ and the coefficients } a\_{i} \left(\sum\_{i=1}^{N\_{a}} \frac{1}{a\_{i}} = 1\right) \text{, for } i = 1, \dots, N\_{a}.$$

$$(2) \quad \text{For } i = 1 \text{ to } N\_a.$$

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

$$d\_j = \mathbf{g}(\mathfrak{m}\_j), \text{ for } j = 1, \dots, N\_{\mathfrak{c}\_\nu} \tag{20}$$

where *<sup>g</sup>*(·) is the nonlinear forward model, *dj* is assignment the model parameters to vector *<sup>m</sup>j* at time zero.

b. For each ensemble member, perturb the observation vector by

$$d\_{\rm nc,j} = d\_{\rm obs} + \sqrt{\alpha\_i} \mathcal{C}\_{\rm D}^{1/2} z\_{d,j}, \text{ for } j = 1, \cdots, N\_{\rm \varepsilon} \tag{21}$$

where *zd*,*j* ∼ <sup>N</sup>0,*INd* .

c. Update the ensemble

$$\mathbf{m}\_{j}^{a} = \mathbf{m}\_{\rangle} + \mathbf{C}\_{\text{MD}} (\mathbf{C}\_{\text{DD}} + \alpha\_{i} \mathbf{C}\_{\text{D}})^{-1} \left(d\_{\text{uc},j} - d\_{\text{j}}\right), \text{ for } j = 1, \cdots, N\_{\text{c}}.\tag{22}$$

$$\begin{array}{ccccc}\text{where} & \mathsf{C} & = & \mathsf{C}\_{DD} + a\_{i}\mathsf{C}\_{D} & = & \begin{bmatrix} \mathsf{C}\_{DD} + a\_{1}\mathsf{C}\_{D} & \mathsf{C}\_{DD} & & \cdots & \mathsf{C}\_{DD} \\ & \mathsf{C}\_{DD} & \mathsf{C}\_{DD} + a\_{2}\mathsf{C}\_{D} & & \cdots & \mathsf{C}\_{DD} \\ & & \vdots & & \ddots & \vdots \\ & & \mathsf{C}\_{DD} & \cdots & & \mathsf{C}\_{DD} + a\_{N\_{d}}\mathsf{C}\_{D} \end{bmatrix}; \\\mathsf{C}\_{MD} = & \frac{1}{N\_{c}-1} \sum\_{j=1}^{N\_{c}} \left( m\_{j} - \overline{m} \right) \left( d\_{j} - \overline{d} \right)^{T}; \ \mathsf{C}\_{DD} = \frac{1}{N\_{c}-1} \sum\_{j=1}^{N\_{c}} \left( d\_{j} - \overline{d} \right) \left( d\_{j} - \overline{d} \right)^{T}; \ \overline{m} \text{ and } \overline{d} \text{ represent the automorphism parameter or annihilate coordinates.} \end{array}$$

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.
