Interpretation of Gas/Water Relative Permeability of Coal Using the Hybrid Bayesian-Assisted History Matching: New Insights

Abstract

The relative permeability of coal to gas and water exerts a profound influence on fluid transport in coal seams in both primary and enhanced coalbed methane (ECBM) recovery processes where multiphase flow occurs. Unsteady-state core-flooding tests interpreted by the Johnson–Bossler–Naumann (JBN) method are commonly used to obtain the relative permeability of coal. However, the JBN method fails to capture multiple gas–water–coal interaction mechanisms, which inevitably results in inaccurate estimations of relative permeability. This paper proposes an improved assisted history matching framework using the Bayesian adaptive direct search (BADS) algorithm to interpret the relative permeability of coal from unsteady-state flooding test data. The validation results show that the BADS algorithm is significantly faster than previous algorithms in terms of convergence speed. The proposed method can accurately reproduce the true relative permeability curves without a presumption of the endpoint saturations given a small end-effect number of <0.56. As a comparison, the routine JBN method produces abnormal interpretation results (with the estimated connate water saturation ≈33% higher than and the endpoint water/gas relative permeability only ≈0.02 of the true value) under comparable conditions. The proposed framework is a promising computationally effective alternative to the JBN method to accurately derive relative permeability relations for gas–water–coal systems with multiple fluid–rock interaction mechanisms.

1. Introduction

The recovery of coalbed methane (CBM) from underground coal seams has multiple benefits. These include the reduction of greenhouse gases released into the atmosphere and the enhancement of underground coal mining safety and in augmenting the supply of natural gas [1]. Most coal seams are initially saturated with water at in situ conditions [2]. Hence, two-phase gas/water flow is frequently encountered in the primary recovery of coalbed methane (CBM) by either routine dewatering or enhanced coalbed methane (ECBM) strategies such as the injection of N₂ and/or CO₂. It is well recognized that the relative permeability of coal exerts a profound effect on multiphase fluid transport behavior [3,4] and controls gas/water production rates [5,6], the cost of produced water disposal [7], and N₂/CO₂ injectivity [8,9].

Experimental methods to obtain the relative permeability curves of coal typically include the unsteady-state [10] core flooding, steady-state [11] core flooding, and microfluidic flow [12] tests. The steady-state method is reputed to produce more accurate and reliable results than the unsteady-state method [13]; however, the former is rarely used for estimating gas/water relative permeability of coal due to its cumbersome nature and experimental difficulties [14]. So far, only Gash [11] and Reznik et al. [13] reported the use of steady-state tests to obtain the gas and water relative permeability of coal. The microfluidic method is capable of visualizing the fluid flow through coal cleats, which however requires an expensive experimental apparatus and complicated operational procedures [15]. As a comparison, the unsteady-state gas-displacing-water flooding test is much easier and more rapid, which is widely used for determining relative permeability of coals [16].

Although the unsteady-state flooding test demonstrates a remarkable superiority over other methods in terms of experimental efficiency and cost-effectiveness, the interpretation of relative permeability from the unsteady displacement data (i.e., pressures and fluid productions) is much more complex [17,18]. To date, various interpretation methods have been proposed to derive relative permeability from unsteady-state experiments, which can be classified into analytical and history matching methods [19]. The analytical methods use relatively simple and straightforward computational schemes and are easy to implement, which therefore are currently most routinely used. Among the analytical methods, the Johnson–Bossler–Neumann (JBN) method [20] is most commonly used in the CBM research community. However, the JBN method was initially developed for water/oil flows in sandstones, and its accuracy remains questionable when applied in a gas–water–coal system that involves more complex mechanisms. The JBN model was initially developed under three fundamental assumptions that (i) capillary end-effects are negligible, (ii) the density of each phase remains constant (i.e., neither phase is compressible), and (iii) piston-like displacement occurs throughout the core sample. The first assumption can be satisfied in reality by using high-speed displacement rates to diminish the relative effect of capillary pressure. However, the latter two assumptions are difficult to accommodate for a gas/water system because (i) gas is extremely compressible with strongly varying density in the displacement process, especially for cases where the pressure gradient along the core sample is high, and (ii) gas and water have a distinctive contrast in viscosities, which may result in viscous fingering during the displacement process that fails the assumption of a piston-like displacement [21,22,23]. Moreover, when an adsorbing gas species (such as methane, nitrogen, and CO₂ that are of interest to CBM recovery processes) is injected into a coal sample, the adsorption of the gas typically results in changes in the coal permeability, which inevitably affects the gas/water production curves that are used to calculate relative permeability. These factors cannot be captured by the JBN method.

Compared with the JBN method, the assisted history matching method, based on numerical simulations, is capable of tracking gas compressibility, gas adsorption, and coal permeability evolution provided that these mechanisms are incorporated in the numerical simulator. The assisted history matching method in essence solves a mathematical optimization problem with an objective function of minimizing the errors between the simulated and experimental production curves by tuning the relative permeability curves with an optimization algorithm(s). The assisted history matching method has been successfully used to derive water/oil relative permeability curves from core-flooding tests [24,25,26]. Shaw et al. [27] used manual history matching to simultaneously derive relative permeability and capillary pressure for gas–water–coal systems. However, to the best knowledge of the authors, no publication has ever reported the application of assisted history matching in deriving the relative permeability (RP) from core-flooding tests for a gas/liquid system.

To date, various algorithms have been applied to assist in history matching, including (i) gradient-based methods, e.g., the Levenberg–Marquardt (LM) [25] and the quasi-Newton [28] methods; (ii) ensemble-based methods, e.g., the ensemble Kalman filter (EnKF) [24] and ensemble smoother with multiple data assimilations (ES-MDA) [29] methods; and (iii) swarm intelligence algorithms [30]. It should be noted that, for these methods, the saturation endpoints for both phases should be assumed to be known prior to conducting history matching, which, if not assumed correctly, may produce unreliable results [31]. For a real gas-displacing-water test, the saturation endpoint for water can be readily determined based on the sample weights before and after gas flooding. However, determination of the saturation endpoint for gas can be extremely difficult without expensive online X-ray computed tomography (X-ray CT) [32] or nuclear magnetic resonance (NMR) imaging [33]. This may significantly hinder the application of assisted history matching in a gas/liquid system. It should also be noted that gradient-based methods are prone to being trapped in local minima if not properly initialized [34]. Moreover, for the ensemble and swarm intelligence methods, a large number of numerical simulation runs are needed in order to achieve a satisfactory result. This makes the estimation of relative permeability a low-efficiency and computationally expensive process.

In summary of previous studies, the JBN method is most frequently used for interpreting relative permeability from unsteady-state displacement experiments, which however is highly questionable with application to the gas/water/coal system due to its inherent theoretical deficiencies. As a comparison, the assisted history matching method is more reliable provided that complex fluid–rock interaction mechanisms are properly accounted for. However, most of the existing assisted history matching methods require a prior assumption or determination of the endpoint saturations, which is difficult to achieve for systems involving compressible phase(s). In addition, existing algorithms are associated with the issue of high computational overhead and in being time-consuming. Addressing these issues will improve the utility of assisted history matching methods in estimating gas/water relative permeability curves in coals—enabling them to be used with both high accuracy and efficiency. This study presents an improved method using a hybrid Bayesian optimization (BO) and mesh adaptive direct search (MADS) algorithm to automatically history-match production profiles to determine relative permeabilities of coals from unsteady-state core-flooding experiments. The improved accuracy of the proposed method is demonstrated on synthetic numerical core-flooding experiments. In addition, the superiority of the method is highlighted by comparison with the routine application of the JBN method.

2. Methods

2.1. Hybridizing BO with MADS

The MADS algorithm is a local optimizer that has been hybridized with global optimization algorithms to solve varying optimization problems related to the petroleum engineering (PE) community [35,36,37]. Although the application of BO within the PE community is only sparingly reported, it has been actively studied in other areas involving expensive cost functions such as in machine learning [38].

2.1.1. Basics of BO

The BO is an efficient framework for solving black box optimization problems that involve objective functions that are computationally expensive to evaluate. A BO typically consists of two consecutive steps (Figure 1). First, a prior over function is built to represent the belief about the objective function. When observations of the objective function are obtained, the posterior distribution over functions can be updated by combining the prior distribution with the observations, which forms a surrogate model on the objective function. Because of its flexibility and tractability, Gaussian process (GP) priors are frequently utilized to represent the assumptive prior distribution on functions [39]. Second, an acquisition function is applied to construct a utility function from the model posterior, which is then maximized to select the next query point at which to evaluate the function. After new function evaluations are conducted, the updated observations are combined with the previous posterior estimation to construct a new posterior distribution over functions. This iterative search process is repeated until an optimal value is obtained.

The covariance and acquisition functions are two key elements that may significantly influence the performance of the Gaussian process. Stationary covariance functions include the automatic relevance determination (ARD), squared exponential (SE), rational quadratic (RQ), Matérn 3/2, and Matérn 5/2 kernels. The acquisition functions that are commonly used include the lower confidence bound (LCB), upper confidence bound (UCB), and expected improvement (EI). In this study, we adopt the recommendation of [40] for the stationary covariance and acquisition functions, which are set to be the ARD RQ and the LCB, respectively.

2.1.2. Basics of MADS

The MADS algorithm is a directional direct search method proposed by [41]. At each iteration, the MADS updates the candidate solutions in two sequential steps, namely the search stage and poll stage (Figure 2). In the search stage, MADS can use strategies such as heuristics and surrogate models (e.g., GPs in the BO process) to generate a finite number of trial points on the mesh for following functional evaluations [42]. The use of a surrogate in this step is generally less costly and can save considerable computational expenditure [43]. The poll stage is performed if the search stage fails to improve the objective function value. It has been theoretically proven that the poll stage ensures the convergence of MADS to local optima for non-smooth functions [44]. If either the search or poll succeeds in finding a mesh point with an improved objective value, the incumbent is updated and the mesh size remains the same or is multiplied by a factor τ > 1. If neither search nor poll are successful, the incumbent does not move and the mesh size is divided by τ. The algorithm proceeds until a stopping criterion is met (e.g., when a budget of evaluations has been exhausted).

2.1.3. Hybridization

The philosophy behind the Bayesian adaptive direct search (BADS) is to utilize the surrogate models (GPs) constructed in the BO process to assist in generating candidate solutions in the search and poll stages of the MADS. Mathematical details on the BADS are given in [40], and the workflow is illustrated in Figure 3, which is briefly introduced as follows for completeness of this paper.

For a D_n-dimension minimization problem ${\mathit{x}}^{*}=\arg {\min }_{x\in {ℝ}^{D_n}}f\left(\mathit{x}\right)$, the BADS method begins function evaluations with an initial variable vector ${\mathit{x}}_0$ and a number of D_n additional points chosen via a space-filling quasi-random Sobol sequence [45] that are forced to be on the mesh grid (Figure 3). The initial mesh size (${\Delta }_0^{mesh}$) and poll size (${\Delta }_0^{poll}$) are set to 2⁻¹⁰ and 1.0 according to [40], respectively. The GP surrogate model is updated every 2D_n to 5D_n function evaluations and whenever the accuracy of the current GP is unreliable by refitting the hyperparameters of GP using the gradient-based fmincon optimizer implemented in the Optimization Toolbox of MATLAB.

At the kth iteration, the candidate solutions in the search stage ${\mathit{x}}_{search}$ are generated by performing a local optimization of the acquisition function in the neighborhood of the incumbent ${\mathit{x}}_k$ using the two-step evolutionary strategy adapted from the covariance matrix adaptation evolutionary strategies (CMA-ES) -based MADS (which is also termed as LTMADS) algorithm [41] and are then scaled to the GP kernel length scales to construct the final polling directions D_k. In this way, the solution candidates ${\mathit{x}}_{poll}$ in the polling stage can be updated according to ${\mathit{x}}_{poll}=\left\{{\mathit{x}}_k+{\Delta }_k^m\upsilon :\upsilon \in {\mathit{D}}_k\right\}$. The polling sets are evaluated according to their ranking given by the acquisition function. If the poll succeeds in sufficiently improving the objective function within three consecutive steps, the incumbent is updated and BADS switches to a new iteration with mesh and poll sizes multiplied by τ = 2; otherwise, the incumbent remains unchanged and BADS switches to a new iteration with the mesh and poll sizes divided by τ = 2. These steps are repeated until a preset maximum number of iterations is met, the algorithm stalls, or the poll size becomes extremely small.

2.2. Representation of Relative Permeability and Capillary Pressure

Power-law models are frequently used to represent relative permeability in history matches due to the small number of tunable parameters. Although Shawn et al. [27] suggest that the simple modified Brooks–Corey-type model [46] is adequate in representing gas and water relative permeability curves of coals, we adopt Chen et al.’s model [2] in this study, which is more versatile than the Brooks–Corey model. Chen et al.’s relative permeability model is developed specifically for coals and incorporates both the match-stick geometry of cleats and its effect on the computation of the capillary pressure. The applicability and accuracy of this model has been confirmed against a large number of relative permeability curves for a variety of coals [47,48]. The general form of Chen et al.’s model is written as follows:

$$k_{rw}=k_{rw}^{*}{S}_{wD}^{\eta +1+2/\lambda }$$

(1)

$$k_{rg}=k_{rg}^{*}{\left(1-S_{wD}\right)}^{\eta }\left(1-S_{wD}^{1+2/\lambda }\right)$$

(2)

$$S_{wD}=\frac{S_w-S_{wc}}{1-S_{wc}-S_{gr}}$$

(3)

where $k_{rw}^{*}$ and $k_{rg}^{*}$ are endpoint relative permeabilities for the water and gas phase, respectively; $S_{wc}$ and $S_{gr}$ are connate water and irreducible gas saturations, respectively; and η and λ are positive fitting exponents that represent the tortuosity and cleat size distribution of the cleat system.

For the unsteady-state flooding test, the displacement rate (or pressure gradient applied between core ends) should be properly set such that the reservoir conditions are appropriately replicated while the capillary end-effect can be neglected [49]. However, it is usually difficult to define the proper values without knowledge of the capillary pressure characteristics in real laboratory tests. To the best knowledge of the authors, experimental data on capillary pressure for gas–liquid–coal systems at high pore pressures and under confined conditions are still lacking due to the difficulty in making such measurements [27]. As such, it is of practical significance to estimate relative permeability and capillary pressure simultaneously from an individual unsteady-state displacement test to eliminate the difficulty of capillary pressure measurement. To represent the capillary pressure curve, which may be simultaneously interpreted with that for relative permeability, the Liu et al. model [50] (developed specifically for coal) was used to represent the capillary pressure curve, which is written as follows:

$$p_c=p_{\mathrm{c},\max }{\left[1-b{S}_{wD}^{\gamma }\right]}^{-\frac{1}{\beta }}$$

(4)

where p_c,max is the maximum pressure on the capillary pressure curve and b, γ, and β are the fitting parameters.

2.3. Implementation of the Assisted History Matching Framework

The workflow for the assisted history matching is shown in Figure 4 and is described as follows:

Step 1: Randomly initialize the control parameters (k^*_rw, k^*_rg, S_wc, S_gr, η, λ, b, β, γ, and p_c,max) shown in Equations (1)–(4) within their respective constraint boundaries that are given in [47,50].

Step 2: Calculate the relative permeability and capillary pressure curves according to Equations (1)–(4).

Step 3: Update the data file for numerical simulation.

Step 4: Call the reservoir simulator to conduct numerical simulation for each solution candidate.

In this study, we set the inlet and outlet pressures as constraints and attempt to simultaneously match the cumulative productions of gas and water at outlet. Thus, the objective function to be minimized is written as follows:

$$f=\frac{{{\displaystyle \sum }}_{i=1}^{N_g}{\left(q_{cal,i}^g-q_{obs,i}^g\right)}^2}{N_g}+\frac{{{\displaystyle \sum }}_{i=1}^{N_w}{\left(q_{cal,i}^w-q_{obs,i}^w\right)}^2}{N_w}$$

(5)

where q is cumulative production; N is the number of points on the cumulative production curve; the subscripts “g” and “w” represent gas and water phase, respectively; and the subscript “cal” and “obs” represent calculated and observed values, respectively.

Step 6: Update the optimal solution, and generate new solutions using the BADS algorithm. In this study, the MATLAB implementation of BADS developed by [40] is used, which is available at https://github.com/lacerbi/bads#reference.

Step 7: Check if the maximum number of function evaluations is met. If not, repeat steps 2 through 6; otherwise, terminate history matching, and output the relative permeability and capillary pressure.

3. Validation of the Framework

Compared to real laboratory tests, a most distinguishing advantage of using synthetic numerical core flooding tests is that the true properties of the core samples are definitively known and can therefore be used to test the accuracy of an interpretation technique [51,52]. Relative to this, numerical methane-displacing-water core flooding tests are constructed, conducted, and used to validate the proposed interpretation framework.

Cylindrical core samples are typically used for laboratory flooding tests (Figure 5a). This geometry can be approximated by a rectangular cuboid in numerical simulations provided that their cross-sectional areas are identical [53,54] with no loss in fidelity (Figure 5b). The tubing network of the experimental equipment exerts a finite influence on gas storage and flow behavior and, therefore, should also be accounted for in the grid model [53,55]. In this study, a rectangular grid model with 50 grid blocks along the x-, y-, and z-directions was constructed for simulating unsteady gas-displacing-water tests with capacitance in the tubing represented by single grid blocks both upstream and downstream. The central 50 grid blocks represented the coal core, whereas 1 grid block on either side represented the upper- and downstream tubing voids. Each grid block had a respective dimension of 0.2, 5.321, and 5.321 cm in the x- (axial), y-, and z- (lateral) directions so that a coal core sample with a respective diameter and length of 5 cm and 10 cm was represented.

The coal core was initially saturated with water, with the basic properties summarized in

Table 1. Coal properties used in the synthetic numerical model.

Property	Value
Cleat porosity, %	3.00
Permeability, mD *	0.1
Rock compressibility, KPa−1	2.5 × 10−5
Sorption strain, fraction	0.009
Pε, KPa	2700
Poison ratio, fraction	0.35
Young’s Modulus, KPa	3.0 × 106
Initial pressure, KPa	101.1
Initial water saturation, %	100
Langmuir volume, m3/t	24.0
Langmuir pressure, KPa	2300
Desorption time, d	0.5
Temperature, K	308.15
Bulk density, g/cm3	1.4

* The permeability value is assumed to be with reference to the averaged inlet and outlet pressures.

. The grid blocks on either end of the grid were penetrated with an injection and production well, respectively, to simulate the injection/production process of the fluids. For the two grid blocks at either end, the reservoir properties were set as follows: (i) Langmuir volume and volumetric strain were zero, and compressibility was extremely small since no gas sorption or volumetric variation occurred in the tubing system; (ii) porosity was assumed to be 0.02% to account for the tubing voids; and (iii) permeability was set to be 100D so that the bottomhole pressures (BHPs) approximated the grid block pressures on either end.

For real unsteady-state flooding tests, relative permeability and capillary pressure are inevitably associated with measurement uncertainties and thus may deviate from the ideal modeled curves. In other words, real relative permeability and/or capillary curves can only be represented with analytical models if noise is added. In this regard, the reference relative permeability and capillary pressure curves are generated by adding noise to the data points generated with the Chen et al. [2] and Liu et al. [50] models, respectively, in order to mimic the real cases (Figure 6). Previous studies conclude that the pressure drop (or fluid injection rate) applied on the core sample exerts a nonnegligible effect on the relative permeability due to the presence of capillary end-effect [56,57]. To examine the robustness and accuracy of the proposed framework, different simulation scenarios with varying pressure drops were investigated. The BHP of the injector was set to 4.0 MPa, which generally approximates the methane desorption pressure of the coal seams in the Qinshui basin. The BHPs of the producer were varied at 3.7, 3.75, 3.8, 3.85, 3.9, and 3.95 MPa to represent pressure drops of 0.3, 0.25, 0.2, 0.15, 0.1, and 0.05 MPa, respectively. The corresponding simulated cumulative gas and water production rates are shown in Figure 7. It should be noted that, for real core flooding tests, applying a large pressure drop may result in viscous fingering due to the relatively large viscosity ratio between gas and liquid phases [58] and, therefore, may affect the relative permeability. However, for the numerical core flooding tests, viscous figuring was not accounted for and the possible differences in the interpreted relative permeability may be attributed solely to the effects of capillary pressure, gas compressibility and adsorption, and the dynamic permeability behavior of coal, which will be discussed in the following section.

4. Results and Discussion

4.1. Convergence Trend of the BADS Algorithm

BADS is a stochastic optimization algorithm; therefore, variations may exist with respect to different independent simulation runs. To diminish the uncertainties in interpretation caused by the stochastic nature of BADS, multiple independent runs were conducted for each of the pressure gradient cases. Figure 8 illustrates the convergence trend of the objective function values for different pressure gradient scenarios. As can be seen, for each scenario, the objective function decreases dramatically during the initial stage (with a number of ≈100 to ≈150 function evaluations) and then converges gradually to local optima up to ≈300 numerical simulation runs. This suggests an alternation in mode from global exploration to local exploitation [5]. It should be noted that the number of function evaluations required for the BADS algorithm to achieve relatively stable convergence is noticeably small compared with previous algorithms relative to the specific problem of relative permeability estimation. For example, Zhang et al. [51], Zhou et al. [25], and Fayazi et al. [59] demonstrated that the EnKF, LM, and genetic algorithm (GA) require respectively ≈1000, ≈950, and ≈500 simulation runs to achieve stable convergence. This highlights the efficiency of the BADS algorithm in assisted history matching.

Also shown in Figure 8, for each pressure difference scenario, the ten independent runs exhibit noticeable variations in the final objective function values after even ≈300 simulation runs. The variations are attributed to the stochastic nature of the iterative process in BADS in the generation of solution candidates, which is a common issue for stochastic optimization algorithms [5]. The averaged final objective function values (even after ≈300 numerical simulations) are also affected by the differential pressure. The final averaged objective function values are ≈0.001, 0.005, and 0.45 for scenarios with pressure drops of >0.15 MPa, 0.1 MPa, and 0.05 MPa, respectively. The variations in the final objective function values indicate a noticeable influence of the pressure drop in the interpretation accuracy, which will be discussed in the following section.

4.2. Effect of Pressure Drop on the Interpretation Accuracy

Figure 9 depicts the interpreted relative permeability with reference to the “true” curves. As shown, the interpreted relative permeability curves from different independent runs are scattered around the true curves, which is possibly due to the random initialization and stochastic nature of the BADS algorithm. Nonetheless, the averaged relative permeability of ten independent runs agree well with the true curves at larger pressure drops (0.2, 0.25, and 0.3 MPa) (Figure 9a–c). For a pressure drop of 0.15 MPa, the interpreted gas relative permeability agrees well with the true permeability, whereas small discrepancies exist for the water relative permeability at water saturations of >0.9 (Figure 9d). For the scenario with a pressure drop of 0.1 MPa, the interpreted gas and water relative permeability curves deviate noticeably from the true behavior especially in the region around the endpoint saturations (Figure 9e). For the scenario with an extreme small pressure drop of 0.05 MPa, significant discrepancies exist between the estimated and true relative permeability curves: the connate water saturation is overestimated, whereas the relative permeabilities for both the gas and water phases are generally underestimated compared with the true ones (Figure 9f). The discrepancies between the interpreted and true connate water saturation for the extreme low pressure drop (0.05 MPa) may be attributed to the dominating role of capillary forces relative to viscous forces that results in a high water saturation remaining at the end of the displacement experiment (Figure 7b). The relatively high residual water saturation misleads the BADS algorithm to derive a relatively high connate water saturation. Once the connate water saturation is overestimated, the relative permeability for both the gas and water phases should be forced to reduce in order to improve the estimation accuracy.

To further investigate the effect of pressure drop on the interpretation accuracy, the dimensionless flow parameter of end-effect number (N_c,end) [60] is included. This is defined as follows:

$$N_{\mathrm{c},\mathrm{end}}=\frac{P_{\mathrm{c},\max }}{\Delta p}$$

(6)

where $\Delta p$ is the pressure drop along the core sample.

For the specific synthetic models in this work, the calculated end-effect number values are shown in Figure 10. As addressed previously, the interpretation accuracy is relatively high provided that the applied pressure drop on the core sample is not lower than ≈0.1 MPa—this corresponds to an N_c,end of ≈0.56. Bentsen [61] suggests that the effect of the capillary term on fluid displacement can be neglected given N_c,end < 0.01 and N_c,end < 0.1 for two-phase systems with mobility ratios of 2 and 10, respectively. For the specific conditions of this study (Table 1), the mobility ratio of gas over water is much higher (≈70) than that of [61]. Nonetheless, by analogy with [61], it should be reasonable to assume a higher critical N_c,end value below which the capillary end-effect can be neglected, specifically for the gas/water system in this study.

Numerous studies have documented the use of significantly larger pressure drops or flow rates than those in real reservoirs in order to minimize the capillary end-effect in core flooding tests and thus to improve the interpretation accuracy of relative permeability [56,57]. However, applying a large pressure gradient tends to increase the possibility of viscous figuring considering the necessarily unfavorably high mobility ratio of gas relative to water [21,22,23]. This may produce relative permeability curves that are not representative of the true reservoir conditions. For the methane–water–coal system, applying larger pressure drops also results in more significant variations in the permeability change and in gas adsorption along the core sample. Therefore, it should be of practical significance to minimize the pressure drop applied on the core sample while ensuring the accuracy of interpretation, which can be readily achieved using the interpretation framework proposed in this study.

Chen et al. [31] argued that endpoint saturations should be assumed prior to conducting assisted history matching in order to ensure the accurate interpretation of results. However, for gas and water systems, the endpoint gas saturation is usually difficult to identify without advanced measurement techniques—such as X-ray CT or NMR relaxation imaging [33]. Moreover, an accurate determination of the true endpoint water saturation usually requires extremely high injection rates and long-duration displacement. From the results above, it is apparent that the proposed framework is capable of accurately reproducing the true relative permeability curves without prior assumption of the endpoint saturations, provided that the displacement pressure is properly assigned—this represents a significant advancement over previous methods.

Figure 11 compares the interpreted capillary pressure curves with the true curve at varying differential pressures. As depicted, the interpreted capillary pressure curves deviate noticeably from the true curve, especially where water saturation is <0.7. Previous studies [34,51] indicate that accurate interpretation of the capillary pressure curve requires that not only the injection/production and pressure data but also time-lapse saturation profiles along the core sample be included in the objective function. However, the acquisition of the time-lapse saturation profiles requires expensive and sophisticated experimental techniques such as real-time X-ray CT scanning, which may not be available for all core flooding tests. Dabbous et al. [62] and Reznik et al. [13] demonstrate that the maximum gas/water capillary pressure for bulk coals, at effective stresses of up to 6.9 MPa, are generally less than 0.25 MPa. Such small capillary pressures have been argued to exert a minor or even negligible effect on fluid flows in naturally fractured formations [63]. As such, it is reasonable to neglect the errors in the interpretation of capillary pressure and to place the accurate estimation of relative permeability as a priority.

4.3. Comparison with the JBN Method

Figure 12 depicts the interpretation results of the routine JBN method for the relatively large pressure drop of 0.3 MPa (corresponding to a large end-effect number value of ≈2.2). As clearly shown, the derived relative permeability curves for both water and gas phases using the JBN method deviate significantly from the true reference data. The derived connate water saturation is approximately 0.738, which is 33% higher than the true values (0.553). The derived water relative permeability values are on the order of magnitude of 10⁻⁵ to 10⁻³, with the endpoint water relative permeability being ≈0.01 of the true value. The derived gas relative permeability is of the order of ≈10⁻², with the endpoint gas relative permeability being ≈0.02 of the true value. Compared with the proposed interpretation framework introduced and benchmarked in this study, the results interpreted using the JBN are, in contrast, highly questionable and unreliable. As noted previously, the JBN method was initially developed to represent incompressible two-phase water and oil flow under the assumption of negligible capillary end-effects. For a real methane-displacing-water test on a coal sample, three discrepancies exist relative to the necessary assumptions for JBN: (i) the capillary end-effects are nonnegligible, (ii) gas compressibility and adsorption effects are significant, and (iii) permeability evolves in coal due to gas injection and as a result of both sorptive-swelling and the pressure sensitive nature of the cleat system. These factors are not accounted for in the JBN method, which therefore result in the discrepancies between the interpreted and true relative permeabilities that we observe.

To ascertain the effect of each of these factors on the JBN interpretation results, three simulation model cases are tuned based on the basic model parameterization (Table 1 and Figure 5 and Figure 6), and the JBN method is applied to interpret the relative permeabilities for each individual case (Figure 13). For all these cases, the pressure drop applied on the core sample is set to be relatively large (0.3 MPa).

(1)

A first case involves only gas compressibility and adsorption by assuming zero capillary pressure and by accommodating no evolution in permeability. As can be seen from Figure 13, for this case, the JBN-interpreted relative permeability curves agree well with the true relations, indicating that gas compressibility and adsorption exert minor effects on the accuracy of interpretation using the JBN algorithm.

(2)

A second case is constructed by adding the capillary pressure while still excluding permeability evolution. Apparent in Figure 13 is that the inclusion of the capillary effect tends to lower the interpreted water relative permeability with reference to the true curve. The interpretation of gas phase relative permeability appears to be less affected compared with that for the water phase over the range of high water saturations; however, the interpreted gas relative permeability is lower than the real magnitude when the water saturation approximates the connate water saturation.

(3)

A third case is constructed by adding the permeability evolution effect while still excluding the impact of capillary pressure but otherwise represents the conditions of the first case. As depicted in Figure 13, for this case, the JBN-interpreted relative permeability curves deviate significantly from the true behavior. These results are actually quite similar to the curves interpreted for the basic model (that incorporate capillary end-effects, gas compressibility and adsorption, and permeability evolution).

In summary, the effect of gas compressibility and adsorption can be neglected in an interpretation using the JBN method. Capillary end-effects exert a nonnegligible but significantly smaller effect on the accuracy of interpretation using the JBN algorithm compared to the impacts of permeability evolution. This latter effect of permeability evolution is the primary source of error causing unreliable interpretation of permeabilities from the JBN method. In other words, the JBN method only yields reliable relative permeability estimates where capillary end-effects and permeability evolution (due to swelling and changes in effective stress) are negligible.

However, swelling and stress-dependent permeability evolution are inherent in fractured coals. This effect may however be reduced by using non-adsorbing gases as the displacing fluid. As stated previously, the impact of the capillary end-effect may be reduced by using a large pressure gradient. To further evaluate the flow regime where the JBN method is valid for cases where permeability evolution is negligible, additional simulation runs are conducted for increased pressure drops (0.4 and 0.5 MPa). The interpreted results are shown in Figure 14. As apparent from the figure, the JBN method gives an accurate estimation for the gas relative permeability at pressure drops of 0.4 and 0.5 MPa. For the water phase, the estimated relative permeability curves at pressure drops of 0.4 and 0.5 MPa also agree with the true values for water saturations lower than the cross-point saturation but are slightly lower than the true values at water saturations higher than this. Overall, the accuracy of interpretation of the JBN algorithm is acceptable for pressure drops larger than 0.4 MPa, which correspond to an N_c,end of ≈3.3. By comparison, the proposed framework introduced in this study is capable of accurately reproducing the relative permeability given an N_c,end of >0.56, which is obviously lower that that required for the JBN method to return reasonable estimates.

5. Conclusions

An improved assisted history matching framework is developed and implemented to improve the fidelity of relative permeability estimates where the impacts of end-cap effects, capillary pressure differentials, and sorptive/stress-induced changes in permeability are present. This new algorithm (BADS) is applied to recover gas/water relative permeability curves of coals from unsteady-state gas-displacing-water tests and to define its improvement over existing methods. Synthetic phase displacement experiments are conducted by numerical simulation and by using varied pressure drops to both validate the proposed framework and to determine its relative improvement. The performance of the framework was compared with the performance of the existing JBN algorithm. We demonstrate that the BADS algorithm is capable of reproducing the relative permeabilities from unsteady-state displacement tests on synthetic samples with relatively high accuracy provided that the end-effect number is higher than 0.56. The BADS algorithm is remarkably efficient, which requires only a small number of numerical simulation runs to achieve stable convergence: ≈100 to 150 runs. This rapid convergence represents a significant improvement over previous algorithms, such as EnKF, LM, and GA, that typically require in excess of 500 numerical simulation runs to converge for the same representative problem. Furthermore, we demonstrate that the JBN method leads to unreliable interpretation results in the presence of capillary end-effects and when permeability evolution as a result of changes in pressure (impacts of sorption/swelling and changes in effective stress) are present and significant. Conversely, the JBN method produces reasonable estimations only where relatively large pressure drops (corresponding to an end-effect number of >3.3) are present and where changes in permeability due to gas injection are negligible. Future work is anticipated regarding application of the proposed framework in interpreting gas/water relative permeability of coal from real core flooding tests so that the effects of varying factors (i.e., pore pressure, displacement rate, etc.) on the relative permeability can be identified. The proposed framework should be rather easily extended for water/oil, gas/oil, and gas/water in shale, sand, and carbonate rocks. As far as we are concerned, the only issue of the proposed method for application in other porous flow problems is the representation of the relative permeability pressure curve. As has been stated previously, for the gas/water/coal system, the relative permeability curve has been proven to be well represented with Chen et al.’s model, which however may be replaced with other representation methods such as the Corey-type or B-spline models [31] for other type of rocks.