Abstract
In recent years, sand production has been frequently observed in offshore weakly consolidated sandstone reservoirs. Permeability changes due to sand migration seriously affect the confidence in well test interpretation, production forecasts, and oilfield development plan schedules. The purpose of this paper is to propose a comprehensive model of coupled sand migration, stress sensitivity, and high viscosity oil and to study the effect of sand production induced permeability zoning on transient pressure behavior by combining discrete boundary and discrete wellbore with the boundary element method. In this two-zone composite model, the reservoir can be divided into the inner zone with the improved permeability due to sand migration and the outer zone with initial reservoir permeability. The multifactor effects of stress-sensitive, highly viscous oil, sand migration, and horizontal well are included in this model. Thus, the seepage equation presents a highly nonlinear and difficult to obtain an accurate analytical solution. In this paper, the boundary element method (BEM) is introduced to separate the boundary and wellbore, and the semianalytical solution of the hybrid model is obtained. The comparative analysis of measured pressure curve fitting from a horizontal well, located in the eastern of the South China Sea, proves that this comprehensive model can be used for pressure transient analysis of the weakly consolidated sandstone reservoir. The flow regime analysis indicates that a two-zone composite system may develop seven flow regimes: the wellbore storage stage, early-time radial stage, first transition stage, inner linear stage, inner pseudoradial flow, transition flow from the inner area to the outer area, and outer pseudoradial flow. Sensitivity analysis indicates that the smaller the sand production radius, the shorter the duration of the transition flow from the inner to the outer zone, which suggests the well is mainly affected by the outer boundary in the later period. The larger the permeability ratio, the higher the pressure curves may move up.
1. Introduction
Sand production is a common phenomenon in weakly consolidated sand reservoirs, which may cause the permeability changes in loose sandstone reservoirs, producing troublesome and costly problems to oil well production. Many researchers have focused on sand controls, and many techniques have been raised on how to control produced sands [1–3]. The sand production process and mechanism were investigated combined with the particle flow model and discrete element method, which indicated that the sand body would collapse when the pressure drop was greater than the critical pressure drop [4]. The influence of several factors on sand production mechanism was studied in unconsolidated sand reservoirs based on which cumulative sand production could be predicted by the gas injection pressure and moisture content [5]. A prediction model for sand production was established for the weak reservoirs of heavy oil flow, where the sanding area can be considered as the plastic zone around the perforated tunnel [6]. The effects of fluid flow rate and external stress on sand production were performed, which indicated sand rate increased gradually when the stress was great than initial sand production stress [7]. The sand-failure criterion based on pressure gradient was built to obtain sand production during the cold heavy-oil production with sand, the results of which indicated the wormhole generation might cause permeability changes [8]. An analytical model describing average sand production rate quantitatively was confirmed by experimental sand production data [9]. A three-dimensional (3D) sand production model coupled multiphase fluid flow, and elastoplastic was built to investigate the transient pressure behavior, which could combine the mechanical failure and fluid erosion [10]. Experiments were conducted to investigate the cause of sand production in weakly consolidated heavy oil sand cores. The model for the effect of gas exsolution phenomena and geomechanics was developed based on the experiment results [11]. A new model for dynamic sand production mechanism was proposed to analyze the effect of field stress and pressure-depletion effect, the result of which indicated sand might produce from wellbore due to either a large pressure gradient or a large porosity gradient [12].
There were a lot of literature reports on sand production mechanism and the prediction of sand production. However, few studies paid attention to the impact of sand production on transient pressure behavior. A comprehensive model to analyze the transient pressure behavior of weakly consolidated offshore sand reservoirs was proposed in this paper, in which two zones were divided based on the permeability changes due to sand production.
At present, two-zone models were mainly introduced to analyze the pressure transient response of unconventional oil and gas reservoirs. An analytical trilinear flow model of fractured horizontal wells was proposed to analyze the pressure transient and production behavior [13]. A simplified comprehensive model with four rectangle region which included the near wellbore engineered area, SRV, the original reservoir area, and the outer area in multistage fractured horizontal well was proposed by Zhao [14]. A new analytical model of natural fractured and hydraulic fractured horizontal well in double porosity tight gas reservoirs in which the typical curves are more suitable for the homogenous and naturally fractured reservoirs was proposed by Xu et al. [15]. Wang et al. built a two-zone composite pressure response model of stimulated reservoir volume (SRV) and unstimulated reservoir volume (un-SRV) in the tight oil formation fracturing vertical well. The SRV was considered as a dual-porosity formation while the un-SRV was modelled as a single porosity formation of stress-sensitive effect [16]. Zhao et al. presented a comprehensive rectangular method to describe the unique percolation mechanism for multiscaled shale fractured horizontal wells [17]. Li et al. developed a dual-porosity mathematical model for fracture horizontal well in tight gas reservoirs, taking the stress-sensitivity effect into account, which can be solved by Laplace transformation, orthogonal transformation, perturbation technique, and numerical inversion [18]. Wu et al. used permeability modulus and dynamic threshold pressure gradient to describe permeability stress-sensitive effect and non-Darcy flow in tight reservoirs [19].
Due to the weak cementation of loose sandstone, the effect of permeability dependence on stress is more pronounced. Many researchers have studied the stress-sensitive effect on the pressure transient behavior [20]. The solution of pressure behavior was obtained from the permeability stress-sensitive formations [21]. The integrated effects of reservoir properties and fluid flow characteristics on well production performance were investigated [22]. The analytical model of the unsteady-state or pseudosteady state flow for pressure transient interpretation was developed in stress-sensitive reservoirs [23].
In loose sandstone formation, the reservoir rocks are poorly cemented, and the formation crude oil has a high viscosity. The high viscosity and non-Newtonian fluid characteristics can be considered as threshold pressure gradient. The previous research results showed that Eq. (1) can be applied to describe the non-Darcy flow [24–26]. where represents the reservoir permeability, mD; is the oil viscosity, mPa·s; is the formation pressure, is the threshold pressure gradient, MPa/m; is the seepage velocity, cm/s; is the radial distance, m.
The previous composite model was used to analyze the pressure behavior in a low permeability reservoir, in which the reservoir was considered as dual porosity and divided into several rectangle zones with the main fracture as the boundary. The bizone composite model was hardly introduced to describe the impact of sand production on pressure curves. Thus, the purpose of this article is to provide a new model integrating the sand production, stress sensitivity effect, and threshold pressure gradient for this poorly consolidated sandstone reservoirs.
The boundary element method based on Green’s function could deal with complex boundaries and reservoir heterogeneity in high accuracy. This method is to replace the domain differential equation with the integral square discrete solution using the Green theorem and so on. The diffusion equation is transformed into an elliptic Poisson equation by Laplace transform of the seepage differential equation, and then the Laplace space solution is obtained using the boundary element square, and finally, the real space solution of the problem is obtained by numerical inversion algorithm. The boundary element theory has been introduced to solve the transient pressure and flow behavior for many years because of high precision. Boundary integral methods of real domain and Laplace domain were compared in unsteady flow aquifers [27]. Gas desorption and multiscale flow effects were considered into a composite model with an arbitrary fractured region, in which a semianalytical transient pressure could be solved combined with Laplace transform and boundary element theory [28]. A discrete fracture network model was proposed to analyze the transient pressure behavior of a composite zone, the solution of which could be derived from the line-source functions and boundary element theory [29]. A Green element method was applied to address the pressure transient behavior of heterogeneous reservoir with discrete fracture networks [30]. The solutions for flow transient behavior in fractured reservoirs were obtained by a robust boundary element method, which could accurately characterize the complex fractures [31]. The solution for arbitrary shapes in composite reservoirs could be obtained through Laplace transform BEM, taking the effects of reservoir heterogeneities and complex boundaries into consideration [32].
There are many two-zone models were built to analyze the pressure transient behavior of unconventional reservoirs. However, the coupled effects of stress sensitivity, heavy oil, and horizontal well are not considered in the previous study. Therefore, the purpose of this paper is to propose a comprehensive model of coupled sand migration, stress sensitivity, and high viscosity oil and to study the effect of sand production induced permeability zoning on transient pressure behavior by combining discrete boundary and discrete wellbore with the boundary element method.
This semianalytical model combines the analytical solution of the bottomhole pressure and the discretization of the well segments, which could deal with complex wellbores such as horizontal wellbores and branch wellbores. To obtain the semianalytic solution of this composite model, the boundary and horizontal wellbore are divided into several segments with the whole composite system coupled.
2. Methodology
2.1. Physical Model
A silting area is formed near the wellbore due to the sand retaining effect of screen completion, which can be simply considered as the skin factor in this model. Besides, more and more sand particles migrate in the formation with the development of oil production. As a result, the reservoir can be divided into two areas from the inside to the outside: the inner zone with the improved permeability due to sand migration and the outer zone with initial formation permeability. Figure 1 is a schematic diagram of the composite model caused by sand migration.
Combined with the characteristics of sand migration in the weakly consolidated sandstone reservoir, a two-zone composite model for horizontal well is proposed. Some assumptions of this new bizone composite model are as follows: (1)The weakly consolidated sandstone reservoir is simplified into a two-zone circular model with radius , and many factors such as stress sensitivity effect, heavy oil and sand migration are considered(2)The horizontal direction is infinite in this new bizonal model with closed top and bottom and initial formation pressure (3)The single-phase slightly-compressed fluid is considered in this model, and the capillary force and gravity are ignored(4)The pressure drop can be neglected at the interface between the inner and outer area
2.2. Mathematical Model
Combined the governing equation with internal boundary condition, external boundary conditions, and initial condition, a two-zone composite mathematical model is obtained as follows.
Flow model of inner zone:
Flow model of outer zone:
Initial condition:
At interface:
Internal boundary condition:
External boundary conditions:
2.3. Solutions for the Bizone Composite Model
2.3.1. Fluid Flow Model in Laplace Domain
The fluid flow equation may be highly nonlinear due to the stress sensitivity coefficient. Perturbation technique Eq. (10) can be introduced to linearize the equation. In general, the zero-order perturbation form can be applied to solve this nonlinear equation [3, 33].
Substituting Eq. (10) into two zone composite mathematical model, we can get the linearized model as follows:
Equation (11) is dealt with the Laplace transformation in order to acquire the pressure distribution quickly. Thus, Eq. (12) can be obtained.
2.3.2. Boundary Element Method
Combined with Green’s function, an equation of simplified form can be obtained using boundary element theory. Thus, pressure of the inner zone can be expressed as:
Like the inner zone, for the outer zone where can be defined as the interior angle between two adjacent elements on the boundary. And then pressure in the inner area can be calculated by Eq. (13), considering the inner area boundary conditions and horizontal well production condition.
However, Eq. (14) is supposed to remove the source terms considering that there is no wellbore in the outer area. Flow equation in the inner area and outer area can refer to the following results using the Green’s function [34]: where is defined as the second zero-order modified Bessel function, and then is expressed as
Besides, there will exist a pressure drop in the inner area when fluid flows into the horizontal wellbore. It is assumed that fluid flow is uniform in this paper so that the outer zone flow equation can be simplified as Eq. (17) using Green’s function: where
2.3.3. Flow Equation Discretization
Pressure and flow rate on the zone boundaries are supposed to be determined. In order to solve Eqs. (13) and (14), the inner boundary should be divided into several boundary elements, which are numbered in a clockwise direction, as is shown in Figure 2.
Considering that flow equations are similar, subscripts of equations are ignored. In order to simplify the derivation process, the discretized form can be expressed as Eq. (22). where is referred to as the th element on the boundary condition.
A new local coordinate system is supposed to be built to deal with the line integrals in Eq. (22). As a result, boundary can be linearized when is applied to each element. Considering the limits of integration between -1 and +1 after coordinate conversion, Eq. (22) can be expressed as follows: where , , and can be expressed as: where variations of and can be derived by interpolation functions. Linear interpolation functions are introduced to solve this problem according to the conclusion that nonlinear boundaries can be transformed by linear, quadratic, or higher-order elements. Thus, the th element can be expressed as Eq. (27) based on the Linear interpolation functions: where and represent the endpoints of the th element pressures. As a result, , for linear elements, Eq. (22) can be expressed as:
Then, the pressure in the inner zone and interface can be obtained by Eq. (28).
In the same way, we can get the pressure in the outer zone by Eq. (28), excluding that line-source terms are supposed to be ignored. Green’s function source can be applied to all points on the inner boundary. As a result, the matrix can be derived as
Appendix A Provides Detailed Expression of A, B, and C. in Order to Solve the Outer Zone Equation, Eq. (29) Can Be Rewritten as
The subscript () and in the coefficient matrices represent the boundary and wellbore of the inner zone. Besides, a fictitious source can be applied to all points of wellbore, and then the pressure of the wellbore segments can be acquired as
The term matrices is provided in Appendix A, and the is an identical matrix. The dimensionless definition in this article is given in Appendix B.
2.3.4. Coupled Solutions
Coupling the boundary element method and the two-zone composite model, the transient flow behavior can be obtained using the continuity condition of pressure and fluid flow on the inner boundary and the external boundary conditions. Therefore, equations of a two-zone model can be acquired by Eq. (31) combined with inner boundary conditions Eq. (7) and the outer boundary conditions Eq. (8). It is worth noting that the Neumann condition is supposed to be applied to the exterior boundary, and then the comprehensive equations may be expressed as:
Therefore, the Gauss elimination method can be used to solve Eq. (32). Based on the Duhamel principle [23, 25], the bottomhole pressure can be expressed as Eq. (33) considering well storage and skin effect
3. Results and Discussion
3.1. Bizone Model Validation
Based on the oil well physical parameters, the pressure and pressure derivative curves were fitted by changing key parameters such as permeability ratio and inner radius. A horizontal well, located in an offshore loose sandstone formation in the eastern of the South China Sea, was converted into a water injection well in March 2018. Well testing was carried out between January 23 and February 10, 2020. Combined with the actual pressure test data, the proposed two-zone composite model (considering stress sensitivity, high viscosity oil, and sand migration) and the conventional composite model caused by permeability change (without considering stress sensitivity and threshold pressure gradient) were used to fit the actual pressure curve. The results are shown in Figure 3.
As shown in Figure 3, in the early stage of pressure testing, both the model proposed in this paper and the conventional composite model caused by permeability changes can fit the actual pressure data. However, in the late stage of pressure testing, the composite model proposed in this paper can better fit the actual pressure data, while the conventional pressure test curve is lower, which further proves that it is necessary to propose a combined model of coupled stress sensitivity, high viscosity oil, and sand migration to study the influence of sand production process on transient pressure behavior of unconsolidated sandstone heavy oil reservoir. The key parameter results are interpreted in Table 1.
3.2. Transient-Pressure Behavior
Based on the well parameters in the offshore loose sandstone reservoir, the parameters of our new bimodel are given in this section. Wellbore storage is 0.0001, skin factor is 1, permeability ratio is 5, stress sensitivity coefficient is 0.01, threshold pressure gradient of dimensionless is 0.01, and sanding production radius of dimensionless is 10. Thus, we can derive typical pressure curves of bizonal composite model, as shown in Figure 4. The flow regime analysis indicates that there may be 7 stages: (I) the wellbore storage stage, (II) early-time radial flow, (III) first transition flow stage, (IV) inner linear flow, (V) inner pseudoradial flow, (VI) transition flow from internal area to external area, and (VII) outer pseudoradial flow. The typical flow stages of early-time radius flow, inner linear flow, and inner zone pseudosteady flow are depicted in Figure 5.
3.3. Sensitivity Analysis of Sand Production
In this two-zone comprehensive model, sand migration is mainly characterized by the change of permeability caused by sand production and the size of the sand production area. Therefore, the analysis of sand migration mainly reflects the influence of permeability ratio and sand production radius on transient pressure behavior. The detailed analysis results are shown as follows.
3.3.1. Effect of Sand Production Radius
Figure 6 indicates that sand production radius mainly affects transition flow from the inner area to the outer area. The smaller the sand production radius, the shorter duration of the transition flow stages from the inner to outer zone. When the sand production radius is smaller, the pseudoradial flow in the inner zone may disappear. It suggests the production well is mainly affected by the outer boundary in the later period, which indicates that the well with a smaller inner zone radius is urgently needed to energy supplement.
3.3.2. Effect of Permeability Ratio
Figure 7 indicates that the permeability ratio mainly affects transition flow stages from inner to outer area and outer pseudoradial flow. The larger the permeability ratio, the higher the pressure curve may rise. This is mainly because the outer zone permeability is relatively lower compared with the inner zone permeability due to the sand migration, resulting in the weak fluid flow in the outer zone.
4. Conclusions
(1)A comprehensive model of coupled sand migration, stress sensitivity, and high viscosity oil is proposed to study the effect of sand production induced permeability zoning on transient pressure behavior(2)The semianalytical solution of this composite model is obtained by the discrete boundary and discrete wellbore with the boundary element method. The flow regimes of bizonal composite model may develop seven stages: the wellbore storage stage, early-time radial stage, first transition stage, inner linear flow, inner pseudoradial flow, transition flow from inner area to the outer area, and outer pseudoradial flow(3)The smaller the inner radius, the shorter duration of transition flow stages between the inner and outer zone. It suggests the production well is mainly affected by the outer boundary in the later period, which provides some guidelines that the well with a smaller inner zone radius is urgently needed to energy supplement. The larger the permeability ratio, the higher the pressure curve may rise. This is mainly because the outer zone permeability is relatively lower compared with the inner zone permeability due to the sand migration, resulting in the weak fluid flow in the outer zone
Appendix
A. Expression of the Coefficients A, B, and C of the Matrix Equation
In the two-zone flow equations (Eq. (25)), the terms in the matrix are expressed as: where the Kronecker delta can be defined as
And
B. Dimensionless Definitions for the New Two Zone Composite Model
For the dimensionless pressure
The dimensionless radial radius is defined as
For the dimensionless wellbore radius
The dimensionless vertical distance is defined as
For the dimensionless horizontal well position
The dimensionless stress sensitivity factor is defined as
The dimensionless threshold pressure gradient is defined as
For the storage coefficient ratio
The permeability ratio is defined as
The dimensionless horizontal length is defined as
For the dimensionless reservoir thickness
Nomenclature
| : | Volume factor, dimensionless |
| : | Dimensionless wellbore storage |
| : | Production rate, m3/d |
| : | Dimensionless radius, dimensionless |
| : | Total compressibility, MPa-1 |
| : | Dimensionless thickness, dimensionless |
| : | Reservoir thickness, m |
| : | First kind modified Bessel function of zero-order |
| : | Second kind modified Bessel function of zero-order |
| : | Production time, h |
| : | Reservoir permeability, mD |
| : | Dimensionless horizontal length, dimensionless |
| : | Inner permeability, mD |
| : | Vertical permeability, mD |
| : | Dimensionless pressure, dimensionless |
| : | -coordinate, m |
| : | Initial reservoir pressure, MPa |
| : | Threshold pressure gradient, MPa/m |
| : | Horizontal well position dimensionless |
| : | Dimensionless threshold pressure gradient |
| : | Wellbore storage, m3/MPa |
| : | Radius, m |
| : | Wellbore radius, m |
| : | Dimensionless wellbore radius, dimensionless |
| : | Natural logarithm, 2.71828 |
| : | Laplace-transform variable, dimensionless |
| : | Skin factor, dimensionless |
| : | First kind modified Bessel function of first-order |
| : | Second kind modified Bessel function of first-order |
| : | Time, dimensionless |
| : | Horizontal length, m |
| : | Dimensionless vertical distance, dimensionless |
| : | Outer permeability, mD |
| : | Horizontal permeability, mD |
| : | Permeability ratio dimensionless |
| : | Current pressure, MPa |
| : | Oil viscosity, mPa·s |
| : | Seepage velocity, cm/s |
| : | Stress sensitivity factor dimensionless. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Additional Points
Highlights. (1) A novel semianalytical method to analyze the effect of sand production on pressure transient behavior for the weakly consolidated sandstone reservoirs. (2) The comprehensive effect of sand migration, stress sensitivity, and Non-Newtonian fluid flow are all incorporated in the model. (3) A workflow to solve the comprehensive model considering the complex wellbore and complex boundary is presented
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to appreciate the financial support by the Research Program of Productivity evaluation of loose sandstone heavy oil reservoir from CNOOC China Limited, Shenzhen Branch, Shenzhen, China. Besides, we would further like to acknowledge that the National Natural Science Foundation of China (No. U1762210 and No. 51774297) provide great financial support in this research.