Multi-Objective History Matching with a Proxy Model for the Characterization of Production Performances at the Shale Gas Reservoir

Abstract

This paper presents a fast, reliable multi-objective history-matching method based on proxy modeling to forecast the production performances of shale gas reservoirs for which all available post-hydraulic-fracturing production data, i.e., the daily gas rate and cumulative-production volume until the given date, are honored. The developed workflow consists of distance-based generalized sensitivity analysis (DGSA) to determine the spatiotemporal-parameter significance, fast marching method (FMM) as a proxy model, and a multi-objective evolutionary algorithm to integrate the dynamic data. The model validation confirms that the FMM is a sound surrogate model working within an error of approximately 2% for the estimated ultimate recovery (EUR), and it is 11 times faster than a full-reservoir simulation. The predictive accuracy on future production after matching 1.5-year production histories is assessed to examine the applicability of the proposed method. The DGSA determines the effective parameters with respect to the gas rate and the cumulative volume, including fracture permeability, fracture half-length, enhanced permeability in the stimulated reservoir volume, and average post-fracturing porosity. A comparison of the prediction accuracy for single-objective optimization shows that the proposed method accurately estimates the recoverable volume as well as the production profiles to within an error of 0.5%, while the single-objective consideration reveals the scale-dependency problem with lesser accuracy. The results of this study are useful to overcome the time-consuming effort of using a multi-objective evolutionary algorithm and full-scale reservoir simulation as well as to conduct a more-realistic prediction of the shale gas reserves and the corresponding production performances.

1. Introduction

The production performances of shale gas reservoirs are quite distinct from those of conventional gas fields in that the former involve a long transient period and are placed under the control of human-made treatment for the entire field life, e.g., in hydraulic fracturing and well stimulation. The typical uncertainty for post-hydraulic-fracturing characterization of gas production is whether the artificial fractures are well-developed and maintained within the shale deposit. In addition to this topological uncertainty, and regardless of whether microseismic monitoring can be used to estimate the stimulated reservoir volume, there is the uncertainty regarding the flow properties, e.g., the reservoir properties within the stimulated reservoir volume as the product of the hydraulic fracturing and the flow characteristics that occur between the fracture and the matrix. The reservoir heterogeneity affects the entire performance, but the lack of available data makes it difficult to obtain the actual values of the stimulated reservoir volume.

For shale gas reservoirs, it is a challenge to reliably and accurately characterize production performance. Widely-used schemes for shale gas reservoirs can be divided into rate transient analysis (RTA) and the time-consuming reservoir simulations. RTA is known as a simple and fast approach with an accuracy that decreases significantly as the reservoir heterogeneity increases [1]. A reservoir simulation can also include both geological features and the nonlinear relationship between the flow parameters to improve the predictability, but a substantial computational effort is needed [2]. Various proxy models have also been studied, including a streamline simulation [3,4,5], an artificial neural network (ANN) [6,7,8], and the fast marching method (FMM) [9,10,11,12,13], to reduce the computational time to characterize unconventional resources. When compared to the other models, FMM has been actively investigated in terms of its production-characteristic analysis for shale gas reservoirs. FMM can estimate the performance using diffusive time of flight (DTOF; the arrival time of pressure front) and can determine key parameters such as drainage volume, gas-production rate, and bottom-hole pressure, all of which can be estimated without the need for time-consuming simulations.

History matching is an inverse-modeling method that is used to generate equiprobable geomodels that show production behavior similar to those of the observed production history. The method satisfies the objective constraints by updating the geomodels, and the reservoir heterogeneity as well as the nonlinear characteristics can therefore be included in the optimized geomodels [14]. To match the nonlinear behavior of the reservoir, the typical history-matching technique is based on single-objective optimization algorithms such as genetic algorithms (GA) [15,16] and ensemble Kalman filters [17,18]. Xie et al. [19] used the FMM to apply GAs to the history matching for shale gas reservoirs, and Leem et al. [20] implemented the FMM to characterize the shale gas reservoirs for which the ensemble Kalman filter had been used. The disadvantage of the use of single-objective optimization algorithms to match the production history is the sensitivity of the optimal solution to the weight factors [21]. Alternatively, for the multi-objective optimization algorithms, several objective functions are covered, and the dynamic time series is integrated. Further, such multi-objective optimization facilitates the creation of sets of optimal solutions, i.e., pareto solutions, instead of merely a single solution [22]. In terms of the reservoir characterization, a number of multi-objective history matching techniques have been applied to quantify the production uncertainties [23,24,25,26]. The representative advantages of this multi-objective history matching are the consideration of scale-different series and the matching of a few dynamic data, whereas the disadvantage is the large computational time required to create the pareto sets. To reduce the cost of the multi-objective optimization process, it is unnecessary to adopt a particular proxy within the model selection framework [27], since any proxy can be used.

For this study, the workflow is developed to reduce the computational effort of the typical full-size numerical simulation within the acceptable error ranges. To accomplish the computational efficiency with acceptable errors, a FMM-based proxy modeling is applied, and its performance is examined in terms of a comparison with the results of a full reservoir simulation. Distance-based generalized sensitivity analysis (DGSA) is implemented to quantify the uncertain parameters for the multi-objective history matching. The multi-objective history matching then covers the available performance metrics at the shale gas reservoir, such as the daily gas rate and the cumulative volume until a given date.

2. Methodology

Figure 1 explains the proposed workflow for the multi-objective history matching for which FMM-based proxy modeling is applied. To efficiently converge to the pareto solutions, the significant parameters that affect the objective functions are used, the parameters of which are determined by the DGSA. Meanwhile, the GA-based multi-objective history matching integrates the well-based production profiles. The theoretical background for each method is as follows.

2.1. FMM for the Proxy Modeling

The FMM evaluates the drainage volume and the gas-production rate by tracking the interfaced propagation and computing the DTOF. The governing FMM equation is the generalized Eikonal equation (Equation (1)) [11,12,13], as follows:

$$F(\overrightarrow{x})\left|\nabla \tau (\overrightarrow{x})\right|=1$$

(1)

$$\mathrm{Initial}\mathrm{condition}:{\tau |}_{F(\overrightarrow{x})=0}=0$$

where $F(\overrightarrow{x})$ is the velocity at each measuring point $\overrightarrow{x}$, and $\tau (\overrightarrow{x})$ is the DTOF. The DTOF ($\tau (\overrightarrow{x})$) is computed from the hydraulic diffusivity that is assigned to each grid, and a smaller value results in a longer DTOF. $F(\overrightarrow{x})$ is always greater than or equal to zero, whereby the pressure is transmitted along only one direction. The pressure diffusion is represented in the Eikonal-equation form, as shown in Equation (2). The hydraulic diffusivity $\alpha (\overrightarrow{x})$ can be determined from the permeability ($k$), porosity ($\phi$), fluid viscosity ($\mu$), and total compressibility ($c_t$), as shown in Equation (3), as follows:

$$\sqrt{\alpha (\overrightarrow{x})}\left|\nabla \tau (\overrightarrow{x})\right|=1$$

(2)

$$\alpha (\overrightarrow{x})=\frac{k}{\phi \mu c_t}$$

(3)

Once the pressure front arrives at the specific grid, grid draining is assumed to commence. The drainage volume at the time, $t$, $V_p(t)$ is the sum of the individual pore volumes within that time contour (Equation (4)), as follows:

$$V_p(t)={\displaystyle {\sum }_{t_i=0}^t\mathrm{Pore}\mathrm{volumes}\mathrm{at}\mathrm{time},t_i}$$

(4)

Equation (5) represents Darcy’s law for the radial flow.

$$q(r,t)=\frac{kA(r)}{\mu }\frac{\partial p(r,t)}{\partial r}$$

(5)

$A=4\pi r^2$ is the surface area for spherical radial flow, $A=2\pi rh$ for cylindrical radial flow, and $A$ = constant for linear flow [12]. By setting the drainage volume as a variable and applying the chain rule, Equation (5) can be expressed in Equation (6), as follows:

$$q(V_p,t)=\frac{k\phi A^2}{\mu }\frac{\partial p(V_p,t)}{\partial V_p}$$

(6)

Now, the Darcy flux ($q(V_p,t)$) in Equation (6) can be approximated according to the multiplication of the well-production rate ($q_{well}$) by the dimensionless flux ($q_D$) along the drainage volume, as given by Equation (7). Equation (8) is derived from the substitution of Equation (7) into Equation (6). The integration of Equation (8) under a constant decrease of the well-flowing pressure ($p_{wf}$) from the initial reservoir pressure ($p_i$) derives Equation (9), as follows:

$$q(V_p,t)\cong q_{well}(t)\cdot q_D(V_p,t)$$

(7)

$$q_{well}(t)\cdot q_D(V_p,t)\cong \frac{k\phi A^2}{\mu }\frac{\partial p(V_p,t)}{\partial V_p}$$

(8)

$$\Delta p=p_i-p_{wf}\cong q_{well}(t)\mu {\displaystyle {\int }_0^{V_p(t)}\frac{q_D(V_p,t)}{k\phi A^2}d{V}_p}$$

(9)

The dimensionless flux $q_D(V_p,t)$ can be formulated using Equation (10), as follows:

$$q_D(V_p,t)=1-\frac{V_p}{V_p(t)}$$

(10)

Equation (11) is the surface-based production rate ($q_{well}$) that is obtained from the substitution of Equation (10) into Equation (9). $B_g$ is the gas-formation-volume factor for the conversion of the reservoir-into-well flow rate into the surface-produced gas rate, as follows:

$$q_{well}(t)\cong \frac{\Delta p}{B_g\mu }\frac{1}{{\displaystyle {\int }_0^{V_p(t)}\frac{1}{k\phi A^2}(1-\frac{V_p}{V_p(t)})d{V}_p}}$$

(11)

2.2. DGSA for the Determination of the Effective Parameters

The DGSA classifies the response variables into a limited set of discrete classes [28]. Unlike conventional sensitivity methodologies, it can be used to explain the effects of the spatiotemporal parameters in the implementation of statistical distributions [29,30]. The fundamental workflow concept is as follows: if the result is sensitive to the parameter, the frequency distribution of the parameter is different in each class, whereas the frequency distribution of the parameter is the same in each class for the parameter-insensitive result.

Figure 2 explains the DGSA process: Figure 2a depicts a distance map and the k-medoid clustering consisting of three clusters, and Figure 2b shows the cumulative-density functions (CDFs) for each cluster. The sensitivity of each parameter is defined as the CDF distance that represents the difference area between the base CDF and each class (Equation (12)), as follows:

$$d_i(k)=A(F_i,F_i(k)),k=1,\dots ,K.$$

(12)

where $d_i(k)$ is the CDF distance of the parameter $i$ for the cluster $k$, $A$ is the function of the numerical calculation of the inter-curve area, $F_i$ is the distribution function of parameter $i$, $F_i(k)$ is the distribution function of parameter $i$ for cluster $k$, and $K$ is the number of clusters. If one or more CDF distances are greater than the given significant level, the parameter can be considered as a sensitive parameter [28]. Figure 3 shows the CDF distances for three parameters, as an example, wherein parameters $i$ and $j$ are sensitive.

2.3. NSGA-II for the Multi-Objective Evolutionary Algorithm

Srinivas and Deb [31] introduced the crowding distance-based non-dominated sorting genetic algorithm (NSGA) to satisfy the multiple constraints, and this was followed by the development of NSGA-II, an improved version of the NSGA, by Deb et al. [22]. The two-principle NSGA-II workflows are as follows: non-dominated sorting and crowding-distance sorting, the latter of which is for diversity preservation (see Figure 4).

At generation $t$, the parent population $P_t$ is used to create the offspring population $Q_t$ for which the genetic operator (selection, crossover, and mutation) is used. The combination of the two populations $R_t$ is classified using non-dominated sorting. The non-dominated sorting ranks the responses in the pareto front. Accordingly, the response $r_1$ dominates the response $r_2$, but only if Equation (13) is satisfied [25], as follows:

$$\forall i\in \{1,\dots ,N\}:r_1^i\le r_2^i\wedge \exists i\in \{1,\dots ,N\}:r_1^i<r_2^i$$

(13)

where $r_j^i$ is the $j^{th}$ response in the direction of the $i^{th}$ objective function, and $N$ is the number of objective functions. Figure 5 depicts two sorting methods in NSGA-II; Figure 5a shows non-dominated sorting and Figure 5b shows crowding-distance sorting. Figure 5a has four non-dominated fronts, and the response here in the higher rank front of $r_1$ is superior to the responses in the lower rank fronts of $r_2$, $r_3$, $r_4$, and $r_5$, but superiority is not evident among the responses in the same rank fronts of $r_3$ and $r_4$.

The likelihood of the selection of the responses in the rank front for the next generation is higher. If the number of the responses from $F_1$ to $F_3$ is greater than $N$, some of the responses from $F_3$ should be rejected by the crowding-distance sorting for the diversity preservation. The crowding distance of the $j^{th}$ response $C_j$ can be defined using Equation (14) [24], as follows:

$$C_j={\displaystyle {\sum }_{i=1}^N\frac{d_j^i}{r_{max}^i-r_{min}^i}}$$

(14)

where $N$ is the number of objective functions, $d_j^i$ is the displacement between two neighboring points with the $j^{th}$ response in the direction of the $i^{th}$ objective function, $r_{max}^i$ is the maximum response in the direction of the $i^{th}$ objective function, and $r_{min}^i$ is the minimum response in the direction of the $i^{th}$ objective function (see Figure 5b). The responses with the lower crowding distance are rejected for diversity preservation, and upon the completion of the crowding-distance sorting, the surviving responses are selected for the next generation $P_{t+1}$.

3. Results and Discussion

3.1. Performance Test for the Validation of the FMM-Based Proxy Modeling

3.1.1. Description of a Synthetic Reservoir and the Testing Method

The illustration in Figure 6 is a synthetic 2D reservoir having been subject to hydraulic fracturing. The reservoir comprises five artificial fractures of different fracture lengths and enhanced permeable zones, assuming a horizontal well and the bi-wing fracture model.

Table 1. Input properties of the synthetic reservoir to validate the FMM-based proxy modeling.

Input Parameters (Unit)	Fixed Value	Uncertain Value †
Reservoir size (x, y, z) (m)	(1200, 600, 165)	–
Grid size (Δx, Δy, Δz) (m)	(3, 3, 3)	–
Fracture half-length (m)	120, 90, 135, 120, 105 (from the left-hand side in Figure 6)	–
Fracture permeability (millidarcy, md)	10	–
Gas viscosity (cp)	0.015	–
Gas saturation (fraction)	0.66	–
Gas formation volume factor (rm3/sm3)	0.012	–
Initial reservoir pressure (MPa)	34.47	–
Well-flowing pressure (MPa)	6.89	–
Reservoir temperature (°C)	104.44	–
Total compressibility (1/kPa)	0.00003	–
Porosity (fraction)	–	Triangular (0.06, 0.08, 0.10)
Matrix permeability (md)	–	Triangular (0.0001, 0.0007, 0.0013)
Enhanced permeability in the stimulated reservoir volume (md)	–	Triangular (0.005, 0.007, 0.009)

^† The reservoir properties are uncertain and statistically selected within the given range. The triangular distribution, Triangular (minimum, mode, maximum), is assumed for the statistical distribution.

summarizes the input parameters and the possible ranges of the uncertain properties. Here, the porosity, matrix permeability, and enhanced permeability in the stimulated reservoir volume, all of which are presented in terms of hydraulic fracturing, are assumed to be the uncertain properties. The reservoir properties have been taken from the major shale gas plays in the USA [32].

The reliability of the FMM-based proxy modeling is investigated using ECL100 (Schlumberger, Houston, USA), a commercial black oil simulator developed by Schlumberger, and its results are compared to those of the FMM-based proxy modeling. For the comparison, the gas production period of the synthetic reservoir is assumed to be 10 years, and 100 realizations are constructed by changing the uncertain parameters, as shown in Table 1. A Monte Carlo simulation is then implemented to select the values of the uncertain parameters for each simulation. The errors in the daily gas-production rates and the cumulative production volume are compared from the results of the full-sized numerical simulation (ECL100) and the FMM-based proxy modeling.

3.1.2. Results of the FMM-Based Proxy Modeling

The performance of the FMM-based proxy modeling is very similar to that of the full-sized simulator. Figure 7 shows a comparison of the production performance for the P50 simulated using FMM-based proxy modeling and the ECL100, as follows: Figure 7a,b shows the daily gas rate and the cumulative volume, respectively. Notably, the production profiles are so similar that FMM-based proxy modeling is a feasible substitute for ECL100. That is, the overall average errors in Figure 7a,b are approximately 0.97% and 0.4%, respectively.

Figure 8 plots the average error for the given time according to 100 realizations, as follows: Figure 8a shows the average error of the daily gas rate between the FMM-based proxy model and ECL100 while Figure 8b shows the average error of the cumulative volume at a given time. In the profile of the daily gas rate, the maximum error observed is approximately 5.8% at the early time, and it is approximately 7.3% for the cumulative trajectories at the same time. The errors decrease with time, and that of the estimated ultimate recovery (EUR) calculated at the end of the gas production is approximately 2%. The above observations are a result of the different amounts of gas released from the matrix into the fracture-governed area that is enhanced by hydraulic fracturing. Furthermore, this is influenced by using the single porosity/single permeability (SPSP) concept for the FMM and considering the dual porosity/dual permeability (DPDP) model for the ECL100. At the early stage, the quantity of drained gas of the DPDP model is larger than that of the SPSP model, and therefore the model selection is influential in terms of the observed differences. Another possible reason is the expected transient-flow governance of the early performance that prevents the proxy modeling from demonstrating this kind of rate-change and making more errors than the ECL100.

The total computation time for ECL100 over 100 realizations is 407.3 min, but that of the FMM-based proxy modeling is 35.5 min. That is, FMM-based proxy modeling is 11 times faster than ECL100. These results confirm the accuracy of the FMM-based proxy model and its sufficient speed for surrogacy to the full-size reservoir simulation, as well as the acceptable accuracy of the ECL100.

3.2. Application: Multi-Objective History Matching with FMM-Based Proxy Modeling

3.2.1. Description of the Gas-Production Data and Assumptions

Figure 9 depicts the daily gas rates at the shale gas field, for which four artificial fractures are known to be generated, and the gas production has been continued for five years without any additional treatment. The daily gas rate starts at approximately 937 Mscf/day, followed by a reduction of approximately 150 Mscf/day at the five-year production stage. The production profile shows the typical trajectories of the shale gas reservoir. During the first stage, it produces the maximum rates, but this rate decreases significantly due to a reduction in the positive effects of hydraulic fracturing.

The objective function for the match production histories consists of the two different production profiles that are defined in Equation (15), with one related to the daily gas rate (refer to Equation (16)) and the other for the cumulative-production volume (see Equation (17)).

$$f=\left(\begin{array}{c}{f}_1\\ f_2\end{array}\right),$$

(15)

$$f_1=\frac{{\displaystyle \sum _{k=1}^n\frac{\left|P_m(k)-P_f(k)\right|}{P_f(k)}}}{n},$$

(16)

$$f_2=\frac{{\displaystyle \sum _{k=1}^n\frac{\left|C_m(k)-C_f(k)\right|}{C_f(k)}}}{n}.$$

(17)

In Equation (16), $P_f(k)$ is the production rate that is calculated by the FMM at day $k$, $P_m(k)$ is the production data at day $k$, and $n$ is the final calculation day. In Equation (17), $C_f(k)$ is the cumulative-production volume that is calculated by the FMM at day $k$, and $C_m(k)$ is the cumulative-production data at day $k$. The qualitative definition of $f_1$ indicates the absolute relative errors of the daily gas rates while that of $f_2$ indicates the absolute relative errors of the cumulative volume.

History matching is carried out using 1.5-year (540-day) data, and the forecast of the production profiles is conducted from 1.5 years to five years (1800 days). The prediction accuracy is evaluated to examine the applicability of the multi-objective history matching in terms of the FMM-based proxy modeling and the DGSA.

3.2.2. Sensitivity Analysis Based on the DGSA

The sensitivity analysis ascertains the parameters that influence the objective functions and quantifies their uncertainty.

Table 2. Input properties for FMM-based modeling and sensitivity analysis.

Input Parameters (Unit)	Fixed Value	Uncertain Value
Reservoir size (x, y, z) (m)	(375, 300, 135)	–
Grid size (Δx, Δy, Δz) (m)	(3, 3, 3)	–
Gas viscosity (cp)	0.02	–
Gas saturation (fraction)	0.6	–
Gas formation volume factor (rm3/sm3)	0.012	–
Initial reservoir pressure (MPa)	10.34	–
Well flowing pressure (MPa)	3.44	–
Reservoir temperature (°C)	37.7	–
Total compressibility (1/kPa)	0.00003	–
Porosity (fraction)		0.005–0.08
Matrix permeability (md)		0.0001–0.0003
Enhanced permeability in the stimulated reservoir volume (md)	–	0.0005–0.005
Fracture permeability (md)	–	0.01–0.1
Fracture half-length (m)	–	15–60
Horizontal enhanced ratio †	–	0.2–0.4
Vertical enhanced ratio ‡	–	0.8–1.6

^† $Horizontalenhancedratio=\frac{horizontallengthofenhancedzone}{fracturelength}$; ^‡ $Verticalenhancedratio=\frac{verticallengthofenhancedzone}{fracturelength}$

includes a summary of the proxy-modeling input properties, with seven parameters that vary while the others remain fixed. The following seven uncertain properties are assumed here: matrix permeability, enhanced permeability in the stimulated reservoir volume, fracture permeability, horizontal enhanced ratio, vertical enhanced ratio, fracture half-length, and porosity. These properties are randomly selected for 200 runs.

The DGSA determines the parameters that influence the objective functions, i.e., $f_1$ and $f_2$ (Equations (16) and (17)). Figure 10a depicts the effects of the uncertain parameters on the objective function $f_1$, and Figure 10b shows the effects for the objective function $f_2$. The enhanced permeability (ke), fracture permeability (kf), and fracture half-length (xf) are the three parameters that affect both objective functions. While the porosity (po) can influence $f_1$, it is less significant for $f_2$.

In brief, four of the parameters, i.e., enhanced permeability, fracture permeability, fracture half-length, and porosity, are selected to change the values of the multi-objective history matching. Although it is known that porosity is not influential regarding $f_1$, it is still included for the history matching to extend the search space.

3.2.3. Multi-Objective History Matching

To examine the applicability of the developed workflow, the following three GA-based single objective optimization comparisons are introduced: $f_1$ (Equation (16)), $f_2$ (Equation (17)), and their arithmetic mean $f_3$ (see Equation (18)). The objective function $f$ for the GA-based multi-objective optimization is defined using Equation (15).

$$f_3=\frac{f_1+f_2}{2}$$

(18)

The FMM-based proxy modeling displaces the reservoir simulation to obtain the production profiles. The changing variables that are used for history matching are determined by the DGSA, and their ranges are defined in Table 2. The other variables, which are not significant regarding the objective functions, were fixed according to their mean values, as follows: the matrix permeability is 0.00002 md, the horizontal enhanced ratio is 0.3, and the vertical enhanced ratio is 1.2. History matching is implemented for a total duration of 540 days (1.5 years) according to the objective functions defined by Equations (15)–(18). The prediction accuracy from 540 days to 1800 days is compared with the true data, as depicted in Figure 9.

The model used for comparison with a result from a single objective function produced only one optimal trajectory while the proposed multi-objective history matching identified four solutions. Figure 11 depicts the errors of the optimal solutions with Figure 11a showing the errors between the optimal solution and the true data during the history-matching period (over 540 days), and Figure 11b describing those that occurred during the prediction period from 540 days to 1800 days. The proposed multi-objective history matching shows a lesser number of errors for both the history matching and the prediction compared to the single-objective optimization. The single-objective optimization shows a scale-dependency problem whereby it leans toward its allocated objective function [23,24,25]. The case of the arithmetic averaged error, $f_3$ is located at the middle point between $f_1$ and $f_2$, as shown in Figure 11a. The four determined solutions of the proposed method are non-dominated and generate the pareto front, whereby the uncertainty can be considered. Figure 11b reveals the importance of considering the uncertainty during history matching. Even though Figure 11a shows a weak difference between the single-objective and the multi-objective optimizations, the forecast performance is still differentiated. The overall errors of the proposed multi-objective history matching are located in an area with smaller errors compared to the single-objective case.

Figure 12 compares the production profiles of the developed workflow with four optimal solutions to that of true data. All obtained solutions are not only compatible, but also reliably predict the future performance. Figure 13 presents box plots of the error data that are calculated at 20 time steps allocated to each optimum solution.

Table 3. The results of the history matching and the prediction for the three comparisons and the proposed model. EUR: estimated ultimate recovery.

Objective Function	During the History-Matching Period		During the Prediction Period		Error of EUR (%)
	Error of Daily Gas Rate (%)	Error of Cumulative Volume (%)	Error of Daily Gas Rate (%)	Error of Cumulative Volume (%)
$f_1$	1.75	2.12	4.83	1.83	2.87
$f_2$	2.58	1.04	4.94	1.15	2.48
$f_3$	3.00	1.88	2.34	1.40	1.35
$f$ (mean value)	2.08	2.34	2.42	0.45	0.43

summarizes the averaged errors of the three comparisons and the proposed method with the true profile. The averaged errors of the multi-objective history matching are less than those of the single-objective comparisons. The optimal solutions that are obtained using the proposed workflow show production-profile errors of around 2%, and this is irrespective of the use of the FMM-based proxy model instead of the time-consuming reservoir simulator for the proposed method. Principally, the four optimum solutions show an error of only 0.43% in terms of the EUR, but the results of the three comparisons show errors of around 2%.

The developed multi-objective history matching is based on the FMM proxy model, and the DGSA forecasts the production performance of the shale gas reservoir by integrating the well-based production performance. The optimal solutions on the pareto front show lesser errors compared to the single-objective optimization, and this is without a scale-dependency problem. This study assumes homogeneity of the rock properties, including permeability, porosity, and initial gas saturation, but to improve the applicability of the proposed method in the field, the reservoir heterogeneity estimated by microseismic and well testing data need to be considered. Validation of the reliable estimation of geomechanical characteristics and a quantitative analysis of the reservoir heterogeneity are still a challenge for hydraulic-fracturing techniques.

4. Conclusions

The authors of this paper developed multi-objective history matching for which the FMM-based proxy model and the DGSA were applied. The FMM proxy model could serve as a surrogate for reservoir simulation by reducing the computation time within the acceptable errors, while the DGSA determined the parameters that have influence in cases with single-objective functions. In this study, the effective parameters were four properties related to hydraulic fracturing: enhanced permeability, fracture permeability, fracture half-length, and porosity. The properties of the stimulated reservoir volume that affected the gas production were, in sequence, the permeability of the enhanced zone, permeability of the hydraulic fractures and the fracture half-length. The multi-objective history matching for which both the daily gas rate and its cumulative volume were considered determined the optimal prediction of the unknown performances without the scale-dependency problem; furthermore, it showed more effective applicability compared with the single-objective optimization through its consideration of the uncertainty of the optimal solutions.

The method proposed in this paper addresses the time-consuming problems of reservoir simulation as well as the scale-dependency of single-objective optimization, whereby the EUR and the gas rates could be predicted in a more-realistic manner. The developed scheme could therefore be applied to evaluate the production performances of a shale gas reservoir and the effects of hydraulic fracturing.