The fracture network formed after volume fracturing of shale gas reservoir is very complex with fractures of multi-scales, resulting in difficult description and simulation of fractures and difficult productivity prediction of shale gas wells. To obtain accurate productivity evaluation results, 3 conditions should be met: (1) reasonable characterization of spatial distribution of the fracture network; (2) accurate simulation of flow process of fluid in the fracture network; and (3) efficient and accurate matching of historical data of gas well production.

Natural fractures in shale are characterized by discrete fracture medium, which is difficult to be simulated by the traditional equivalent continuous model^[1,2]. The discrete fracture network (DFN) model describes fractures as fracture sheets with different scales and shapes, and uses the stochastic simulation method, ant body tracking algorithm and well-seismic data to obtain actual fracture network^[3,4,5]. But limited by the stochastic simulation method, this model has large uncertainties in the key characterization parameters. Since the geometric configuration of natural fractures has fractal characteristic^[6,7,8] (i.e., the fracture characteristic parameter is the statistical probability field of fractal dimension), compared with the uniform random theory that the DFN model founded on, the fractal geometry theory uses the fractal dimension to describe stochastic, irregular but statistically self-similar fracture system^[7]. The generated fractal discrete fracture model (FDFM) can more accurately describe the heterogeneity and complex geometries of natural fractures. Fracture systems with the same fractal parameter may have different distribution patterns and geometric shapes but almost the same reservoir flowing capacity, i.e., the same equivalent permeability of the fracture network^[8,9], and the corresponding dynamic curves of pressure are also similar^[10,11]. The self-similarity and scale invariance of natural fractures make the fractal parameters able to control the overall distribution of natural fractures, this ensures the retainment of key mathematical information while improving the characterization accuracy of natural fractures.

Restricted by the limitations of local grid refinement^[12], equivalent continuous model^[13,14] and perpendicular bisection^[15,16] techniques in characteristic scale, grid subdivision, simulation accuracy and computational efficiency, the embedded discrete fracture model (EDFM) has emerged in recent years. EDFM can avoid the complexity of perpendicular bisection and non-convergence of the associated algorithm^[17,18], and ensure the flow simulation accuracy between fractures and between fracture and matrix^[19,20], so it has obvious technical advantages in efficient simulation of productivity of any complex fracture network. Moreover, the history matching based on the productivity model is a necessary link for parameter inversion and productivity prediction. The automatic history matching processes based on ordinary optimization algorithms such as Newton optimization algorithm, particle swarm optimization algorithm and genetic algorithm are cumbersome, time-consuming and labor consuming, making the most representative matching results difficult to seek^[21]. In contrast, the auxiliary history matching technique based on intelligent algorithm can fully consider the strong correlation between the main control factors of productivity, improve the history matching efficiency and reduce the interference of human factors, but because of the complexity of simulation itself, the research and development of appropriate matching algorithm are very challenging.

This paper presents a productivity evaluation method from fracture modeling, dynamic simulation, history matching and productivity prediction, including: (1) Using fractal geometry theory to build discrete natural fracture characterization model. (2) Based on the generated 3D natural/artificial fracture model, combined with the new generation of EDFM, accurately simulating the production performance of shale gas wells in the whole life cycle. (3) Performing machine learning through automatic sampling with the help of intelligent algorithm, realizing effective inversion of key reservoir and fracture parameters, and quantifying risk of fracture parameters and productivity.

Fractal theory is divided into 2D fractal and 3D fractal, which have their specific parameters respectively. The fractal dimension of fractures can be obtained by outcrop, imaging logging and core analysis. The parameters such as central point, opening, length, strike and distribution of fractures can control the cluster distribution of fracture system. At present, 2D fractal fracture system is universally adopted in the world, and 3D fracture adopts the disc model, but this is not consistent with the actual fracture morphology^[9].

The Euclidean dimension corresponded by 3D space is 3. The number of fractures with fracture length between l and l+dl is denoted asn(I,L)dl, which meets^[6]:

α_3D determines the number of fractures in 3D space, which is independent of the model scale and needs to be calculated. According to the Piggott equation^[22], the relation between 2D and 3D density constants is obtained as follows:

Usually, the parameter α_3D is introduced to characterize the ratio of long fractures to short fractures, which is named power law length distribution index by Darcel et al.^[22] When α_3D=3, the long fractures predominate, when α_3D=$\infty$, the fractures are consistent and approach the minimum l_min in length. Its relationship with D_l3D meets the following relation^[3]:

At the same time, the relationship between 2D fractal dimension and 3D fractal dimension meets^[6]:

In Eq. (4), the range of D_c is 1-2, and the range of D_l is 1- $\infty$. Through integration on Eq. (1), the total number of fractures with length greater than l_min in the 3D model is obtained^[9]:

Eq. (5) is also suitable for the calculation of number of 2D fractal fractures. The distribution of geometric central point of fracture in 2D plane space can be calculated using fractal dimension D_c and pairwise correction function $C_{2}(r)$ ^[23]:

It should be emphasized that the multiple information superposition algorithm suitable for the generation of geometric parameters of 2D fractal fracture is also suitable for 3D fracture model, and the algorithm is mainly used in the calculation of regional distribution probability of natural fractures. The 2D model is decomposed to multiple subregions after repeated iterations (7-9 iterations), the probability of natural fractures occurring in each region is constrained by specific equation^{[6, 22]}, and the total probability of all subregions is not necessarily equal to 1 (otherwise, it is a completely random distribution). The iteration times of multiple information superposition meet^[9]:

Based on Python language, a 3D fractal fracture generation module has been developed in this study, which can generate fracture attribute parameters. Fig. 1 shows the 3D distribution of central points of natural fractures generated by multiple information superposition process; points in Fig. 1 represent the spatial positions of central points of natural fractures, and the specific parameters given are: D_c3D=2.8, D_l3D=3.5, L=1000 m, α=0.5, l_min=50 m, and N_iter=6.

The 3D fractal model presented by Kim et al.^[9] assumed that the fracture was disc shape and the simulation area was a box with equal length, width and height, which is inconsistent with the actual fracture morphology; therefore, the equivalent area method is used to correct the disc model to a rectangular model. Because the length of natural fracture follows the power law distribution, rather than the normal distribution, the power law distribution is used to simulate the length of natural fractures^[11]:

α_3D and D_c3D control the distribution characteristics of fracture system in 3D space: (1) α_3D determines the ratio of the number of long fractures to that of short fractures (values of 3-$\infty$). When α_3D=3.0, long fractures predominate in the model. As α_3D increases, the proportion of long fractures in the model decreases, and when α_3D approaches infinity, all fractures in the model are consistent in length (equal to the minimum length l_min). (2) D_c3D determines the heterogeneity of fracture distribution (in the range of 2.0 to 3.0). When D_c3D=2, the fractures are unevenly distributed and occur in clusters with many blank areas; as D_c3D increases, the fracture distribution tends to be uniform. (3) Under extreme conditions of α_3D=$\infty$ and D_c3D=3.0, the fractures are consistent in length and uniform in distribution, and this case is equivalent to a continuous model.

Fig. 2 shows the statistical results of lengths, strikes and dip angles of natural fractures under two groups of fracture conditions. It can be seen that the fracture lengths basically exhibit skewed normal distribution, whereas the azimuth and dip angle of natural fractures tally with Gaussian normal distribution.

When the 3D fracture is described as a rectangle, the corresponding fracture height is determined by the product of the number between 0 and 1 randomly generated and the main length of the natural fractures. By using the calculation results of geometric parameters such as central point, length, height, strike and dip angle of fractures, taking the central point as the basic point, the spatial morphology of fractures is determined based on the strikes and angles of fractures, the geometric size of fractures is determined based on the length and height of fractures, and the spatial distribution map of 3D fractal fracture geometry is obtained finally (Fig. 3).

The 3D fractal discrete fracture model is modified based on 2D model algorithm. For the convenience of intuitive analysis, 2D model was taken as an example to verify the fracture generation algorithm. The parameters controlling the 2D fractal discrete fracture model include D_c, the fractal dimension of 2D fracture distribution; D_l, the fractal dimension of 2D fracture length distribution; α, the 2D fracture density constant; and the median of fracture azimuth. In the example, the size of the generated area was set at 10 m× 10 m, the minimum fracture length l_min=0.05 m, and the rest required fractal parameters are listed in Table 1. The parameters were input, the fractal discrete fracture network was obtained through the fractal fracture generation model (Fig. 4), with 9326 fractures in total of the model.

Based on the distribution of generated fractures in Fig. 4, the distributions of azimuth, length and position of fractures were counted (Fig. 5); by using Eq. (5) under 2D conditions (in the form of $\frac{N\left( L \right)}{{{L}^{{{D}_{\text{c}}}}}}=\frac{\alpha }{{{D}_{\text{l}}}}l_{\min }^{-{{D}_{\text{l}}}}$) and Eq. (6), data fitting was used to obtain fractal parameters of fractures (Table 1). The simulation results show that the fitting results of D_c, D_l and median azimuth of fractures are basically the same as the input values. Although the calculation result of α is somewhat different from the input value, it is still within the acceptable error range, indicating that the fractal fracture generation algorithm in this paper is reliable, that is, the reliable 2D fractal discrete fracture generation algorithm can be used to obtain accurate 3D fractal fracture network model^{[9, 11]}.

On the generated 3D fractal natural fracture system, the staged fracturing horizontal well was set as the basic well type to characterize the post-fracturing complex fracture system with coupled natural fractures and artificial fractures (Fig. 6), and the fracture system was combined with the numerical simulator.

The black oil was used to conduct the flow simulation. Considering the geological and development characteris

tics of shale gas reservoir, the following assumptions were made for the model.

(1) The reservoir was homogeneous, horizontal and isopachous, and isotropic in planar matrix permeability; as bedding in shale affects the vertical flow of gas, the vertical permeability of the shale was set at 0.1 times of the plane permeability.

(2) The gas reservoir adopted 3D single pore numerical model, which included 3 different scales of spaces: matrix, natural fracture and artificial fracture. The matrix system was described by regular orthogonal grids, the fracture system (natural and artificial fractures) was described by virtual grids generated by EDFM. The greatest advantage of this method is that it considers fluid transfer characteristics between different types of pores.

(3) In order to simplify the gas desorption process, the instantaneous desorption model was used to describe the desorption process of adsorbed gas, and the BET (Brunauer-Emmett-Teller) multilayer isothermal adsorption model was selected as the adsorption model:

where, C and n are dimensionless experimental constants. When n=1, the BET multilayer model is simplified into the Langmuir isothermal adsorption model:

where, v_L=v_m, p_L=p_o/C.

(4) Considering that the flowback of fracturing fluid, the fluid flow in the fractures in the course of production is the gas-water two-phase flow, so the Corey relative permeability model was employed:

where, K_rwo, K_rgo, N_w and N_g are experimental fitting values in Eqs. (11)-(12). The model adopted the non-equilibrium initialization method, i.e., the water saturation was higher in the fractures and lower than the bound water saturation in the matrix.

(5) Considering the stress sensitivity of permeability in artificial fracture system, the variation of permeability with pressure meets the attenuation equation:

(6) Considering fluid gravity and capillary pressure between phases, in this study, the gas and water phases were regarded as 2 components without mass transfer. Based on the Darcy's law, the volume flow velocity of component j_gw (j_gw=g represents gas phase, j_gw=w represents water phase) meets the following equation:

The flow behaviors in fractures and matrix are characterized by K and $K_{r,j_{gw}}$ separately.

(7) By combining Eqs. (9)-(14) with the mass conservation equation, the pressure control equations of gas and water phases are obtained as follows:

The orthogonal grid was used to conduct finite difference processing on the mathematical model (Eq. (15)). In this study, the independently developed embedded discrete fracture model (EDFM) technique was used to process the fracture connection^[24]. On the basis of orthogonal grids, EDFM can quickly simulate complex fracture systems by using equivalent processing method while maintaining high calculation efficiency of regular grids.

In order to keep the consistency of fracture volume before and after processing, the porosity of virtual fracture grid should meet the following equation:

The core of EDFM technique is the transmissibility calculation of non-neighboring connection pairs (NNC), which is mainly used to deal with the flow exchange between grids which are neighboring in physical model but not neighboring in computing domain. The general calculation equation of T_NNC, the transmissibility of non-neighboring grid connection pairs, is:

where, when the fracture is connected with the matrix, d_NNC is the average distance from the matrix block to the fracture plane, and K_NNC is the matrix permeability; when two fractures are connected, d_NNC is the normal distance between the fracture elements, and K_NNC is the average permeability of the fractures.

The reliability and computational efficiency of the EDFM model were verified by comparing the calculation results of local grid refinement (LGR), perpendicular bisection (PEBI) and embedded discrete fracture models (EDFM). Considering the limitation of LGR method, a regular staged fracturing horizontal well was taken as an example. The gas reservoir size was 1700 m×400 m×20 m, the grid had a single layer vertically, and the fracture was assumed to run through the reservoir with a fracture aperture of 0.01 m and a porosity of 35%. The gas well produced single-phase gas for 20 years at the constant bottomhole pressure of 1.5 MPa. The other basic parameters are listed in Table 2.

In LGR and EDFM, the number of basic grids was 170×40×1. In LGR model, the logarithmic refinement principle was adopted to split the grid in the direction perpendicular to the fracture plane, and uniform refinement grid was adopted in the direction of the fracture plane. In PEBI model, the grids near the fracture were refined radially and locally, and constructed according to the flow field around the fracture. In EDFM model, fractures were connected with the original basic grid in the form of virtual grid, and the real attributes of fractures were described by explicit equations, with no need of equivalent processing.

Fig. 7 shows the production performance data of the gas well simulated by LGR, PEBI and EDFM, and the results of EDFM and PEBI are highly consistent. Limited by the grid subdivision degree, the result of LGR has certain error, which is still in an acceptable range. It should be emphasized that the accuracy of LGR simulation results is related to the grid subdivision degree, the denser the grid near the fracture, the higher the calculation accuracy is, but the calculation amount of the model also increase greatly at the same time^[19].

Table 3 compares the computational efficiency of EDFM, PEBI and LGR. LGR needs a large number of fracture grids to describe the flow characteristics around the fracture. Due to the limitations of orthogonal grid itself, this method can’t depict complex fracture network.

Moreover, when the reservoir permeability is low, more refined grids are needed, which greatly increases the calculation amount. Although PEBI has the smallest number of grids on the whole, the number of grids needed to depict fractures is still large. When the complex fracture network is described, more refined grids are needed to calculate the fluid exchange between fracture intersections, which seriously impairs the PEBI computational efficiency. In contrast, EDFM can effectively reduce the number of fracture grids while ensuring the calculation accuracy, and thus significantly save simulation cost and improve simulation efficiency. Especially, EDFM has more obvious calculation advantage in simulating fracture network with complex configuration.

The automatic history matching of gas reservoir is essentially a method to obtain the minimum value of objective function. By making use of the powerful ability of EDFM in handling complex fracture network, combined with Markov chain-Monte Carlo (MCMC) machine learning algorithm, intelligent history matching module (NN-MCMC)^[25] based on neural network (NN) has been developed in this study. This new method has reached efficient and accurate evaluation of complex fracture system (including key parameters such as effective fracture height, fracture length, conductivity, cluster efficiency and stress sensitivity of fractures).

Fig. 8 shows the intelligent history matching process. Firstly, the measured data of gas wells such as gas/water production rate and wellhead pressure are sorted out. Secondly, the ranges of unknown parameters are set based on the engineering/geological data. The sampling method in the orthogonal experimental method is used to conduct primary iteration sampling on the parameter combinations to obtain 50 samples. The modular input file and the EDFM model are used to synchronously simulate 50 samples to obtain the simulation results, and then, 2 output parameters such as bottomhole flowing pressure and gas/water ratio or water production rate are selected to conduct global error calculation. It should be emphasized that the data at a specific time point (named index number of data point) should be selected to compare the real data with the simulation output, e.g., Eq. (18):

After obtaining the results of the initial cycle, the artificial intelligence-neural network is used to build a related surrogate model for 50 input variables and output variables. Furthermore, the Markov chain-Monte Carlo sampling method is used to conduct systematic sampling on the input variables^[25] to ensure that the inversion combination is representative and the global error value can be optimized. After completing the sampling, thousands of sample combinations can be obtained, and the trained artificial intelligence-neural network surrogate model is used to conduct simulation and prediction of the representative sample combinations. Then, Eq. (18) is used to calculate the global errors of thousands of representative samples, and 25 samples with the smallest global error are selected as the input values of the next iteration to conduct simulation operation. Repeat this step until the prediction of the surrogate model converges or reaches the maximum number of iterations.

By means of repeated iterations, this method can ensure the accuracy of the surrogate model improves with the increase of the number of iterations. After running all iterations, the errors of the objective functions (i.e., the two selected output parameters) need to be calculated and compared to screen the history matching solution. The calculation method of the objective function is similar to Eq. (18), and the difference is that the error is calculated separately for each matching parameter, rather than summing the matching parameters to obtain the error of objective function. By setting the screening threshold of each objective function, multiple history matching solutions meeting the error requirements are screened out. And finally, statistical analysis is conducted, and long-term and representative productivity prediction results are simulated.

Horizontal well A in the shale gas reservoir of the Zhaotong region was taken as an example. In this well block, the natural fractures are mostly in near SN, NW and NE strikes, with dip angle averaging 80°. The measured data of fracture extension lengths in the outcrop area show that most of the fractures are less than 100 m in length; therefore, in the example, according to the simulation requirements, fracture slits were set at 0-100 m in length, about 50 m averagely, and 2:1 in length/height ratio. The fractal parameters of natural fractures in the control range of the well are listed in Table 4. In this well block, the shale has a porosity of 5.5%-6.5%, and an original matrix permeability of about (1.0-10.0)×10^-8 μm². Other geological and engineering parameters are listed in Table 5. The well was put into production in August 2015 after staged fracturing and flowback.

The fracturing simulation software was used to simulate the hydraulic fracturing of the well. The statistical results of hydraulic fracture parameters are shown in Fig. 9, among which, the hydraulic fractures are 49-712 m long, 321 m on average; and 2-107 m high, 58 m on average; the propped fractures are 4-606 m long, 231 m on average; and 0-105 m high, 32 m on average.

As the difference of two horizontal principal stresses is large in Southern Sichuan Basin, and the simulated artificial fractures are generally straight, which is similar to the standard staged fracturing horizontal well; therefore, the average value of fracture model parameters was used to set the geometric parameters of artificial fractures.

Thus, the productivity evaluation model of the shale gas well composed of matrix, natural fracture and artificial fracture was built. Some model parameters were modified to match the production history, especially match the development indexes such as daily gas production rate, daily water production rate and bottomhole flowing pressure. Based on dynamic production data, the automatic history matching and inversion technique were used to effectively carry out the post-fracturing effect evaluation to reduce the uncertainty of core parameters. Table 6 lists the upper and lower limits of parameter values participating in the matching.

In this example, 12 steps of automatic iterations (10 million samples per iteration) were carried out, among which 325 gas reservoir models were selected to conduct history matching. The total error decreased with the increase of iterations. Setting the global error ε<50%, 75 sets of history matching solutions were picked out, and the matching effects are shown in Fig. 10.

The combinations of uncertain parameters involved in automatic history matching can be represented by parallel coordinate diagram, and Fig. 11 shows the simulation results of all 325 times of automatic history matching. The abscissa represents every uncertain parameter, the ordinate represents the range of each uncertain parameter, the red line represents the combination of the optimal solution, the orange lines represent the combination of all optimized history matching solutions (75 times), and the grey line represents the combination of non-history matching solutions (250 times).

P₅₀ values of matched parameters from statistics are listed in Table 7. After history matching, the cumulative water production and cumulative gas production of the gas well in 20 years were predicted. Based on the obtained interpretation parameters, the constant bottomhole flowing pressure mode was used to predict the gas well productivity (at the bottomhole pressure of 1.5 MPa), and the results are shown in Fig. 12. The productivity prediction curves in this figure are the results under combinations of different model parameters, and the curve cluster contains the actual performance curve. Fig. 13 shows the statistics on the prediction results of ultimate produced fluid in 20 years. It can be seen that in the P₁₀-P₉₀ confidence domain, the cumulative gas production of the gas well is (1.05-1.22)×10⁸ m³, the optimal solution is 1.15×10⁸ m³, the cumulative water production of the gas well is (0.42-0.81)×10⁴ m³, and the optimal solution is 0.72×10⁴ m³. The production data analysis method combining a variety of production decline models is an important means for productivity evaluation of shale gas well, and the cumulative gas production and cumulative water production of single well obtained using this method were 1.01×10⁸ m³ and 0.81×10⁴ m³ respectively^[26]. The interpretation results of above two methods are highly similar, confirming that the prediction results of this method have high confidence level.

The 3D natural fracture generation method based on multiple information superposition algorithm of fractal theory can use fractal parameters to control the overall distribution of fracture network, also can characterize the complex multiscale fracture system in shale after fracturing combined with artificial fractures.

The embedded discrete fracture model, with the advantages of small number of fracture grids and short computation time, can flexibly characterize the properties of natural and artificial fractures, and accurately simulate the exchange process of fluid in matrix and fracture while effectively reducing the calculation amount.

The combination of embedded discrete fracture model and intelligent history matching algorithm can reduce the uncertainty of calculation of unknown parameters such as fracture and reservoir, realize effective inversion of key reservoir and fracture parameters and quantitative prediction of gas well productivity. Practical application has verified that the prediction results of the integrated shale gas well productivity evaluation model have high confidence level.

A_NNC—contact area between connection pairs, m²;

c—Constant;

C—relevant constant of net heat consumption, dimensionless;

C₂(r)—correction function;

c_ag—proportion of adsorbed gas in pore volume, dimensionless;

d_NNC—distance between connection pairs, m;

D—formation depth, m;

D_c—2D fractal dimension;

D_c3D—3D fractal dimension, dimensionless;

D_l3D—fractal dimension of 3D fracture length, dimensionless;

D_l—fractal dimension of 2D fracture length, dimensionless;

F—random number between 0 and 1, dimensionless;

i—matching parameter No.;

j—matching data No.;

j_gw—fluid No. (j_gw=g represents gas phase, j_gw=w represents water phase);

K—matrix permeability, 10^-3 μm²;

K_i—permeability under initial formation pressure, 10^-3 μm²;

K_NNC—connection permeability, 10^-3 μm²;

$K_{r,j_{gw}}$—relative permeability of fluid;

K_rw, K_rg—relative permeability of water and gas;

K_rwo, K_rgo—relative permeability endpoint values of water phase and gas phase, dimensionless;

l—Fracture length, m;

l_ratio—side ratio of region to subregion, dimensionless, with value of 2;

l_min, l_max—the minimum and maximum natural fracture length, m;

L—model scale, m;

n—number of adsorption layers;

n(I,L)—number of fractures with length of l-l+dl in unit scale, fractures/m;

N—total number of fractures;

N_w, N_g—relative permeability indexes of water phase and gas phase, dimensionless;

N(L)—total number of fractures with length greater than l_min;

N_iter—iteration times;

N_p(r)—number of fracture pairs with distance between centers of 2 fractures less than r;

p—pore pressure, MPa;

p_i—initial formation pressure, MPa;

p_L—Langmuir gas pressure, MPa;

p_o—saturation pressure of gas phase, MPa;

P—number of matching parameters;

P₁₀, P₅₀, P₉₀—values at the cumulative probability of 10%, 50% and 90%;

Q—quantity of matching data;

q_g—gas production rate, m³/s;

q_w—water production rate, m³/s;

r—main length of natural fracture, m;

S_g, S_w—gas phase and water phase saturations, %;

S_gr—gas phase endpoint saturation, %;

S_seg—fracture plane area, m²;

S_wirr—water phase endpoint saturation, %;

t—time, s;

T_NNC—transmissibility of non-neighboring grid connection pairs, 10^-3 μm²•m;

u—fluid flow rate, m/s;

v_ag—gas volume adsorbed by unit mass of rock, m³/kg;

v_L—Langmuir volume, m³/kg;

v_m—the maximum adsorbed gas volume per unit mass of rock, m³/kg;

V_b—microelement volume, m³;

V_f—fracture volume, m³;

w_f—fracture width, m;

w_i,j—weight value of the ith data point corresponded by the jth parameter, dimensionless;

x, y, z—coordinate axis, m;

X—function argument;

x_model—simulated production index, dimensionless;

x_history—actual production index, dimensionless;

α—2D fracture density constant, dimensionless;

α_3D—3D fracture density constant, dimensionless;

$\Gamma \left( X \right)$—gamma function;

γ—permeability attenuation constant, MPa^-1;

$r_{j_{gw}}$—fluid gravity, kg/(m•s)²;

μ—fluid viscosity, mPa·s;

u—fluid flow rate, m/s;

u_g, u_w—flow rate of gas phase and water phase, m/s;

ρ_g, ρ_ag, ρ_s, ρ_w—density of gas, adsorbed gas, rock and water, kg/m³;

ϕ—porosity, %;

ϕ_f—fracture grid porosity, %;

ε—global error, dimensionless.