Investigation of the Gas Pressure Field and Production Rate for Two Typical Proppants: Small-Sized Continuous and Large-Sized Discontinuous Proppants

Abstract

The small-sized proppant is widely used in the traditional hydraulic fracturing reservoir stimulation, but the theoretical research shows that the large-sized proppant can greatly improve the fracture permeability. Although many scholars have proposed the method of using large-sized proppant, the characteristics of the pressure field change and gas drainage rate during reservoir development are still unclear after pumping large-sized proppant. In this article, to study the change regularity of the reservoir pressure field, two representative proppants are selected: the small-sized proppant and large-sized proppants. First, the Navier–Stokes equations are solved using the numerical simulation method, and the characteristics of the reservoir pressure field are finely reproduced. Then, the production rate is discussed to reveal the huge potential of large-sized proppant for the natural gas development. The results show that under the same conditions, if the particle size of the proppant increases by 5 times, then the reservoir permeability will increase by 27 times approximately, and the single well production efficiency will increase by 17~19 times in the first 3000 days. In addition, a new quantitative model is proposed to evaluate the permeability magnification of the fracture and reservoir when adopting the large-sized proppant. This study further confirms that the method of large-sized proppant proposed by the author in the earlier stage has great potential. This study is helpful for the researchers and engineers to better understand the evolution regularity of the reservoir pressure field and gas production rate in the process of oil and gas exploitation after using the large-sized proppant.

1. Introduction

The development and utilization of natural gas (shale gas, coalbed methane, etc.) is of great significance to the sustainable development of the world [1]. Under the influence of in situ stress, the reservoir porosity and permeability are generally low [2,3]. At present, the low reservoir permeability is the main factor restricting the development efficiency of natural gas. Hydraulic fracturing is a key technology to improve reservoir permeability [4]. Hydraulic fracturing refers to the technology that uses the high-pressure fluid, such as water or liquid nitrogen, to fracture the rock formation to form a fracture network, thereby improving the permeability of the rock formation [5]. The high pressure forces the rock rupture to form fractures and eventually form artificial fractures in the formation [6,7,8]. At present, hydraulic fracturing technology has been widely used in the field of reservoir stimulation and permeability enhancement of coal, shale, and tight sandstone [9,10,11,12,13].

The reservoir stimulation mainly includes three steps. For the first step, the rock is fractured using the hydraulic fracturing method. For the second step, proppant is pumped into the fracture network. Third, fracturing fluid backflow. Finally, the proppant is retained in the fracture network to support the fracture and prevent its closure. The fracture has good conductivity under the support of proppant, so the permeability of the reservoir will increase after hydraulic fracturing stimulation [14].

Although the fracturing techniques are used, the reservoir permeability increase is still very limited in some situations, and the support effect of the proppant tends to weaken over time [15]. Proppant plays an important role in the reservoir stimulation and permeability enhancement [16]. How to improve the reservoir permeability from the perspective of proppant is the key research direction.

In the past decades, quartz sand and ceramsite were mostly used as proppants [17]. The shapes of proppants were round, oval, cylindrical, and so on [18,19]. The shape of the proppants had a certain impact on fracture permeability, but the impact was limited. The size of the proppant has an obvious influence on the crack permeability. According to the traditional permeability evaluation theory, for example, Osiptsov’s model k/R² = 0.204φ^4.58 [20], the larger the proppant particle size, the greater the fracture permeability under the premise of constant porosity. This conclusion can also be confirmed by the other permeability models [21,22]. At present, the commonly used proppant type is mainly between 8 and 140 mesh, and the corresponding particle size range is 2.36–0.106 mm [17,23,24]. Usually, these small-sized proppants exist between the fractures in the form of continuous packet, as shown in the left part of Figure 1.

Gillard’s research shows that the use of discontinuous proppant distribution can significantly improve the fracture permeability [25]. The discontinuous proppant mentioned by Gillard means that the proppant is stacked into discrete packets, and there are channels of a large size between the proppant packets, as shown in the right part of Figure 1. Hou et al. verified Gillard’s view with the indoor experiments [26]. Hou et al. found that the proppant in the form of discrete packages can significantly improve the fracture conductivity [26]. Similar conclusions were given by Medvedev et al. (2013) [27]. Medvedev et al. (2013) [27] found that the permeability of a fracture containing large-sized proppant was much higher than that of proppant packs. Li et al. [13] drew on Gillard’s research idea and proposed a layout scheme of large-sized proppant, as shown in Figure 2. Li’s research shows that the large-sized proppant can significantly improve fracture permeability. Essentially, the discrete proppant package mentioned by Gillard is equivalent to the large cylindrical proppant mentioned by Li et al. [13].

Although Gillard and Li pointed out that the large-sized proppant can improve the fracture permeability, the improvement level has not been quantitatively evaluated yet. Furthermore, the research on the reservoir pressure field and production efficiency is still lacking after using the large-sized proppant. Compared with the traditional scheme of the continuous small-sized proppant package, the dynamic change characteristics of the reservoir pressure and gas production rate under the large-sized proppant scheme still lack understanding, which needs further research. This paper focuses on this topic.

2. Equations and Methods

2.1. Conservation Equations

There is a fracture network in the reservoir after the hydraulic fracturing of rock or shale. Fluid (gas and water) migration in the porous media follows the laws of mass conservation and momentum conservation [28]. The mass equations of the gas phase and the water phase are:

$$-\nabla (\frac{v_g}{B_g})+q_{sg}+q_{smfg}=\frac{\partial }{\partial t}(\frac{{\phi }_f{s}_g}{B_g})$$

(1)

$$-\nabla (\frac{v_w}{B_w})+q_{sw}=\frac{\partial }{\partial t}(\frac{{\phi }_f{s}_w}{B_w})$$

(2)

where v_g is the gas velocity, m/s; B_g is the gas volume coefficient, dimensionless; q_sg is the source term of formation gas, 1/s; q_smfg represents the gas source term that diffuses from the matrix into the fracture, 1/s; φ_f is the average fracture porosity of the reservoir, dimensionless; s_g is the gas saturation, dimensionless; v_w is the water phase velocity, dimensionless; B_w is the volume coefficient of the water phase, dimensionless; q_sw is the formation water source, 1/s; and s_w is the gas saturation, dimensionless. It is generally believed that the migration process of fluid from the matrix to the fracture only includes the diffusion process of the gas phase and the diffusion process of the anhydrous phase.

For the gas and water phases, the Darcy equations are written as:

$$v_g=-\frac{k_{\infty }{k}_{rg}}{{\mu }_g}(\nabla P_g-{\rho }_{g}g\nabla h)$$

(3)

$$v_w=-\frac{k_{\infty }{k}_{rw}}{{\mu }_w}(\nabla P_w-{\rho }_{w}g\nabla h)$$

(4)

where k_∞ is permeability, m²; k_rg is the relative permeability of gas, dimensionless; P_g is the gas pressure, Pa; μ_g is the gas viscosity, Pa·s; ρ_g is the gas density, kg/m³; k_rw is the relative permeability of the water phase, dimensionless; P_w is the water phase pressure, Pa; μ_w is the viscosity of the water phase, Pa·s; ρ_w is the density of the water phase, kg/m³; and h is the reservoir depth, m. Combine Equations (1) and (3) to obtain the gas phase pressure control equation, as follows:

$$\nabla \left[\frac{k_{\infty }{k}_{rg}}{B_g{\mu }_g}(\nabla P_g-{\rho }_{g}g\nabla h)\right]+q_{sg}+q_{smfg}=\frac{\partial }{\partial t}(\frac{{\phi }_f{s}_g}{B_g})$$

(5)

Combine Equations (2) and (4) to obtain the pressure control equation of the water phase, as follows:

$$\nabla \left[\frac{k_{\infty }{k}_{rw}}{B_w{\mu }_w}(\nabla P_w-{\rho }_{w}g\nabla h)\right]+q_{sw}=\frac{\partial }{\partial t}(\frac{{\phi }_f{s}_w}{B_w})$$

(6)

Supplement the auxiliary equation, and the gas–liquid saturation meets the following relationship:

$$s_g+s_w=1$$

(7)

The capillary pressure meets

$$P_{cgw}=P_g-P_w$$

(8)

So far, there are 4 unknown variables (P_g, P_w, s_g, s_w) and 4 Equations (5)–(8). Hence, the equation can be solved.

2.2. Numerical Computing Method

First, Equations (5) and (6) are deformed and then added to obtain the gas–liquid two-phase pressure control equation, as follows:

$$B_g\nabla \left[{\beta }_g\nabla (P_g-{\rho }_{g}gh)\right]+B_w\nabla \left[{\beta }_w\nabla (P_w-{\rho }_{w}gh)\right]+B_g(q_{sg}+q_{smfg})+B_w{q}_{sw}={\phi }_f{C}_t\frac{\partial P_g}{\partial t}$$

(9)

where ${\beta }_g=\frac{k_{\infty }{k}_{rg}}{B_g{\mu }_g}$, ${\beta }_w=\frac{k_{\infty }{k}_{rw}}{B_w{\mu }_w}$. The definition of the total compression coefficient is $C_t=C_f+C_g{s}_g+C_w{s}_w$. For Equation (9), the common discretization methods include explicit discretization and implicit discretization. In this paper, the implicit discretization is adopted, and Equation (9) can be discretized into

$$a_m{P}_{i,j,k}^{n+1}+a_w{P}_{i-1}^{n+1}+a_e{P}_{i+1}^{n+1}+a_s{P}_{j-1}^{n+1}+a_n{P}_{j+1}^{n+1}+a_t{P}_{k-1}^{n+1}+a_b{P}_{k+1}^{n+1}=B_{i,j,k}$$

(10)

where

$$\begin{array}{l}{a}_w=B_{g,i,j,k}{T}_{g,i-1/2}+B_{w,i,j,k}{T}_{w,i-1/2}\\ a_e=B_{g,i,j,k}{T}_{g,i+1/2}+B_{w,i,j,k}{T}_{w,i+1/2}\\ a_s=B_{g,i,j,k}{T}_{g,j-1/2}+B_{w,i,j,k}{T}_{w,j-1/2}\\ a_n=B_{g,i,j,k}{T}_{g,j+1/2}+B_{w,i,j,k}{T}_{w,j+1/2}\end{array}$$

$$\begin{array}{l}{a}_t=B_{g,i,j,k}{T}_{g,k-1/2}+B_{w,i,j,k}{T}_{w,k-1/2}\\ a_b=B_{g,i,j,k}{T}_{g,k+1/2}+B_{w,i,j,k}{T}_{w,k+1/2}\end{array}$$

$$a_m=-\left[a_w+a_e+a_s+a_n+a_t+a_b+{{(\mathrm{V}}_B{\phi }_f{C}_t)}_{i,j,k}^n/\mathsf{\Delta }t\right]$$

$$B_{i,j,k}=-\left[B_{g,i,j,k}{(-\mathrm{G}}_g{+\mathrm{Q}}_{sg,i,j,k}{+\mathrm{Q}}_{smfg,i,j,k})+B_{w,i,j,k}{(-\mathrm{G}}_w{+\mathrm{Q}}_{sw,i,j,k})+P_{g,i,j,k}^n{{(\mathrm{V}}_B{\phi }_f{C}_t)}_{i,j,k}^n/\mathsf{\Delta }t\right]$$

$${\mathrm{G}}_g=T_{g,j+1/2}{\rho }_{g}g(h_{j+1}-h_j) - T_{g,j-1/2}{\rho }_{g}g(h_j-h_{j-1})$$

$$\begin{array}{l}{\mathrm{G}}_w=T_{w,j+1/2}{\rho }_{w}g(h_{j+1}-h_j)-T_{w,j-1/2}{\rho }_{w}g(h_j-h_{j-1})\\ \hspace{1em} +T_{w,i+1/2}(P_{cgw,i+1}-P_{cgw,i})-T_{w,i-1/2}(P_{cgw,i}-P_{cgw,i-1})\\ \hspace{1em} +T_{w,j+1/2}(P_{cgw,j+1}-P_{cgw,j})-T_{w,j-1/2}(P_{cgw,j}-P_{cgw,j-1})\\ \hspace{1em} +T_{w,k+1/2}(P_{cgw,k+1}-P_{cgw,k})-T_{w,k-1/2}(P_{cgw,k}-P_{cgw,k-1})\end{array}$$

$$\begin{array}{l}{\mathrm{Q}}_{sg,i,j,k}=V_B{q}_{sg,i,j,k}\\ {\mathrm{Q}}_{smfg,i,j,k}=V_B{q}_{smfg,i,j,k}\\ {\mathrm{Q}}_{sw,i,j,k}=V_B{q}_{sw,i,j,k}\end{array}$$

where $T_{g,j+1/2}$ is given as follows:

$$T_{g,j+1/2}=\frac{\mathsf{\Delta }{y}_{j+1/2}\mathsf{\Delta }{z}_{j+1/2}}{\mathsf{\Delta }{x}_{j+1/2}}{(\frac{k_g{k}_{rg}}{B_g{\mu }_g})}_{j+1/2}$$

The coefficients on other nodes are deduced in the same way. Without special instructions, the above coefficients are all calculated with physical quantities at time n. For the numerical calculation of pressure in Equation (10), because it is a 7-diagonal large sparse matrix equation system, it can be solved using the iterative method. The commonly used iterative methods are Jacobin iteration, Gauss Seidel iteration, SOR iteration, etc. These numerical methods can be found in the classical book [29].

Both the coal and shale reservoirs contain rich organic matter, which has the characteristics of adsorbing methane gas. A large part of gas is stored in organic matter in the form of adsorbed gas. For example, in shale, adsorbed gas accounts for 20~80% of the total gas storage, and in coal rock, this proportion is higher (usually more than 85%). In the middle and late stages of production, the adsorbed gas has a significant impact on the gas production rate. Therefore, it is necessary to describe the adsorption and desorption of the adsorbed gas. The Langmuir adsorption model is the most frequently used and simplest model to characterize the gas adsorption and desorption in the organic matter.

Generally, q_smfg represents the matrix diffusion source term; that is, because the gas pressure in the fracture decreases, the gas in the matrix is desorbed from the matrix surface and diffused into the fracture. The existence of q_smfg can continuously supplement the gas in the fracture, so it is called the source term, whose unit is 1/s, and its definition is

$$q_{smfg}=-F_G\frac{d{\overline{V}}_m}{dt}$$

(11)

where

$$\frac{d{\overline{V}}_m}{dt}=\frac{1}{\tau }(V_E-{\overline{V}}_m)$$

(12)

where ${\overline{V}}_m$ represents the gas concentration in the matrix, dimensionless; V_E is the gas concentration in the fracture, dimensionless; and F_G is the matrix shape factor, which can be taken from 2 to 6. τ is the desorption time. The above formula shows that the greater the concentration difference, the greater the source term. According to the Langmuir desorption model, there are the following relationships

$${\overline{V}}_m=\frac{V_L{P}_{mg}}{P_L+P_{mg}}$$

(13)

$$V_E=\frac{V_L{P}_{fg}}{P_L+P_{fg}}$$

(14)

where the gas concentration ${\overline{V}}_m$ is related to the gas pressure P_mg in the matrix, and the gas concentration V_E is related to the gas pressure P_fg in the fracture where V_L is the Langmuir volume constant. V_L represents the maximum adsorption capacity (saturated adsorption capacity) of gas in unit volume matrix. P_L is the Langmuir pressure constant. If the free gas in the matrix is considered, the free gas term shall be supplemented on the right side of the above equation including φ_m/B_mg and φ_m/B_fg, where B_mg and B_fg represent gas volume coefficients in matrix and fracture, respectively.

In this paper, the Gauss Seidel iterative method is used to solve the 7-diagonal sparse matrix equations, and the iterative convergence must be guaranteed at each time level. Therefore, it is necessary to test the iterative convergence. By comparing the residuals, it is found that the simulation results obtained by dropping the residuals by three orders of magnitude and four orders of magnitude are consistent, which indicates that the simulation on this time layer has converged when dropping the residuals by three orders of magnitude, and the next iteration can be performed.

2.3. Boundary Condition

In the simulation of oil and gas reservoir production, the frequently used boundary conditions include the constant pressure boundary, constant pressure gradient boundary, constant flow rate boundary, and symmetrical boundary. Generally, the pressure gradient boundary condition is used for the outer boundary. The symmetrical boundary condition is used at the symmetrical location. The constant pressure or constant flow boundary condition is used at the wellbore location. In this paper, the zero-pressure gradient boundary condition is adopted for the outer boundary, the symmetric boundary condition is used for the symmetric position, and the constant pressure boundary condition is used for the borehole wall. For the numerical calculation of pressure in Equation (10), because it is a 7-diagonal large sparse matrix equation system, it can be solved using an iterative method. The commonly used iterative methods are Jacobin iteration, Gauss–Seidel iteration, SOR iteration, etc. In this paper, we used the Gauss–Seidel iteration to ensure the stability and convergence of calculation when adopting the static pressure boundary condition.

3. Validation

3.1. Basic Assumptions

For the reservoir simulation, there are many flow models, such as the double medium model and the triple medium model [30,31]. The most representative triple medium model is the trilinear seepage model. This paper introduces the traditional triple media model theory by taking the trilinear seepage model as an example. The purpose of this introduction is that this paper develops a numerical simulation program based on the trilinear model theory, which can be used for the reservoir simulation and gas production assessment. The premise assumption adopted by the trilinear seepage model is that after the large-scale hydraulic fracturing, as shown in Figure 3, the reservoir in the near well area is divided into several cuboid areas due to the existence of artificial fractures, which are called the inner stimulation region (referred to as the inner region). The region outside the stimulation area is called the original reservoir region (referred to as the outer region) because it is not subject to the stimulation effect of artificial fracturing. The specific assumptions of the trilinear flow model include [13,30,31,32]:

(1) The three regions of the reservoir are isothermal, homogeneous, uniform in thickness, and closed in outer boundary;

(2) The gas flow in the reservoir matrix and fractures is single-phase and unsteady;

(3) The adsorption desorption process of gas in the matrix obeys the Langmuir theory, and the diffusion process satisfies Fick’s first law;

(4) The inner zone contains matrix and micro fractures (including natural fractures and secondary fractures caused by fracturing), and the matrix and micro fracture networks are regarded as an equivalent medium;

(5) The outer zone contains matrix and microcrack (only natural fractures), and the matrix and microcrack networks are treated as an equivalent medium;

(6) The artificial fracture area does not contain a matrix, and there are no desorption or adsorption processes. The gas flow in the artificial fracture region conforms to the Darcy seepage, and the gas flow in the outer and inner regions consider the slip effect (i.e., non-Darcy seepage);

(7) Gravity and capillary force are not considered;

(8) The gas in the outer zone first flows into the inner zone, then the gas in the inner zone flows into the artificial fracture, and finally the artificial fracture gas flows into the wellbore.

3.2. Case Validation

In this paper, based on the above model equations and numerical calculation methods, the one-dimensional, two-dimensional, and three-dimensional calculation programs are developed. The reliability of the developed program has been examined and verified by comparing the simulation results with the field production data; for example, by comparing with the Barnett Basin gas production data. The detailed validation can be found in previous research [13]. Here, we give an example to further validate the reliability of the present computation program.

The detailed reservoir parameters are listed in

Table 1. Barnett shale reservoir and well parameters.

Parameter	Value	Parameter	Value
Reservoir half length (m)	550	Initial gas pressure (MPa)	20.34
Reservoir half width (m)	145	Well pressure (MPa)	3.69
Reservoir thickness (m)	90	Langmuir pressure (MPa)	4.48
Artificial fracture half length (m)	70	Langmuir volume	7.2
Artificial fracture width (m)	0.003	Matrix porosity (inner and outer zones)	0.04
Artificial fracture thickness (m)	90	Fracture porosity (inner and outer zones)	0.04
Artificial fracture distance (m)	30.5	Artificial fracture porosity	0.5
Artificial fracture number	28	Permeability in outer zone (mD)	4 × 10−5
Well half length (m)	426.7	Permeability in inner zone (mD)	1 × 10−4
Reservoir temperature (K)	352	Permeability in artificial fracture (mD)	1800

. The main reservoir parameters include the reservoir length, width, thickness, averaged matrix porosity, fracture porosity, initial gas pressure, and so on. Furthermore, the adsorption and desorption are considered in the present simulation. The adsorption and desorption parameters are listed in Table 1.

It is important to conduct the grid independence test. In order to accurately calculate the flow field information, the grid near the artificial fracture should be fine, and the grid in the area far away from the artificial fracture can be sparse, as shown in Figure 4. To ensure the accurate calculation of the gas production rate, it should be ensured that the first grid point near the well should fall in the linear area of pressure distribution; that is, it must be ensured that dp/dy does not depend on the grid selected for calculation. The grid distribution is shown in Figure 4. The X direction grid dependency test shows that the calculation results are consistent between N_x × N_y = 581 × 41 and N_x × N_y = 871 × 41. There is a slight difference (Figure 5) between N_x × N_y = 581 × 41 and N_x × N_y = 291 × 41. Therefore, to ensure that the calculation results do not depend on the number of grids, the number of grids used in the X direction is Nx = 581. Similarly, the grid dependency test is conducted in the Y direction, and it is found that the calculation results are consistent between N_x × N_y = 581 × 41 and N_x × N_y = 581 × 61. There is a slight difference (Figure 5) between N_x × N_y = 581 × 41 and N_x × N_y = 581 × 21. Therefore, we finally select Nx × Ny = 581 × 41 as the computational grid. The minimum grid spacing corresponding to the X direction is dx_min = 0.003 m, and the minimum grid spacing corresponding to the Y direction is dy_min = 0.5 m. In fact, we did not solve the region of the artificial fracture. The vicinity of the artificial fracture was defined as the boundary with constant well pressure. As a result, there is not any grid in the vicinity of the artificial fracture.

For the selection of time step, in this paper, the implicit scheme is used to solve the problem. Compared with the explicit scheme, the implicit scheme can adopt a larger time step. But, even so, the time step size is still affected by the minimum grid size. The smaller the grid size, the smaller the calculation time step size should be. When the time step size is too large, it is difficult to ensure the convergence, which leads to more iterations required for each time step to ensure convergence, and the more time it takes. In order to minimize the total calculation time, it is necessary to balance the relationship between the time step size and the grid size. Finally, the time step selected in the calculation process of the trilinear seepage model is dt = 30 s.

The software platform developed in this paper is used to calculate the gas reservoir productivity with Barnett Basin reservoir data as the input parameters. With the increase in drainage time, the reservoir pressure gradually decreases, and the range of the low-pressure area gradually increases, as shown in Figure 6. With the increase in drainage time, the pressure gradient of the reservoir gas near wells and artificial fractures gradually decreases, which leads to the gradual decrease in the gas production rate. Figure 7 verifies this conclusion. Figure 7 also shows that the gas reservoir productivity calculated based on the present software platform is consistent with the actual field production data, and it is also consistent with the simulation results provided by Yu and Wu [30,32]. Therefore, the simulation program developed by us is reliable.

The dashed line is the previous simulation data, and the diamond is the field data.

3.3. Parameters and Details

In the process of oil and gas production, it is very important to evaluate the production rate. For a large number of studies in the past, the production evaluation formula commonly used is as follows [33]:

$$P_r=A\cdot 2\pi L_w\frac{k_{rg}\sqrt{k_y{k}_z}}{B_g{\mu }_g}\mathsf{\Delta }P$$

(15)

where A = $\frac{1}{\mathit{ln}(r_e/r_w)+S-0.75}$ is related to the equivalent supply radius r_e, well radius r_w, and skin factor S. Usually, A = 0.5~1, and we choose A = 0.7 in the present study.

It is worth noting that for the horizontal well, L_w is the discrete grid length dx in the above formula, and then the production is accumulated along the well length. For vertical wells, L_w is directly used in the above equation, without integration. Output is recorded only once. In fact, the above two methods are essentially the same, both of which are used to calculate the total production on the total length of the well.

Fracture permeability is related to many factors, such as fracture opening and proppant size. Generally, the larger the fracture opening, the higher its permeability. When the proppant size is constant, the greater the porosity in the fracture, the greater the fracture permeability. When the porosity in the fracture is constant, the larger the proppant size, the greater the fracture permeability. For the small-size continuous proppant, there are many models for evaluating fracture permeability, such as Osiptsov’s (2017) [20] model $k/r_e^2=0.204{\phi }^{4.58}$. The detailed models for evaluating fracture permeability are summarized in

Table 2. Main permeability models.

Author	Model	Annotation
Carman-Kozeny (1937) [34]	$\frac{k}{r_e^2}=\frac{1}{9{\kappa }_0}\frac{{\phi }^3}{{\left(1-\phi \right)}^2}$	${\kappa }_0$ is the Kozeny constant, taking a value between 4.17 and 5 [18]
Gebart (1992) [35]	$\frac{k}{r_e^2}=C{\left(\sqrt{\frac{1-{\phi }_c}{1-\phi }}-1\right)}^{5/2}$	C$\mathrm{is} \mathrm{a} \mathrm{geometric} \mathrm{factor}, C=16/9\sqrt{2}\mathsf{\pi }$ and ${\phi }_c$ = 1 − π/4 for a square array
Coelho (1997) [36]	$\frac{k}{r_e^2}=0.117{\phi }^{4.57}$	/
Koponent et al. (1998) [37]	$\frac{k}{r_e^2}=\frac{5.55}{e^{10.1(1-\phi )}-1}$	/
Clague et al. (2000) [38]	$\frac{k}{r_e^2}=b_1{\left(\sqrt{\frac{1-{\phi }_c}{1-\phi }}-1\right)}^2·e^{b_2\left(1-\phi \right)}$	b1 and b2 are curve-fitting parameters. For example, b1 = 17.01 and b2 = −5.861 for the sandstone
Nabovati et al. (2009) [21]	$\frac{k}{r_e^2}=C_1{\left(\sqrt{\frac{1-{\phi }_c}{1-\phi }}-1\right)}^{C_2}$	C1 = 0.491, C2 = 2.31
Osiptsov (2017) [20]	$\frac{k}{r_e^2}=0.204{\phi }^{4.58}$	/
Ezzatabadipour and Zahedi (2018) [22]	$\frac{k}{r_e^2}=b_1{\left(\sqrt{\frac{1-{\phi }_c}{1-\phi }}-1\right)}^n·e^{b_2\left(1-\phi \right)}$	For example, n = 1.194, b1= 8.664, and b2 = −2.229 × 10−14 for the sandstone
Li and Huang (2022) [13]	$\frac{k}{D_e^2}=\alpha e^{\beta \phi }\left(1-e^{-s_p/s_f}\right)$	$\alpha =9.18\times {10}^{-5}{\left(\frac{S_p}{S_f}\right)}^{-1.35}$$\beta =7.39{\left(\frac{S_p}{S_f}\right)}^{0.0412}$

Annotation: r_e is the equivalent radius of proppants. φ is the porosity of porous medium or proppant packs. φ_c is the critical porosity below which there is no permeating flow (the percolation threshold) [21]. $D_e=\sqrt{Dh}$ is the equivalent diameter of the large-sized proppant.

. These models are mainly applicable in the case of the small-sized continuous proppant packet. For the problem of the large-sized discontinuous proppant, the traditional permeability model is no longer applicable. For this reason, Li et al. (2022) [13] proposed a new permeability evaluation model $k/D_e^2=\alpha e^{\beta \phi }(1-e^{-s_p/s_f})$, which is suitable for the large-sized proppant fracture system. Here, $D_e^{}=\sqrt{Dh}$ is the equivalent diameter of the large-sized proppant.

In this paper, we mainly discuss the effect of small-sized and large-sized proppant on the permeability. So, it is necessary to compute the fracture permeability for the above two cases. For the small-sized proppant, Osiptsov’s model (2017) [20] is used here. For the large-sized proppant, Li’s model (2022) [13] is adopted. These two models are suitable for the small-sized and large-sized cases, respectively.

Another important problem is how to calculate reservoir permeability from fracture permeability. Generally, the greater the fracture density, the better the reservoir permeability. This paper focuses on the influence of proppant size on the reservoir pressure field and productivity. Therefore, it assumes that the crack density is a fixed value for different proppant types. Therefore, if the fracture permeability is increased by n times due to the change in proppant size, it can be considered that the overall reservoir permeability is also increased by n times.

In this study, the large size and small size are relative. Small size mainly refers to the fact that the size of the proppant is much smaller than the crack opening; that is, the proppant exists in the crack in the form of a continuous package. A large-size proppant means that the proppant size is equivalent to the crack opening, and the proppant exists in the crack in the form of dispersion. Therefore, large size and small size are concepts relative to the crack opening. For example, the proppant with a diameter of 1 mm can be considered a small-size proppant (assuming a crack opening of 10 mm), and the proppant with a diameter of 1 mm can also be considered a large-size proppant (assuming a crack opening of 1 mm). Therefore, the large size and small size are relative concepts. In the following discussion, it is assumed that the crack opening is 1 mm, so the proppant of 0.1 mm can be regarded as the small-size proppant, and the proppant of 1 mm can be regarded as the large-size proppant.

4. Results and Discussion

4.1. Gas Pressure Field and Production Rate around a Horizontal Well

According to the reservoir characteristics, the drilling and exploitation forms widely used at present include the vertical exploitation well and horizontal exploitation well. This paper mainly discusses these two forms. In addition, there will be multi-well production. Different drainage layouts have an important influence on the reservoir pressure field and production rate. This section first discusses the problem of the natural gas production rate under the condition of a single horizontal production well.

The main reservoir parameters include the reservoir length, width, thickness, averaged matrix porosity, fracture porosity, initial gas pressure, and so on. Furthermore, the well length and well pressure are another two key parameters that affect the production rate of gas. The above parameters are given in

Table 3. Simulation parameters of reservoir including horizontal well.

Parameter	Value	Parameter	Value
Reservoir length (m)	100	Initial gas pressure (MPa)	10
Reservoir width (m)	100	Well pressure (MPa)	1
Reservoir thickness (m)	100	Langmuir pressure (MPa)	4.48
Drainage length of well (m)	50	Langmuir volume (m3/m3)	7.2
Average matrix porosity	0.01	Reservoir temperature (K)	352

. Most of these parameters are chosen upon referring to the previous references [13,30,31,32]. The three-dimensional case is conducted based on the above numerical method.

After the hydraulic fracturing, the reservoir forms a fracture network, as shown in Figure 8. In the fracture network, proppants exist, and these proppants support the fracture to avoid fracture closure under the action of ground stress. Without losing generality, it is assumed that the average fracture aperture is 1 mm and the fracture porosity is 0.5. For case 1, the fracture is continuously filled by the small-sized proppants, and its permeability can be assessed using Osiptsov’s model [20]. The results show that the fracture permeability is 8.53 × 10⁴ mD for case 1. For case 2, the fracture is discontinuously filled by the large-sized proppants, and its permeability can be assessed using Li’s model [13]. The results show that the fracture permeability is 2.34 × 10⁶ mD for case 2. The fracture permeability of case 2 is 27.43 times that of case 1. Assuming that the average permeability of reservoir 1 is 0.001 mD, it can be considered that the average permeability of reservoir 2 is 0.02743 mD, as shown in

Table 4. The proppant size and the corresponding permeability for two cases.

Case 1: Small-Sized Proppant	Value	Case 2: Large-Sized Proppant	Value
Size de (mm)	0.2	Size De (mm)	1
Fracture permeability k (mD)	8.53 × 104	Fracture permeability k (mD)	2.34 × 106
Reservoir permeability kav (mD)	0.001	Reservoir permeability kav (mD)	0.02743

Remarks: The permeability of the fracture containing small-sized proppant is derived from Osiptsov’s model [20], where d_e = 2 r_e. The permeability of fracture containing large-sized proppant is obtained from Li’s model [13]. Suppose that the average permeability of the fracturing reservoir is 0.001 mD for the continuous small-size proppant reservoir. The fracture width is 1 mm. The single fracture porosity is 0.5.

.

Before the discussion, it is vital to assess the grid independency. Four kinds of grids are tested, including the coarse grid, N_x × N_y = 41 × 41, medium grid, N_x × N_y = 81 × 81, fine grid, N_x × N_y = 121 × 121, and the finest grid, N_x × N_y = 161 × 161. The grid dependency test shows that the calculation results are consistent between the medium grid and the fine grid. There are only slight differences (Figure 9) between the coarse grid and the medium grid. Therefore, to ensure that the calculation results do not depend on the number of grids, the medium grid is adopted here for the following simulation and discussion.

Another important issue is the choice of time step. If the selection is too small, the calculation efficiency will be reduced. If it is too large, the convergence of the calculation results is difficult to ensure. In this paper, we choose two steps to test, including dt = 0.1 d and dt = 0.01 d. As shown in Figure 10, the results show that the calculation results given by the two different time steps are consistent, which indicates that the convergence solutions can be given by the two kinds of time steps mentioned above. Considering the computational efficiency, dt = 0.1 d is chosen as the final time step in the following simulations.

According to the above reservoir parameters, well parameters, and calculation parameters, the flow process is simulated in the reservoir fracture systems containing small-sized proppant and large-sized proppant, respectively. Partial simulation results are shown in Figure 11. It is seen that the gas pressure in the reservoir decreases gradually with the increase in time (Figure 11 and Figure 12). Nearest the wellbore area, the gas pressure is the lowest. The farther away from the wellbore, the higher the gas pressure. Comparing the reservoirs with different proppant types, it is found that the pressure of reservoirs with large-sized proppant decreases faster, and the corresponding natural gas production rate is higher, which is confirmed by Figure 13. With the increase in production time, the production efficiency decreases gradually, but the gas production rate of the large proppant reservoir is higher.

4.2. Gas Pressure Field and Production Rate around Single and Multiple Vertical Wells

For the vertical well mining, this paper mainly considers the single well, double wells, three wells, and four wells. As shown in Figure 14, the reservoir area range considered in the simulation process is L_x = 200 m, L_y = 200 m. The distance between the wells is 50 m. Detailed reservoir parameters are listed in

Table 5. Simulation parameters of reservoir including vertical well.

Parameter	Value	Parameter	Value
Reservoir length (m)	200	Initial gas pressure (MPa)	10
Reservoir width (m)	200	Well pressure (MPa)	1
Reservoir thickness (m)	10	Langmuir pressure (MPa)	4.48
Drainage length of well (m)	10	Langmuir volume (m3/m3)	7.2
Average matrix porosity	0.01	Well distance (m)	50
Average fracture porosity of reservoir	0.034	Reservoir temperature (K)	352

. The reservoir thickness is 10 m. The initial reservoir pressure is 10 MPa, and the well hole pressure is 1 MPa. The desorption and adsorption are also considered where the Langmuir pressure is 4.48 MPa and the Langmuir volume is 7.2.

Similarly to the horizontal well, assuming that the permeability of the reservoir with small-size proppant is 0.001 mD, it can be deduced that the permeability of the reservoir with large-size proppant is 0.027 mD. The detailed permeability parameters are listed in

Table 6. The proppant size and the corresponding permeability in two cases.

Case 1: Small-Sized Proppant	Value	Case 2: Large-Sized Proppant	Value
Size de (mm)	0.2	Size De (mm)	1
Fracture permeability k (mD)	8.53 × 104	Fracture permeability k (mD)	2.34 × 106
Reservoir permeability kav (mD)	0.001	Reservoir permeability kav (mD)	0.02743

Where d_e = 2r_e.

.

The results are shown in Figure 15. Case 1 represents the reservoir with small-sized proppant and case 2 represents the reservoir with large-sized proppant. It is seen that the average reservoir pressure decreases from 10 MPa to 9.8 MPa for case 1, where the decrease is only 2%. The average reservoir pressure decreases from 10 MPa to 6.7 MPa for case 2, where the decrease is up to 33%. The production curves indicate that the gas production rate is in the range of 50~100 m³/d for case 1. However, the gas production rate is only in the range of 3~6 m³/d for case 2. These results show that the fracture propped by the large-sized proppant can remarkably increase the gas production rate by enhancing the reservoir permeability.

Figure 16 clearly shows the variation in the reservoir pressure during the first 3000 days. Compared to the case 1, the continuous proppant pack model, case 2 provides a faster drop in gas pressure in the reservoir. The change in the pressure field corresponds to the gas production rate curves.

For the reservoir with double wells, the simulation results are shown in Figure 17. Case 1 represents the reservoir with small-sized proppant and case 2 represents the reservoir with large-sized proppant. The results show that the average reservoir pressure decreases from 10 MPa to 9.55 MPa for case 1, where the decrease is only 4.5%. The average reservoir pressure decreases from 10 MPa to 5.3 MPa for case 2, where the decrease is up to 47%. The production curves indicate that the gas production rate is in the range of 80~200 m³/d for case 1. However, the gas production rate is only in the range of 6~12 m³/d for case 2. These results show that the fracture propped by the large-sized proppant can remarkably increase the gas production rate by enhancing the reservoir permeability.

Compared to the single well, the double wells plan gives a faster pressure drop and higher production rate. Figure 18 clearly shows the variation in the reservoir pressure during the first 3000 days. Compared to case 1 of the continuous proppant pack model, case 2 provides a faster pressure drop in the reservoir. The change in the pressure field corresponds to the gas production rate curves. In conclusion, the large-sized proppant plan can provide a more efficient drainage rate.

For the reservoir with three wells, the simulation results are exhibited in Figure 19 and Figure 20. Case 1 represents the reservoir containing small-size proppant and case 2 represents the reservoir containing large-size proppant. The results show that during the first 3000 days, the average reservoir pressure decreases from 10 MPa to 9.3 MPa for case 1, where the decrease is only 7%. Meanwhile, the average reservoir pressure decreases from 10 MPa to 4.4 MPa for case 2, where the decrease is up to 56%. The production curves indicate that the gas production rate is in the range of 100~300 m³/d for case 1. However, the gas production rate is only in the range of 8~16 m³/d for case 2. These results further confirm that the fracture propped by the large-size proppant can remarkably increase the gas production rate by enhancing the reservoir permeability.

Compared to the single well and double wells, the three wells plan provides a faster pressure drop and higher production rate. Figure 20 clearly shows the variation in the reservoir pressure during the first 3000 days. Compared to case 1 of the continuous proppant pack model, case 2 provides a faster pressure drop in the reservoir. The change in pressure field corresponds to the gas production rate curves. The faster the pressure drops, the higher the gas production rate. In conclusion, the discontinuous large-size proppant plan can provide a more efficient drainage rate.

For the reservoir with four wells, the simulation results are exhibited in Figure 21 and Figure 22. Case 1 represents the reservoir containing the continuous small-size proppant and case 2 represents the reservoir containing the discontinuous large-size proppant. The results show that during the first 3000 days, the average reservoir pressure decreases from 10 MPa to 9 MPa for case 1, where the decrease is only 10%. Meanwhile, the average reservoir pressure decreases from 10 MPa to 3.8 MPa for case 2, where the decrease is up to 62%. The production curves indicate that the gas production rate is in the range of 100~400 m³/d for case 1. However, the gas production rate is only in the range of 10~20 m³/d for case 2. These results further confirm that the fracture propped by the large-size proppant can remarkably increase the gas production rate by enhancing the reservoir permeability.

Compared to the single well, double wells, and three wells, the four wells plan provides the fastest pressure drop and highest production rate. Figure 22 clearly shows the variation in the reservoir pressure during the first 3000 days. Compared to case 1, case 2 provides a faster pressure drop in the reservoir. The pressure field variation affects the gas production rate curves. The faster the pressure drops, the higher the gas production rate. In conclusion, the discontinuous large-size proppant plan can provide a more efficient drainage rate.

The number of wells has an important influence on reservoir pressure and gas production rate. Figure 23 shows that with the increase in the number of wells, the average gas pressure of the reservoir decreases, and the gas production rate increases. Comparing case 1 and case 2, it can be clearly seen that the gas pressure of reservoir 2 is lower and the gas production rate is higher. These results indicate that the discontinuous large-sized proppant can give a higher production rate when the other reservoir parameters are the same between reservoir 1 and reservoir 2.

Under the condition of same porosity, the permeability of fracturing reservoir 1 containing large-sized proppant is higher, causing the gas production rate to be higher. Taking a single well as an example, the natural gas production rate is increased by 19.24 times (from 5.546 m³/d to 106.7 m³/d) at the initial stage of production (t = 1 d) after the use of large-sized proppant. In the later period of exploitation (t = 3000 d), the natural gas exploitation rate is increased by 17.01 times (from 3.037 m³/d to 51.65 m³/d) when adopting the large-sized proppant. Under the same conditions, increasing the number of production wells can significantly improve the drainage efficiency of natural gas. However, the increased times of the production efficiency will gradually decrease with time. Taking the reservoir containing large-sized proppant as an example, the production rate of natural gas will increase by four times (from 106.7 m³/d to 426.7 m³/d) at the initial stage of production (t = 1 d). In the later period of exploitation (t = 3000 d), the natural gas exploitation efficiency is increased by 2.12 times (from 51.65 m³/d to 106.37 m³/d) when the well number increases from one to four.

4.3. Discussion

It is valuable to evaluate the enhancement of fracture permeability when adopting the large-size proppant instead of the small-sized proppant. For the fracture containing the large-sized cylindrical proppant, the definition of porosity is

$$\phi =1-\frac{\pi D^2}{4L_p^2}$$

(16)

where D is the proppant diameter and L_p is the fracture element length. Then, the ratio of the wetted area is written as

$$\frac{S_p}{S_f}=\frac{\pi Dh}{2L_p^2}=\frac{2h}{D}(1-\phi )$$

(17)

where S_p is the wetted area of the proppant, S_f is the fracture element area, and h is the proppant height. Considering a particular situation of D = h, the ratio of the wetted area can be further written as

$$\frac{S_p}{S_f}=2(1-\phi )$$

(18)

Then, the permeability proposed by Li et al. (2022) [13] can be rewritten as

$$\begin{array}{l}k=D^2\alpha e^{\beta \phi }(1-e^{-2(1-\phi )})\\ \alpha =9.18\times {10}^{-5}{(2-2\phi )}^{-1.35}\\ \beta =7.39{(2-2\phi )}^{0.0412}\end{array}$$

(19)

For the fracture containing the small-size proppant, the fracture permeability can be assessed using the power-law model proposed by Osiptsov (2017) [20], as follows:

$$k=0.204r_e^2{\phi }^{4.58}$$

(20)

where r_e is the equivalent radius of the small-sized proppant. Based on the analysis, we obtain the fracture permeability magnification when increasing the proppant size. The assessment model is given as

$$k_r=20{\left(\frac{D}{d}\right)}^2\frac{\alpha e^{\beta \phi }\left(1-e^{-2\left(1-\phi \right)}\right)}{{\phi }^{4.58}}$$

(21)

where d = 2r_e is the equivalent diameter of the small-sized proppant and D is the diameter of the large-sized proppant. It is noteworthy that Formula (21) is based on the premise of D = h. Nonetheless, the Formula (21) can be used to assess the permeability magnification of the fracture.

The detailed magnification curves are exhibited in Figure 24. It is seen that when D/h is constant, the permeability magnification increases with the increase in porosity. When the porosity is constant, the permeability magnification also increases with D/h.

In fact, Formula (21) can also be approximately used to assess the permeability enhancement of the reservoir. Objectively, the reservoir after fracturing contains many fractures. The magnification of fracture permeability can be mapped to the enhancement of reservoir permeability. Therefore, Formula (21) can be chosen as an approximate model to assess the potential of permeability enhancement of the reservoir when adopting the large-sized proppant.

5. Conclusions

This paper analyzes the influence of two typical proppant types on fracture and reservoir permeability and then discusses the differences in reservoir pressure field and production rate between reservoir 1, containing small-sized proppant, and reservoir 2, containing large-sized proppant. The main findings are concluded as follows.

(1) On the premise of the same porosity, the permeability of the fracturing reservoir 1 containing large proppant is higher, which makes the average pressure of the reservoir drop faster, and the natural gas production efficiency is higher. Taking a single well as an example, the natural gas production efficiency is increased by 19.24 times (from 5.546 m³/d to 106.7 m³/d) at the initial stage of production (t = 1 d) after the use of large-sized proppant. In the later period of exploitation (t = 3000 d), the natural gas exploitation efficiency is increased by 17.01 times (from 3.037 m³/d to 51.65 m³/d) when adopting the large-size proppant.

(2) Under the same conditions, increasing the number of production wells can significantly accelerate the reduction of reservoir pressure and significantly improve the efficiency of natural gas production. However, the increased times of the production efficiency will gradually decrease with time. Taking the reservoir containing large-size proppant as an example, the production efficiency of natural gas will increase by four times (from 106.7 m³/d to 426.7 m³/d) at the initial stage of production (t = 1 d). In the later period of exploitation (t = 3000 d), the natural gas exploitation efficiency is increased by 2.12 times (from 51.65 m³/d to 106.37 m³/d) when the well number increases from one to four.

(3) A new evaluation model is proposed to quantitatively evaluate the permeability magnification of the fracture and reservoir when adopting the large-sized proppant. The present model shows that the fracture permeability can be increased 27~240 times in the porosity range of 0.5 ≤ φ ≤ 0.9 when using the discontinuous proppant. The present finding further confirms that there is a huge potential for applying the discontinuous large-sized proppant to enhance the reservoir permeability.