Abstract
The channels may be formed in the unconsolidated sands reservoir due to formation failure during high-pressure water injection or frac-packing. Based on the continuum mechanics, a mathematical model has been established to simulate the formation process of big channels in unconsolidated sands reservoir during fluid injection. The model considers the effect of reservoir heterogeneity, solid particles erosion, and deposition. The dynamic formation process of channels around the borehole and its influencing factors are analyzed by this model. The results indicate that the seepage erosion plays a very significant role in the formation of the channels during fluid injection for the unconsolidated sands with extremely low strength. The formation of the channels is closely related to the duration of fluid injection, injection pressure, reservoir heterogeneity, formation plugging, and critical fluid velocity. The long channels are more likely to form as injection time increases. Higher injection pressure will lead to higher flow rate, thus eroding the solid particles and forming big channels. The increase of the rock strength will enhance the value of critical fluid velocity, which makes it difficult for the occurrence of erosional channelization. The near-wellbore damage of the formation will decrease the flow rate, and the preferential flow channels are less likely to be induced under the same injection pressure when compared with undamaged formation. In addition, we also found that the reservoir heterogeneity is essential to the formation of preferential flow channels. The channels are especially prone to be formed in the regions with high porosity and permeability at the initial time. The study can provide a theoretical reference for the optimal design of high-pressure water injection or frac-packing operation in the unconsolidated sands reservoir.
1. Introduction
There are a large number of unconsolidated sands reservoirs discovered in the Bohai Bay of China [1, 2], including Bozhong, Penglai, and Shengli oil fields. They are the major oil production areas in China. During long-term exploitation, several problems have arisen in these reservoirs, such as sand production and injectivity decline. To mitigate sand production, frac-packing has been widely used in these unconsolidated sand reservoirs [3], which can form high-permeability channels near the production well without causing the reduction of oil production too much. Although the porosity and permeability of the unconsolidated sands are high [4, 5], the poor water quality (such as suspended particles, bacteria, oil particles) can still lead to formation plugging and permeability reduction near the injection wells during produced water reinjection (PWRI), which causes the injectivity decline of reservoir [6]. For this reason, the fracturing injection regimes are commonly used to enhance injectivity. The fracture morphology will greatly affect the sweep efficiency of the waterflooding reservoir. Therefore, the evaluation of fracturing injection is very essential to optimize the waterflooding schedule. It can be seen that both frac-packing and PWRI involve channelization problem in the unconsolidated sands, so an in-depth understanding of channel formation during fluid injection in unconsolidated sands is essential for the high-efficient development of these reservoirs.
The strength of the unconsolidated sands is very low, and most of the uniaxial compressive strength is lower than 1-2 MPa [7, 8]. Sometimes there is even no cohesion, and only friction exists between solid particles. It is for this reason that the formation of channels in the unconsolidated sands during fluid injection is quite unlike the hard rocks with high strength [9]. For hard rocks, the hydraulic fractures as defined in linear elastic fracture mechanics are used to represent the channels [10]. However, for the unconsolidated sands, high-permeability channels instead of fractures are formed [11, 12]. This is because the large driving pressure gradient owing to the fluid injection leads to internal erosion in the unconsolidated sands, which is similar to the seepage erosion phenomena in soils [13–15]. The internal erosion causes the detachment and transport of the particles. As a result, the preferential flow paths within the formation are formed [16]. Based on experiments, some concepts such as viscous fingering, fluidization, and channelization had been proposed for unconsolidated sands. Johnsen et al. [17] investigated the displacement patterns in granular materials confined in a radial Hele-Shaw cell during central air injection. As injection pressure increases, the grains start to move, and fingers were observed in the granular materials. Huang et al. [18] conducted similar experiments by injecting an aqueous glycerine solution into the loose sands compacted in the Hele-Shaw cell. In their experiments, they observed that some channels which are similar to ramified finger are formed as the fluid viscosity and the injection velocity increase. Mahadevan et al. [19] used a Hele-Shaw cell filled with glass beads of different sizes to observe the occurrence of erosion and channelization during water injection. It is observed that only the smaller beads can be eroded and transport through the voids between the larger beads. Gago et al. [20] investigated the process of channel formations in soft sand confined into a radial Hele-Shaw cell. They found that the channel formations are greatly affected by fluid injection pressure. For low fluid injection pressures, it behaves as a solid porous material while for high enough injection pressures, grain rearrangement occurs.
Besides experiments, many researchers also focus on different numerical methods of modelling the channelization in unconsolidated sands. Most of them used the concept of internal erosion to interpret the formation of channels. One typical way is based on the continuum theory. The interaction between fluid flow and solid particles is used to model the process of internal erosion. For example, Stavropoulou et al. [21] proposed a mathematical model which accounts for the coupling between the fluid flow and mechanical damage to model the sand erosion. In this model, porosity is taken as the coupling parameters, and rock elasticity and cohesion become lower as porosity increases. They found that the failure of rock near the wellbore is closely related to sand erosion. Another attempt to model the channel formation in the unconsolidated sands was provided by Mahadevan et al. [19]. They proposed a continuum model that involves three phases (immobile solid phase, mobile solid phases, and fluid phase) and used the finite volume method to model the erosional channelization induced by fluid flow in a saturated porous medium. In their model, the internal erosion rate of the granular medium was calculated by assuming that the grain will be dislodged when the pressure gradient exceeds a threshold. Based on the model proposed by Mahadevan et al., Ameen and Taleghani [22] have studied the effect of injection rate, fluid viscosity, and rock properties on channel formation during fluid injection. Yerro et al. [23] attempted to use the material point method to model internal erosion process in bimodal internally unstable soils, which assumes that erosion rate is proportional to the velocity difference between fluid and solids. They performed a numerical test consisting of a soil column subjected to a vertical water flow and observed that the fine grains can be eroded and are able to move freely through the stable skeleton of the soil.
Instead of using the continuum theory, some researchers use the discrete element method (DEM) to model channel formation. For the DEM, Newton’s second law of motion is applied to determine the position of each particles under the driven of drag force of the pore fluid. For example, Harshani et al. [24] attempted to use the coupled DEM-LBM (Lattice Boltzmann Method) scheme to model the onset of erosion in granular materials at a microscale when the fluid pass through a solid matrix in an opposing direction to gravity. In his method, the DEM was used to model particles movement, while the LBM was used to model fluid flow at a mesoscale, and the coupled DEM-LBM scheme is effective to study the internal erosion phenomena. The results showed that the onset of erosion is mainly controlled by hydraulic gradient; Zhang et al. [25] used a discrete element method code PFC2D to model the particles movement during fluid injection. Sun et al. [26] investigated the mechanisms of opening-mode fracture initiation in granular media by the discrete element method and computational fluid dynamics. However, due to the expensive cost of calculation, the discrete element method can only be used to study the microscopic mechanism of channelization in a small scale.
At present, little researches focus on the channel formation near the wellbore caused by water injection or fracturing. To evaluate the effect of fluid injection on formation near the wellbore, a numerical simulation of the channel formation process in the unconsolidated sands during water injection was performed based on the principle of conservation of mass in continuum mechanics. The model considered the planar heterogeneity of reservoir properties. The erosion and sedimentation rate equations are used to describe the erosion and deposition of fine particles. The dynamic evolution process of the channels near the wellbore and its influencing factors are analyzed, which can give a theoretical reference for the optimal design of high-pressure water injection and frac-packing operation.
2. Mathematical Model for the Formation of Channels in Unconsolidated Sands
When the fluid flows through the pores of unconsolidated sands, it will produce hydraulic drag forces on the solid particles. These hydraulic drag forces mainly include volume forces and viscous forces. Meanwhile, some resistance on the particles also exists, such as cementing forces, contact forces, and adsorption forces. When the drag force of the fluid on the sand particles is greater than the resistance, the solid particles begin to move. For the hard rock, the solid particles are difficult to be detached because the cementing force is high. However, due to the low strength of the unconsolidated sands, fine particles are prone to be dislodged and migrate through the pore throat between the coarse particle skeleton. As a result, the seepage erosion occurs and leads to the change of physical properties and the formation of high-permeability channel in the unconsolidated sands [27].
The channels in the unconsolidated sands are the zones with high porosity and permeability. The evolution of porosity and permeability is closely related to the forces applied on the sand particles. However, due to the heterogeneity of the formation and the complexity of the pore structure, it is difficult to simulate the formation of channels at a microscale. In this paper, we rely on the continuous mechanics to perform the numerical simulation of formation of channels at a macroscale. Channelization is mainly the result of formation-porosity evolution driven by fluid flow as a function of space and time. The change of porosity in unconsolidated sands during fluid injection is used as an index to analyze the formation of channels.
For simplicity, the following three assumptions are required: (1) The solid particles and pore fluid are incompressible; (2) the pore fluid flow obeys Darcy’s law; and (3) the suspended particles flow together with the pore fluid, and both of them have the same velocity.
The governing equations mainly include mass conservation equations, fluid flow equation, porosity-permeability relation, erosion rate equation, particle deposition rate equation, and spatial porosity distribution equation of the reservoir.
2.1. Mass Conservation Equations
Because the unconsolidated sands is a porous medium, it can be divided into three phases, that is, immobile solid phase, mobile solid phase, and pore fluid phase when the porous fluid contains migrated solid particles [19]. Their volume fractions are expressed as , , and (the porosity of the rock is ), and the relationship between them satisfies
Based on the assumption of incompressibility of pore fluid and solid, mass conservation for the individual phases implies that [19] where is the erosion rate, is the deposition rate, and is the time. and denote the true velocity of particles and pore fluid, respectively.
2.2. Fluid Flow Equation
The pore fluid flow is described by Darcy’s law: where is the porosity, is the true pore fluid velocity, is the permeability, is the viscosity of pore fluid, and is the pore fluid pressure.
Thus, according to the continuity equation, the pore fluid flow equation for the incompressible fluid can be written as
2.3. Kozeny-Carman Relationship
The relationship between porosity and permeability is described by the classic Kozeny-Carman equation [28]: where is the average diameter of solid particles and is a constant that can be back-calculated from the actual permeability.
2.4. Erosion Rate Equation
The formation of high-permeability channels in unconsolidated sands is essentially caused by the erosion and migration of solid particles. The previous experiments have shown that there is a close relation between the erosion rate of solid particles and the fluid flow velocity. Erosion only occurs when the fluid flow velocity becomes larger than a critical threshold. Based on the erosion rate equation proposed by Shen et al. [29], a modified erosion rate equation is given by where is the erosion rate coefficient, is the stable porosity, and is the critical fluid flow velocity.
Equation (6) shows that there exists a critical fluid velocity for the estimate of detachment of sand particles. When , seepage erosion does not occur. Otherwise, the erosion rate increases linearly as the fluid velocity increases. We can also see that the erosion rate decreases with the increase of porosity. This is because the hydraulic gradient decreases as the pore throat diameter increases as a result of the increase of porosity. In addition, since only fine particles will migrate through the pore between the coarse particles [30], we add a stable porosity in the erosion rate equation which represents the maximum attainable porosity.
2.5. Deposition Rate Equation
When the solid particles move in the porous media, they may accumulate near the pore throat and become immobile particles again. This refers to the problem of deposition of solid particles in porous media. In practice, the deposition rate of solid particles is related to the concentration of solid particles, the size and shape of the pores, the particle size, the nature of the fluid, the flow rate, and so on [31]. For simplicity, we consider a relatively simple model to estimate the deposition rate, which assumes that the deposition of particles is a function of their concentration: where is the deposition rate coefficient.
2.6. Spatial Porosity Distribution Equation
The unconsolidated sand reservoirs are all heterogeneous. Fluid tends to flow through the areas with high porosity and high permeability, which may cause the local hydraulic scour and formation of high-permeability channels. The heterogeneity of the reservoir is very complex and can hardly be estimated accurately. Here, we use the Gaussian distribution to describe the porosity variation in the reservoir, which indicates that the probability of distribution is getting larger as the porosity approaches to the average value. The probability density function of porosity can be expressed as where is the standard deviation of porosity and is the average porosity.
3. Finite Difference Solution of Mathematical Model
Because of the heterogeneity and coupling between fluid and solids, it is not possible to obtain analytical solutions of these equations. The finite difference method is chosen for the analysis. For each time step, these equations should be solved simultaneously. First, the pore pressure at time is calculated by fluid flow equation, and then the fluid velocity can be determined by the calculated pressure. Subsequently, the erosion and deposition rate are calculated. , , and are then obtained by the mass conservation equations. At the end of each current time step, the porosity and permeability of the formation are updated for the calculation of the next time step.
The wellbore boundary has the maximum fluid velocity and is most sensitive to the initialization of the channel. Therefore, in order to accurately simulate the evolution of porosity near the wellbore, the mesh in the near-wellbore zone needs to be refined. Here, we use the geometric progression mesh to divide the domain where is the mesh size ratio and and are the radius at and , respectively.
To solve the model, the radial coordinate needs to be converted to the equidistant rectangular coordinate . Let , and we have where is the number of mesh.
Substituting Equation (10) into Equation (4), we obtain the pore fluid flow equation in terms of : where is the angle in the theta direction.
The explicit finite difference method is used to solve Equations (2) and (11) numerically. The corresponding finite difference discretization scheme is given in the Appendix. They can be solved by the Gauss-Seidel iterative method.
4. Model Validation
In order to verify the numerical model, porosity distribution of beads mixture obtained from the numerical simulation was compared with experimental results obtained by Mahadevan et al. [19]. Mahadevan et al. have conducted physical experiments in a chamber with the aim of demonstrating erosion and channelization in a saturated porous medium. The experiments were carried out in a vertical quasi-two-dimensional chamber based on a Hele-Shaw cell filled with mixture of glass beads. The mixture of beads with different radii was used: 60% of the initial volume has large beads with a diameter of 4 mm, and 24% has small beads with a diameter of 0.7 mm. Therefore, the corresponding maximum attainable absolute porosity is 0.4 since only fine beads can be eroded. The porosity of the beads mixture can be approximately described by the Gaussian distribution function. During the experiments, different inlet flow rates are applied uniformly to the lower boundary of the experimental domain to cause the erosion of beads. In the numerical simulation, the following parameters were used: (1) height is 30 cm, width is 30 cm, (2) inlet flow rates , (3) outer pressure is 0.1 MPa, (4) initial average porosity , standard deviation , (5) water viscosity , (6) average diameter of beads , , (7) stable porosity , (8) critical fluid flow velocity , and (9) deposition rate coefficient .
In the numerical simulation, the mesh size is 0.1 mm. The time step is 1 second. Figure 1 shows the comparison of probability density of normalized porosity after steady state is achieved between experimental and simulation results.
In Figure 1, the absolute porosity is normalized by the maximum attainable absolute porosity 0.4. It indicates that the numerical simulation result fits well with the experiments. Initially, the distribution is broadly peaked at 0.4. As time increases, both the numerical simulation and experiment show that the distribution broadens further due to seepage erosion and a second peak that represents completely eroded channels is developed at . The differences of porosity in some places between experimental and simulation results are caused by uncertainty of some parameters, especially the initial porosity distribution of the beads mixture.
5. Numerical Simulation of Channels Formation and Its Influencing Factors
A 2D model is performed to simulate the channelization during fluid injection. Since the formation of channels is due to the increase of porosity in some areas, the porosity changes are used to characterize the formation of channels in unconsolidated sands. The effects of injection pressure, critical fluid velocity, reservoir heterogeneity, and formation plugging near the wellbore on the formation of high-permeability channels are analyzed. In the numerical simulation, the borehole diameter is set to be 215.9 mm, and the outer boundary radius is 10 m. The number of mesh in radial and circumferential direction is 80 and 40, respectively. In order to obtain the evolution of porosity and permeability distribution near the wellbore more accurately, the nonuniform geometric progression mesh is used to divide the whole domain. Figure 2 shows the geometry and the geometric progression finite difference mesh used in the solution of the mathematical model. An enlarged view of mesh around the wellbore is shown at the upper right of the figure. The ratio of radial dimension between the adjacent elements is 1.05, and a total of 8000 meshes are used.
The fluid viscosity is 1 mPa s. The value of critical fluid velocity is . is set to be s-1. and . It is assumed that at the initial time, which means that the grains are immobile. The initial average porosity is 0.26, and the standard deviation is set to be 0.05.
If the flow boundary condition is specified on the inner boundary, the fluid flow velocity should be different for different nodes on the wellbore perimeter since the porosity and permeability are nonuniform around the wellbore. However, it is difficult to determine the accurate fluid flow velocity for different nodes. And inaccurate velocity will greatly influence the flow pattern, then affect the final channel. For this reason, the pressure boundary condition is chosen to be specified because the pressure of all nodes on the wellbore perimeter are the same in actual situation. For the outer boundary, a fixed pressure which is equal to the initial pore pressure is specified. In addition, it is assumed that the mobile particles will not be trapped at the outer boundary.
6. Results and Discussion
Let the injection pressure be 22 MPa, and calculate the change of formation porosity near the injection well after 1000, 5000, 10000, and 20000 s. When the porosity of some area is getting larger, it means that it is more likely to form high-permeability channel in such areas. Figure 3 shows the porosity distribution near the well after different injection times. The channel near the wellbore is shown in an enlarged view at the upper right of the figure ().
(a)
(b)
(c)
(d)
Figure 3 indicates that some large radial channels (red areas) appear around the wellbore for the given injection pressure. It can be seen that the area with increased porosity becomes larger and larger with time. At the beginning, the length of the channels increases more rapidly. However, after 10000 s, the length of channels only exhibits a slight increase since the pressure gradient decrease gradually when channels go far from the wellbore. It can be seen more clearly from the pore pressure distribution in Figure 4. It means that internal erosion rate is getting smaller with the growth of the channel. Figure 4 also shows that the slope of the pore pressure curve near the wellbore began to flatten as time increases. This is because the existence of channels increases the permeability of the formation near the wellbore. As a result, the fluid injectivity of the formation is enhanced. If the area with a porosity greater than 0.5 is regarded as the channel, the length of the longest channel at are 0.22, 0.58, 0.63, and 0.82 m, respectively. It can be seen that the length of the channel in the formation around the wellbore is getting larger and larger as the time increases.
To observe the initiation and evolution of the channels near the wellbore, the porosity distribution profiles surrounding the wellbore at at different times are shown in Figure 5.
It can be seen from Figure 5 that the porosity at the initial time () obeys the Gaussian distribution. After the fluid is injected, the porosity gradually increases with time since the value of the seepage velocity at is greater than the value of critical fluid velocity. 11 porosity peaks with value greater than 0.5 occurs, which indicates that 11 large channels are formed near the wellbore. These channels are the preferential flow paths generated by fluid injection, and the distribution of these dominant channels is coincident with the high-porosity area at the initial time. It shows that the channels are prone to be formed along the initial high-porosity areas in unconsolidated sands where porosity is larger than 0.26.
6.1. Effect of Injection Pressure on the Formation of Channels
The formation porosity distributions near the wellbore for different injection pressures (16, 18, 20, and 21 MPa) when are shown in Figure 6. The channels near the wellbore are shown in an enlarged view at the upper right corner ().
(a)
(b)
(c)
(d)
When the injection pressure is 16 MPa, there is no obvious change in porosity near the wellbore, indicating that no obvious channels are formed. When the injection pressure is increased to 18 MPa, several high-porosity channels begin to appear near the wellbore due to the increase of the seepage velocity. As the injection pressure continues to increase to 20 MPa and 21 MPa, both the number and length of the channels increase. This is because the pressure gradient increases with the injection pressure. As a result, the solid particles are more likely to be dislodged and form the high-porosity channels.
6.2. Effect of Critical Fluid Velocity on the Formation of Channels
The critical fluid velocity refers to the minimum seepage velocity required for the onset of solid particles erosion. It depends strongly on the strength and in situ stress of unconsolidated sands. Higher strength and in situ strength lead to larger critical fluid velocity. Its magnitude may affect the formation of channels. Figure 7 shows the porosity distribution near the wellbore for different values of critical fluid velocity (, , , and ) after 20000 s when the injection pressure is 20 MPa.
(a)
(b)
(c)
(d)
As can be seen from Figure 7, the length of the channels is getting shorter as the critical fluid velocity increases. In addition, the number of the channels also decreases with the increase of the critical fluid velocity. This is because the increase of critical fluid velocity makes the fluid flow velocity required for the erosion of solid particles to enhance. Thus, it is becoming increasingly difficult to form the channels. If the critical fluid velocity is greater than the highest fluid flow velocity in the formation, large channels will not be formed.
6.3. Effect of Heterogeneity on the Formation of Channels
The heterogeneity of porosity will affect the pore fluid flow field in the unconsolidated sands [32]. Thus, it may have significant influence on the formation of channels. In the foregoing examples, we have obtained the porosity distribution around the wellbore in the heterogeneous formation. For comparison, we also calculate the porosity distribution near the wellbore for the homogeneous formation, which is shown in Figure 8. In this example, the porosity is set to be 0.26, and the injection pressure is 18 MPa. To observe the evolution of porosity near the wellbore and the formation of channels more clearly, we take porosity distribution profiles when at , which are shown in Figure 9.
(a)
(b)
The comparison analysis of Figures 3, 8, and 9 shows that the heterogeneity of formation properties is a necessary condition for the formation of channels. For the homogeneous formations, fluid injection will only cause uniform erosion. As a result, it increases the porosity uniformly around the wellbore. It is impossible to form a channel along a certain direction (see Figures 8 and 9; the porosity is uniformly distributed around the wellbore after different times; the porosity is all 0.26 at the initial time and increases to 0.41 at and 0.48 at ). This is because the erosion rate around the wellbore is the same for homogeneous formations. On the contrary, when the porosity obeys the Gaussian distribution, the particle erosion rate differs at different area. That is, nonuniform erosion occurs in the formation. The fluid flow velocity is higher in the area with a higher initial porosity. As a result, the channels are prone to be formed in such area (Figure 5).
6.4. Effect of Plugging on the Formation of Channels
The formation near the wellbore is easily blocked under the pollution of various working fluids (drilling fluid, completion fluid, etc.) or long-term water injection conditions, which may cause the decrease of porosity and permeability near the wellbore. To assess the effect of the influence of plugging near the wellbore on the formation of channels, the porosity of the formation within 0.5 m of the borehole is reduced uniformly (the average porosity of the plugging zone is set to 0.2, 0.15, and 0.1, respectively; the standard deviation of porosity is 0.03). Figure 10 shows the porosity distribution around the wellbore after 40000 s when the injection pressure is 20 MPa.
(a)
(b)
(c)
If the average formation porosity near the wellbore reduces to 0.2, multiple short channels are formed near the wellbore after 40000 s during fluid injection. However, it does not penetrate through the plugging zone. When the average porosity of the plugging zone is reduced to 0.15, the channels can still be formed, but the length of channels is shorter, and the number of channels is less. When the average porosity is reduced to 0.1, it is difficult to form obvious channel due to serious blockage near the wellbore. This is because the flow velocity is lower than the critical fluid velocity, and the skeleton particles cannot be dislodged. It implies that the bottomhole pressure needs to be enhanced to form the channel by increasing the fluid flow velocity. To study the effect of increasing injection pressure on the formation of channels in the plugging zone, the injection pressure is increased to 21.5 MPa. Figure 11 shows the porosity distribution near the wellbore after 40000 s.
As can be seen in Figure 11, when the injection pressure increases to 21.5 MPa, two channels that penetrate through the plugging zone are formed. The results illustrate that plugging could be mitigated by enhancing the injection pressure. Compared to the unblocked formation, the injection pressure required to form the channels in the blocked formation is more higher. If we consider this injection pressure as the fracture initiation pressure, it is concluded that the plugging zone near the wellbore can enhance the fracture initiation pressure.
7. Conclusions
A numerical model for simulating the channel formation process in the unconsolidated sand process during fluid injection was performed based on the principle of conservation of mass in continuum mechanics. Sensitivity analyses are conducted for this model to understand how different parameters may affect the creation and propagation of channels. The main conclusions are as follows: (1)The unconsolidated sands are cohesionless, and the strength is very low. The fracture patterns in the unconsolidated sands are distinctly different from fractures in the hard formation during fluid injection. The seepage erosion caused by hydraulic scour plays a leading role in the formation of channels in unconsolidated sands(2)Channelization might occur when the fluid-fluid velocity becomes locally larger than a critical threshold, and small grains are dislodged and move away. So the channelization depends largely on the injection pressure and the critical fluid velocity. As the injection pressure increases, the solid particles are more likely to be dislodged and form the high-porosity channels. However, due to the increasing difficulty of particle detachment, the length and number of the channels decrease with the increase of the critical fluid velocity. Below the critical injection rate, which depends on the initial porosity distribution in the medium, fluid will flow into the formation with limited damage near the injector, but no erosion or channelization occurs(3)The plugging of the formation near the wellbore will cause the decrease of seepage velocity. It is more difficult to induce channels under the same injection pressure. However, the simulation results show that the plugging zone can be penetrated by enhancing the injection pressure, which indicates that water injection at high pressure can contribute to the plug removal(4)The heterogeneity of formation properties is a necessary condition for the onset of channelization. Different initial porosity leads to different final-porosity distribution, and the channels are prone to be formed in area with initial higher porosity.
The currently established mathematical model for the formation of channels does not consider the effect of non-Darcy flow on particle seepage erosion. In addition, the effect of in situ stress on erosion is also not considered. To simulate the channelization in unconsolidated sands more realistically, further study is needed to establish a mathematical model including the effect of non-Darcy seepage and in situ stress.
Appendix
Finite Difference Scheme for Solving the Equations
The explicit finite difference method is used to discretize Equations (2) and (11). For the two-order spatial derivative terms, the central difference is used. Thus, the finite difference scheme of Equation (11) can be written as where the superscript denotes the time and the subscripts and denote the parameter at and .
Using the fluid flow equation, combined with the Kozeny-Carman relation, erosion rate equation, and the deposition rate equation, the finite difference schemes of the mass conservation equations are
Nomenclature
| : | Mesh size ratio |
| : | A constant that can be back-calculated from the actual permeability |
| : | Deposition rate (1/s) |
| : | Average particle diameter (m) |
| : | Erosion rate (1/s) |
| : | Permeability (m2) |
| : | Pore fluid pressure (Pa) |
| : | Radius at and (m) |
| : | Standard deviation of porosity |
| : | Time (s) |
| : | True velocity of particles and pore fluid (m/s) |
| : | Seepage velocity (m/s) |
| : | Critical fluid flow velocity (m/s) |
| : | Deposition rate coefficient (1/s) |
| : | Theta coordinate (rad) |
| : | Erosion rate coefficient (1/m) |
| : | Dynamic viscosity of fluid (Pa s) |
| : | Porosity |
| : | Average porosity |
| : | Volume fraction of mobile particles, immobile particles, and pore fluid |
| : | Stable porosity. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by supported by the Hainan Provincial Natural Science Foundation of China (No. 2018CXTD346) and National Key Research and Development Program of China (No. 2019YFC0312301).