Nonlinear Moving Boundary Model of Low-Permeability Reservoir

: At present, there are two main methods for solving oil and gas seepage equations: analytical and numerical methods. In most cases, it is difﬁcult to ﬁnd the analytical solution, and the numerical solution process is complex with limited accuracy. Based on the mass conservation equation and the steady-state sequential substitution method, the moving boundary nonlinear equations of radial ﬂow under different outer boundary conditions are derived. The quasi-Newton method is used to solve the nonlinear equations. The solutions of the nonlinear equations with an inﬁnite outer boundary, constant pressure outer boundary and closed outer boundary are compared with the analytical solutions. The calculation results show that it is reliable to solve the oil-gas seepage equation with the moving boundary nonlinear equation. To deal with the difﬁculty in solving analytical solutions for low-permeability reservoirs and numerical solutions of moving boundaries, a quasi-linear model and a nonlinear moving boundary model were proposed based on the characteristics of low-permeability reservoirs. The production decline curve chart of the quasi-linear model and the recovery factor calculation chart were drawn, and the sweep radius calculation formula was also established. The research results can provide a theoretical reference for the policy-making of development technology in low-permeability reservoirs.


Introduction
Classical Darcy's law [1] has been widely used to describe seepage mechanics and numerical simulation for a long time. However, in recent years, the experimental monitoring data of pore seepage in tight oil and gas reservoirs, such as shale oil and gas, low-permeability oil and gas reservoirs and weak permeable layers, confirm that the fluid seepage in low-permeability tight media does not follow Darcy's law, and the seepage velocity and pressure gradient reflected from the seepage characteristic curve is not linear. The study of nonlinear seepage in tight media has important guiding significance and application value for the effective development of tight oil and gas reservoirs.
Bear [2] proposed the concept of starting pressure gradient and gave the non-Darcy seepage equation of low-permeability formation during the process of studying the seepage mechanism of low-permeability formation through a laboratory test with a one-dimensional core in 1972. Prada [3] and Civan [3,4] found that the linear fitting curve of seepage velocity and pressure gradient experimental data had obviously deviated from the origin in the experimental study of rock fluid flow in different low-permeability porous media, and they developed the quasi-linear motion equation by modifying the Darcy equation.
When J ≤ J 0 , q = 0; and J > J 0 , fluid flow occurs in the porous media.
where λ is the start-up pressure gradient, MPa/cm; τ o is the ultimate shear stress, N; K is the core air permeability, 10 −3 µm 2 ; φ is the porosity. For the convenience of calculation, many experts and scholars have conducted laboratory experiments to obtain an empirical formula for calculating the starting pressure gradient. Chengyuan Lu [8] selected sandstone reservoir cores from 10 oilfields in the Shengli Oilfield in 2002 to represent low-permeability sandstone reservoirs in the Shengli Oilfield. The start-up pressure gradient of single-phase oil phase percolation and the air permeability of the sample were measured by the capillary balance method, and the empirical formula for the start-up pressure gradient was established by the linear regression method, which would be a better way to meet the requirements of engineering calculations.
The commercial reservoir numerical simulator represented by Eclipse is developed based on the linear Darcy flow law, which cannot accurately describe the nonlinear flow characteristics existing in the developing process of low-permeability reservoirs. Even with the simplest quasi-linear non-Darcy flow model, the common reservoir numerical simulators can produce significant errors in the optimization design of reservoir development parameters, such as well pattern and spacing and recovery factors. For low-permeability reservoirs, when the oil and water seepage occurs in the formation during production, the area where the formation pressure gradient is greater than the start-up pressure gradient is the pressure drop sweep area; the area where the formation pressure gradient is less than the starting pressure gradient is the non-sweep area; and the formation pressure in this area remains as the original formation pressure. On the boundary between the two regions, the formation pressure gradient is equal to the starting pressure gradient. As the seepage process in the reservoir is unstable, the pressure gradient in the low-permeability formation also changes, resulting in the boundary between the two regions moving with time and space. The mathematical expression that characterizes the moving boundary in mathematical modeling is called the moving boundary condition. The model based on the quasi-linear seepage equation is actually a moving boundary problem. In 1981, H. Pascal [9] first studied the one-dimensional unsteady seepage flow boundary model in porous media with a start-up pressure gradient, solved the model by an approximate analytical method and numerical method, and then analyzed the influence of the start-up pressure gradient on pressure and flow distribution in the flow system. In 1992, WuY.S. [10,11] used the law of flow continuity and conservation of mass to obtain the approximate analytical solution of the control partial differential equation in the sense of an average integral using the integration method and established the corresponding reservoir test interpretation method. In 2008, Feng [12] studied the unstable radial seepage model in a low-permeability gas reservoir considering the start-up pressure gradient through the numerical approximate Green function method; the results show that the moving boundary can characterize the control radius of a single well. Chen [13] used pore scale network modeling to study the seepage and displacement of porous media with the starting pressure gradient or yield stress. Lu J. [14,15] presented a boundary-dominated flow model in radial shape geometry and developed the analytical solutions. Yun [16,17] proposed a fractal model of starting pressure gradient in the flow process of non-Newtonian Bingham fluid in porous media. In 2010, Xie K. H [18] obtained a general approximate analytical solution for a one-dimensional clay soil consolidation moving boundary model considering the starting pressure gradient. Liu et al. [19] discussed analytical and numerical solutions of a semi-infinite low permeability reservoir under one-dimensional flow with a moving boundary and TPG. However, their solutions, which are processed by similarity transformations, are complicated and difficult, and the influence of the TPGs is not obvious. Wang Xiao-dong [20] proposed the analytical solutions under fixed-borehole-rate and fixed-borehole-pressure by applying Green's function. Two transcendental equations for the moving boundary model were obtained and solved by using the Newton-Raphson iteration. Huang Y [21] presented a new unsteady production decline model of dual-porosity medium (matrix and fracture) and composite gas reservoir, proposed a complex mathematical solution, and discussed the influence of different reservoir parameters on the production decline regularity of composite reservoirs.
The approximate analytical solution, numerical solution, Bundle of capillary tube modeling (BCTM), direct pore scale modeling (DPSM) and Pore Network Modeling (PNM) are the main research methods and tools; however, the exact analytical solution of a porous media seepage flow boundary model with the starting pressure gradient, the proof of the existence of the solution and the moving law of the moving boundary are rarely studied. Various nonlinear seepage models based on laboratory experiments and theoretical derivation enrich the seepage theory of low-permeability reservoirs, but most models greatly increase the difficulty of calculation and are difficult to use in engineering. From the perspective of engineering calculation and calculation error, the quasi-linear model can better simulate the seepage mechanism of low-permeability reservoirs, and it is the mainstream of seepage theory research on a low-permeability reservoir.
When the steady-state sequential replacement method is used to solve the unstable seepage problems, it often regards each instantaneous state of the unstable seepage process as stable, and this can transform the unstable seepage partial differential equation into an ordinary differential equation without time. Numerous researchers are inclined to apply it to solve the flow radius of different reservoir models.
The quasi-linear seepage model of a low-permeability reservoir is used within the moving boundary, and the fluid flow does not occur outside the moving boundary. Based on this idea, this paper deduces and establishes the nonlinear moving boundary equation based on the steady-state sequential replacement method, which can be accurately and simply used to solve the quasi-linear model of a low-permeability reservoir. By solving the model, the calculation chart of recovery factor, sweep radius and other parameters of low permeability reservoirs under different starting pressure gradients can be obtained. In our model, only a nonlinear formula needs to be solved at each time-step, which achieves the same result with low-permeability reservoir numerical simulation. Therefore, smaller errors according to the analytic solution and fewer calculations compared to a routine numerical simulation method are the two maximum advantages of our model. This research provides a new idea for nonhomogeneous partial differential equations as well. The results can provide a theoretical reference for low-permeability reservoir development technology policy formulation and numerical simulation of low-permeability reservoirs.  Basic assumptions are as follows: 1 The reservoir is isotropic homogeneous and has a circular shape with equal thickness, closed top and bottom, and uniform original formation pressure; 2 The fluid is single-phase and micro-compressible, seepage flow is under constant temperature without physical and chemical changes and satisfies Darcy's law; 3 Permeability, porosity, compressibility, and viscosity, etc., are all constants and do not change with time or pressure; 4 The influence of gravity is not considered, and the pressure gradient is small.

Moving Boundary Nonlinear Equation
Considering a constant pressure production well in the center of a circular homogeneous area, the seepage equation is as follows: Constant pressure inner boundary conditions: Infinite boundary: lim Stable boundary: Closed boundary:

Moving Boundary Nonlinear Equation
Steady state sequential replacement method [22]: when solving a series of problems related to the unstable seepage, each instantaneous state of the unstable seepage process is often regarded as stable. Considering the time period t 1 →t 1 + ∆t, the pressure sweep radius is from r 1 →r 2 (see Figure 2). Steady state pressure equation: The formation pressure P 1 distribution at time t is obtained by an integral solution: The formation pressure P 2 distribution at time t + ∆t can be obtained similarly.
Dimensionless bottom-hole production calculation formula: Among them: The same is shown below. With the increase in r iD , the yield decreased gradually. Mass conservation equation of the stratum in ∆t time: Among them: Meanwhile, Substituting (14) and (15) into (13),we can obtain: Equation (16) gave the initial value of the pressure sweep radius.
(2) Considering the time period t 1 →t 1 + ∆t, the pressure sweep radius is from r 1 →r 2 (require r 2 < r e ).
The formation pressure P 2 distribution at time t 1 →t 1 + ∆t can be obtained similarly.
Dimensionless bottom-hole production calculation formula: The same is shown below.
With the increase in r iD , the yield decreased gradually. Mass conservation equation of the stratum in ∆t time: Among them, Substituting (20) and (21) into (19),we can obtain: Simplify and then obtain the infinite boundary ripple radius r 1D , r 2D calculation formula: Equation (23) is the iterative formula of r 2D , which can be solved by the quasi-Newton method [16].
(3) After the pressure reaches the outer boundary, consider the t 2: +∆t time period, and the pressure at the outer boundary changes from P i1 →P i2 (see Figure 3). Formation pressure P 1 distribution at time t 2: : Distribution of formation pressure P 2 at time t 2: +∆t: The yield calculation formula becomes: The mass conservation equation after the pressure wave reaches the boundary: Substitute (24)-(26) into (27) to solve P i2 (the initial value of P i1 is P i ): Equation (28) gave the outer dimensionless boundary pressure. Equation (26) gave the dimensionless wellbore production. They comprised a nonlinear solution of the moving boundary nonlinear equation.
For the constant pressure outer boundary, if the sweep radius is greater than the outer boundary radius, the sweep radius is taken as the outer boundary radius.
For the closed boundary, if the sweep radius is greater than the outer boundary radius, the outer boundary radius will not change, but the pressure at the outer boundary will change. The oil reservoir is isotropic and homogeneous with circular, uniform thickness, closed top and bottom, and uniform original formation pressure; 2 The fluid is single-phase and micro-compressible, and the seepage flow is under constant temperature without physical and chemical changes and satisfies Darcy's law; permeability, porosity, compressibility, and viscosity, etc., are all constants and do not change with time or pressure; 3 The influence of gravity is ignored, and the pressure gradient is small.
Considering a constant pressure production well in the center of a circular homogeneous area, the flow equation of the flow zone in a low-permeability reservoir is as follows: 1 Internal boundary conditions of constant pressure: Infinite outer boundary: lim Constant pressure outer boundary: Closed outer boundary:

Moving Boundary Nonlinear Equation
The steady-state pressure equation in the flow zone: Integral solution to obtain the formation pressure distribution: p = c 1 ln r + Gr + c 2 (35) c 1 , c 2 is the integral constant. Formation pressure P 1 within the range of moving boundary r 1 at time t: Formation pressure P 2 in the range of moving boundary r 2 at time t 2: +∆t: Dimensionless bottom-hole production calculation formula: Among them: G D is the dimensionless starting pressure gradient, dimensionless, As r increases, the yield gradually decreases. Formation quality conservation equation in ∆t time: Of which: Substituting (40) and (41) into (39), we obtain: After the pressure reaches the outer boundary, consider the time period t 2 →t 2 :+∆t, and the pressure at the outer boundary starts from P i1 →P i2 .
Distribution of formation pressure P 1 at time t 2 : Distribution of formation pressure P 2 at time t 2 :+∆t: The yield calculation formula becomes: When the pressure wave reaches the boundary, the mass conservation equation can be expressed as: Substitute Equations (44)-(46) into Equation (47) to solve for p i2D (p i1D initial value is p i ).
Equations (42) Similarly, for the external boundary with constant pressure, if the sweep radius is larger than the outer boundary radius, the sweep radius is taken as the outer boundary radius. For the closed boundary, if the sweep radius is greater than the outer boundary radius, the outer boundary radius will not change, but the pressure at the outer boundary will change. The nonlinear solution calculating process of the low-permeability reservoir is presented below.
If r 1D f r eD , taking p i2D = 0(outer boundary pressure is pi). Solving q D by Equation (46), calculation terminated. 3

If infinite boundary
Step2: Solving r 2D by Equation (42).Taking r 2D replace r 1D , solving q D by Equation (38). Repeat step2.  The analytical solution could be solved by Stehfest's numerical inversion [23]. Comparing the analytical solution with nonlinear solution, the relative error is calculated by the following formula. δ = q as − q ns q as × 100% (52)

Results and Error Analysis
Among them: q as is an analytical solution, calculated by Equations (49)-(52); q ns is a nonlinear solution; δ is relative error. Tables 1-3 show the comparison results of the nonlinear solution and the analytical solution under different external boundary conditions. The error of the nonlinear solution and the analytical solution is very small, which fully meets the calculation accuracy requirements.

Propagation Mechanism of Moving Boundary
The moving boundary propagation mechanism is a difficult point in the study of low-permeability reservoirs, and the pressure sweep radius is the core indicator of lowpermeability reservoir development technology policy formulation. At present, mainstream commercial numerical simulation software uses the Darcy flow method to simulate the development of low-permeability reservoirs. Many experts and scholars adopt a step-bystep numerical simulation method, which firstly determines the flowable grid iteratively, and then performs the simulation calculation on the flowable grid. When the mainstream nonuniform grid is used to perform step-by-step simulation calculations, the calculation error will be significant if the grid is too large.
The model in this study can accurately calculate the low-permeability reservoir sweep radius, and the pressure sweep radius with different start-up pressure gradients was calculated. A convenient formula for calculating the swept radius under different starting pressure gradients must be established, and it can provide theoretical guidance for the formulation of low-permeability reservoir development technology policies.
Different models are used for fitting. When a power function is used for fitting ( Figure 5), the correlation coefficient is the greatest (R = 0.999). Therefore, it is reasonable to use the power function model to fit. Calculate the curve of the sweep radius versus time under different starting pressure gradients G D , and, respectively, fit a and b under different starting pressure gradients. The coefficient a and exponent b linear fitting curves were drawn under different starting pressure gradients (Figures 6 and 7). The fitting results indicate that the linear model can reach a higher fitting accuracy (Table 4, R exceeds 0.998), and the linear model can be used for an approximate calculation.   By substituting the power function, the calculation formula of the sweep radius under constant flow pressure conditions of a low-permeability reservoir can be obtained: The sweep radius is a function of a dimensionless start-up pressure gradient and dimensionless time. The propagation law of the sweep radius in a low-permeability reservoir can be determined through dimensionless formula conversion, which can provide theoretical guidance for well spacing and well pattern design in a low-permeability reservoir.

Analysis of Production Decline
The production decline curve reflects the stable production situation of the reservoir. Generally, the production decline curve generated by statistical field data is compared with the theoretical curve chart to evaluate the reservoir development effect. According to the established nonlinear calculation model of the outer boundary of a low-permeability reservoir, the production decline curve chart of a low-permeability reservoir with different outer boundary conditions under constant flow pressure can be obtained.
For infinite reservoirs, the boundary cannot be detected within the calculation time, so the key factor affecting the production decline curve is only the dimensionless start-up pressure gradient. For infinite reservoirs, the larger the start-up pressure gradient is, the faster the production will decline (Figure 8). For closed outer boundary reservoirs, both dimensionless start-up pressure gradient and outer boundary radius will affect the shape of the production decline curve (Figures 9 and 10). The starting pressure gradient affects the slope of the curve. The larger the dimensionless starting pressure gradient is, the faster the output will decline; the smaller the outer boundary radius is, the faster the output will be depleted.  For constant pressure boundary reservoirs, the larger the starting pressure gradient is, the faster the production will decline. After the final production is stable, the larger the starting pressure gradient, the lower the production (Figures 11 and 12).

Recovery Factor of Low-Permeability Reservoir
The recovery factor determines the economic and technical limits of low-permeability reservoirs and is the core indicator of the low-permeability reservoir development program.
At present, empirical formula methods, field statistics methods, numerical simulation methods, and theoretical derivation methods, etc., are mainly used for calculating the recovery factor of low-permeability oil reservoirs. In this study, the production curve chart of closed reservoirs under different starting pressure gradients has been established. Therefore, the cumulative production can be obtained by numerical calculation methods, and the recovery factor of low-permeability reservoirs can be calculated using the cumulative production. According to the closed-boundary production decline curve, the cumulative production is approximated by the numerical integration method.
Among them, t i and q i are output decline curve data. Then, the cumulative production and geology reserves can be used to calculate the recovery factor and the oil recovery factor under different outer boundary radii and start-up pressure gradients and draw the low-permeability reservoir recovery factor map ( Figure 13). For a given starting pressure gradient and well spacing, the oil recovery factor can be queried from the chart.

Conclusions
(1) Based on the steady-state sequential substitution method of the seepage mechanics theory, the moving boundary nonlinear equations of the infinite outer boundary, constant pressure outer boundary and closed outer boundary in production with constant flow pressure of a homogeneous reservoir under Darcy's seepage condition are derived, respectively. The nonlinear equation is solved by the quasi-Newton method, and the obtained nonlinear solution is compared with the analytical solution.
The comparison shows that the method has high accuracy and can be used to solve the problem of unstable seepage. The research in this paper provides a third method for solving the oil and gas seepage equation, which is different from the analytical solution and the numerical solution and provides a new method for the solution of linear partial differential equations.