A Research on the Effect of Heterogeneities on Sandstone Matrix Acidizing Performance

Matrix acidizing is one of the common methods to enhance production in sandstone reservoirs. Conventional acidizing designs generally neglected the effect of heterogeneities of mineral and flow field distributions both in areal and vertical directions and assumed that the acid front propagates with a piston-like style. However, sandstone formations inevitably have small-scale heterogeneities of minerals and flow properties that may give rise to acid propagation in a manner much different from what is predicted based on homogeneous assumptions. In this paper, we conduct a research to numerically investigate how the heterogeneities affect acidizing performance under in situ conditions. Firstly, a heterogeneity model is built for mineral and porosity distributions by using the semivariogram model of geological statistics, based on which we generate spatially correlated porosity and mineral distributions. Next, a model of radial acid flooding is developed based on mass balance and the chemical reactions between the acids and minerals occurring during the acidizing process. The model is numerically solved to investigate the permeability response, acid distributions, precipitate distributions, and the effect of the heterogeneities on acidizing. The results show that the heterogeneities both in areal and vertical directions have a significant effect on acidizing. The flow field heterogeneities have a more serious impact than the mineral heterogeneities. In a plane, strong porosity heterogeneity can give rise to acid fingering and even channeling, which make the acid penetration distance longer than the homogeneous cases. The secondary precipitate has a significant effect when fast-reacting mineral content is high. Vertically, several-fold permeability contrast creates the acid break through the high-perm zone leaving the low-perm zone understimulated. Both flow field and mineral heterogeneities make it possible to create high-permeability channels during the acidizing process and to obtain a longer acid penetration distance.


Introduction
Drilling, completion, or some other well operations inevitably damage formations, which may seriously decrease the well productivity in sandstone reservoirs. Matrix acidizing is a common method to remove the formation damage and recover well productivity. In acidizing, the acid flows into the porous media, reacts with the minerals, and increase formation porosity as well as permeability. Many researchers conducted studies on acid flooding in experiments. One of the unignorable risks in acidizing in sandstone is the secondary precipitate. Crowe [1], Almond et al. [2], and Underdown et al. [3] investigated the conditions that produce precipitate.
They found that silica gel is produced at a slow rate due to the slow reaction at room temperature. When the temperature is higher than 55°C, the fast reaction will generate a lot of precipitate and damage the formation seriously.
Due to the small size of the core sample and the complexity of parallel core flooding, field-scale simulation and design necessitates numerical modeling. Sandstone acidizing models include the capillary model, the microscopic model, the kinetic model, the lumped parameter model, the distributed parameter model, the one-acid two-mineral model, the two-acid three-mineral model, and the generalized geochemistry model [4][5][6][7][8][9][10], most of which made the assumption of homogenous formation. The two-acid three-mineral model developed by Bryant [4] takes into account the effect of precipitates and H 2 SiF 6 and can simulate well the physical phenomenon of acidizing. At present, the most sophisticated model [9] is the generalized geochemistry model. The model accounts for the reaction of multiple kinds of minerals and the acid and can analyze the effect of mineral content, reactant, product, and reaction kinetics. They found that dissolution and precipitation do not occur simultaneously. Li et al. [11] built a linear acid flooding model to simulate linear core flooding at an experimental scale with the consideration of the heterogeneities of minerals and flow fields. Li et al. [12] built a radial acidizing model and considered the effect of the temperature field on the reaction, but they did not consider the effect of the heterogeneities. Leong et al. [13,14] did a modeling of the acidizing of HBF 4 in linear cores with COMSOL.
In this paper, we modeled radial acid flooding based on a two-acid, three-mineral model and considered the heterogeneities of minerals and flow fields in both planar direction and vertical directions. Based on the model, we did extensive numerical simulations to analyze the effect of heterogeneities on acidizing performance. Also considered are multilayer flooding and secondary precipitates.
(1) Cylindrical coordinate system is used for the radial flooding (3) Hydrofluoric acid concentration distribution equation (4) Fluosilicic acid concentration distribution equation where ∅ is the porosity; C 1 is the hydrofluoric acid concentration; C 2 is fluosilicic acid concentration; u is the velocity vector (u r is the velocity in the radial direction, and u θ is the velocity in the circumferential direction); V j is the mineral volume fraction; j is the fast-reacting mineral, slow-reacting mineral, and silica gel; MW i is the mole weight of acid, S * j is the specific surface area of minerals, i is the hydrofluoric acid or fluosilicic acid; E f ,i,j is the reaction rate constant; α is the reaction order; ρ is the density of minerals; and β i,j is the gravimetric dissolving power.  p r, θ, t = 0 = p r , where V 0 1 , V 0 2 , and V 0 3 are initial mineral 1, mineral 2, and mineral 3 distribution, respectively. ∅ 0 r,θ is initial porosity distribution. p r is reservoir pressure.
The boundary conditions are as follows: p r e , θ = p e , where Q inj is pump rate. C 0 1 and C 0 2 are the injection concentrations of acids 1 and 2, respectively. p e is the pressure of the outer boundary.

Porosity and Mineral Distributions
3.1. Porosity Distribution. Change of porosity directly affects the permeability of the sandstone formation, thus affecting the distribution and reaction rate of acid fluid. Therefore, it is very important to generate the porosity distribution of a real formation for the accuracy of the whole model. In terms of geostatistics, porosity distribution shows characteristics of directionality [15][16][17]; that is, they are spatially correlated instead of completely random. A semivariogram model [8,18], which is an important tool in geostatistics to describe spatially correlated distributions, is used to generate porosity distribution. We use the geostatistical software GSLIB [19], which incorporates a semivariogram model for spatial correlation, to generate spatially correlated numbers, with which we generate porosity distribution as ϕ = 1, ϕ ≥ 1 ϕ 0 + ϕ 0 c vĜ l x , l y , 0 01 < ϕ < 1 0 01, ϕ ≤ 0 01 , 10

Geofluids
whereĜ l x , l y is the spatially correlated number output from GSLIB and l x and l y represent the dimensionless correlation lengths of x direction and y direction, respectively. c v is the coefficient of variation of porosity. Bigger c v means strong inhomogeneity; otherwise, it means slight inhomogeneity. ϕ 0 is the average porosity. Figures 2 and 3 show the porosity distributions for two layers with ϕ = 0 15. The figures indicate that porosity distributions vary a lot with different geostatistic parameters so as to affect flow fields significantly.

3.2.
Pore-Permeability Correlation. Acid-rock reaction increases porosity, resulting in a change of permeability, which in turn affects acid flow. We use the following model to describe the relationships [20] between permeability and local porosity.
where β is a constant which is dependent on the structure of the medium.

Mineral Distribution.
Using the same method for porosity distribution, we generate the spatial association distribution function which is defined as follows:   where M1 and M2 are the average volume fractions of the fast-reacting mineral and slow-reacting mineral, respectively. Figure 4 is an example of the distributions of the volume fraction (V j ) of fast-and slow-reacting minerals.

Numerical Solution
The equations are discretized with the finite volume method and solved sequentially. Firstly we generate initial porosity and mineral distributions. Then, we solve equations (2) and (3) to obtain the pressure and velocity fields. Next, we solve equations (4) and (5) to get the acid concentration distribution. Finally, we update the mineral content and porosity distributions based the acid-rock reaction (equations (6) and (7)).

Result and Discussion
In this section, extensive numerical simulations are conducted to analyze acid distribution patterns as well as the effect of heterogeneities on acidizing performance with the parameters in Table 1. And acid-rock reaction parameters are determined through references [21][22][23].

Comparison of Homogenous and Heterogeneous Cases.
The parameters used in the simulation are listed in Table 2.
Other parameters are kept the same for the two cases. Figure 5 shows acid concentration and mineral content distribution in the radial direction. Since they are different in different directions for the heterogeneous cases, only 0 degree (relative to the black line in Figure 5, for example) is chosen for displaying. VF_M1, VF_M2, and VF_M3 stand for fast-reacting mineral, slow-reacting mineral, and silica gel precipitate content. Compared to the homogeneous case, the curves of the heterogeneous case fluctuate apparently, the acid front is longer, and precipitation is deeper. Since the fast-reacting mineral has a much higher reaction rate than the slow-reacting mineral, most the fast-reacting mineral is spent, but a small amount the slow-reacting mineral is dissolved in the live acid zone.   6 Geofluids gives a regular shape of distribution. Meanwhile, the heterogeneous case gives rise to an irregular distribution shape. The porosity distribution is strongly correlated in the horizontal direction. The horizontal direction is the leading flow direction, so acid penetrates deeper in this direction. The acid fingers through the high-perm zone, leaving the acid front to propagate nonuniformly. Some low-perm zones are bypassed by the acid. For the homogeneous case, from the wellbore to the outer zone, there exist zones of high perm, low perm, and original formation perm. Near the wellbore zone, the permeability is enhanced due to rock dissolution by the acid. In the live acid zone, most of the fast-reacting mineral is dissolved. The low-perm zone is cause by the secondary precipitate of spent acid. Silica gel precipitation due the reaction between H 2 SiF 6 and the fast-reacting mineral is inevitable [24], which is unfavorable for sandstone acidizing. For the heterogeneous case, the reduced perm zone is not distinctly apparently. The precipitate does exist, but the drastic porosity variation from one place to another makes it hard to identify. The significant difference of the penetration distance between the homogeneous and heterogeneous cases shows just how important it is to consider heterogeneity in field treatments. The acid penetration distance is much longer than the one predicted by the homogenous model, especially when there exists major flow channels.

Effect of Vertical Heterogeneity.
A well often consists of multilayers with much different properties. For this kind of well, an important aspect in acidizing is acid placement. Here, the difference of the properties of the layers is called vertical heterogeneity. Vertical heterogeneity considers the difference of permeability and the mineral mainly in this simulation.

Effect of Vertical Permeability
Heterogeneity. The vertical permeability difference of layers can be characterized by the permeability ratio, which is defined as the ratio of the average permeability of layer A to layer B as follows: where α β is the permeability ratio and k avgA and k avgB represent the average permeabilities of layer A and layer B, respectively. We did simulations for α β = 2, 5, and 10, respectively. Parameters are listed in Table 3. Except for the different porosity and permeability distribution values, other parameters are kept the same. Results after 90-minute acid injection are shown as follows. As shown in Figure 11, vertical permeability heterogeneity has a notable effect on the acid distribution of the vertical profile. The higher the permeability ratio, the more seriously acid intake differs. When ∝ β = 10, the acid entering into the low-perm zone is negligible compared to that entering into the high-perm zone. Except for the permeability difference which causes the discrepancy of acid intake, the permeability increase due to rock dissolution also magnifies the difference. The high-permeability zone takes more acid and raises permeability faster, so it accepts more and more acid. This demonstrates the necessity of acid diversion in acidizing multilayer formations. Figure 12 shows fast-reacting mineral dissolution in two layers. In line with the acid distribution, the dissolution is greater in the higher perm layer, and the dissolutions in the lower perm layer decreases with the increase of the permeability ratio. Figure 13 shows the distribution of precipitation. In the higher perm layer, the precipitate spreads wider because  Figure 10: Silica precipitation distribution after acidizing for homogeneous and heterogeneous initial distributions. 7 Geofluids of more spent acid generated by the acid entering into the layer. Precipitate will decrease the permeability; more precipitate in the higher perm layer will counteract the permeability increase caused by rock dissolution, but due to areal heterogeneity, the acid can channel through the lowered perm zone.   8 Geofluids Figures 14-16 show the results after acidizing for 90 minutes. Since the two layers have identical initial porosity and permeability distribution, the dissolution geometry is similar, but the extent of layer A is notably larger than layer B. Layer B has a much higher fast-reacting content than layer A. Due to the relative fast reaction rate of the fast-reacting mineral, most of the fast-reacting mineral is dissolved in the area covered by the live acid. The acid front propagates much slower in layer B. Figure 17 shows the pump rate allo-cation with time for the two layers. The initial rate is the same. As acidizing proceeds, the rate of layer A increases and the rate of layer B decreases, since the total rate is fixed. There are two reasons for this phenomenon. The first is that the acid front moves faster in layer A than layer B. The second is more serious precipitation in layer B due to the higher content of the fast-reacting mineral. Even though the vertical heterogeneity of the mineral content is obvious, its effect is less remarkable than the effect of permeability.  9 Geofluids 5.5. Effect of Areal Heterogeneity. Areal heterogeneity can be expressed by the coefficient of variation of porosity (c v ). The higher the c v , the more serious the heterogeneity is. c v = 0 means homogeneous case. Figures 18-20 show the results for c v = 0 1, 0 3, 1. The dimensionless correlation lengths are 0.08 and 0.0001 in the x and y directions. c v = 0 1 means very weak heterogeneity. The porosity distribution is very close to the homogeneous case. c v = 1 means very strong heterogeneity, which generates some leading flow paths similar to the channels in the reservoir. The more serious the heterogeneity, the more irregular the geometry of acid coverage is. The case of c v = 1 makes the acid finger through some highperm paths, making the acid penetration distance much longer than the other two cases. Another areal heterogeneity is natural fractures. For example, there is a fracture connecting to the wellbore as shown in Figure 21. The acid mainly flows through the fracture to reach a deeper distance. These indicate that in strong heterogeneous reservoirs, the acid  10 Geofluids penetration distance may be much deeper than the one predicted by the conventional homogeneous model.

Conclusion
In this paper, we built a sandstone matrix acidizing model accounting for areal and vertical heterogeneities and the effect of secondary precipitates. Based on the model, extensive numerical simulation was conducted to investigate the effect of heterogeneity on acidizing performance. The study reached the following conclusions: (1) The heterogeneities both in areal and vertical directions have a significant effect on acidizing     Figure 19: HF distribution after acidizing in two layers.  12 Geofluids performance. The areal heterogeneity influences the acid dissolution pattern and live acid penetration distance. Heterogeneous cases give longer live acid penetration distance than the homogeneous cases. The vertical heterogeneities determine the vertical acid intake profile (2) Both flow field and mineral heterogeneities make it possible to create high-permeability channels in acidizing and to obtain longer acid penetration distance. The heterogeneity of a flow field has a more serious effect on the acid dissolution pattern than the mineral heterogeneities (3) Serious vertical heterogeneities make an unfavorable acid intake profile, with the high-permeability zone taking excessive acid and the low-permeability zone taking insufficient acid. In multilayer formations, heterogeneity creates acid fingering instead of a piston-like pattern, which indicates the necessity of an acid diversion (4) The secondary precipitate has a significant effect on acidizing performance when fast-reacting mineral content is high Nomenclature C 1 : Hydrofluoric acid concentration C 2 : Fluosilicic acid concentration C 0 1 : Injection concentration of acid 1 C 0 2 : Injection concentration of acid 2 c v : Coefficient of variation of porosity E f ,i,j : Reaction rate constant G l x , l y : Spatially correlated number output from GSLIB k avgA : Average permeability of Layer A k avgB : Average permeability of Layer A l x , l y : Dimensionless correlation lengths of x direction and y direction, respectively MW i : Mole weight of acid M1, M2: Average volume fractions of fast-reacting mineral and slow-reacting mineral, respectively Q inj : Pump rate p r : Reservoir pressure p e : Pressure at outer boundary S * j : Specific surface area of minerals u: Velocity vector u r : Velocity in the radial direction u θ : Velocity in the circumferential direction V j : Mineral volume fraction V 0 1 , V 0 2 , V 0 3 : Initial distribution of the volume fractions of mineral 1, mineral 2, and mineral 3, respectively α: Reaction order α β : Permeability ratio β i,j : Gravimetric dissolving power β: Constant dependent on the structure of the medium ∅: Porosity ∅ 0 r,θ : Initial porosity distribution ϕ 0 : Average porosity ρ : Density.

Data Availability
The numerical simulation data used to support the findings of this study are available from the corresponding author upon request.