#### Abstract

Low-velocity non-Darcy flow can be described by using the threshold pressure gradient (TPG) in low-permeability porous media. The existence of the TPG yields a moving boundary so that fluid starts to flow inside this boundary when the pressure gradient overcomes the viscous forces, and beyond this boundary, there will be no flow. A mathematical model of considering the TPG is developed to describe the flow mechanism in multiple-porosity media. By defining new dimensionless variables, the nonlinear mathematical model can be solved analytically. This new approach has been validated with several approximate formulas and numerical tools. The diffusion of the moving boundary varying with time is analyzed in detail in multiple-porosity media, and then the effect of the moving boundary on pressure transient response is investigated and compared with that of the traditional three boundary types (closed boundary, infinite-pressure boundary, and constant-pressure boundary). Sensitivity analysis is conducted to study the effect of the TPG on pressure and pressure derivative curves and rate decline curves for single-porosity media, dual-porosity media, and triple-porosity media, respectively. The results show that the moving boundary exerts a significant influence on reservoir performance at a relatively early time, unlike the other three boundary types, and only a boundary-dominated effect at the late time. The larger the threshold pressure gradient, the smaller the diffusion distance of the moving boundary and the rate of this well at a given dimensionless time. At the same time, the pressure transient response exhibits a higher upward trend because of a larger TPG. All behavior response might be explained by more pressure drop consumed in low-permeability reservoirs. The finding is helpful to understand the performance of low-permeability multiple-porosity media and guide the reasonable development of low-permeability reservoirs.

#### 1. Introduction

For flow in porous media, Darcy’s law is a fundamental relationship to describe mathematically the fluid flow through porous media. It is well known Darcy flow assumes that the flow velocity-pressure gradient relationship is a straight line which passes through the origin. Recently, many non-Darcy flow phenomena are found in low-permeability media [1, 2], tight carbonate oil and gas reservoirs [3], shale oil and gas reservoirs [4, 5], and clay soil [6]. Non-Darcy flow indicates that the flow velocity-pressure gradient relationship is not a straight line, also called nonlinear flow in porous media. Kutilek [7] summarized 12 different types of non-Darcy flow, which can be divided into two categories: high-velocity non-Darcy flow and low-velocity non-Darcy flow. At present, the low-velocity non-Darcy flow phenomenon has aroused wide attention, especially in tight oil and gas reservoirs.

The mechanism of non-Darcy flow is very complex because of the complicated pore structure and solid-liquid interaction. In recent years, many laboratory experiments have been conducted to study the low-velocity non-Darcy flow of oil and gas in low-permeability cores. Prada and Civan [8] pointed out that fluids can flow through porous media only if the fluid force is sufficient to overcome the frictional effects, and therefore, Darcy’s law should be corrected for the effect of the threshold pressure gradient (TPG). Hao et al. [1] discovered that different types of fluids gave a different TPG versus permeability power function and that the TPG for the two-phase oil and water was greater than that for the single-phase flow by a core displacement experiment. Zeng et al. [2] used 23 natural cores taken from an ultralow-permeability reservoir in Western China to perform a laboratory investigation and concluded that the flow characteristic in ultralow-permeability cores does not follow Darcy’s law and minimum threshold pressure gradient and pseudothreshold pressure gradient exist widely. Yao and Ge [9] conducted a number of core tests on low-permeability reservoir rocks to determine the non-Darcy flow equation and critical Reynolds number. When the Reynolds number (Re) is less than the critical Reynolds number, the non-Darcy seepage index obtained from core flow tests was 0.3987; in contrast, when the Re is larger than the critical Reynolds number, the non-Darcy seepage index was 0, i.e., Darcy flow region. According to the above-mentioned experiment results, the whole flow velocity-pressure gradient is a combination of a straight line and a concave curve. Xiong et al. [10] pointed out that three zones exist in the low-permeability reservoirs: dead oil zone, nonlinear seepage zone, and pseudolinear seepage zone. When the pressure gradient in a zone is smaller than the minimum threshold pressure gradient, the zone is called the dead oil zone; when the pressure gradient in a zone is smaller than the maximum pressure gradient but larger than the minimum threshold pressure gradient, the zone is called the nonlinear seepage zone; when the pressure gradient in a zone is larger than the maximum pressure gradient, the zone is called the pseudolinear seepage zone.

Some mathematical models have been proposed to describe the low-velocity non-Darcy flow, which can be separated into three categories: the quasilinear flow model based on the pseudothreshold pressure gradient [11], the two-parameter flow model based on the boundary layer theory [12], and the nonlinear flow model based on the maximum, minimum, and pseudothreshold pressure gradients [13]. A major challenge in the calculation of the nonlinear mathematical model is to solve accurately the diffusion of the moving boundary. Because of using the TPG to describe the low-velocity non-Darcy flow, a moving boundary front would appear [14]. As a result, the mathematical model becomes a nonlinear equation that many researchers try to solve using different methods. In one-dimensional flow problems in porous media with a moving boundary, it is similar to the classical Stenfan problem of heat conduction. Pascal et al. [15, 16] investigated the effect of the threshold gradient on the pressure and flow rate distributions in the nonsteady flow through porous media by the integral method. Song et al. [17] calculated control radius, control distance, control area, and control coefficient for ultralow-permeability reservoirs based on ellipse flow theory. Liu et al. [18, 19] used similarity transformation to obtain exact analytical solutions. Wang et al. [20] applied Green’s function to obtain exact analytical solutions. In the well test, some approximation methods were used, such as the integral method [21], steady-state successive approach [22], Green’s function method [23], and advanced mathematical methods [24]. In addition, for tight carbonate reservoirs, dual-porosity and triple-porosity models with a moving boundary have been developed [25–32]. Numerical methods, such as the finite difference method [33–36], also provide a powerful tool to calculate moving boundary problems. A nonlinear seepage numerical simulation software has been developed to study oil and water two-phase flow in low permeability [37, 38].

The aim of this article is to present the results of the diffusion of the moving boundary varying with time and focus on effects of the TPG on single well production performance in multiple-porosity media based on low-velocity non-Darcy flow. We firstly develop mathematical models considering the interaction between the TPG and the moving boundary. Secondly, the location of the moving boundary front is determined. Thirdly, sensitive studies of pressure and pressure derivate curves are analyzed. Finally, rate decline curves are plotted.

#### 2. Mathematical Model

##### 2.1. Physical Model

The carbonate reservoir in Pre-Caspian Basin is a typical example of reservoirs composed of matrix, fracture, and vug systems. A pore-type classification system has been developed and is shown in Figure 1 [39]. This triangular diagram is closely linked to sedimentological fabric and pore types. But it is not directly linked to flow properties. In the oil and gas reservoir engineering field, a multiple-porosity media classification system is most widely used as follows: single-porosity, dual-porosity, and triple-porosity media. Table 1 gives the core pictures and thins of single-porosity, dual-porosity, and triple-porosity media. Single-porosity media are mainly composed of the matrix system including pore and throat, dual-porosity media consist of matrix and fracture systems, and triple-porosity media are mainly composed of matrix, fracture, and vug systems.

Figure 2 shows the physical model scheme of multiple-porosity media. Fractures and vugs are uniformly distributed in a reservoir, and the media directly connected with a wellbore are the fracture system for dual-porosity and triple-porosity media. Interporosity flows from the matrix system to the fracture system and from the vug system to the fracture system are assumed.

**(a)**

**(b)**

**(c)**

Physical model assumptions are as follows:(1)A single well production located a center in a cylindrical reservoir by a constant rate or a constant pressure; the external boundary of the reservoir may be closed boundary, infinite-pressure boundary, constant-pressure boundary, and moving boundary proposed in this article(2)Homogeneous and isotropic reservoirs with constant physical properties, such as permeability, porosity, reservoir thickness, oil viscosity, and compressibility(3)Isothermal and non-Darcy flows, ignoring gravity and capillary forces(4)Wellbore storage and skin effect are considered

##### 2.2. Non-Darcy Flow

Figure 3 shows the typical velocity-pressure gradient relationship curves for Darcy flow and non-Darcy flow. It is clear that the Darcy flow curve is a straight line going through the origin, while the non-Darcy flow curve is a concave curve at the small pressure gradient. If the short segment of the concave curve is ignored for engineering calculation, the flow curve can be simplified into a straight line not going through the origin, with an intercept of TPG.

The most common form of the corrected Darcy’s law for describing the low-velocity non-Darcy flow is

Initially, the reservoir pressure is constant everywhere. With the flowing time increase, the following conditions at the moving boundary front are used:

Note that is controlled by the TPG, and the fluid that is beyond the moving boundary cannot flow.

##### 2.3. Mathematical Model

The dimensionless variables are defined as follows.

Dimensionless pressure:

Dimensionless wellbore pressure:

Dimensionless rate:

Dimensionless time:

Dimensionless radius:

Dimensionless radius of the external boundary:

Dimensionless radius of the moving boundary:

Dimensionless wellbore storage coefficient:

Dimensionless threshold pressure gradient:

With the above dimensionless parameters, the radial transient flow governing equations for multiple-porosity media can be written in the Laplace domain as

Initial condition:

Inner boundary condition:

The outer boundary conditions are as follows.

Infinite-pressure boundary:

Closed boundary:

Constant-pressure boundary:

Moving boundary:where

Substituting the initial condition from equation (13) into equation (12), we obtain

The general solution to equation (20) can be given aswhere

Considering the inner boundary condition from equation (14) and four different outer boundary conditions from equations (15)–(18), the final solution can be obtained for all outer boundary conditions, which is shown in Table 2.

Especially, substituting the dimensionless pressure solutions of the moving boundary into equation (18), the movement equation of the moving boundary can be given as

According to Wronskian’s formula for modified Bessel functions,

Reducing equation (24), we get

Equations (21) and (25) are solved using the iterative method in the Laplace domain. Then, the dimensionless pressure in the real domain can be obtained by Stehfest’s numerical inversion [40]. The calculation procedure is shown in Figure 4.

Applying the Duhamel theorem in the Laplace domain, the dimensionless rate can be given aswhere

#### 3. Comparison and Verification

##### 3.1. Comparison of Our Moving Boundary with Approximate Formulas

Because of the nonlinearity of the mathematical model with the TPG, the exact analytical solutions for moving boundary problems are not easy to obtain. Several moving boundary front formulas have been proposed by some researchers to describe approximately the behavior of the moving boundary varying with time [41]. The following five formulas were compared with the solution from this work.

K-C:

Liu C-Q:

Wu Y-S:

Song F-Q:

Wang X-D:

Comparison of the results is presented in Figure 5 when the TPG is equal to 0.01. It can be seen that the calculated result presented in this study is larger than that of Wu Y-S but less than other four results at the large dimensionless time. However, our result is the biggest compared to other five solutions at the small dimensionless time. In our model, the interaction of the TPG and moving boundary was considered, which indicated the pressure did not instantaneously arrive at the infinite boundary. Unlike the other solutions based on various assumptions and solved by the approximate method, our solution was strictly solved on the basis of the derived mathematical model. So it is believed that the proposed result in this work is found to be reasonably accurate and with a high degree of precision.

##### 3.2. Comparison of Our Method with Numerical Simulation

In order to validate further the results, the numerical simulation was compared with our solution. A vertical well was located in the center of a reservoir with a reservoir radius of 100 m. The production time is 30 days with a constant oil rate of 3.5 m^{3} per day. Table 3 shows the basic reservoir and fluid property data. Two cases were considered: without TPG and with a TPG of 65 kPa/m assumed to exist in the low-permeability reservoir.

Figure 6 shows the comparison of formation pressure from our method with that from the numerical method. The calculated pressure in our model is in excellent consistency with the numerical method regardless of TPG existence. As the drainage radius increases, the formation pressure increases rapidly and then becomes relatively flat. Moreover, the pressure curve without TPG is higher than that with a TPG of 65 kPa/m. The results indicate a more difficult fluid flow when the TPG exists. In other words, a low-permeability reservoir needs more energy to produce. Because of the effects of grid size and gravity, the results of the numerical method are slightly smaller compared to our solutions.

Figure 7 compares the pressure distribution with and without TPG by the numerical simulation. The figure indicates that the existence of the TPG increases the pressure decline near the wellbore and that at the same distance, the pressure gradient when considering the TPG will be larger compared to that of the case without TPG. This is easy to explain that maintaining the same constant rate, additional reservoir energy will be consumed because of the existence of the TPG. Figure 8 shows moving boundary diffusion with time. Obviously, the moving boundary is gradually expanded with time. That is to say, the outside of the moving boundary will keep the initial formation pressure of 25 MPa. Only inside the moving boundary, the fluid can flow. Unlike the case without TPG, its pressure quickly reaches the outermost boundary.

**(a)**

**(b)**

#### 4. Moving Boundary Analysis

For dimensionless threshold pressure gradients equal to the values of 0.001, 0.01, and 0.1, the results of the moving boundary at different dimensionless time are plotted for single-porosity media in Figure 9. It is observed that, at a given dimensionless time *t*_{D}, the smaller the dimensionless threshold pressure gradient, the larger the dimensionless radius of the moving boundary *r*_{fD}. A smaller TPG means a smaller resistance to flow. Thus, controlling radius *r*_{fD} of the well is larger, meaning that the moving boundary front is further away from the wellbore. Moreover, the change of the distance of the moving boundary with time is almost parallel at the middle and late time.

For dual-porosity media, the effect of the TPG on the moving boundary is presented in Figure 10. Because there is an interporosity flow from the matrix system to the fracture system in dual-porosity media, these curves exhibit flat straight line periods at the middle time. This means that the moving boundary stops temporarily when the interporosity flow occurs. It is a similar trend that as the TPG increases, the curve of the moving boundary is larger.

The effect of the TPG on the moving boundary is also shown for triple-porosity media in Figure 11. The curve of the moving boundary exhibits two relative flat straight line periods at the middle time. This is because the first flat straight line period is related to the interporosity flow from the vug system to the fracture system, and the second flat straight line period corresponding to the interporosity flow is between the matrix system and the fracture system. The moving boundary shows a temporary stop when these two interporosity flows occur. Similarly, these curves of the moving boundary are almost parallel.

#### 5. Pressure Transient Analysis

##### 5.1. Effects of Moving Boundary on Pressure Transient Response

Figure 12 shows the pressure transient response under four different types of outer boundaries in single-porosity media. We can see that the pressure derivative curve with the moving boundary has an upward trend from a relatively early time compared to the other three boundaries. In addition, the slope of pressure derivative for the moving boundary is smaller than that for the closed boundary. This result shows that, unlike traditional boundary types (closed boundary, infinite-pressure boundary, and constant-pressure boundary) only influencing the flow period at the latter stage, the moving boundary can occur from the middle to the latter stage. From the above study, we have understood that the existence of the TPG in low-permeability media leads to the propagation of the moving boundary. Therefore, the existence of the moving boundary means more pressure loss and thus more difficulty for the flow.

For dual-porosity media, the pressure transient response under four different types of boundaries is presented in Figure 13. As we all know, the pressure derivative has a characteristic V-shaped intermediate segment in dual-porosity media, which suggests an interporosity flow from the matrix system to the fracture system. If the moving boundary is considered, the derivative curve will be a deviation after a V-shaped interporosity flow. The latter part is similar to the single-porosity media.

Similarly, the pressure derivative has a characteristic VV-shaped intermediate segment in triple-porosity media presented in Figure 14. The first V-shape represents an interporosity flow from the vug system to the fracture system. The second V-shape means an interporosity flow from the matrix system to the fracture system. The moving boundary will lead to an upward trend from an intermediate time compared to the other three boundary types.

Generally, the non-Darcy transient flow with the moving boundary and TPG shows the typical feature that pressure derivative curves exhibit a straight line of slope that is larger than zero after a wellbore flow regime of early time.

##### 5.2. Effects of TPG on Pressure Transient Response

In order to investigate the effect of the TPG on the pressure transient response for multiple-porosity media, Figures 15–17 present pressure transient response under different TPG values. Figure 15 shows the pressure and pressure derivate curves under a TPG of 0.001, 0.01, and 0.1, respectively, for single-porosity media, which demonstrates that a larger TPG leads to a higher upward trend of pressure and pressure derivate curves. Moreover, a larger TPG causes a deviation earlier. This can be explained from the pressure loss view. A lower permeability means a larger TPG and a more difficult fluid flow. So the flow needs to overcome more resistance. At the same time, we can also see that different pressure derivative curves are almost parallel at the latter stage. This suggests that the flow speed changes as a whole when the value of the TPG is different.

Similarly, as the TPG increases, pressure derivate curves move upward and higher in dual-porosity media, as shown in Figure 16. The V-shaped intermediate segment also moves upward. It is easy to observe that the pressure curve is a flat straight line period when the interporosity flow between the matrix system and the fracture system occurs. When the TPG is smaller than 0.1, both of the pressure curves are almost the same. However, the pressure curve jumps a step when the TPG is equal to 0.1. This reduction order of magnitude for permeability needs consumption of more energy.

Now it is easy to understand the behavior of pressure and pressure derivative for triple-porosity media plotted in Figure 17. The effect of the TPG on pressure and derivative curves shows that a larger TPG value leads to a larger pressure change and a larger pressure derivative. It can be explained that a larger TPG causes a smaller flow speed, thus leading to a larger pressure loss and larger pressure derivative. We can also observe that the first V-shape of pressure derivative curves is corresponding to the flat straight line period of pressure curves and that the second V-shape of pressure derivative curves relates to the transition period of pressure curves, which are consistent with the propagation of the moving boundary.

#### 6. Rate Transient Analysis

Under the constant-pressure production condition, dimensionless rate decline curves are plotted in Figures 18–20 based on the previous mathematical model. Figure 18 shows the effect of the TPG on rate decline curves in single-porosity media. As the dimensionless time increases, the dimensionless rate declines rapidly except the case without TPG at the late time. It can also be found that, at a given dimensionless time, the larger the dimensionless threshold pressure gradient is, the smaller the dimensionless rate is. This result is a totally consistent trend with the above-mentioned curves of the moving boundary and pressure transient response for single-porosity media.

It is also easy to observe that the dimensionless rate decline curve has a flat straight line at the middle time for dual-porosity media shown in Figure 19. This means that the rate of a well can be kept at a relatively stable status when the interporosity flow between the matrix system and the fracture system occurs. Similarly, the dimensionless rate decline curve has two relatively flat straight line periods for triple-porosity media at the middle time presented in Figure 20. This is the same reason that two interporosity flows start to play a role. These results prove further that the TPG indeed affects the production and fluid flow. The non-Darcy transient flow with the TPG should be emphasized.

#### 7. Conclusions

A single-phase flow mathematical model incorporating the TPG was developed to describe the non-Darcy flow behavior in multiple-porosity media. The Laplace transform and Stehfest numerical inversion based on dimensionless variables were applied to study the moving boundary, pressure transient response, and rate transient decline for the low-permeability porous media. Sensitivity analysis was conducted to investigate the effect of the TPG on the moving boundary, pressure transient response, and rate transient decline for single-porosity media, dual-porosity media, and triple-porosity media, respectively. The following conclusions can be drawn:(1)The moving boundary varies with time because of the existence of the TPG. Unlike traditional boundary types (closed boundary, infinite-pressure boundary, and constant-pressure boundary), they are immobile and only flow dominated at the late time. The moving boundary diffuses slowly corresponding to larger TPGs.(2)It is shown that the larger the dimensionless TPG is, the more the pressure drop is required to maintain a constant rate production. The curves of dimensionless pressure derivative have an upward characteristic at a relatively early time compared to the other three outer boundaries.(3)The dimensionless rate under the constant-pressure production declines more severely because of larger TPGs. This result is related to the more pressure drop required because of larger TPGs.(4)For dual-porosity and triple-porosity media, there exists a similar characteristic of the flat straight line on the moving boundary and dimensionless rate decline curves and the V-shape on dimensionless pressure derivative response when interporosity flow occurs.

#### Nomenclature

: | Flow velocity (10^{−3} m/s) |

μ: | Fluid viscosity (mPa·s) |

k: | Permeability (10^{−3} μm^{2}) |

: | Pressure gradient (MPa/m) |

G: | Threshold pressure gradient (MPa/m) |

P: | Formation pressure (MPa) |

P_{i}: | Initial formation pressure (MPa) |

r: | Radial cylindrical coordinate (m) |

t: | Time (h) |

: | Moving boundary (m) |

P_{j}: | Formation pressure (MPa) |

h: | Reservoir thickness (m) |

q: | Rate (m^{3}/d) |

B: | Formation volume factor (m^{3}/m^{3}) |

: | Wellbore pressure (MPa) |

φ: | Porosity (fraction) |

c_{t}: | Total compressibility (1/MPa) |

: | Well radius (m) |

r_{e}: | Outer boundary distance (m) |

C: | Wellbore storage coefficient (m^{3}/MPa) |

u: | Laplace transform variable |

ω: | Fluid storage capacitance coefficient |

λ: | Interporosity flow factor |

S: | Skin factor (fraction). |

#### Special Functions

: | The first-type derivative Bessel function of order j (j = 1, 2) |

: | The second-type derivative Bessel function of order j (j = 1, 2). |

#### Subscripts

D: | Dimensionless |

m: | Matrix system |

f: | Fracture system |

: | Vug system |

j: | Number. |

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.