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

Authors: Zhenghe Yan; Chong Cao; Mingying Xie; Pin Jia; Guotao Liu; Xihui Wei; Linsong Cheng; Zhiyong Tu
Journal: Geofluids
 (2021)

---

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


---

## Body

## 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].
(1)v=−kμ⋅∂p∂r1−λ∂p/∂r∂p∂r>λ,0∂p∂r≤λ,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 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. Figure1 is a schematic diagram of the composite model caused by sand migration.Figure 1
Physical model of composite zone.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 radiusr, 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 pressurepi(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:(2)1rD∂∂rDrD∂p1D∂rD−α∂p1D∂rD2+λDrD+LD2∂2p1D∂zD2−αD∂p1D∂zD2=hDLD2eαp1D∂p1D∂tD0≤rD≤rmD.Flow model of outer zone:(3)1rD∂∂rDrD∂p2D∂rD−αD∂p2D∂rD2+λDrD+LD2∂2p2D∂zD2−αD∂p2D∂rD2=hDLD2MωeαDp2D∂p2D∂tDrmD≤rD≤r2D.Initial condition:(4)p1DtD=0=p2DtD=0=0.At interface:(5)p1DtD=0=p2DtD=0=0,(6)∂p1D∂rDrD=rmD=1M12∂p2D∂rDrD=rmD.Internal boundary condition:(7)limεD→0limr0→0∫zwD−εD2zwD+εD2rDe−αDP1D∂p1D∂rD+λDdzwD=−12,zD−zwD≤εD2.External boundary conditions:(8)p2DrD→∞=0,(9)∂p1D∂zDzD=0,1=∂p2D∂zDzD=0,1=0.


### 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].
(10)pDrD,tD=−1αDln1−αDηDrD,tD.Substituting Eq. (10) into two zone composite mathematical model, we can get the linearized model as follows:
(11)1rD∂∂rDrD∂η1D∂rD+λDrD+LD2∂2η1D∂zD2=hDLD2∂η1D∂tD0≤rD≤rmD,1rD∂∂rDrD∂η2D∂rD+λDrD+LD2∂2η2D∂zD2=hDLD2Mω∂η2D∂tDrmD≤rD≤r2D.Equation (11) is dealt with the Laplace transformation in order to acquire the pressure distribution quickly. Thus, Eq. (12) can be obtained.
(12)1rD∂∂rDrD∂η¯1D∂rD+λDsrD+LD2∂2η¯1D∂zD2=hDLD2s∂η¯1D∂tD0≤rD≤rmD,1rD∂∂rDrD∂η¯2D∂rD+λDsrD+LD2∂2η¯2D∂zD2=hDLD2Mωs∂η¯2D∂tDrmD≤rD≤r2D.


#### 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:(13)θ1p¯1DxD,yD;s=∫Γ1DG¯1D∂p¯1D∂nD−p¯1D∂G¯1D∂nDdΓ+π∫Ω1Dq¯wDG¯1DdΩ.Like the inner zone, for the outer zone(14)θ2p¯2DxD,yD;s=∫Γ1D+Γ2DG¯2D∂p¯2D∂nD−p¯2D∂G¯2D∂nDdΓ,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]:
(15)G¯1DxD,yD,XD,YD;s=−K0rDs,G¯2DxD,yD,XD,YD;s=−K0rDsMFs,where K0 is defined as the second zero-order modified Bessel function, and then rD is expressed as
(16)rD=xD−XD2+yD−YD2.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:
(17)π∫Ω1Dq¯wDG¯1DdΩ=π∑j=1NFq¯wDjS¯wDjxD,yD,xwDj,ywDj,θwj;s,where
(18)S¯wDjxD,yD,xwDj,ywDj,θwj;s=12ΔLwDj∫−11K0rDsdξ,(19)rD=xD−xDj2+yD−yDj2,(20)xDj=xwDj+ξ12ΔLwDjcosθwj,(21)yDj=ywDj+ξ12ΔLwDjsinθwj.


#### 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.Figure 2
Schematic diagram of circular boundary linearization.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).
(22)θp¯DxD,yD;s=∑L=1Nb∫ΔΓDLG¯D∂p¯D∂nDdΓ−∑L=1Nb∫ΔΓDLp¯D∂G¯D∂nDdΓ+π∑j=1Nwq¯wDjS¯wDj,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:
(23)θp¯DxD,yD;s=∑L=1Nb∫−11G¯D∂p¯D∂nDJLηLdηL−∑L=1Nb∫−11p¯D∂G¯D∂nDJLηLdηL+π∑j=1NFq¯wDjS¯wDj,where JLηL, G¯D, and p¯D can be expressed as:
(24)JLηL=dXDdηL2+dYDdηL21/2,(25)G¯D=G¯DxD,yD,XDηL,YDηL;s,(26)p¯D=XDηL,YDηL;s,where variations of p¯D and ∂p¯D/∂nD 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:
(27)FLηL=p¯DLf1ηL+p¯DL+1f2ηL,where p¯DL and p¯DL+1 represent the endpoints of the Lth element pressures. As a result, JLηL=XDL+1−XDL2+YDL+1−YDL21/2/2, for linear elements, Eq. (22) can be expressed as:
(28)θp¯DxD,yD;s=∑L=1NbJL∂p¯DL∂nD∫‐11G¯DηLf1ηLdηL+∂p¯DL+1∂nD∫‐11G¯DηLf2ηLdηL−∑L=1NbJLp¯DL∫‐11∂G¯DηL∂nDf1ηLdηL+p¯DL+1∫‐11∂G¯DηL∂nDf2ηLdηL+π∑j=1NFq¯wDjS¯wDj.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
(29)Ap¯D+B∂p¯D/∂nD+Cq¯wD=0AppendixA Provides Detailed Expression of A, B, and C. in Order to Solve the Outer Zone Equation, Eq. (29) Can Be Rewritten as
(30)A1,11p¯D1,1+B1,11∂p¯D/∂nD1,1+Cw1q¯wD=0.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
(31)E1,1Fp¯D1,1+B1,1F∂p¯D/∂nD1,1+CwFq¯wD+Hp¯wD=0.The term matricesE is provided in Appendix A, and the H 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:
(32)A2,22A2,12−MB2,12A2,21A2,11−MB2,11A1,11B1,11CF1E1,1FB1,1FCwFRHT−bTbb×p¯D2,2p¯D1,1∂p¯D/∂nD1,1q¯wDp¯wDp¯wD=−πsB2,22qeD−πsB2,21qeDr.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
(33)η¯wD=sη¯1D+Ss+CDs2sη¯1D+S.


## 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. Figure1 is a schematic diagram of the composite model caused by sand migration.Figure 1
Physical model of composite zone.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 radiusr, 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 pressurepi(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:(2)1rD∂∂rDrD∂p1D∂rD−α∂p1D∂rD2+λDrD+LD2∂2p1D∂zD2−αD∂p1D∂zD2=hDLD2eαp1D∂p1D∂tD0≤rD≤rmD.Flow model of outer zone:(3)1rD∂∂rDrD∂p2D∂rD−αD∂p2D∂rD2+λDrD+LD2∂2p2D∂zD2−αD∂p2D∂rD2=hDLD2MωeαDp2D∂p2D∂tDrmD≤rD≤r2D.Initial condition:(4)p1DtD=0=p2DtD=0=0.At interface:(5)p1DtD=0=p2DtD=0=0,(6)∂p1D∂rDrD=rmD=1M12∂p2D∂rDrD=rmD.Internal boundary condition:(7)limεD→0limr0→0∫zwD−εD2zwD+εD2rDe−αDP1D∂p1D∂rD+λDdzwD=−12,zD−zwD≤εD2.External boundary conditions:(8)p2DrD→∞=0,(9)∂p1D∂zDzD=0,1=∂p2D∂zDzD=0,1=0.


## 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].
(10)pDrD,tD=−1αDln1−αDηDrD,tD.Substituting Eq. (10) into two zone composite mathematical model, we can get the linearized model as follows:
(11)1rD∂∂rDrD∂η1D∂rD+λDrD+LD2∂2η1D∂zD2=hDLD2∂η1D∂tD0≤rD≤rmD,1rD∂∂rDrD∂η2D∂rD+λDrD+LD2∂2η2D∂zD2=hDLD2Mω∂η2D∂tDrmD≤rD≤r2D.Equation (11) is dealt with the Laplace transformation in order to acquire the pressure distribution quickly. Thus, Eq. (12) can be obtained.
(12)1rD∂∂rDrD∂η¯1D∂rD+λDsrD+LD2∂2η¯1D∂zD2=hDLD2s∂η¯1D∂tD0≤rD≤rmD,1rD∂∂rDrD∂η¯2D∂rD+λDsrD+LD2∂2η¯2D∂zD2=hDLD2Mωs∂η¯2D∂tDrmD≤rD≤r2D.


### 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:(13)θ1p¯1DxD,yD;s=∫Γ1DG¯1D∂p¯1D∂nD−p¯1D∂G¯1D∂nDdΓ+π∫Ω1Dq¯wDG¯1DdΩ.Like the inner zone, for the outer zone(14)θ2p¯2DxD,yD;s=∫Γ1D+Γ2DG¯2D∂p¯2D∂nD−p¯2D∂G¯2D∂nDdΓ,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]:
(15)G¯1DxD,yD,XD,YD;s=−K0rDs,G¯2DxD,yD,XD,YD;s=−K0rDsMFs,where K0 is defined as the second zero-order modified Bessel function, and then rD is expressed as
(16)rD=xD−XD2+yD−YD2.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:
(17)π∫Ω1Dq¯wDG¯1DdΩ=π∑j=1NFq¯wDjS¯wDjxD,yD,xwDj,ywDj,θwj;s,where
(18)S¯wDjxD,yD,xwDj,ywDj,θwj;s=12ΔLwDj∫−11K0rDsdξ,(19)rD=xD−xDj2+yD−yDj2,(20)xDj=xwDj+ξ12ΔLwDjcosθwj,(21)yDj=ywDj+ξ12ΔLwDjsinθwj.


### 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.Figure 2
Schematic diagram of circular boundary linearization.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).
(22)θp¯DxD,yD;s=∑L=1Nb∫ΔΓDLG¯D∂p¯D∂nDdΓ−∑L=1Nb∫ΔΓDLp¯D∂G¯D∂nDdΓ+π∑j=1Nwq¯wDjS¯wDj,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:
(23)θp¯DxD,yD;s=∑L=1Nb∫−11G¯D∂p¯D∂nDJLηLdηL−∑L=1Nb∫−11p¯D∂G¯D∂nDJLηLdηL+π∑j=1NFq¯wDjS¯wDj,where JLηL, G¯D, and p¯D can be expressed as:
(24)JLηL=dXDdηL2+dYDdηL21/2,(25)G¯D=G¯DxD,yD,XDηL,YDηL;s,(26)p¯D=XDηL,YDηL;s,where variations of p¯D and ∂p¯D/∂nD 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:
(27)FLηL=p¯DLf1ηL+p¯DL+1f2ηL,where p¯DL and p¯DL+1 represent the endpoints of the Lth element pressures. As a result, JLηL=XDL+1−XDL2+YDL+1−YDL21/2/2, for linear elements, Eq. (22) can be expressed as:
(28)θp¯DxD,yD;s=∑L=1NbJL∂p¯DL∂nD∫‐11G¯DηLf1ηLdηL+∂p¯DL+1∂nD∫‐11G¯DηLf2ηLdηL−∑L=1NbJLp¯DL∫‐11∂G¯DηL∂nDf1ηLdηL+p¯DL+1∫‐11∂G¯DηL∂nDf2ηLdηL+π∑j=1NFq¯wDjS¯wDj.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
(29)Ap¯D+B∂p¯D/∂nD+Cq¯wD=0AppendixA Provides Detailed Expression of A, B, and C. in Order to Solve the Outer Zone Equation, Eq. (29) Can Be Rewritten as
(30)A1,11p¯D1,1+B1,11∂p¯D/∂nD1,1+Cw1q¯wD=0.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
(31)E1,1Fp¯D1,1+B1,1F∂p¯D/∂nD1,1+CwFq¯wD+Hp¯wD=0.The term matricesE is provided in Appendix A, and the H 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:
(32)A2,22A2,12−MB2,12A2,21A2,11−MB2,11A1,11B1,11CF1E1,1FB1,1FCwFRHT−bTbb×p¯D2,2p¯D1,1∂p¯D/∂nD1,1q¯wDp¯wDp¯wD=−πsB2,22qeD−πsB2,21qeDr.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
(33)η¯wD=sη¯1D+Ss+CDs2sη¯1D+S.


## 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].
(10)pDrD,tD=−1αDln1−αDηDrD,tD.Substituting Eq. (10) into two zone composite mathematical model, we can get the linearized model as follows:
(11)1rD∂∂rDrD∂η1D∂rD+λDrD+LD2∂2η1D∂zD2=hDLD2∂η1D∂tD0≤rD≤rmD,1rD∂∂rDrD∂η2D∂rD+λDrD+LD2∂2η2D∂zD2=hDLD2Mω∂η2D∂tDrmD≤rD≤r2D.Equation (11) is dealt with the Laplace transformation in order to acquire the pressure distribution quickly. Thus, Eq. (12) can be obtained.
(12)1rD∂∂rDrD∂η¯1D∂rD+λDsrD+LD2∂2η¯1D∂zD2=hDLD2s∂η¯1D∂tD0≤rD≤rmD,1rD∂∂rDrD∂η¯2D∂rD+λDsrD+LD2∂2η¯2D∂zD2=hDLD2Mωs∂η¯2D∂tDrmD≤rD≤r2D.


## 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:(13)θ1p¯1DxD,yD;s=∫Γ1DG¯1D∂p¯1D∂nD−p¯1D∂G¯1D∂nDdΓ+π∫Ω1Dq¯wDG¯1DdΩ.Like the inner zone, for the outer zone(14)θ2p¯2DxD,yD;s=∫Γ1D+Γ2DG¯2D∂p¯2D∂nD−p¯2D∂G¯2D∂nDdΓ,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]:
(15)G¯1DxD,yD,XD,YD;s=−K0rDs,G¯2DxD,yD,XD,YD;s=−K0rDsMFs,where K0 is defined as the second zero-order modified Bessel function, and then rD is expressed as
(16)rD=xD−XD2+yD−YD2.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:
(17)π∫Ω1Dq¯wDG¯1DdΩ=π∑j=1NFq¯wDjS¯wDjxD,yD,xwDj,ywDj,θwj;s,where
(18)S¯wDjxD,yD,xwDj,ywDj,θwj;s=12ΔLwDj∫−11K0rDsdξ,(19)rD=xD−xDj2+yD−yDj2,(20)xDj=xwDj+ξ12ΔLwDjcosθwj,(21)yDj=ywDj+ξ12ΔLwDjsinθwj.


## 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.Figure 2
Schematic diagram of circular boundary linearization.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).
(22)θp¯DxD,yD;s=∑L=1Nb∫ΔΓDLG¯D∂p¯D∂nDdΓ−∑L=1Nb∫ΔΓDLp¯D∂G¯D∂nDdΓ+π∑j=1Nwq¯wDjS¯wDj,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:
(23)θp¯DxD,yD;s=∑L=1Nb∫−11G¯D∂p¯D∂nDJLηLdηL−∑L=1Nb∫−11p¯D∂G¯D∂nDJLηLdηL+π∑j=1NFq¯wDjS¯wDj,where JLηL, G¯D, and p¯D can be expressed as:
(24)JLηL=dXDdηL2+dYDdηL21/2,(25)G¯D=G¯DxD,yD,XDηL,YDηL;s,(26)p¯D=XDηL,YDηL;s,where variations of p¯D and ∂p¯D/∂nD 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:
(27)FLηL=p¯DLf1ηL+p¯DL+1f2ηL,where p¯DL and p¯DL+1 represent the endpoints of the Lth element pressures. As a result, JLηL=XDL+1−XDL2+YDL+1−YDL21/2/2, for linear elements, Eq. (22) can be expressed as:
(28)θp¯DxD,yD;s=∑L=1NbJL∂p¯DL∂nD∫‐11G¯DηLf1ηLdηL+∂p¯DL+1∂nD∫‐11G¯DηLf2ηLdηL−∑L=1NbJLp¯DL∫‐11∂G¯DηL∂nDf1ηLdηL+p¯DL+1∫‐11∂G¯DηL∂nDf2ηLdηL+π∑j=1NFq¯wDjS¯wDj.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
(29)Ap¯D+B∂p¯D/∂nD+Cq¯wD=0AppendixA Provides Detailed Expression of A, B, and C. in Order to Solve the Outer Zone Equation, Eq. (29) Can Be Rewritten as
(30)A1,11p¯D1,1+B1,11∂p¯D/∂nD1,1+Cw1q¯wD=0.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
(31)E1,1Fp¯D1,1+B1,1F∂p¯D/∂nD1,1+CwFq¯wD+Hp¯wD=0.The term matricesE is provided in Appendix A, and the H 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:
(32)A2,22A2,12−MB2,12A2,21A2,11−MB2,11A1,11B1,11CF1E1,1FB1,1FCwFRHT−bTbb×p¯D2,2p¯D1,1∂p¯D/∂nD1,1q¯wDp¯wDp¯wD=−πsB2,22qeD−πsB2,21qeDr.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
(33)η¯wD=sη¯1D+Ss+CDs2sη¯1D+S.


## 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 Figure3.Figure 3
Field test data and history matching curves by conventional model and our new bizone model.As shown in Figure3, 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.Table 1
Interpretation results of field test data.
ParametersValuesParametersValuesWellbore storage, m3/MPa594Skin factor, dimensionless-8.05Stress sensitivity coefficient, dimensionless0.01Inner radius, m67Threshold pressure gradient, MPa/m0.1Reservoir pressure, MPa11.5Inner permeability, mD124Outer permeability, mD300Vertical permeability, mD30Permeability ratio0.41


### 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 storageCD is 0.0001, skin factor S is 1, permeability ratio M12 is 5, stress sensitivity coefficient αD is 0.01, threshold pressure gradient of dimensionless λD is 0.01, and sanding production radius of dimensionless rD 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.Figure 4
Typical pressure curves of bizonal composite model.Figure 5
Schemata of typical flow regimes.


### 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

Figure6 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 6
Influence of internal radius on typical pressure curves.


#### 3.3.2. Effect of Permeability Ratio

Figure7 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.Figure 7
Influence of permeability ratio on typical pressure curves.


## 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 Figure3.Figure 3
Field test data and history matching curves by conventional model and our new bizone model.As shown in Figure3, 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.Table 1
Interpretation results of field test data.
ParametersValuesParametersValuesWellbore storage, m3/MPa594Skin factor, dimensionless-8.05Stress sensitivity coefficient, dimensionless0.01Inner radius, m67Threshold pressure gradient, MPa/m0.1Reservoir pressure, MPa11.5Inner permeability, mD124Outer permeability, mD300Vertical permeability, mD30Permeability ratio0.41


## 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 storageCD is 0.0001, skin factor S is 1, permeability ratio M12 is 5, stress sensitivity coefficient αD is 0.01, threshold pressure gradient of dimensionless λD is 0.01, and sanding production radius of dimensionless rD 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.Figure 4
Typical pressure curves of bizonal composite model.Figure 5
Schemata of typical flow regimes.


## 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

Figure6 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 6
Influence of internal radius on typical pressure curves.


### 3.3.2. Effect of Permeability Ratio

Figure7 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.Figure 7
Influence of permeability ratio on typical pressure curves.


## 3.3.1. Effect of Sand Production Radius

Figure6 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 6
Influence of internal radius on typical pressure curves.


## 3.3.2. Effect of Permeability Ratio

Figure7 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.Figure 7
Influence of permeability ratio on typical pressure curves.


## 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



---

*Source: 6658875-2021-02-23.xml*