Analysis of Pressure and Production Transient Characteristics of Composite Reservoir with Moving Boundary

: The mathematical model of composite reservoir has been widely used in well test analysis. In the process of oil recovery, due to the injection or replacement of the displacement agent, the model boundary can be moved. At present, the mathematical model of a composite reservoir with a moving boundary is less frequently studied and cannot meet industrial demand. In this paper, a mathematical model of a composite reservoir with a moving boundary is developed, with consideration of wellbore storage and skin effects. The characteristics of pressure transient in moving boundary composite reservoir are studied, and the influences of parameters, such as initial boundary radius, moving boundary velocity, skin factor, wellbore storage coefficient, diffusion coefficient ratio, and mobility ratio on pressure and production, are analyzed. The moving boundary effects are noticeable mainly in the middle and late production stages. The proposed model provides a novel theoretical basis for well test analysis in these types of reservoirs.


Introduction
A composite reservoir model is widely used to model reservoirs with distinct regional properties variations. Among others, these variations may be caused by changes in formation physical properties, foreign fluids intrusion, and changes in formation fluid characteristics. These two different areas are separated by a discontinuous interface. Actual production cases often need to be modeled as a composite reservoir, such as in the cases of reservoir flooding (water, gas, chemical, steam) and formation stimulation. In addition, if the formation is stimulated or damaged, it can also be considered a composite reservoir. Therefore, the composite reservoir model has been widely investigated in literature [1][2][3][4][5][6][7][8][9][10][11]. However, there is very little research on composite reservoirs with a moving boundary. There are many cases in which a moving boundary is formed in the actual production process. For example, piston displacement or steam flooding may cause boundary movement. It can be seen that studying the moving boundary has far-reaching effects on the theory of composite reservoirs and industrial demand.
Current pressure transient analysis is an important method for analyzing reservoir characteristics. Most investigators only consider closed, constant pressure, and infinite acting boundaries in their composite reservoir models. These assumptions may lead to biased interpretation if moving boundaries do actually exist. Therefore, it is particularly important to analyze the transient pressure response of composite reservoir with moving boundary. Many papers and methods have been published on the flow of fluids in porous media [12][13][14][15][16][17][18], but the most common method used in petroleum engineering is still pressure-transient analysis.
Research activities on composite reservoir models were increasing in 1960s. Many models have been built and many methods have been applied to analyze the pressure characteristics of composite reservoir. Analytical solutions for media with the same center but different permeability have already appeared in the materials heat shrinkage literature published by Carslaw and Jaeger [19]. Subsequently, some research projects in the literature used composite reservoir models to study the pressure characteristics of reservoirs. In 1960, Hurst proposed a point-sink solution for infinitely large composite reservoirs and described many of the effects of oilfield disturbances involving fluid flow, volumetric behavior, and stratigraphic differences [20]. In 1961, Lucks and Guerrero studied the pressure characteristics of composite reservoir and calculated their pressure distribution. They compared the results with a uniform reservoir and found that the permeability and the size of the inner zone of the composite reservoir under certain conditions can be determined by the pressure drop curves [1]. In 1966, Carter proposed a method to describe the pressure behavior of a finite circular composite reservoir, and gave an analytical solution to the problem of micro-compressible fluid flow [2]. Other studies on pressure drop include Closmann and Ratliff [21], Bixel and Van Poollen [22], and Merrill et al. [23].
Among the research of the above scholars, the well test curves of composite reservoir are particularly important [24]. Scholars have done a lot of research on the effects of various factors on the well test curves. In 1985, Satman studied the interference problem in a composite reservoir and proposed a new solution to analyze the pressure response of composite reservoir observation wells, taking into account the effects of wellbore storage and skin [25]. Based on this, Olarewaju and Lee established an analytical model for the two-zone radial composite reservoir system and proposed a well test curves analysis method. Important parameters such as permeability and wellbore storage coefficient can be determined and production can be predicted [26]. Boussalem et al. studied the effect of mobility on pressure and pressure derivatives in closed composite reservoirs [27]. The influence of each parameter on the well test curve is significant, and the influence of different boundary conditions on the well test curve cannot be ignored. In 1996, Issaka summarized the analytical solutions for the pressure characteristics of composite reservoir under different flow geometries [28]. Then, Issaka and Ambashha proposed a new design equation related to the specific flow patterns observed in the test of spherical and linear composite reservoir, and studied the pressure derivative curves [29]. Based on the composite reservoir in the second zone, Ambastha et al. studied the pressure characteristics of the composite reservoir in the three zones and further improved the well testing theory of the composite reservoir [30].
Jordan and Mattar compared the pressure characteristics of composite reservoir with two-layer reservoirs [31]. They concluded that their stress characteristics are similar during transient flow, but if the model is not suitable, there will be large deviations at later production time. At present, a well testing theory has been developed, and the mathematical models of composite reservoir are also very mature, but the analysis of pressure transients and transient production of composite reservoir with moving boundary remains understudied [39]. Khatib and Noaman proposed a mathematical model to describe transient pressure behavior in composite reservoiraquifer system [40]. They used the finite difference method to solve equations in Laplace space and estimated the position of the moving front using an iterative procedure. However, the accuracy of the solutions obtained by the finite difference method is lower than the precision of the analytical solution, and they do not consider the influence of the initial interface radius and the diffusion coefficient ratio on the pressure and production. These aspects also require more in-depth research. The purpose of this paper is to study the pressure and production characteristics of composite reservoir with moving boundary. First, we used a mathematical model to accurately characterize composite reservoir, including governing equations, initial conditions, and boundary conditions. Second, in order to study the influence of various reservoir parameters on pressure and production, we carried out parameter sensitivity analysis and summarized its influence law. Based on the Darcy's Law, we established a mathematical model. We transformed the model into Laplace space [41], and found that the transformed model is a Bessel function [42]. Therefore, we obtained the analytical solution of the transformed model. Compared with previous studies, the accuracy of our model's solution is higher [40]. The Stehfest algorithm is then used to invert the value of the solution into realtime space [43]. This paper consists of five parts, arranged as follows: Section 1 is the introduction part, which briefly introduces the research content, analyzes the literature related to this research, and summarizes the previous research results in this field. Section 2 is the description part of the model and introduces the basis for building mathematical model, which is based on the Darcy's Law. Section 3 is the solution part of the model. Section 4 discusses the results of the second and third sections, mainly analyzing the pressure characteristics and the influence of various factors on pressure and output. Section 5 summarizes this paper based on research content.

Physical Model and Assumptions
Moving boundary composite reservoir means that the radius of the composite reservoir interface may move under certain conditions. Figure 1 is a schematic diagram of the moving boundary composite reservoir, and Figure 2 is a schematic diagram of boundary movement. It can be seen from Figure 1 that the composite reservoir is divided into two areas: the inner zone (indicated by 1) and the outer zone (indicated by 2). It is assumed here in that the interface between the inner and outer zones is a moving boundary that moves outward at a certain speed, as shown in Figure 2. Based on the fixed boundary, the well test interpretation model of the moving boundary composite reservoir was established, and the numerical solution of the model was obtained. At the same time, the parameter sensitivity analysis was carried out to obtain the influence of each important parameter on the well test curves.
The following assumptions were made before establishing a well test interpretation model for a moving boundary composite reservoir: (1) The reservoir is homogeneous, horizontal, uniform in thickness, and isotropic; (2) The original formation pressure (indicated by Pi) is evenly distributed, and the production rate of well is fixed after the well is opened; (3) The formation rocks and fluids are all slightly compressible; (4) The formation fluid conforms to the Darcy seepage law during seepage; (5) Wellbore storage and skin effects are accounted; (6) Gravity and capillary forces are ignored; (7) At the interface between the inner and outer zones, there is no flow loss and the formation pressure is not abrupt, and the interface moves with the propagation of mass wave; (8) The outer boundary condition is an infinite outer boundary; (9) The inner boundary expands to the outer region over time.

Dimensionless Variables Definition
In order to make the equation easy to solve, the following dimensionless variables are defined: where is the dimensionless pressure; i = 1,2 (when i = 1, it denotes moving boundary composite reservoir inner zone; when i = 2, it denotes moving boundary composite reservoir outer zone); is dimensionless time; is the dimensionless radius of the formation centered on the well. 1 is the dimensionless inner radius; is the dimensionless wellbore storage factor; 1 is the inner zone permeability (unit: μm 2 ); h is the oil layer thickness (unit: m); is the original formation pressure (unit: MPa); ( , ) is the pressure at a certain point in the formation(unit: MPa); 1 is ground flow (unit: m 3 /d ); 1 is fluid viscosity(unit: mPa • s ); B is the volume factor (unit: m 3 /m 3 ); ∅ 1 is porosity; 1 is the total system compressibility (unit: MPa −1 ); is the wellbore radius (unit: m); r is the distance from the center of the well (unit: m); 1 is the inner radius (unit: m); S is the skin factor (dimensionless).

Governing Equations
Based on the above nine assumptions combined with the general three-dimensional seepage differential equation, the dimensionless well test model of the composite boundary reservoir can be obtained: Comprehensive differential equations are shown in Equations (5)-(6): Inner area: Outer area: where = 2 /∅ 2 2 2 1 /∅ 1 1 1 denotes the diffusion coefficient ratio, dimensionless.

Initial and Boundary Conditions
Initial conditions, internal boundary conditions and interface conditions of the dimensionless well test model of the composite boundary reservoir are shown as following: Initial conditions: 1 ( , 0) = 2 ( , 0) = 0 (7) Internal boundary conditions: Interface conditions: where denotes outer zone mobility ratio and denotes mobility.
If it is an infinite outer boundary, Equation (11) can be obtained as follows:

Model Solution
Performing a Laplace transform on the Equations (5)- (11) to obtain the following formula: Comprehensive differential equations are shown in Equations (12)- (13).
Inner area: where ̅ and are image functions in a pull space and z is the Laplace parameter. Outer area: Initial conditions: Internal boundary conditions: Interface conditions: Outer boundary conditions: ̅ 2 (∞, ) = 0 (18) The governing equations (Equations (12)-(13)) are the Bessel-type ordinary differential equations, and the general solution of Equations (12)-(13) are solved as follows: where A,B,C,D are the coefficient; 0 , 1 are the first type of virtual zero-order and first-order Bessel functions; 0 , 1 are the second type of virtual zero-order and first-order Bessel functions. In Equation (20), when → ∞, there must be C=0, thus: Bringing Equation (19) and Equation (21) into Equations (15)-(17) produces: Thus, the bottom hole pressure solution is expressed in Equation (28) as follows: where ̅ is the dimensionless downhole pressure solution without considering the influence of wellbore storage and skin effects under pull space.
Then, introducing the wellbore storage and skin effects and applying the Duhamel's principle, the downhole pressure solution considering the influence of wellbore storage and skin effects is shown in Equation (29): where ̅ is the dimensionless downhole pressure solution considering the influence of wellbore storage and skin effects under pull space.
Assuming t=0, 1 = , when t=0, 1 = − it is assumed that the interface moves with the propagation of the mass wave and satisfies = √ ℎ , where c is the coefficient, c=1 when the piston is displaced, and c is the water rising rate when the piston is not displaced. Before the mass wave reaches the boundary, the interface position does not change. After the mass wave reaches the boundary, the interface begins to move. Suppose that at a certain moment 1 , the mass wave reaches the boundary (where 1 is affected by the parameters of the formation, and the value can be analyzed), the following formula can be obtained: When: where a is initial interface radius (unit: m); Q is production (unit: m 3 /d); t is production time (unit: d); 1 is the propagation time of the mass wave at the interface (unit: d).
When time is turned into a dimensionless form Equation (32) is obtained: Equations ( and define f as the moving speed of the moving boundary, then, the following formula can be obtained: According to the actual production situation of the oil field, the value of f in Equation (35) can be analyzed.
According to the definition of dimensionless variables, the relationship between dimensionless production at the bottom hole pressure and the dimensionless bottom hole flow at constant production can be derived: where ̅̅̅ is the dimensionless production at the bottom hole pressure under the pull space.
The dimensionless production at the bottom of the well can be obtained from Equation (36). The Stehfest numerical inversion method is used to obtain a solution for the true space production, and the parameter sensitivity analysis is performed by MATLAB programming (MATLAB R2017b, MathWorks, America). It can be seen from Figure 3 that the dimensionless pressure curves (solid line) and the derivative curves (dashed line) of the moving boundary composite reservoir can be divided into five flow states.

Flow Regimes Analysis
The first stage is the pure wellbore storage flow period, which appears as a straight line with a slope of 1.
The second stage is the transition period of the pure wellbore storage flow to the inner zone radial flow. The skin and wellbore storage coefficients affect the duration of the phase and the height of the peak.
The third stage is the inner zone radial flow period. This stage is affected by the interface radius between the inner and outer zones.
The fourth stage is the transition period of the radial flow of the inner zone to the radial flow of the moving boundary.
The last stage is the radial flow period of the moving boundary. Since the boundary is always moving, the derivative curves are not a horizontal line but an upturned curve. This stage can be approximated as a radial flow, and the movement of the boundary affects it. As shown in Figure 4, the larger the value of , the longer the duration of the pure wellbore storage phase, which appears as a shift in the pressure and pressure derivative curves to the right. This is because when the value of is large, the inner zone is greatly affected by the wellbore storage effect, and the radial flow phase is delayed. When the wellbore storage coefficient is relatively large ( > 1 ), the inner zone radial flow period and transition flow period will disappear. Conversely, the smaller the value of , the less affected the inner zone is by the wellbore storage. The time to reach the radial flow period of the inner zone is shortened, and the radial flow period of the inner zone is prolonged. When the mass wave propagates to the moving boundary, the influence of the wellbore storage effect will disappear, thus, the pressure derivative curves of the infinite radial flow period in the later outer zone will tend to be consistent. As shown in Figure 5, the larger the value of S, the larger the peak of the pressure and pressure derivative curves. Conversely, the smaller the value of S, the smaller the peak of the pressure and pressure derivative curves. This is due to the stimulation and formation damage, which has caused the permeability of the formation in the near-well zone to change, thus, creating additional resistance. The larger the skin factor, the greater the additional resistance generated; thus, the pressure drop in the near-well zone is greater, and the larger the opening is on the curves. Similarly, if the skin factor is too large and the inner zone radius is small, the inner zone radial flow period may also disappear. Since the skin factor only affects the near-well zone, the flow state of the fluid during the transition period and the radial flow period of the moving boundary is not affected, and the pressure derivative curves coincide at these stages. As shown in Figure 6, the pressure curves will be upturned in the radial flow period of the moving boundary. If the mobility of the inner zone is larger than the outer zone, then the fourth stage (transition period) will tilt down. On the other hand, if the inner zone's mobility is smaller than the outer zone, the fourth stage will be lifted upward. The closer the mobility ratio value is to one, the smoother the transition phase of the curve. This is because the closer the mobility ratio is to one, the smaller the difference between the inner and outer zones, and the shorter the duration of the transition phase. The larger the value of M, the larger the upturn. This is because the greater the mobility ratio, the longer the duration of the transition phase, the more obvious the effect of the moving boundary, and the larger the upturn on the curves. As shown in Figure 7, σ mainly affects the transition phase of the radial flow from the inner zone to the radial flow of the moving boundary. The larger the value of σ, the more gentle the pressure derivative curves. The smaller the value of σ, the more obvious the "concave" produced during the transition phase. The main reason is that when the fluid ratio in the inner and outer zone and the viscosity of the reservoir fluid are constant, the larger the value of σ, the faster the mass wave propagates and spreads at the interface, and the curves shows that the pressure derivative curves is gentler. When the mass wave propagates to the moving boundary, the above influence will gradually disappear, thus, the pressure derivative curves of the infinite radial flow stage in the outer zone will be more consistent again in the later stage of development. As shown in Figure 8, as the radius of the interface becomes larger, the transition from the radial flow of the inner zone to the radial flow of the moving boundary on the pressure derivative curves moves to the right. This shows that the duration of the radial flow phase of the inner zone becomes longer, and the time during which the mass wave propagates to the outer zone is also delayed. If the interface radius between the inner and outer regions is relatively large, the fluid in the inner region will gradually enter the radial flow period, where the derivative value is 0.5 in the dimensionless coordinates. At this point, the straight line segment will appear on the single logarithmic plot. However, if the radius at the interface between the inner and outer regions is relatively small, the radial flow period of the fluid may disappear and it will go directly to the fourth stage. As shown in Figure 9, as the value of f increases, the duration of the radial flow period of the moving boundary becomes longer. Therefore, the influence of the production boundary in the middle and late time production is greater. However, when the moving speed of the moving boundary is low, the pressure derivative curve becomes smoother in advance. The influence caused by the moving boundary disappears quickly. It can be seen from Figure 3 that the pressure derivative curves of the moving boundary and the fixed boundary differ greatly in the middle and late time of production. Affected by the boundary, the pressure derivative and the pressure curves rise. It can be seen that the moving boundary causes pressure loss, which is not conducive to production. As can be seen from Figure 10, the effect of the skin factor on the production is relatively large at the initial stage. The production at the initial moment decreases as the skin factor increases. However, after a long period of production, the effect of the skin on production will become smaller and smaller until it disappears. A large S value represents a reduction of the permeability in the nearwell bore zone, or formation damage. For example, the blockage caused by the migration of rock particles in the reservoir or the incompatibility of the drilling fluid and the formation fluid makes the seepage resistance in the near-well zone become larger and the wells' production becomes lower. The implementation of stimulation measures will reduce the seepage resistance of the near-well belt, that is, S becomes smaller or even becomes a negative number with large absolute values, and the production will also increase significantly. Figure 11. Effect of mobility ratio on transient production for moving boundary composite reservoir.

Effect of Mobility Ratio
As can be seen from Figure 11, when M=1, the rate of decline in production is constant, because the formation is not damaged or stimulated. When M>1, the mobility in the outer region becomes higher, the pressure drop becomes smaller, and the rate of decline in production becomes slower.
When M<1, the mobility in the outer region becomes lower, the pressure drop becomes larger, and the rate of decline in production becomes faster. In addition, it can be concluded from Figure 11 that in the final stage, the effect of the internal regional properties on production is gradually reduced until it disappears. As shown in Figure 12, the value of σ has little effect on the early production decline curves. After a certain period of production, as σ becomes larger, the rate of decline in production becomes faster, because the faster the diffusion rate, the greater the pressure drop. Moreover, in the fluid to infinity, the influence of the internal area gradually disappears. As can be seen from Figure 13, the value of a has little effect on production in the early stage. However, as the production time becomes longer, the influence of a on production gradually increases. This is related to the ratio of mobility in the inner and outer zones. If the mobility of the inner region is smaller than the mobility of the outer region, as the value of a becomes larger, the rate of decrease in production is faster. Because the mobility of the inner zone is small, the pressure drop is large. Therefore, as the value of a becomes larger, the rate of production decrease becomes faster. In the final stage, the influence of physical parameters such as skin factor and mobility ratio on production gradually disappeared.

Conclusions
In this paper, the seepage mechanism of the moving boundary composite reservoir and the solution method of the well test model are studied. A well test interpretation model for moving boundary composite reservoirs was established and the model was solved. Then, the pressure, pressure derivative, and production characteristic curve were plotted. Finally, the parameter sensitivity analysis is carried out, and the following conclusions are drawn: (1) The composite reservoir consists mainly of two zones: the inner zone and the outer zone. The pressure and pressure derivative curves include a total of five flow stages, including pure wellbore storage flow period, transition period of the pure wellbore storage flow to the inner zone radial flow, inner zone radial flow period, transition period of the radial flow of the inner zone to the radial flow of the moving boundary, and radial flow period of the moving boundary. (2) The wellbore storage coefficient increases, and the duration of the pure wellbore storage phase also increases. The larger the skin factor, the higher the peak of the pure wellbore storage flow to the inner zone radial flow transition phase. The greater the mobility ratio in the inner and outer zones, the higher the radial flow phase in the outer zone and the higher the upturn in the radial flow section of the moving boundary. The smaller the diffusion coefficient, the deeper the "concave" of the inner zone radial flow to the outer zone radial flow transition phase. The larger the initial interface radius, the longer the inner radial flow segment lasts. As the moving speed of the boundary increases, the duration of the radial flow period of the moving boundary becomes longer. Affected by the boundary, the pressure derivative and the pressure curves rise. It can be seen that the moving boundary causes pressure loss, which is not conducive to production. (3) For composite reservoir with moving boundary, as the skin factor increases, the initial production gets lower and lower. As the ratio of internal and external flow increases, the rate of decline is faster. The diffusion coefficient has little effect on the early production decline curves, but after a certain period of production, the production decreases rapidly as the diffusion coefficient increases. The interface radius has little effect on early production, but its influence increases as production time increases. If the mobility of the inner zone is greater than the mobility of the outer zone, the production decreases as the radius becomes larger, whereas the production increases as the radius becomes larger.
Author Contributions: Z.S., X.Y. and Y.J. and conceived and designed the numerical investigation; Z.S., X.Y. and S.S. created the numerical model; X.Y. and M.W. summarized the thermal recovery evaluation system; X.Y. and S.S. examined the accuracy of the proposed model; M.W. and X.Y. analyzed the result; Z.S. and X.Y. contributed analysis tools; X.Y. and Y.J. wrote the manuscript; Z.S., X.Y. and Y.J. make a great contribution to the revision of manuscript. All authors have read and agreed to the published version of the manuscript. involved in this project, and also wish to thank the journal editors and the reviewers whose constructive comments improved the quality of this paper greatly.