Pressure Behavior Analysis of Permeability Changes Due to Sand Production in Offshore Loose Sandstone Reservoirs Using Boundary-Element Method

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 twozone 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.


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][2][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 stresssensitive 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][25][26].
where k represents the reservoir permeability, mD; μ is the oil viscosity, mPa·s; p is the formation pressure, λ is the threshold pressure gradient, MPa/m; v is the seepage velocity, cm/s; r 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 2 Geofluids 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.

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 r, and many factors such as stress sensitivity effect, heavy oil and sand migration are considered Flow model of outer zone: Initial condition: At interface: Internal boundary condition: Un-sanding zone Wellbore zone

Geofluids
External boundary conditions:  (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.
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 K 0 is defined as the second zero-order modified Bessel function, and then r D 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 x Dj = x wDj + ξ 1 2 ΔL wDj cos θ wj , ð20Þ

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.

Geofluids
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 ΔΓ DL is referred to as the Lth 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 ΔΓ D 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 jJ L ðη L Þj, G D , and p D can be expressed as: where variations of p D and ∂ p D /∂n D 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 Lth element can be expressed as Eq. (27) based on the Linear interpolation functions: where p DL and p DL+1 represent the endpoints of the Lth element pressures. As a result, jJ L ðη L Þj = ½ðX DL+1 − X DL Þ 2 + ðY DL+1 − Y DL Þ 2 1/2 /2, for linear elements, Eq. (22) can be expressed as: Then, the pressure in the inner zone and interface Γ 1 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 The subscript (1, 1) and w in the coefficient matrices represent the boundary Γ 1 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 E is provided in Appendix A, and the H is an identical matrix. The dimensionless definition in this article is given in Appendix B.

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  6 Geofluids 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.

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 C D is 0.0001, skin factor S is 1, permeability ratio M 12 is 5, stress sensitivity coefficient α D is 0.01, threshold pressure gradient of dimensionless λ D is 0.01, and sanding production radius of dimensionless r D is 10. Thus, we can derive typical pressure curves of bizonal composite model, as shown in 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. 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.

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   (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 S wDK,j = S wDj x DK , y DK , x wDj , y wDj ; s , ðA:6Þ The dimensionless stress sensitivity factor is defined as α D = 1:842 × 10 −3 q sc μB k 1 h α: ðB:6Þ The dimensionless threshold pressure gradient is defined as λ D = k 1 hrr w λ 1:842 × 10 −3 q sc μB : ðB:7Þ For the storage coefficient ratio The permeability ratio is defined as The dimensionless horizontal length is defined as